引言
GEE,google earth engine,是一个地理数据的云计算平台,可以快速、批量处理数据。
我想统计一下自己的家乡,四川省乐至县,这十几年的各乡镇土地覆盖变化情况。考虑到本地计算需要进行繁琐的数据处理,我这里使用GEE进行计算。
研究思路
我将研究过程分为四个大的步骤,前三个步骤在GEE中进行,最后一个步骤是在本地对excel数据进行处理。
数据筛选与预处理
这一个步骤主要是对影像数据进行土地覆盖数据加载、研究区加载以及数据标准化。
逐年影像加载
我在GEE中找到了几个土地覆盖数据:ESA数据与esri数据,10米分辨率,但只有2020年。modis数据,有几十年的数据,但分辨率为500米,太低。
这里我采用的数据是刘小平老师制作的全球2000年—2015年30m分辨率逐年土地覆盖数据。
确定好数据之后,按照年份顺序加载到GEE中:
//加载土地覆盖数据 该数据是2000-2015年的土地覆盖数据,参考文献:刘小平等.全球2000年—2015年30m分辨率逐年土地覆盖制图
var imageCollection = ee.ImageCollection("users/xxc/GLC_2000_2015");
//加载逐年土地覆盖影像
for(var year=0; year<16; year++){
var year_LEZHI=ee.Number(2000).add(year)
//根据名称波段选择土地覆盖影像
var year_band=ee.String("land_cover_").cat(ee.Number(year_LEZHI))
var Lezhi_image = imageCollection.select(year_band)
//配色方案设定,这里用的开源ee.palette,github地址为:https://github.com/gee-community/ee-palettes
var palettes = require('users/gena/packages:palettes');
var palette = palettes.colorbrewer.Dark2[3,4,5,6,7,8];
//加载逐年的土地覆盖影像
if (year>=10){Map.addLayer(Lezhi_image,{min: 10, max: 100, palette: palette},'20' +year +"年土地覆盖影像");}
if (year<10){Map.addLayer(Lezhi_image,{min: 10, max: 100, palette: palette},'200' +year +"年土地覆盖影像");}
}
使用上述方法,我们可以看到乐至县的乡镇矢量数据和逐年土地覆盖影像:
//导出土地覆盖影像数据
Export.image.toDrive({
image: Lezhi_image,
scale: 30,
region:LeZhi,
description:'Lezhi_landcover_2000-2015'
});
同时我们也可以导出乐至县的土地覆盖数据,这个数据下载后一共有16个波段,每一个波段是一个年份,单个波段的像素值是一种地物:
转为collection
在影像预处理步骤,我们需要将乐至县包含15年土地覆盖波段的影像转为影像合集,即由image转为collection,原先的band转为image。
//将image转为collection
var collection=ee.ImageCollection.fromImages([
Lezhi_image.select("land_cover_2000"),
Lezhi_image.select("land_cover_2001"),
Lezhi_image.select("land_cover_2002"),
Lezhi_image.select("land_cover_2003"),
Lezhi_image.select("land_cover_2004"),
Lezhi_image.select("land_cover_2005"),
Lezhi_image.select("land_cover_2006"),
Lezhi_image.select("land_cover_2007"),
Lezhi_image.select("land_cover_2008"),
Lezhi_image.select("land_cover_2009"),
Lezhi_image.select("land_cover_2010"),
Lezhi_image.select("land_cover_2011"),
Lezhi_image.select("land_cover_2012"),
Lezhi_image.select("land_cover_2013"),
Lezhi_image.select("land_cover_2014"),
Lezhi_image.select("land_cover_2015"),
])
单个要素单年份统计
在进行多年多要素统计之前,我们需要写几个函数实现单个要素某一年的土地覆盖统计功能。
Mask掉其他类别
为了待会能够使用reduce工具统计土地覆盖变化情况,我们需要将每一年的影像按照10个地物类别分割为10景影像,并在每一景影像中只保留一种地物。最后将这十景影像转为collection,并return,以便统计。
//将一年影像按照类别拆分为10个单波段的image
var splitdataset = function(img){
//只有值相同的像素点才保留,并将添加一个波段的属性val,val的值为类别
//10耕地、20林地、30草地、40灌木、50湿地、60水体、70苔原、80不透水面、90裸地、100永久冰雪
var class1 = img.eq(10).selfMask().set('val', 1);
var class2 = img.eq(20).selfMask().set('val', 2);
var class3 = img.eq(30).selfMask().set('val', 3);
var class4 = img.eq(40).selfMask().set('val', 4);
var class5 = img.eq(50).selfMask().set('val', 5);
var class6 = img.eq(60).selfMask().set('val', 6);
var class7 = img.eq(70).selfMask().set('val', 7);
var class8 = img.eq(80).selfMask().set('val', 8);
var class9 = img.eq(90).selfMask().set('val', 9);
var class10 =img.eq(100).selfMask().set('val', 10);
//返回一个collection
return ee.ImageCollection.fromImages([
class1, class2, class3, class4, class5, class6, class7, class8, class9, class10]);
};
reduce统计地物类别
在每一个要素中,选择一个年份的土地覆盖影像dataset,将dataset进行mask,获得的结果是一个collection,又将collection通过toBands函数转为image,并使用reduce函数进行image的像素统计。
//是针对每个要素进行多年统计的函数
var single_feature_static = function(feature) {
//addBands是针对每个要素进行一年统计的函数
var single_feature_single_year = function(year, feat){
//地物类别数量
var totalClass = 10;
//选择年份
year = ee.Number(year).toInt()
//拼接波段名称 例如2000年就是:land_cover_2000
var actual_year = ee.Number(2000).add(year)
//选择指定年份的波段
dataset = ee.Image(collection.toList(16).get(ee.Number(year)))
//调用上面的splitdataset函数,并返回带有10种地物image的collection
nd = splitdataset(dataset);
//将待统计的矢量要素转为feature
feat = ee.Feature(feat);
//循环地物类别,统计每种地物面积
for(var a=0; a<totalClass; a++){
var showList = nd.filterMetadata('val', 'equals', a+1).map(function(img){
return img.clip(feat);
})
var combineVal = showList.toBands().reduceRegion({
reducer: ee.Reducer.sum(),
geometry: feature.geometry(),
scale: 30
});
}
通过上述步骤,获得的每一个地物在每一年的像素个数字典combineVal,使用getNumber函数的作用是获取字典中该属性的像元个数。则统计的面积=像素个数*单个像元面积900米²。
//特定的波段的名称拼接
var cls = ee.String(ee.Number(a)).cat(ee.String("_land_cover_").cat(ee.Number(actual_year)));
//统计特定地物的面积,getNumber函数的作用是获取字典中该属性的value
var area = ee.Number(combineVal.getNumber(cls).multiply(900));
//待会导出表的名称
var name = ee.String(ee.Number(actual_year)).cat(ee.String("_Class_")).cat(ee.Number(a+1))
//给每种地物的面积赋值:10耕地、20林地、30草地、40灌木、50湿地、60水体、70苔原、80不透水面、90裸地、100永久冰雪
if ((a+1)===1){ var name_1 = name;var area_1 = area;}
if ((a+1)===2){ var name_2 = name;var area_2 = area;}
if ((a+1)===3){ var name_3 = name;var area_3 = area;}
if ((a+1)===4){ var name_4 = name;var area_4 = area;}
if ((a+1)===5){ var name_5 = name;var area_5 = area;}
if ((a+1)===6){ var name_6 = name;var area_6 = area;}
if ((a+1)===7){ var name_7 = name;var area_7 = area;}
if ((a+1)===8){ var name_8 = name;var area_8 = area;}
if ((a+1)===9){ var name_9 = name;var area_9 = area; }
if ((a+1)===10){var name_10 = name;var area_10 = area; }
Set添加属性
上述函数,获得了每个乡镇在每一年每一种地物的名称与面积,使用set函数,将其设置到每个乡镇要素的属性中:
//将属性值添加到矢量feat中
return feat.set(name_1, area_1, name_2, area_2, name_3, area_3, name_4, area_4, name_5, area_5,
name_6, area_6, name_7, area_7, name_8, area_8, name_9, area_9, name_10, area_10)
多要素多年份统计
List循环每一年
//代表年份的列表,一共15年,1代表2001年,15代表2015年
var list = ee.List.sequence(0, 15);
//newfeat是通过一个迭代,针对每个要素进行多年统计
var newfeat = ee.Feature(list.iterate(single_feature_single_year, feature));
return newfeat;
Map循环每一个feature
//是针对每个要素进行多年统计的函数
var single_feature_static = function(feature) {...}
//针对多个要素进行多年统计
var result = LeZhi.map(single_feature_static);
导出结果
//导出统计结果数据
Export.table.toDrive({
collection: result,
description:'Lezhi_landcover_statics'
});
至此我们的数据已经统计完毕,待数据下载后,我们就能获得数据表格。将导出的结果数据CSV重新编码为:UTF-8-BOM,然后放进excel中即可查看。导出的结果统计表包含每个乡镇2000年-2015年每一种地物的面积,非常方面土地覆盖变化的分析:
数据统计与分析
首先我想看看乐至县的各个乡镇的耕地变化情况,结果出乎我的意料,原本以为会有很多弃耕的耕地,但是结果数据表示,各乡镇的耕地面积变化不大。只有天池镇的耕地面积下降。
我们专注于看不透水面的面积:
地物8为不透水面,可以看到乐至县这十几个乡镇除了天池镇,其他乡镇基本没有扩张。
天池镇的建筑面积也由2000年的3000亩发展到2015年的10000亩,这也是天池镇耕地面积下降的原因。
回忆一下,我老家是在宝林镇,但10多年前和现在的宝林镇基本上没有变化。乡镇的资源基本都往县城(即天池镇)倾斜,现在乡镇的小孩基本都在县城里读书。一方面是不断进步的小县城,另一方面是不断凋零的乡镇与农村。
其他
读者可以把数据源换为感兴趣的数据源(比如10米、300米、500米的土地覆盖数据),也可以扩展研究区,比如研究长时间序列的四川省各个县城的土地覆盖变化情况,这个代码都能满足读者需求。
全部源代码
链接:https://code.earthengine.google.com/f8028d085dda6f8e098847e50355c2b0
/*
函数作用:计算长时间序列的多个要素的土地覆盖面积信息
hi替换为读者的研究区域
注:如果使用其他土地覆盖影像计算,需要保证collection、band和年份的对应关系
*/
//加载土地覆盖数据 该数据是2000-2015年的土地覆盖数据,参考文献:刘小平等.全球2000年—2015年30m分辨率逐年土地覆盖制图
var imageCollection = ee.ImageCollection("users/xxc/GLC_2000_2015");
//加载逐年土地覆盖影像
for(var year=0; year<16; year++){
var year_LEZHI=ee.Number(2000).add(year)
//根据名称波段选择土地覆盖影像
var year_band=ee.String("land_cover_").cat(ee.Number(year_LEZHI))
var Lezhi_image = imageCollection.select(year_band)
//配色方案设定,这里用的开源ee.palette,github地址为:https://github.com/gee-community/ee-palettes
var palettes = require('users/gena/packages:palettes');
var palette = palettes.colorbrewer.Dark2[3,4,5,6,7,8];
//加载逐年的土地覆盖影像
if (year>=10){Map.addLayer(Lezhi_image,{min: 10, max: 100, palette: palette},'20' +year +"年土地覆盖影像");}
if (year<10){Map.addLayer(Lezhi_image,{min: 10, max: 100, palette: palette},'200' +year +"年土地覆盖影像");}
}
//加载研究区域
Map.addLayer(LeZhi, {color: 'FF0000'}, "乐至县");
Map.centerObject(LeZhi, 9);
//裁剪研究区域
var Lezhi_image=imageCollection.filterBounds(LeZhi).max().clip(LeZhi)
//将image转为collection
var collection=ee.ImageCollection.fromImages([
Lezhi_image.select("land_cover_2000"),
Lezhi_image.select("land_cover_2001"),
Lezhi_image.select("land_cover_2002"),
Lezhi_image.select("land_cover_2003"),
Lezhi_image.select("land_cover_2004"),
Lezhi_image.select("land_cover_2005"),
Lezhi_image.select("land_cover_2006"),
Lezhi_image.select("land_cover_2007"),
Lezhi_image.select("land_cover_2008"),
Lezhi_image.select("land_cover_2009"),
Lezhi_image.select("land_cover_2010"),
Lezhi_image.select("land_cover_2011"),
Lezhi_image.select("land_cover_2012"),
Lezhi_image.select("land_cover_2013"),
Lezhi_image.select("land_cover_2014"),
Lezhi_image.select("land_cover_2015"),
])
//将一年影像按照类别拆分为10个单波段的image
var splitdataset = function(img){
//只有值相同的像素点才保留,并将添加一个波段的属性val,val的值为类别
//10耕地、20林地、30草地、40灌木、50湿地、60水体、70苔原、80不透水面、90裸地、100永久冰雪
var class1 = img.eq(10).selfMask().set('val', 1);
var class2 = img.eq(20).selfMask().set('val', 2);
var class3 = img.eq(30).selfMask().set('val', 3);
var class4 = img.eq(40).selfMask().set('val', 4);
var class5 = img.eq(50).selfMask().set('val', 5);
var class6 = img.eq(60).selfMask().set('val', 6);
var class7 = img.eq(70).selfMask().set('val', 7);
var class8 = img.eq(80).selfMask().set('val', 8);
var class9 = img.eq(90).selfMask().set('val', 9);
var class10 =img.eq(100).selfMask().set('val', 10);
//返回一个collection
return ee.ImageCollection.fromImages([
class1, class2, class3, class4, class5, class6, class7, class8, class9, class10]);
};
//初始化变量
var dataset
var feat
var year
var nd
//代表年份的列表,一共15年,1代表2001年,15代表2015年
var list = ee.List.sequence(0, 15);
//是针对每个要素进行多年统计的函数
var single_feature_static = function(feature) {
//addBands是针对每个要素进行一年统计的函数
var single_feature_single_year = function(year, feat){
//地物类别数量
var totalClass = 10;
//选择年份
year = ee.Number(year).toInt()
//拼接波段名称 例如2000年就是:land_cover_2000
var actual_year = ee.Number(2000).add(year)
//选择指定年份的波段
dataset = ee.Image(collection.toList(16).get(ee.Number(year)))
//调用上面的splitdataset函数,并返回带有10种地物image的collection
nd = splitdataset(dataset);
//将待统计的矢量要素转为feature
feat = ee.Feature(feat);
//循环地物类别,统计每种地物面积
for(var a=0; a<totalClass; a++){
var showList = nd.filterMetadata('val', 'equals', a+1).map(function(img){
return img.clip(feat);
})
var combineVal = showList.toBands().reduceRegion({
reducer: ee.Reducer.sum(),
geometry: feature.geometry(),
scale: 30
});
//特定的波段的名称拼接
var cls = ee.String(ee.Number(a)).cat(ee.String("_land_cover_").cat(ee.Number(actual_year)));
//统计特定地物的面积,getNumber函数的作用是获取字典中该属性的value
var area = ee.Number(combineVal.getNumber(cls).multiply(900));
//待会导出表的名称
var name = ee.String(ee.Number(actual_year)).cat(ee.String("_Class_")).cat(ee.Number(a+1))
//给每种地物的面积赋值:10耕地、20林地、30草地、40灌木、50湿地、60水体、70苔原、80不透水面、90裸地、100永久冰雪
if ((a+1)===1){ var name_1 = name;var area_1 = area;}
if ((a+1)===2){ var name_2 = name;var area_2 = area;}
if ((a+1)===3){ var name_3 = name;var area_3 = area;}
if ((a+1)===4){ var name_4 = name;var area_4 = area;}
if ((a+1)===5){ var name_5 = name;var area_5 = area;}
if ((a+1)===6){ var name_6 = name;var area_6 = area;}
if ((a+1)===7){ var name_7 = name;var area_7 = area;}
if ((a+1)===8){ var name_8 = name;var area_8 = area;}
if ((a+1)===9){ var name_9 = name;var area_9 = area; }
if ((a+1)===10){var name_10 = name;var area_10 = area; }
}
//将属性值添加到矢量feat中
return feat.set(name_1, area_1, name_2, area_2, name_3, area_3, name_4, area_4, name_5, area_5,
name_6, area_6, name_7, area_7, name_8, area_8, name_9, area_9, name_10, area_10)
};
//newfeat是通过一个迭代,针对每个要素进行多年统计
var newfeat = ee.Feature(list.iterate(single_feature_single_year, feature));
return newfeat;
};
//针对多个要素进行多年统计
var result = LeZhi.map(single_feature_static);
print(result)
//导出影像数据
Export.image.toDrive({
image: Lezhi_image,
scale: 30,
region:LeZhi,
description:'Lezhi_landcover_2000-2015'
});
//导出统计结果数据
Export.table.toDrive({
collection: result,
description:'Lezhi_landcover_statics'
});
参考:
许晓聪, 李冰洁, 刘小平,等. 全球2000年—2015年30 m分辨率逐年土地覆盖制图[J]. 遥感学报, 2021, 25(9):21.
https://gis.stackexchange.com/questions/390948/resolving-dictionary-does-not-contain-key-error-applying-function-to-feature-c
https://medium.com/@srivastava.saumyata/computing-lulc-zonal-statistics-using-modis-on-google-earth-engine-7b97a408cf24