打开APP
userphoto
未登录

开通VIP，畅享免费电子书等14项超值服

开通VIP

首页

好书

留言交流

下载APP

联系客服
【好文推荐】统计多个研究区多年土地覆盖变化教程
userphoto 关注

GEEer成长日记 2022-04-10 12:00

编者荐语：

锐哥YYDS！

以下文章来源于锐多宝的地理空间 ，作者锐多宝

锐多宝的地理空间想写一千篇博客的技术宅男

Part1引言

GEE，google earth engine，是一个地理数据的云计算平台，可以快速、批量处理数据。

我想统计一下自己的家乡，四川省乐至县，这十几年的各乡镇土地覆盖变化情况。考虑到本地计算需要进行繁琐的数据处理，我这里使用GEE进行计算。

乐至县各乡镇区划图

Part2研究思路

我将研究过程分为四个大的步骤，前三个步骤在GEE中进行，最后一个步骤是在本地对excel数据进行处理。

统计多年土地覆盖流程图

Part3数据筛选与预处理

这一个步骤主要是对影像数据进行土地覆盖数据加载、研究区加载以及数据标准化。

1逐年影像加载

我在GEE中找到了几个土地覆盖数据：ESA数据与esri数据，10米分辨率，但只有2020年。modis数据，有几十年的数据，但分辨率为500米，太低。

这里我采用的数据是刘小平老师制作的全球2000年—2015年30m分辨率逐年土地覆盖数据。

AGLC-2000-2015

确定好数据之后，按照年份顺序加载到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个波段，每一个波段是一个年份，单个波段的像素值是一种地物：

下载后的本地影像

2转为collection

在影像预处理步骤，我们需要将乐至县包含15年土地覆盖波段的影像转为影像合集，即由image转为collection，原先的band转为image。

image转为collection
//将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"),
  ])

Part4单个要素单年份统计

在进行多年多要素统计之前，我们需要写几个函数实现单个要素某一年的土地覆盖统计功能。

3Mask掉其他类别

为了待会能够使用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]);
};

4reduce统计地物类别

在每一个要素中，选择一个年份的土地覆盖影像dataset，将dataset进行mask，获得的结果是一个collection，又将collection通过toBands函数转为image，并使用reduce函数进行image的像素统计。

reduce统计地物类别示意图
//是针对每个要素进行多年统计的函数
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; }

5Set添加属性

上述函数，获得了每个乡镇在每一年每一种地物的名称与面积，使用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)

Part5多要素多年份统计

6List循环每一年

 
//代表年份的列表,一共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;

7Map循环每一个feature

//是针对每个要素进行多年统计的函数
var single_feature_static = function(feature) {...}
//针对多个要素进行多年统计
var result = LeZhi.map(single_feature_static);

8导出结果

//导出统计结果数据
Export.table.toDrive({
  collection: result,
  description:'Lezhi_landcover_statics'
});

至此我们的数据已经统计完毕，待数据下载后，我们就能获得数据表格。将导出的结果数据CSV重新编码为：UTF-8-BOM，然后放进excel中即可查看。导出的结果统计表包含每个乡镇2000年-2015年每一种地物的面积，非常方面土地覆盖变化的分析：

1649158900166.png

Part6数据统计与分析

首先我想看看乐至县的各个乡镇的耕地变化情况，结果出乎我的意料，原本以为会有很多弃耕的耕地，但是结果数据表示，各乡镇的耕地面积变化不大。只有天池镇的耕地面积下降。

乐至县各个区域的耕地面积统计

我们专注于看不透水面的面积：

乐至县各个区域的不透水面统计

地物8为不透水面，可以看到乐至县这十几个乡镇除了天池镇，其他乡镇基本没有扩张。

乐至县城区面积发展趋势

天池镇的建筑面积也由2000年的3000亩发展到2015年的10000亩，这也是天池镇耕地面积下降的原因。

回忆一下，我老家是在宝林镇，但10多年前和现在的宝林镇基本上没有变化。乡镇的资源基本都往县城（即天池镇）倾斜，现在乡镇的小孩基本都在县城里读书。一方面是不断进步的小县城，另一方面是不断凋零的乡镇与农村。

Part7其他

读者可以把数据源换为感兴趣的数据源（比如10米、300米、500米的土地覆盖数据），也可以扩展研究区，比如研究长时间序列的四川省各个县城的土地覆盖变化情况，这个代码都能满足读者需求。

Part8全部源代码

链接：https://code.earthengine.google.com/f8028d085dda6f8e098847e50355c2b0

/*
函数作用：计算长时间序列的多个要素的土地覆盖面积信息
作者：锐多宝的地理空间
公众号：锐多宝的地理空间
使用方法：将研究区域LeZhi替换为读者的研究区域
注：如果使用其他土地覆盖影像计算，需要保证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

喜欢此内容的人还喜欢

GEEer成长日记

不喜欢

确定

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。
打开APP，阅读全文并永久保存 查看更多类似文章
猜你喜欢
类似文章
【热】打开小程序，算一算2024你的财运
如何实现分享网站文章到微信朋友圈时显示指定缩略图或LOGO
PHP函数库範例(1)
类WebOS（添加了主界面，及相关功能代码）
curl 模拟登录微信公众平台带验证码
经验分享：10个简单实用的 jQuery 代码片段
JQuery中each()的使用方法说明
更多类似文章 >>
生活服务
热点新闻
留言交流
回顶部
联系我们
分享 收藏 导长图 关注 下载文章
绑定账号成功
后续可登录账号畅享VIP特权！
如果VIP功能使用有故障，
可点击这里联系客服！

联系客服