C32 农业环境（Part 7 与人类活动相关的应用）
本期导读
介绍如何使用GEE中的数据集和功能来绘制大规模的农业地图。示例：绘制美国中西部作物类型，美国中西部是世界上主要的产粮区之一。也可绘制其他农业特征的能力，如产量和管理实践。
主要内容
- 使用耕地数据层(CDL)作为标签对美国中西部某县的作物类型进行分类。
- 利用ee.Reducer.linearRegression对作物时间序列拟合二阶调和回归，提取调和系数。
- 利用叶绿素植被指数(GCVI)进行作物类型分类。
- 训练随机森林根据调和系数对作物类型进行分类。
- 将训练好的随机森林分类器应用于研究区域并评估其性能。
section1 ：调出研究区所有的Landsat数据
section2 ：给Landsat影像 Add Bands 为谐波分析（Harmonic Regression）做准备
谐波分析（Harmonic Regression）是常用的云缺失填补方法， 它是一种频域时序列分析，从原始时间序列数据的正弦和余弦函数重构周期/季节性波。对遥感特征指数使用谐波分析，有利于反映特征的周期性。
section3 ：对每一个像素拟合谐波回归
section4 ：训练和评估一个随机森林分类器
//========================section 1=====================================================
// 定义研究区域。
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');
// 导入陆地卫星图像。
var landsat7 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2');
var landsat8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2');
//定义了一些函数，使处理陆地卫星数据变得更容易
// 波段重命名（这里是对所有波段，也可先筛选波段再重命名）
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'
]);
}
// 云掩膜、阴影和其他不需要的特征的函数。
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'));
}
// 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]);
}
// 定义时间。
var start_date = '2020-01-01';
var end_date = '2020-12-31';
// 在研究区域的开始日期和结束日期之间调取陆地卫星7号和8号的图像。
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和8收集。
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);
// 可视化一个点位的 时间序列GCVI
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);
// Get crop type dataset.
var cdl = ee.Image('USDA/NASS/CDL/2020').select(['cropland']);
Map.addLayer(cdl.clip(geometry), {}, 'CDL 2020');
//========================section 2 =====================================================
// 定义函数：为图像添加时间波段.
function addTimeUnit(image, refdate) {
var date = image.date();
var dyear = date.difference(refdate, 'year');
var t = image.select(0).multiply(0).add(dyear).select([0], ['t']) .float();
var imageplus = image.addBands(t);
return imageplus;
}
// 定义函数addHarmonics： 为图像添加谐波分析参数omega
function addHarmonics(image, omega, refdate) {
image = addTimeUnit(image, refdate);
var timeRadians = image.select('t').multiply(2 * Math.PI * omega);
var timeRadians2 = image.select('t').multiply(4 * Math.PI * omega);
return image
.addBands(timeRadians.cos().rename('cos'))
.addBands(timeRadians.sin().rename('sin'))
.addBands(timeRadians2.cos().rename('cos2'))
.addBands(timeRadians2.sin().rename('sin2'))
.addBands(timeRadians.divide(timeRadians)
.rename('constant'));
}
// 对影像应用addHarmonics函数
var omega = 1;
var landsatPlus = landsat.map(
function(image) {
return addHarmonics(image, omega, start_date);
});
print('Landsat collection with harmonic basis: ', landsatPlus);
//========================section 3 对每一个像素拟合谐波回归=====================================================
//In the previous section, we added the harmonic basis (cosine, sine, and constant
//terms) to each image in our Landsat collection.
//前面部分，我们给每一幅影像添加了谐波回归参数（sin,cos，和常数项）
//现在，我们可以用这个基作为自变量进行线性回归，Landsat波段和GCVI作为因变量。
//继续定义几个函数，我们完成这个回归。
//定义函数arrayimgHarmonicRegr,对图像运行线性回归
function arrayimgHarmonicRegr(harmonicColl, dependent, independents) {
independents = ee.List(independents);
dependent = ee.String(dependent);
var regression = harmonicColl .select(independents.add(dependent))
.reduce(ee.Reducer.linearRegression(independents.length(), 1));
return regression;
}
//定义函数imageHarmonicRegr,提取和重命名回归系数.
function imageHarmonicRegr(harmonicColl, dependent, independents) {
var hregr = arrayimgHarmonicRegr(harmonicColl, dependent, independents);
independents = ee.List(independents);
dependent = ee.String(dependent);
var newNames = independents.map(function(b) {
return dependent.cat(ee.String('_')).cat(ee.String( b));
});
var imgCoeffs = hregr.select('coefficients')
.arrayProject([0])
.arrayFlatten([independents])
.select(independents, newNames);
return imgCoeffs;
}
//定义函数getHarmonicCoeffs,应用imageHarmonicRegr并创建一个多波段图像。
function getHarmonicCoeffs(harmonicColl, bands, independents) {
var coefficients = ee.ImageCollection.fromImages(bands.map(
function(band) {
return imageHarmonicRegr(harmonicColl, band, independents);
}));
return coefficients.toBands();
}
//调用函数getHarmonicCoeffs.
var bands = ['NIR', 'SWIR1', 'SWIR2', 'GCVI'];
var independents = ee.List(['constant', 'cos', 'sin', 'cos2', 'sin2']);
var harmonics = getHarmonicCoeffs(landsatPlus, bands, independents);
harmonics = harmonics.clip(geometry);
harmonics = harmonics.multiply(10000).toInt32();
//计算拟合值。
var gcviHarmonicCoefficients = harmonics
.select([
'3_GCVI_constant', '3_GCVI_cos',
'3_GCVI_sin', '3_GCVI_cos2', '3_GCVI_sin2'
])
.divide(10000);
var fittedHarmonic = landsatPlus.map(function(image) {
return image.addBands(
image.select(independents)
.multiply(gcviHarmonicCoefficients)
.reduce('sum')
.rename('fitted')
);
});
//在图表中可视化拟合的谐波。
var harmonicsChart = ui.Chart.image.series( fittedHarmonic.select(
['fitted', 'GCVI']), point, ee.Reducer.mean(), 30)
.setSeriesNames(['GCVI', 'Fitted'])
.setOptions({
title: 'Landsat GCVI time series and fitted harmonic regression values',
lineWidth: 1,
pointSize: 3,
});
print(harmonicsChart);
//将CDL作为band添加到谐波中。
var harmonicsPlus = ee.Image.cat([harmonics, cdl]);
// Export image to asset.
var filename = 'McLean_County_harmonics';
Export.image.toAsset({
image: harmonicsPlus,
description: filename,
assetId: 'your_asset_path_here/' + filename,
dimensions: null,
region: region,
scale: 30,
maxPixels: 1e12
});
//在地图上显示谐波系数。
var visImage = ee.Image.cat([
harmonics.select('3_GCVI_cos').divide(7000).add(0.6),
harmonics.select('3_GCVI_sin').divide(7000).add(0.5),
harmonics.select('3_GCVI_constant').divide(7000).subtract( 0.6)
]);
Map.addLayer(visImage, {
min: -0.5,
max: 0.5
}, 'Harmonic coefficient false color');
//========================section 4 训练和评估一个随机森林分类器=====================================================
// Define a random forest classifier.
var rf = ee.Classifier.smileRandomForest({
numberOfTrees: 50,
minLeafPopulation: 10,
seed: 0
});
// Get harmonic coefficient band names.
var bands = harmonicsPlus.bandNames();
bands = bands.remove('cropland').remove('system:index');
// Transform CDL into a 3-class band and add to harmonics.
var cornSoyOther = harmonicsPlus.select('cropland').eq(1)
.add(harmonicsPlus.select('cropland').eq(5).multiply(2));
var dataset = ee.Image.cat([harmonicsPlus.select(bands), cornSoyOther]);
// Sample training points.
var train_points = dataset.sample(geometry, 30, null, null, 100, 0);
print('Training points', train_points);
// Train the model!
var model = rf.train(train_points, 'cropland', bands);
var trainCM = model.confusionMatrix();
print('Training error matrix: ', trainCM);
print('Training overall accuracy: ', trainCM.accuracy());
// Sample test points and apply the model to them.
var test_points = dataset.sample(geometry, 30, null, null, 50, 1);
var tested = test_points.classify(model);
// Compute the confusion matrix and accuracy on the test set.
var testAccuracy = tested.errorMatrix('cropland', 'classification');
print('Test error matrix: ', testAccuracy);
print('Test overall accuracy: ', testAccuracy.accuracy());
// Apply the model to the entire study region.
var regionClassified = harmonicsPlus.select(bands).classify(model);
var predPalette = ['gray', 'yellow', 'green'];
Map.addLayer(regionClassified, {
min: 0,
max: 2,
palette: predPalette
}, 'Classifier prediction');
// Visualize agreements/disagreements between prediction and CDL.
Map.addLayer(regionClassified.eq(cornSoyOther), {
min: 0,
max: 1
}, 'Classifier agreement');
