《gee大神之路》月下醉吟

14、十四 在

　　在本实习中，我们将使用 Earth Engine 提取 Landsat 图像并对美国中西部的农作物类型进行分类。美国是世界上最大的玉米生产国和第二大大豆生产国；因此，美国保持高产量对于全球粮食安全和价格稳定至关重要。作物类型测绘是利用卫星图像估算农业产量的重要先决条件。我们将展示如何通过拟合调和回归并提取系数来从时间序列中获取特征。然后我们将使用随机森林对作物类型进行分类，其中“地面实况”标签将来自美国农业部的农田数据层（CDL）。CDL 本身就是使用卫星图像作为输入和基于调查的地面实况作为训练标签的分类器的产物（USDA 2021）。

　　在中西部，玉米和大豆占主导地位。我们将绘制伊利诺伊州麦克莱恩县的农作物类型地图。该县是美国玉米和大豆产量最多的县，分别为62.9和1,930万蒲式耳。让我们首先定义研究区域并将其可视化（图 A1.1.1）。我们从美国人口普查局的 TIGER 数据集中获取县边界，该数据集可以在 Earth Engine 数据目录中找到。我们通过名称并使用伊利诺伊州的FIPS 从 TIGER 中提取麦克莱恩县特征。

　　// Define study area.
　　var TIGER=ee.FeatureCollection(\'TIGER/2018/Counties\');
　　var region=ee.Feature(TIGER
　　.filter(ee.Filter.eq(\'STATEFP\', \'17\'))
　　.filter(ee.Filter.eq(\'NAME\', \'McLean\'))
　　.first());
　　var geometry=region.geometry();
　　Map.centerObject(region);
　　Map.addLayer(region,{
　　\'color\': \'red\'
　　}, \'McLean County\');

　　接下来，我们导入 Landsat 7和8图像，特别是，我们使用2级数据，其中包含经过大气校正的表面反射率和陆地表面温度。虽然大气顶层数据也将为许多农业应用带来良好的结果，但我们更喜欢在可用时针对大气条件进行校正的数据，因为大气变化只是输入中的又一个噪声源。

　　// Import Landsat imagery.
　　var landsat7=ee.ImageCollection(\'LANDSAT/LE07/C02/T1_L2\');
　　var landsat8=ee.ImageCollection(\'LANDSAT/LC08/C02/T1_L2\');

　　我们定义了一些函数，使处理 Landsat 数据变得更加容易。前两个函数renameL7 和renameL8将 Landsat 7或Landsat 8图像的波段重命名为更直观的名称。例如，我们没有将 Landsat 7图像的第一个波段称为“B1”，而是将其重命名为“BLUE”，因为波段1捕获可见光谱蓝色部分的光。重命名还可以更轻松地同时处理 Landsat 7和8图像，因为它们的传感器（ETM+ 和 OLI）的波段编号不同。

　　// Functions to rename Landsat 7 and 8 images.
　　function renameL7(img){
　　return img.rename([\'BLUE\', \'GREEN\', \'RED\', \'NIR\', \'SWIR1\',
　　\'SWIR2\', \'TEMP1\', \'ATMOS_OPACITY\', \'QA_CLOUD\',
　　\'ATRAN\', \'CDIST\',
　　\'DRAD\', \'EMIS\', \'EMSD\', \'QA\', \'TRAD\', \'URAD\',
　　\'QA_PIXEL\',
　　\'QA_RADSAT\'
　　]);
　　}

　　function renameL8(img){
　　return img.rename([\'AEROS\', \'BLUE\', \'GREEN\', \'RED\', \'NIR\',
　　\'SWIR1\',
　　\'SWIR2\', \'TEMP1\', \'QA_AEROSOL\', \'ATRAN\', \'CDIST\',
　　\'DRAD\', \'EMIS\',
　　\'EMSD\', \'QA\', \'TRAD\', \'URAD\', \'QA_PIXEL\', \'QA_RADSAT\'
　　]);

　　}

　　Landsat 7 Level-2产品有七个表面反射率和温度带，而 Landsat 8 Level-2产品有八个。两者都有许多其他波段用于图像质量、大气条件等。在本章中，我们将主要关注近红外 (NIR)、短波红外（SWIR1 和 SWIR2）和像素质量 (QA_PIXEL)乐队。在电磁波谱的 NIR、SWIR1 和 SWIR2 波段中可以看到植被类型之间的许多差异。在使用光学图像时，像素质量带对于遮蔽云层非常重要。下面，我们定义两个函数：addMask 函数将QA_PIXEL 位掩码转换为多个遮罩层，以及maskQAClear 函数可从每个图像中删除所有不清晰的像素。

　　// Functions to mask out clouds, shadows, and other unwanted features.
　　function addMask(img){
　　// Bit 0: Fill
　　// Bit 1: Dilated Cloud
　　// Bit 2: Cirrus (high confidence) (L8) or unused (L7)
　　// Bit 3: Cloud
　　// Bit 4: Cloud Shadow
　　// Bit 5: Snow
　　// Bit 6: Clear
　　// 0: Cloud or Dilated Cloud bits are set
　　// 1: Cloud and Dilated Cloud bits are not set
　　// Bit 7: Water
　　var clear=img.select(\'QA_PIXEL\').bitwiseAnd(64).neq(0);
　　clear=clear.updateMask(clear).rename([\'pxqa_clear\']);

　　var water=img.select(\'QA_PIXEL\').bitwiseAnd(128).neq(0);
　　water=water.updateMask(water).rename([\'pxqa_water\']);

　　var cloud_shadow=img.select(\'QA_PIXEL\').bitwiseAnd(16).neq(0);
　　cloud_shadow=cloud_shadow.updateMask(cloud_shadow).rename([
　　\'pxqa_cloudshadow\'
　　]);

　　var snow=img.select(\'QA_PIXEL\').bitwiseAnd(32).neq(0);
　　snow=snow.updateMask(snow).rename([\'pxqa_snow\']);

　　var masks=ee.Image.cat([
　　clear, water, cloud_shadow, snow
　　]);

　　return img.addBands(masks);
　　}

　　function maskQAClear(img){
　　return img.updateMask(img.select(\'pxqa_clear\'));
　　}

　　除了 ETM+ 和 OLI 感知的原始波段外，植被指数 (VI) 也可以帮助区分不同的植被类型。之前的工作发现绿色叶绿素植被指数 (GCVI) 对于区分中西部的玉米和大豆特别有用（Wang 等人，2019）。我们会将其作为一个波段添加到每个 Landsat图像中。

　　/ Function to add GCVI as a band.
　　function addVIs(img){
　　var gcvi=img.expression(\'(nir / green) - 1\',{
　　nir: img.select(\'NIR\'),
　　green: img.select(\'GREEN\')
　　}).select([0], [\'GCVI\']);

　　return ee.Image.cat([img, gcvi]);
　　}

　　现在我们准备为我们的研究区域提取 Landsat 7和8图像并应用上述函数。我们将访问2020年1月1日至12月31日期间与麦克莱恩县相交的所有图像。由于不同的作物有不同的生长特征，因此获得时间序列图像来区分作物何时生长、衰老和收获是有价值的。各个时间点的光谱反射率差异也可以提供信息。

　　// Define study time period.
　　var start_date=\'2020-01-01\';
　　var end_date=\'2020-12-31\';

　　// Pull Landsat 7 and 8 imagery over the study area between start and end dates.
　　var landsat7coll=landsat7
　　.filterBounds(geometry)
　　.filterDate(start_date, end_date)
　　.map(renameL7);

　　var landsat8coll=landsat8
　　.filterDate(start_date, end_date)
　　.filterBounds(geometry)
　　.map(renameL8);
　　接下来，我们合并 Landsat 7 和 Landsat 8 集合，遮盖云层，并将 GCVI 添加为波段。
　　// Merge Landsat 7 and 8 collections.
　　var landsat=landsat7coll.merge(landsat8coll)
　　.sort(\'system:time_start\');

　　// Mask out non-clear pixels, add VIs and time variables.
　　landsat=landsat.map(addMask)
　　.map(maskQAClear)
　　.map(addVIs);

　　我们可以使用ui根据 CDL 的说法，图表对象可可视化特定点处的综合 Landsat GCVI 时间序列，在本例中该点位于玉米田中。

　　// Visualize GCVI time series at one location.
　　var point=ee.Geometry.Point([-88.81417685576481,
　　40.579804398254005
　　]);
　　var landsatChart=ui.Chart.image.series(landsat.select(\'GCVI\'),
　　point)
　　.setChartType(\'ScatterChart\')
　　.setOptions({
　　title: \'Landsat GCVI time series\',
　　lineWidth: 1,
　　pointSize: 3,
　　});
　　print(landsatChart);

　　在伊利诺伊州麦克莱恩县的一个玉米点可视化的陆地卫星 GCVI 时间序列。GCVI 越高，叶片叶绿素含量越高；农作物在六月开始变绿，在九月开始衰老。

　　最后，我们还看一下将用于训练和评估分类器的作物类型数据集。我们加载2020年的 CDL 图像并选择包含作物类型信息的图层。

　　// Get crop type dataset.
　　var cdl=ee.Image(\'USDA/NASS/CDL/2020\').select([\'cropland\']);
　　Map.addLayer(cdl.clip(geometry),{}, \'CDL 2020\');

　　在 CDL 中，玉米田以黄色显示，大豆田以绿色显示。
　　伊利诺伊州麦克莱恩县上空的 CDL 可视化。玉米田显示为黄色，而大豆田显示为深绿色。其他土地覆盖类型包括灰色的城市地区、青色的落叶林和蓝色的水域。