问题背景
本次实践主要是对论文《Time-series land cover mapping and urban expansion analysis using OpenStreetMap data and remote sensing big data: A case study of Guangdong-Hong Kong-Macao Greater Bay Area, China》的代码复现（DOI: DOI），现在共享出来供大家学习。
整体流程如下图所示：
数据预处理
在GEE平台上编写代码，引入相关数据，选择武汉市作为研究区，数据预处理代码如下所示：
var Landsat5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
var Landsat7 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
var Landsat8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
// var table = table1;
var wuhan = table;
/** function definition **/
function applyScaleFactorsL89(image) {
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
return image.addBands(opticalBands, null, true)
.addBands(thermalBands, null, true);
}
function applyScaleFactorsL457(image) {
var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
var thermalBand = image.select('ST_B6').multiply(0.00341802).add(149.0);
return image.addBands(opticalBands, null, true)
.addBands(thermalBand, null, true);
}
function cloudmaskL89(image) {
// Bits 3 and 5 are cloud shadow and cloud, respectively.
var cloudShadowBitMask = (1 << 4);
var cloudsBitMask = (1 << 3);
// Get the pixel QA band.
var qa = image.select('QA_PIXEL');
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return image.updateMask(mask);
}
function cloudMaskL457(image) {
var qa = image.select('QA_PIXEL');
// 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 << 3)
.and(qa.bitwiseAnd(1 << 9))
.or(qa.bitwiseAnd(1 << 4));
// Remove edge pixels that don't occur in all bands
var mask2 = image.mask().reduce(ee.Reducer.min());
return image.updateMask(cloud.not()).updateMask(mask2);
}
/** bands rename function **/
function bandRenameL89(image) {
var blue = image.select(['SR_B2']).rename('blue');
var green = image.select(['SR_B3']).rename('green');
var red = image.select(['SR_B4']).rename('red');
var nir = image.select(['SR_B5']).rename('nir');
var swir1 = image.select(['SR_B6']).rename('swir1');
var swir2 = image.select(['SR_B7']).rename('swir2');
var new_image = blue.addBands([green, red, nir, swir1, swir2]);
return new_image;
}
function bandRenameL457(image) {
var blue = image.select(['SR_B1']).rename('blue');
var green = image.select(['SR_B2']).rename('green');
var red = image.select(['SR_B3']).rename('red');
var nir = image.select(['SR_B4']).rename('nir');
var swir1 = image.select(['SR_B5']).rename('swir1');
var swir2 = image.select(['SR_B7']).rename('swir2');
var new_image = blue.addBands([green, red, nir, swir1, swir2]);
return new_image;
}
/** caculate the NDVI, NDWI, NDBI **/
function NDVI(img) {
var nir = img.select("nir");
var red = img.select("red");
var ndvi = img.expression(
"(nir - red)/(nir + red)",
{
"nir": nir,
"red": red
}
).rename("NDVI");
return ndvi;
}
function NDWI(img) {
var green = img.select("green");
var nir = img.select("nir");
var ndwi = img.expression(
"(green - nir)/(green + nir)",
{
"green": green,
"nir": nir
}
).rename("NDWI");
return ndwi;
}
function NDBI(img){
var nir = img.select("nir");
var swir = img.select("swir1");
var ndbi = img.expression(
"(swir - nir)/(nir + swir)",
{
"nir": nir,
"swir": swir
}
).rename("NDBI");
return ndbi;
}
function MNDWI(img) {
var green = img.select("green");
var swir1 = img.select("swir1");
var mndwi = img.expression(
"(green - swir1)/(green + swir1)",
{
"green": green,
"swir1": swir1
}
);
return mndwi;
}
function glcm(img) {
var nir = img.select("nir").toUint16();
return nir.glcmTexture();
}
/** generate the bi-temporal image**/
function generateImage(start_date, end_date, cloudCover, wuhan){
// landsat 8
var l8_col = Landsat8.filterBounds(wuhan)
.filterDate(start_date, end_date)
.filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
.map(applyScaleFactorsL89)
.map(cloudmaskL89)
.map(bandRenameL89);
print('landsat8', l8_col.size())
// landsat 7
var l7_col = Landsat7.filterBounds(wuhan)
.filterDate(start_date, end_date)
.filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
.map(applyScaleFactorsL457)
.map(cloudMaskL457)
.map(bandRenameL457)
;
print('landsat7', l7_col.size())
// landsat 5
var l5_col = Landsat5.filterBounds(wuhan)
.filterDate(start_date, end_date)
.filter(ee.Filter.lt('CLOUD_COVER', cloudCover))
.map(applyScaleFactorsL457)
.map(cloudMaskL457)
.map(bandRenameL457)
;
print('landsat5', l5_col.size())
/** composite the image **/
var l5 = l5_col
var l7 = l7_col.mean()
var l8 = l8_col
print('l5', l5)
print('l7', l7)
print('l8', l8)
/** merge image **/
var image = l8.merge(l7).merge(l5);
print('image', image)
// print("final image count", image.size(), image)
/** start caculate **/
var final_image = image.mean().clip(wuhan);
var image_ndvi = NDVI(final_image)
var image_ndwi= NDWI(final_image)
var image_ndbi = NDBI(final_image)
var optim_image = final_image.addBands(image_ndvi.select("NDVI")).addBands(image_ndwi.select("NDWI"))
.addBands(image_ndbi.select("NDBI")).clip(wuhan)
// print("optim image", optim_image)
return optim_image
}
上述代码主要针对Landsat卫星数据进行处理。
主成分提取
var year = 2021;
for(var i = year; i<=year; i++){
var year_name = i;
var start_date = (year_name) + '-01-01';
var end_date = (year_name + 1) + '-01-01';
var cloudCover = 20
/** IMAGE HANDLE **/
var t1 = generateImage(start_date, end_date, cloudCover, table)
print("---------------------")
var start_date = (year_name + 1) + '-01-01';
var end_date = (year_name + 2) + '-01-01';
var t2 = generateImage(start_date, end_date, cloudCover, table)
/** PCA Analysis **/
// print("PCA Image is:", t2)
// var image = t1;
var region = table;
// var bandNames = image.bandNames();
var getNewBandNames = function(image, prefix)
{
var seq = ee.List.sequence(1, image.bandNames().length());
return seq.map(function(b) {return ee.String(prefix).cat(ee.Number(b).int());});
};
var getPrincipalComponents = function(image, table)
{
var arrays = image.toArray();
var covar = arrays.reduceRegion({
reducer: ee.Reducer.centeredCovariance(),
geometry:table,
scale: 30,
maxPixels: 1e13});
var covarArray = ee.Array(covar.get('array'));
var eigens = covarArray.eigen();
var eigenValues = eigens.slice(1, 0, 1);
var eigenVectors = eigens.slice(1, 1);
var arrayImage = arrays.toArray(1);
var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
var sdImage = ee.Image(eigenValues.sqrt()).arrayProject([0]).arrayFlatten([getNewBandNames(image, 'sd')]);
return principalComponents.arrayProject([0]).arrayFlatten([getNewBandNames(image, 'pc')]).divide(sdImage);
};
//pcImage中各波段名叫pc1、pc2
var pcImage1 = getPrincipalComponents(t1, table);
var pcImage2 = getPrincipalComponents(t2, table);
//取主成分分析的第一主成分pc1
// var pc1 = pcImage.select([pcImage.bandNames().get(0).getInfo()]);
// pc1 = pc1.float();
// print("pcImage:", pcImage)
var imgA = pcImage1.select(["pc1", "pc2", "pc3"])
var imgB = pcImage2.select(["pc1", "pc2", "pc3"])
// print("imgA:", imgA)
// 区域就是设置的roi，分辨率30米
// Export.image.toAsset({
// image: imgA,
// description: 'pca'+i,
// assetId: 'Article/pca'+i,
// scale: 30,
// maxPixels:1e13
// });
}
变化检测挑选变化点和不变点
/** CVA **/
print("start CVA analysis")
var roi = table;
// CVA没有亲和变换的一致性，因此需要对影像进行均值标准化，以消除辐射不一致性
// 1.计算均值、标准差
var meanA = imgA.reduceRegion({
reducer:ee.Reducer.mean(),
geometry:roi,
scale:30,
maxPixels:1e13});
var meanB = imgB.reduceRegion({
reducer:ee.Reducer.mean(),
geometry:roi,
scale:30,
maxPixels:1e13});
var stdA = imgA.reduceRegion({
reducer:ee.Reducer.stdDev(),
geometry:roi,
scale:30,
maxPixels:1e13});
var stdB = imgB.reduceRegion({
reducer:ee.Reducer.stdDev(),
geometry:roi,
scale:30,
maxPixels:1e13});
// 2.归一化
function get_mutibands_image(band_1, band_2, band_3) {
var image1 = ee.Image.constant(ee.Number(band_1)).rename(band1_name);
var image2 = ee.Image.constant(ee.Number(band_2)).rename(band2_name);
var image3 = ee.Image.constant(ee.Number(band_3)).rename(band3_name);
return image1.addBands(image2).addBands(image3);
}
// 设定波段组合，Sentinel-2数据中真彩色:432,假彩色:843
var band_list = ee.List(["pc1", "pc2", "pc3"]);
var band1_name = ee.String(band_list.get(0));
var band2_name = ee.String(band_list.get(1));
var band3_name = ee.String(band_list.get(2));
// print(meanA.get(band1_name))
var imgA_mean = get_mutibands_image(meanA.get(band1_name), meanA.get(band2_name), meanA.get(band3_name));
var imgB_mean = get_mutibands_image(meanB.get(band1_name), meanB.get(band2_name), meanB.get(band3_name));
var imgA_std = get_mutibands_image(stdA.get(band1_name), stdA.get(band2_name), stdA.get(band3_name));
var imgB_std = get_mutibands_image(stdB.get(band1_name), stdB.get(band2_name), stdB.get(band3_name));
var normal_imgA = imgA.subtract(imgA_mean).divide(imgA_std);
var normal_imgB = imgB.subtract(imgB_mean).divide(imgB_std);
// Map.addLayer(normal_imgA, {min: 0, max: 1}, "normal_imgA");
// Map.addLayer(normal_imgB, {min: 0, max: 1}, "normal_imgB");
// 3.CVA
function cva(normal_imgA, normal_imgB) {
var img_diff = normal_imgA.subtract(normal_imgB);
var square_img_diff = img_diff.pow(2);
var sum_image = square_img_diff.expression("B1+B2+B3", {
B1: square_img_diff.select(band1_name),
B2: square_img_diff.select(band2_name),
B3: square_img_diff.select(band3_name),
});
var l2_norm = sum_image.rename("sum_image").sqrt();
return l2_norm;
}
var l2_norm = cva(normal_imgA, normal_imgB);
// 4.OTSU
function otsu(histogram) {
var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
var size = means.length().get([0]);
var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
var mean = sum.divide(total);
var indices = ee.List.sequence(1, size);
var bss = indices.map(function(i) {
var aCounts = counts.slice(0, 0, i);
var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
var aMeans = means.slice(0, 0, i);
var aMean = aMeans.multiply(aCounts)
.reduce(ee.Reducer.sum(), [0]).get([0])
.divide(aCount);
var bCount = total.subtract(aCount);
var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
bCount.multiply(bMean.subtract(mean).pow(2)));
});
return means.sort(bss).get([-1]);
}
var histogram = l2_norm.reduceRegion({
reducer: ee.Reducer.histogram(),
geometry: roi,
scale: 30,
maxPixels: 1e13,
});
var threshold = otsu(histogram.get("sum_image"));
print(threshold,'分割阈值');
// 结果展示
var mask = l2_norm.gte(threshold).clip(roi);
var change_map = mask.rename("change_map");
// print("mask", mask)
Map.addLayer(change_map.clip(roi),{min:0,max:1,palette:['#000000','#ffffff']},'change_result');
// Export the change_map
Export.image.toAsset({
image: change_map.clip(roi),
description: 'change_map',
assetId: 'Article/change_map',
scale: 30,
maxPixels:1e13
});
var change_map = change_map.clip(roi)
print("extract value to point")
根据变化检测结果选点
//1.Copland
var cropland = change_map.sampleRegions({
collection: Cropland,
properties: ["fclass"],
scale: 30,
geometries:true,
tileScale:16
});
print("Cropland: ", Cropland)
var unchanged_cropland = cropland.filter(ee.Filter.eq('change_map', 0));
//2.Forest
var forest = change_map.sampleRegions({
collection: Forest,
properties: ["fclass"],
scale: 30,
geometries:true,
tileScale:16
});
var unchanged_forest = forest.filter(ee.Filter.eq('change_map', 0));
//3.Isurface
var isurface = change_map.sampleRegions({
collection: Isurface,
properties: ["fclass"],
scale: 30,
geometries:true,
tileScale:16
});
var unchanged_isurface = isurface.filter(ee.Filter.eq('change_map', 0));
//4.Water
var water = change_map.sampleRegions({
collection: Water,
properties: ["fclass"],
scale: 30,
geometries:true,
tileScale:16
});
var unchanged_water = water.filter(ee.Filter.eq('change_map', 0));
print(unchanged_water.size())
训练随机森林并测试
//==================================================================================
// now, I get the unchanged points, then I can extract the landsat and other images
// 1.SRTM
var dataset = ee.Image('USGS/SRTMGL1_003').clip(roi);
var elevation = dataset.select('elevation');
var slope = ee.Terrain.slope(elevation);
var e = elevation.select("elevation").clip(roi)
var s = slope.select("slope").clip(roi)
//2.1 SAR_1, pre image
var sentinel_1 = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(wuhan).filterDate(year+'-01-01', year+'-12-20').median().clip(roi);
//2.2 SAR_2, post image
var sentinel_2 = ee.ImageCollection('COPERNICUS/S1_GRD').filterBounds(wuhan).filterDate((year+1)+'-01-01', (year+1)+'-12-20').median().clip(roi);
//3.GLCM
function glcm(img) {
var nir = img.select("nir").toUint16();
var gt = nir.glcmTexture();
var savg = gt.select('nir_savg').rename('savg')
var gvar = gt.select('nir_var').rename('var')
var idm = gt.select('nir_idm').rename('idm')
var con = gt.select('nir_contrast').rename('con')
var diss = gt.select('nir_diss').rename('diss')
var ent = gt.select('nir_ent').rename('ent')
var asm = gt.select('nir_asm').rename('asm')
var corr = gt.select('nir_corr').rename('corr')
return savg.addBands(gvar).addBands(idm)
.addBands(con).addBands(diss).addBands(ent).addBands(asm).addBands(corr)
}
var glcm_1 = glcm(t1)
var glcm_2 = glcm(t2)
print("glcm: ",glcm_1)
//4.concat the image
var t1 = t1.addBands(e).addBands(s).addBands(sentinel_1).addBands(glcm_1)
var t2 = t2.addBands(e).addBands(s).addBands(sentinel_2).addBands(glcm_2)
var t1 = t1.clip(roi)
var t2 = t2.clip(roi)
print("t1", t1)
print("t2", t2)
//----------------------------------------------------------
//----------------------------------------------------------
// start classifier, first do 2022, then pre years
var past_training = unchanged_cropland.merge(unchanged_forest).merge(unchanged_isurface).merge(unchanged_water);
print("2021 training:", past_training.limit(10))
var training = Cropland.merge(Forest).merge(Isurface).merge(Water);
print("2022 training: ",training.limit(10));
var training = past_training;
print("training: ",training.size());
var training = training.remap(
["cropland", "forest", "Impervious Surface", "water"], // 原始类别标签（0表示无效像素）
[0, 1, 2, 3], // 映射后的数字（0留给无效像素）
"fclass" // 字段
);
var trainingData = training.randomColumn('random')
var sample_training = trainingData.filter(ee.Filter.lte("random", 0.8));
var sample_validate = trainingData.filter(ee.Filter.gt("random", 0.8));
print("sample_training size is: ", sample_training.size())
print("sample_validate size is: ", sample_validate.size())
// 利用样本点拾取特征值用于模型训练和验证
var training = t1.clip(roi).sampleRegions({
collection: sample_training,
properties: ["fclass"],
scale: 30,
tileScale:16
});
var validation = t1.clip(roi).sampleRegions({
collection: sample_validate,
properties: ["fclass"],
scale: 30,
tileScale:16
});
print("validation: ", validation)
var classifier = ee.Classifier.smileRandomForest(50)
.train({
features: training,
classProperty: 'fclass',
inputProperties: t1.bandNames()
});
var dict = classifier.explain();
var variable_importance = ee.Feature(null, ee.Dictionary(dict).get('importance'));
// 混淆矩阵法
var validated = validation.classify(classifier);
// 混淆矩阵
var testAccuracy = validated.errorMatrix('fclass', 'classification');
// 总体分类精度
var accuracy = testAccuracy.accuracy();
// 用户分类精度
var userAccuracy = testAccuracy.consumersAccuracy();
// 生产者精度
var producersAccuracy = testAccuracy.producersAccuracy();
// Kappa系数
var kappa = testAccuracy.kappa();
print('混淆矩阵:', testAccuracy);//
print('用户分类精度:', userAccuracy);//用户分类精度
print('生产者精度:', producersAccuracy);//生产者精度
print('总体分类精度:', accuracy);//总体分类精度
print('Kappa:', kappa);
现在已经可以进行土地覆盖制图，一些数据预处理的流程尚未提及，例如拓扑检查等操作。应该会在空闲时间写一下教程。
301
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
