GEE学习之土地覆盖产品统计

时而支楞时而摆烂的小刘

已于 2023-04-07 13:39:28 修改

阅读量924

点赞数 3

文章标签：学习 googlecloud

于 2023-04-07 11:39:01 首次发布

本文链接：https://blog.csdn.net/qq_50291870/article/details/130008105

版权

想通过长时序土地覆盖产品，调查某地土地覆盖类型变化，并通过GIF图查看。

加载数据集

var CLCDCol = ee.ImageCollection("users/GEEStudyRoom/CLCDImageCollection");
//可以通过slice来筛选需要的年份
//var LC_col = ee.ImageCollection(CLCDCol.toList(32).slice(12))

这里介绍一下数据集：基于Google Earth Engine上335,709景Landsat数据，武汉大学黄昕老师团队制作了中国年度土地覆盖数据集（annual China Land Cover Dataset, CLCD），包含1985—2019中国逐年土地覆盖信息。基于5,463个目视解译样本，CLCD的总体准确率达80%。

详细信息：武汉大学遥感院黄昕&amp李家艺团队

添加时间信息并绘制时序曲线

var add_time = function(image){
  var name = ee.String(image.get('system:index')).slice(9)
  var starttime = ee.Date.fromYMD(ee.Number.parse(name),1,1)
  var endtime = ee.Date.fromYMD(ee.Number.parse(name),12,31)
  var image = image.set('system:time_start',starttime.millis()).set('system:time_end',endtime.millis())
  return image

}
var col = CLCDCol.map(add_time)

//以下是绘图代码，一般只需要更改palette和names
function addLegend(palette, names) {
 //图例的底层Panel
 var legend = ui.Panel({
   style: {
     position: 'bottom-left',
     padding: '5px 10px'
   }
 });
 //图例标题
 var title = ui.Label({
   value: '图例',
   style: {
     fontWeight: 'bold',
     color: "red",
     fontSize: '16px'
   }
 });
 legend.add(title);

 //添加每一列图例颜色以及说明
 var addLegendLabel = function(color, name) {
       var showColor = ui.Label({
         style: {
           backgroundColor: '#' + color,
           padding: '8px',
           margin: '0 0 4px 0'
         }
       });


       var desc = ui.Label({
         value: name,
         style: {margin: '0 0 4px 8px'}
       });
     //颜色和说明是水平放置
       return ui.Panel({
         widgets: [showColor, desc],
         layout: ui.Panel.Layout.Flow('horizontal')
       });
 };

 //添加所有的图例列表
 for (var i = 0; i < palette.length; i++) {
   var label = addLegendLabel(palette[i], names[i]);
   legend.add(label);
 }  

 Map.add(legend);
}


var Landcover = col.select('b1');
//计算面积
var LandcoverArea = Landcover.map(function(image){
  var area = ee.Image.pixelArea().addBands(image).reduceRegion({
        reducer: ee.Reducer.sum().group({
          groupField: 1, 
          groupName: "b1"
        }), 
        geometry: table, 
        scale: 30,
        maxPixels:10e15,
      });
  var image = image.set({"area":area});
  return image;
});
print('LandcoverArea',LandcoverArea)
// 把字典面积数据转换成feature
var area = LandcoverArea.reduceColumns(ee.Reducer.toList(), ['area'])
                            .get('list');

var time = LandcoverArea.reduceColumns(ee.Reducer.toList(), ['system:time_start'])
                            .get('list');

var feature = ee.List(time).zip(area).map(function(timeArea){
  var time = ee.Date(ee.List(timeArea).get(0)).get("year");
  var area = ee.List(timeArea).get(1);
      area = ee.Dictionary(area).get("groups");
  var feature = ee.List(area).map(function(dict){
      dict = ee.Dictionary(dict);
  var value = dict.values();
    return ee.Feature(null, { 
      time: time,
      land:value.get(0),
      area:value.get(1)});
    });
  return feature;
});

feature = feature.flatten();
var feature = ee.FeatureCollection(feature);
//print(feature,"feature");

var palette = ['66CD00', '006400', '00FF7F', '7FFF00', '1C86EE', 'FFFFFF', '8B7500','8C8C8C', '698B69']
var names = [ 'cropland',
              'forest',
              'shrub',
              'grassland',
              'water',
              'Snow/Ice',
              'Barren',
              'Impervious',
              'Wetland',]

// 面积统计（有时候也需要改一改）
var chart = ui.Chart.feature.groups(feature, "time", "area", "land")
.setSeriesNames({1:'cropland',
              2:'forest',
              3:'shrub',
              4:'grassland',
              5:'water',
              6:'Snow/Ice',
              7:'Barren',
              8:'Impervious',
              9:'Wetland',
              })//1Cropland、2Forest、3Shrub、4Grassland、5Water、6Sonw/Ice、7Barren、8Impervious、9Wetland
    .setOptions({
      title: 'Landcover area',
      hAxis: {'title': 'year'},
      vAxis: {'title': 'area (m^2)'},
      pointSize: 3,
      colors:['#66CD00', '#006400', '#7FFF00', '#1C86EE',  '#8B7500','#8C8C8C', '#698B69']
    })
print(chart)

// 地图显示（最大最小值根据分类更改）
var landcoverVisualization = {palette: palette, min: 1, max: 9, format: 'png'};
//展示一下第一张（1985年的）
Map.addLayer(Landcover.first(), landcoverVisualization, 'landcover');

//添加图例
addLegend(palette, names);

//研究区的边界
/*var outline = ee.Image().byte().paint({
  featureCollection: table,
  color: 1,
  width: 3
});
Map.addLayer(outline, {palette: ['blue']}, 'ROI');*/
/*

其中绘图函数基本不用更改，一般只修改palette和names，下面addlayer也是。计算面积那一部分代码也是借鉴的，基本不用更改。

得到研究区（自行设置，代码部分的table是我设置的研究区）的土地利用类型变化。

绘制GIF图

// 可视化参数
var args = {
  crs: 'EPSG:3857',  // Maps Mercator
  dimensions: '500',
  region: rect,
  min: 1,
  max: 9,//2000
  palette: ['66CD00', '006400', '00FF7F', '7FFF00', '1C86EE', 'FFFFFF', '8B7500','8C8C8C', '698B69'],
  framesPerSecond: 2, //自行设置帧数
};

// 制作动图，并且加载到地图上
var thumb = ui.Thumbnail({
//为“image”指定一个集合会使图像序列产生动画效果
  image: Landcover,
  params: args,
  style: {
    position: 'bottom-right', //位置
    width: '640px' //宽度
  }});
Map.add(thumb);

print(Landcover.getVideoThumbURL(args)) //获得链接

这里要注意，动图绘制一定是规则边界。所以要先画出一个矩形，也就是上面args里面的region，我命名为rect。

分区绘制

这里我沿用了之前用到的笨笨的方法。得到不同分区的土地利用变化时序曲线。其实和前面代码是相同的，只是把绘图区域变成了roi（不同的子区域）。不多赘述。

var point = table.aggregate_array('NAME')

for(var i = 0; i<5 ; i++){
  var roi = table.filterMetadata('NAME','equals',point.get(i));

  var LandcoverArea = Landcover.map(function(image){
  var area = ee.Image.pixelArea().addBands(image).reduceRegion({
        reducer: ee.Reducer.sum().group({
          groupField: 1, 
          groupName: "b1"
        }), 
        geometry:roi, 
        scale: 30,
        maxPixels:10e15,
      });
  var image = image.set({"area":area});
  return image;
  });
  print('LandcoverArea',LandcoverArea)
  // 把字典面积数据转换成feature
  var area = LandcoverArea.reduceColumns(ee.Reducer.toList(), ['area'])
                            .get('list');

  var time = LandcoverArea.reduceColumns(ee.Reducer.toList(), ['system:time_start'])
                            .get('list');

  var feature = ee.List(time).zip(area).map(function(timeArea){
    var time = ee.Date(ee.List(timeArea).get(0)).get("year");
    var area = ee.List(timeArea).get(1);
        area = ee.Dictionary(area).get("groups");
    var feature = ee.List(area).map(function(dict){
        dict = ee.Dictionary(dict);
    var value = dict.values();
      return ee.Feature(null, { 
        time: time,
        land:value.get(0),
        area:value.get(1)});
      });
    return feature;
  });
  feature = feature.flatten();

  var feature = ee.FeatureCollection(feature);
  // 面积统计
  var chart = ui.Chart.feature.groups(feature, "time", "area", "land")
  .setSeriesNames({1:'cropland',
              2:'forest',
              3:'shrub',
              4:'grassland',
              5:'water',
              6:'Snow/Ice',
              7:'Barren',
              8:'Impervious',
              9:'Wetland',
              })//1Cropland、2Forest、3Shrub、4Grassland、5Water、6Sonw/Ice、7Barren、8Impervious、9Wetland
    .setOptions({
      title: 'Landcover area',
      hAxis: {'title': 'year'},
      vAxis: {'title': 'area (m^2)'},
      pointSize: 3,
      colors:['#66CD00', '#006400',  '#7FFF00', '#1C86EE', '#8B7500','#8C8C8C', '#698B69']
      })
  print(chart)
}