GEE23:基于植被物候实现农作物分类

2024-03-31 02:28

本文主要是介绍GEE23:基于植被物候实现农作物分类,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

地物分类

  • 1. 写在前面
  • 2. 北京作物分类

1. 写在前面

  今天分享一个有意思的文章,用于进行农作物分类。文章提出了一个灵活的物候辅助监督水稻(PSPR)制图框架。主要是通过提取植被物候,并自动对物候数据进行采样,获得足够多的样本点,再使用随机森林等机器学习方法进行分类。这种方法有效解决了样本量不足或者样本位置不够精确的问题,并且分类结构相较于之前的方法更高。我认为这是一种比较有意思的文章,当然这种方法还可以用到其他植被类型分类中。

灵活的物候辅助监督水稻(PSPR)制图框架:
在这里插入图片描述

2. 北京作物分类

   首先,使用Landsat5、7、8数据获取植被物候信息,再提取随机采样点。
  各植被类型的物候:
在这里插入图片描述

/**************************使用PSPR方法生成随机采样点*************************/
// PSPR: 灵活的物候辅助监督水稻制图框架
//设置研究区位置: 北京
var table = ee.FeatureCollection("users/cduthes1991/boundry/China_province_2019");
var BJGrid4 = table.filter(ee.Filter.eq('provinces','beijing'));
var roi = BJGrid4;
Map.addLayer(roi,{'color':'grey'},'roi',false);
Map.centerObject(roi,7);var GridTest = GridRegion(BJGrid4.geometry(),6,6).filterBounds(roi);
print("GridTest size:",GridTest.size());
var color = {'color':'0000FF','fillColor':'FF000000'};
Map.addLayer(GridTest.style(color),null,'GridTest');/********* cropland mask**************************************************/
//这个数据为中科院的30LUCC数据,其中11和12分别表示水田和旱地
var CAS_LULC_2018 = ee.Image("users/chengkangmk/Global-LULC-China/LULC_CAS/CAS_LULC_30m_2018");
var cropland = CAS_LULC_2018.eq(11).or(CAS_LULC_2018.eq(12)).clip(roi);
Map.addLayer(cropland.randomVisualizer(),null,'cropland',false);/*******************************自定义函数*********************************/
// Landsat 4, 5 and 7 去云
function rmL457Cloud(image) {var qa = image.select('pixel_qa');// If the cloud bit (5) is set and the cloud confidence (7) is high// or the cloud shadow bit is set (3), then it's a bad pixel.var cloud = qa.bitwiseAnd(1 << 5).and(qa.bitwiseAnd(1 << 7)).or(qa.bitwiseAnd(1 << 3));// Remove edge pixels that don't occur in all bandsvar mask2 = image.mask().reduce(ee.Reducer.min());var mask3 = image.select("B1").gt(2000);return image.updateMask(cloud.not()).updateMask(mask2).updateMask(mask3.not()).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start",'system:time_end']);
}// Landsat-8 去云
function rmL8Cloud(image) { var cloudShadowBitMask = (1 << 3); var cloudsBitMask = (1 << 5); var qa = image.select('pixel_qa'); var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) .and(qa.bitwiseAnd(cloudsBitMask).eq(0)); var mask2 = image.select("B2").gt(2000);                 return image.updateMask(mask).updateMask(mask2.not()).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start",'system:time_end']);
}// Sentinel-2 去云
function rmS2cloud(image) {var qa = image.select('QA60');// Bits 10 and 11 are clouds and cirrus, respectively.var cloudBitMask = 1 << 10;var cirrusBitMask = 1 << 11;// Both flags should be set to zero, indicating clear conditions.var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));var mask2 = image.select("B2").lte(2000);return image.updateMask(mask).updateMask(mask2).toDouble().divide(1e4).copyProperties(image).copyProperties(image, ["system:time_start", "system:time_end"]);
}// 计算相关指数,包括NDVI、LSWI(植被水分含量指数)、EVI
function addIndex(image){// original bandsvar blue = image.select('blue'); var red = image.select('red');var green  = image.select('green');var nir = image.select('nir');var swir1 = image.select('swir1');var ndvi = image.normalizedDifference(["nir", "red"]).rename("NDVI").toDouble();var lswi = image.normalizedDifference(["nir", "swir1"]).rename("LSWI").toDouble();var lswi2ndvi = lswi.subtract(ndvi).rename("LSWI2NDVI").toDouble();var evi = image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {'NIR': nir,'RED': red,'BLUE':blue}).rename("EVI");return image.addBands(ndvi).addBands(lswi).addBands(lswi2ndvi).addBands(evi);
}// 为影像的特定波段指定名称
var LC8_BANDS = ['B2',   'B3',    'B4',  'B5',  'B6',    'B7']; //Landsat 8
var LC7_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7']; //Landsat 7
var LC5_BANDS = ['B1',   'B2',    'B3',  'B4',  'B5',    'B7']; //Llandsat 5
var S2_BANDS  = ['B2',   'B3',    'B4',  'B8',  'B11',   'B12']; // Sentinel-2
var STD_NAMES = ['blue', 'green', 'red', 'nir', 'swir1', 'swir2'];// 定义时间变量、LSWI阈值、采样点数量
var year = '2015';
var th_lswi2NDVI = 0;
var pointNum = 1000;// 导入数据
// landsat 8
var l8Col = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR').map(rmL8Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01')// .filter(ee.Filter.lte('CLOUD_COVER',10)).select(LC8_BANDS, STD_NAMES); 
// landsat 7 
var l7Col = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR').map(rmL457Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01').select(LC7_BANDS, STD_NAMES); 
// landsat 5  
var l5Col = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR').map(rmL457Cloud).filterBounds(roi).filterDate(year+'-04-01',year+'-11-01').select(LC5_BANDS, STD_NAMES); var L578COl = ee.ImageCollection(l8Col.merge(l7Col).merge(l5Col)).sort("system:time_start").map(addIndex);
print("L578COl: ", L578COl);//****************************************************************************
//***************************** 基于LSWI的水稻制图****************************
//****************************************************************************
var palette1 = {min: 0, max: 1.0, palette: ['000000', '00FF00']};// 1.水淹阶段,定义植被的时间窗口
var imgCol_flood = L578COl.filterDate(year+'-05-01',year+'-06-30') .filterBounds(roi);// 使用 LSWI-NDVI 提取水稻
// "rice_0_lswi": 表示初步产生水稻的时间
var rice_0_lswi =  imgCol_flood.select("LSWI2NDVI").map(function(image){return image.select("LSWI2NDVI").gt(th_lswi2NDVI);  
});// 从rice_0_lswi图像集合中获取每个像素位置上的LSWI指数最大值(找到最大的 LSWI 值)。
// LSWI是一种用于检测土地表面水体的指数,其数值通常与水体的存在程度成正比。在稻田中,LSWI 的最大值对应于稻田灌溉时水体充足的情况。因此这样可以最大程度的保证水稻位置的获取。
var rice_0_lswi = rice_0_lswi.reduce(ee.Reducer.max()).clip(roi).updateMask(cropland);
Map.addLayer(rice_0_lswi, palette1, "rice_0_lswi candidates",false);//2.生长阶段,水稻幼苗被移植到稻田中,直到水稻成熟
var imgCol_growth =  L578COl.filterDate(year+'-07-01',year+'-10-31').filterBounds(roi);// 不同阶段图像
var PalettePanel = {bands:["swir1","nir","red"],min:0,max:0.3};var imgCol_flood_qmosaic = imgCol_flood.qualityMosaic('LSWI2NDVI').clip(roi);
Map.addLayer(imgCol_flood_qmosaic, PalettePanel, 'imgCol_flood_qmosaic');var imgCol_growth_qmosaic = imgCol_growth.qualityMosaic('NDVI').clip(roi);
Map.addLayer(imgCol_growth_qmosaic, PalettePanel, 'imgCol_growth_qmosaic', false);// 使用 0 值来填充rice_0_lswi的区域
var rice_0_map = rice_0_lswi.unmask(0).clip(roi);
Map.addLayer(rice_0_map.selfMask(), {"palette":'#FF0000'}, "rice_0_map",false);//****************************************************************************
//*********************基于CCVS 方法提取水稻物候******************************
//****************************************************************************
var visParam = {min: -0.2,max: 0.8,palette: 'FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400,' +'3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301'
};// 使用LSWI_max, LSWI_min, NDVI_max, NDVI_min来推导RCLN,其中RCLN表示LSWI相对于与NDVI的变化幅度
var RCLE = imgCol_growth_qmosaic.select("LSWI").subtract(imgCol_flood_qmosaic.select('LSWI')).abs().divide(imgCol_growth_qmosaic.select("EVI").subtract(imgCol_flood_qmosaic.select('EVI')).abs()).rename("RCLE");
Map.addLayer(RCLE, visParam,'RCLE', false);// 创建一个布尔类型的图像,其中大于 0.1 的像素被设置为 true,小于或等于 0.1 的像素被设置为 false
var LSIW_min = imgCol_flood.select('LSWI').min().gt(0.1).clip(roi).updateMask(cropland);
Map.addLayer(LSIW_min.selfMask(), {"palette":'#FF0000'}, "LSIW_min",false);// 
var RCLE_rice = RCLE.updateMask(RCLE.gt(0)).updateMask(LSIW_min).lt(0.6).unmask(0).clip(roi); 
Map.addLayer(RCLE_rice.selfMask(), {palette: 'green'}, "CCVS_rice",false);var riceCombine = RCLE_rice.add(rice_0_map);// 将以上两种方法得到的数据图像进行合并求交
var stableMask = riceCombine.eq(0).or(riceCombine.eq(2));
var LULU_mutual = riceCombine.updateMask(stableMask).where(riceCombine.eq(2),1).rename("RiceMap").updateMask(cropland);
Map.addLayer(LULU_mutual,palette1,'LULU_mutual',false);var riceMapCol = GridTest.toList(GridTest.size()).map(function(ROIFea){// set study areavar roi = ee.FeatureCollection([ROIFea]);var ricemap = riceTrainData(roi);return ee.FeatureCollection(ricemap);
});
var samplePoint = ee.FeatureCollection(riceMapCol).flatten();
Map.addLayer(samplePoint,null,'samplePoint',false);// check the generated sample data
print("samplePoint:",samplePoint.limit(10));var ricePoint_1 = samplePoint.filter(ee.Filter.eq('landcover',1));
Map.addLayer(ricePoint_1,{'color':'#FFA500'},'ricePoint_1');
print("ricePoint_1 size",ricePoint_1.size());var NonricePoint_1 = samplePoint.filter(ee.Filter.eq('landcover',0));
print("NonricePoint_1 size",NonricePoint_1.size());Export.table.toAsset(samplePoint,"PSPRSampleGeneration"+year,"PSPRSampleGeneration"+year)function riceTrainData(roiRegion){var roi = roiRegion;var LULU_mutual2 = LULU_mutual.clip(roi);/***********************************************************************Define neighboor function* and generate the samples***********************************************************************/function neighFun(img,kernalRadius,roi){var kernel = ee.Kernel.square(kernalRadius,'pixels',false);var kernelArea = (ee.Number(kernalRadius).multiply(2).add(1)).pow(2);var imgNeibor = ee.Image(img).convolve(kernel).eq(kernelArea).set("system:footprint",roi.geometry());return img.updateMask(imgNeibor);}var samplePoint = ee.List([]);for(var i=0;i<=1;i++){var class_Num = i;var class_i_mask = neighFun(LULU_mutual2.eq(class_Num),1,roi);var class_i = LULU_mutual2.updateMask(class_i_mask);samplePoint = samplePoint.add(class_i);}var samplePoint = ee.ImageCollection(samplePoint).mosaic().rename("landcover").updateMask(cropland);var pointSample =  samplePoint.stratifiedSample({numPoints:pointNum,classBand:"landcover",region:roi.geometry(),scale:30,seed:0,// tileScale:8,geometries:true});return pointSample;
}/**************************************************************************
generate the grid
***************************************************************************/
function generateGrid(xmin, ymin, xmax, ymax, dx, dy) {var xx = ee.List.sequence(xmin, ee.Number(xmax).subtract(0.0001), dx);var yy = ee.List.sequence(ymin, ee.Number(ymax).subtract(0.0001), dy);var cells = xx.map(function(x) {return yy.map(function(y) {var x1 = ee.Number(x);var x2 = ee.Number(x).add(ee.Number(dx));var y1 = ee.Number(y);var y2 = ee.Number(y).add(ee.Number(dy));var coords = ee.List([x1, y1, x2, y2]);var rect = ee.Algorithms.GeometryConstructors.Rectangle(coords); return ee.Feature(rect);});}).flatten(); return ee.FeatureCollection(cells);
}function GridRegion(roiRegion,xBlock,yBlock){//roiRegion: area of interest in the form of geometry// compute the coordinatesvar bounds = roiRegion.bounds();var coords = ee.List(bounds.coordinates().get(0));var xmin = ee.List(coords.get(0)).get(0);var ymin = ee.List(coords.get(0)).get(1);var xmax = ee.List(coords.get(2)).get(0);var ymax = ee.List(coords.get(2)).get(1);var dx = (ee.Number(xmax).subtract(xmin)).divide(xBlock); //4var dy = (ee.Number(ymax).subtract(ymin)).divide(yBlock);var grid = generateGrid(xmin, ymin, xmax, ymax, dx, dy);  grid = grid.filterBounds(roiRegion); return grid;
}

植被水分含量指数(LSWI) 在其公式中使用了NIR和SWIR通道,SWIR波段对植被含水量的变化较敏感:
在这里插入图片描述
在这里插入图片描述

结果展示:
在这里插入图片描述

这篇关于GEE23:基于植被物候实现农作物分类的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/863181

相关文章

Springboot处理跨域的实现方式(附Demo)

《Springboot处理跨域的实现方式(附Demo)》:本文主要介绍Springboot处理跨域的实现方式(附Demo),具有很好的参考价值,希望对大家有所帮助,如有错误或未考虑完全的地方,望不... 目录Springboot处理跨域的方式1. 基本知识2. @CrossOrigin3. 全局跨域设置4.

Spring Boot 3.4.3 基于 Spring WebFlux 实现 SSE 功能(代码示例)

《SpringBoot3.4.3基于SpringWebFlux实现SSE功能(代码示例)》SpringBoot3.4.3结合SpringWebFlux实现SSE功能,为实时数据推送提供... 目录1. SSE 简介1.1 什么是 SSE?1.2 SSE 的优点1.3 适用场景2. Spring WebFlu

基于SpringBoot实现文件秒传功能

《基于SpringBoot实现文件秒传功能》在开发Web应用时,文件上传是一个常见需求,然而,当用户需要上传大文件或相同文件多次时,会造成带宽浪费和服务器存储冗余,此时可以使用文件秒传技术通过识别重复... 目录前言文件秒传原理代码实现1. 创建项目基础结构2. 创建上传存储代码3. 创建Result类4.

SpringBoot日志配置SLF4J和Logback的方法实现

《SpringBoot日志配置SLF4J和Logback的方法实现》日志记录是不可或缺的一部分,本文主要介绍了SpringBoot日志配置SLF4J和Logback的方法实现,文中通过示例代码介绍的非... 目录一、前言二、案例一:初识日志三、案例二:使用Lombok输出日志四、案例三:配置Logback一

Python如何使用__slots__实现节省内存和性能优化

《Python如何使用__slots__实现节省内存和性能优化》你有想过,一个小小的__slots__能让你的Python类内存消耗直接减半吗,没错,今天咱们要聊的就是这个让人眼前一亮的技巧,感兴趣的... 目录背景:内存吃得满满的类__slots__:你的内存管理小助手举个大概的例子:看看效果如何?1.

Python+PyQt5实现多屏幕协同播放功能

《Python+PyQt5实现多屏幕协同播放功能》在现代会议展示、数字广告、展览展示等场景中,多屏幕协同播放已成为刚需,下面我们就来看看如何利用Python和PyQt5开发一套功能强大的跨屏播控系统吧... 目录一、项目概述:突破传统播放限制二、核心技术解析2.1 多屏管理机制2.2 播放引擎设计2.3 专

Python实现无痛修改第三方库源码的方法详解

《Python实现无痛修改第三方库源码的方法详解》很多时候,我们下载的第三方库是不会有需求不满足的情况,但也有极少的情况,第三方库没有兼顾到需求,本文将介绍几个修改源码的操作,大家可以根据需求进行选择... 目录需求不符合模拟示例 1. 修改源文件2. 继承修改3. 猴子补丁4. 追踪局部变量需求不符合很

idea中创建新类时自动添加注释的实现

《idea中创建新类时自动添加注释的实现》在每次使用idea创建一个新类时,过了一段时间发现看不懂这个类是用来干嘛的,为了解决这个问题,我们可以设置在创建一个新类时自动添加注释,帮助我们理解这个类的用... 目录前言:详细操作:步骤一:点击上方的 文件(File),点击&nbmyHIgsp;设置(Setti

SpringBoot实现MD5加盐算法的示例代码

《SpringBoot实现MD5加盐算法的示例代码》加盐算法是一种用于增强密码安全性的技术,本文主要介绍了SpringBoot实现MD5加盐算法的示例代码,文中通过示例代码介绍的非常详细,对大家的学习... 目录一、什么是加盐算法二、如何实现加盐算法2.1 加盐算法代码实现2.2 注册页面中进行密码加盐2.

MySQL大表数据的分区与分库分表的实现

《MySQL大表数据的分区与分库分表的实现》数据库的分区和分库分表是两种常用的技术方案,本文主要介绍了MySQL大表数据的分区与分库分表的实现,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有... 目录1. mysql大表数据的分区1.1 什么是分区?1.2 分区的类型1.3 分区的优点1.4 分