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 bands
  var 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 bands
  var 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 area
  var 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 coordinates
  var 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); //4
  var dy = (ee.Number(ymax).subtract(ymin)).divide(yBlock);
  
  var grid = generateGrid(xmin, ymin, xmax, ymax, dx, dy);  
  grid = grid.filterBounds(roiRegion); 
  
  return grid;
}