GEE必须会教程——一文教会你研究区域分块计算与合并

还记得吗，上期小编教了你如何实现研究区域分块处理。有友友可能会问：为什么要分块？一整块算不是挺好的吗？emmmm，用过的人都知道，GEE好比幼年期的“紫妍”，要用“丹药”来养的，除非你是萧炎小友，到哪都有人给爆金币。很显然，大家都不是，你可别嫌麻烦，不分块就是算不出来啊啊啊！所以分块是专治GEE各种算力不服详情可见（http://t.csdnimg.cn/ALH3v），好了，不扯那么多了，今天我们在此基础上聊聊如何在分块上进行计算并且批量导出计算结果。

一、学习目标

1、巩固分块处理的代码基础

2、学会在分块基础上进行计算和批量下载

3、学会利用R语言进行合并

二、 案例应用

我们今天想做的是计算2020年夏季西藏自治区的植被覆盖度（FVC），FVC的估算需要使用归一化植被指数（NDVI数据），这里我们用到的产品是MODIS的相关数据产品，该产品的分辨率是25 0m。植被覆盖度的计算和分级可以查看小编的以往文章（http://t.csdnimg.cn/OAmgx）。

我们在做相关研究时，每次做到青藏地区，总会发现出现算力不足的情况，因为这个地方真的是太大了啊啊啊。比如这次，小编计算时产生的问题就是计算超时。这个时候，就要首先对西藏自治区进行分块处理了，处理完了再进行计算。

三、华丽上代码

Ⅰ 导入研究区域

//导入西藏自治区作为研究区域
var Province = ee.FeatureCollection("users/hesuixinya511/Province");
var roi = Province.filterMetadata("简称","equals","藏");
Map.centerObject(roi,6);
Map.addLayer(roi,{"color":"red"},"Tibet");

 此处，我们通过省份数据，按照属性筛选，得到西藏自治区的行政边界。

Ⅱ 研究区域分块处理

//定义切分研究区的函数，将研究区分为n*m块
function gridSplit (roi, n, m) {
var bounds = roi.geometry().bounds().getInfo().coordinates[0];
//为了确保边缘像元包括在范围内，提高容错率，将范围向四周扩大0.1
var xmin = bounds[0][0]-0.1;
var xmax = bounds[1][0]+0.1;
var ymin = bounds[0][1]-0.1;
var ymax = bounds[2][1]+0.1;
print(xmin,ymin,xmax,ymax);
// 计算分块的单位长度
var dx = (xmax-xmin)/n;
var dy = (ymax-ymin)/m;
//循环生成分块
var xx = ee.List.sequence(xmin, xmax, dx);
var yy = ee.List.sequence(ymin, ymax, dy);
//为了简化代码，这里生成了一些超出研究区范围的分块
var rects = xx.map(function(i){
  return yy.map(function(j){
  var x1 = ee.Number(i);
  var x2 = ee.Number(i).add(ee.Number(dx));
  var y1 = ee.Number(j);
  var y2 = ee.Number(j).add(ee.Number(dy));
  var coords = ee.List([x1, y1, x2, y2]);
  var rect = ee.Algorithms.GeometryConstructors.Rectangle(coords);
  return ee.Feature(rect);
  });
}).flatten();
//筛选出和研究区相交的分块
var rects_col = ee.FeatureCollection(rects).filterBounds(roi.geometry());
//计算总共有多少个分块
var GridNum =rects_col.size().getInfo();
print('GridNum: ', GridNum);
//为每个分块赋予编号
var idList=ee.List.sequence(0, GridNum-1);
var grid = ee.FeatureCollection(idList.map(function(i) {
  return ee.Feature(rects_col.toList(rects_col.size()).get(i)).set("grid_id",ee.Number(i).add(1));
}));
return grid;
}
//将研究区分为4*4块，其中有2块被筛掉，因此共分为14个小区
var subroi = gridSplit(roi, 4, 4);
//可视化所有小区
Map.addLayer(subroi,{}, 'subroi_all');
print(subroi,'subroi_list');
var roi5=subroi.filterMetadata('grid_id','equals', 15);
Map.addLayer(roi5,{}, 'subroi_5');

 此处对研究区域进行分块处理，主要搭建了gridSplit函数，将研究区域分成了n*m块，最后经过筛选得到了14块小区。

Ⅲ 分块计算FVC并批量导出

// 将FeatureCollection转换为列表
var list = subroi.toList(subroi.size());
print(list);
var Vis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};
// 定义一个函数来处理和导出每个地块
var processAndExport = function(index) {
  // 从列表中获取地块
  var feature = ee.Feature(list.get(index));
  function FVC_Cal_MOD (img){
    var dict = img.reduceRegion({
                      reducer: ee.Reducer.percentile([5, 95]),
                      geometry: feature.geometry(),
                      scale: 500,
                      maxPixels: 1e13,
                      tileScale: 16
                    });
    var _p5 = ee.Number(dict.get("NDVI_p5"));
    var _p95 = ee.Number(dict.get("NDVI_p95"));
    var imgFVC = ((img.subtract(_p5)).divide(_p95.subtract(_p5))).float();
    var FVC=((imgFVC.lt(0)).multiply(0))
    .add(((imgFVC.gt(0)).and(imgFVC.lte(1))).multiply(imgFVC))
    .add((imgFVC.gt(1)).multiply(1));
    return FVC.rename("FVC"); 
  }
  function addFVCbands(img){
    var FVC=FVC_Cal_MOD (img);
    return img.addBands(FVC);
  }
  //MODIS data
  var Dataset=ee.ImageCollection("MODIS/006/MOD13Q1")
                 .filterBounds(feature.geometry())
                 .filter(ee.Filter.calendarRange(6, 8, 'month'))
                 .filter(ee.Filter.calendarRange(2020, 2020, 'year'))
                 .select('NDVI')
                 .map(addFVCbands)
                 .mean()
                 .select('FVC')
                 .clip(feature.geometry())
  var exportParams = {
    image: Dataset.select('FVC'),
    description: 'Export_' + i,
    scale: 100,
    region: feature.geometry()
  };
   Map.addLayer(Dataset.select('FVC'), Vis, 'FVC'+i);
  // 导出图像
  Export.image.toDrive(exportParams);
};

此处设计了一个主函数用于对每个小块计算FVC，并在最后进行批量导出。该函数共包括两个方面，一个是根据NDVI计算FVC，一个是定义导出影像的参数函数。

// 使用客户端JavaScript循环遍历列表
var n = subroi.size().getInfo();
for (var i = 0; i < n; i++) {
  processAndExport(i);
}

最后，写一个for循环即可实现影像的批量导出。不过因为一共有14张影像，正常情况下我们需要在“task”栏中，一个一个点，经过小编一番资料的查阅，发现了一个好东西。来来来，我们打开商店，安装一个扩展程序——Open Earth Engine extension。

有友友可能会问，这个软件有什么用呢？看安装之后的效果：

欸，发现顶端多了一行欸“RUN ALL”，漂亮！！这样点一次就可以全部下载喽~，是不是超级无敌有用呢！大家还不快快用起来吧。

Ⅳ 结果拼接处理

下载好的结果毕竟是一块一块的，怎么对下载好的结果进行拼接，变成一张大图呢，别急，接下来，我们上科研小霸王——R语言。

library(raster)

# 设置文件夹路径
folder_path <- "F:\\result\\"

# 读取文件夹中所有的TIFF文件
tiff_files <- list.files(folder_path, pattern = "\\.tif$", full.names = TRUE)

# 初始化一个空的RasterLayer列表
rasters <- list()

# 逐个读取TIFF文件
for (file in tiff_files) {
  # 读取TIFF文件为RasterLayer
  r <- raster(file)
  # 将处理后的RasterLayer添加到列表中
  rasters[[length(rasters) + 1]] <- r
}

# 合并所有RasterLayer为一个单一的RasterLayer
merged_raster <- do.call(merge, rasters)

# 将合并后的RasterLayer保存为新的TIFF文件
plot(merged_raster)
writeRaster(merged_raster, filename="F:\\result\\merge.tif")

在ArcMap里面看看效果吧

 不过，上述代码仍然存在几个问题：

①分块的代码并不是严格按照西藏自治区的边界进行的，因此我们还需要再次进行裁剪。

②合并出来的结果，块与块之间分割明显。

因此代码仍然有许多优化的空间，各位友友有没有更好的解决办法，有的话私信小编噢！！

好了，今天的分享到这里就结束了，更多内容欢迎同步关注小编的公众号——梧桐GIS，咱们下期再会啦！！