（文章中代码全部由ChatGPT生成）

在GIS的分析中，有一个十分常见的问题：在很多的实操中经常需要我们根据栅格数据生成一组与每一个像元完全重叠的格网【注意是完全重叠而不仅仅是与每个像元等大，以往生成等大的格网时每个网格经常会与对应的栅格像元发生错位，这是这个问题的关键和难点所在】。有了这个格网能够帮助实现分区统计等各种后续操作目的。被这个问题困扰了好久，在网上也一直没有找到能完全实现该目的的教程。如果使用ArcGIS操作会相对容易一些，大概思路如下：

1.获取栅格的四至坐标和分辨率（在栅格属性中即可查看）；

2.以四至坐标为范围，以分辨率作为宽度（单位用度或者米与栅格属性中的分辨率保持一致，ArcGIS同时支持以度或者米为单位创建渔网，这里为后面GEE埋个雷），使用“创建渔网”工具即可生成符合要求的格网。

大概的思路还是很简单的，在ArcGIS上也相对比较容易实现。但是在GEE当中原来一直没将其实现出来，网上也几乎找不到能够有帮助的代码。这次试着在GEE上真正实现一下，实现过程中踩了很多雷，主要是因为GEE中生成格网时如果要用buffer的方法必须以米为单位，而常见的栅格数据的坐标系都是地理坐标系WGS84，因此需要涉及度和米的转换。下面是实现的全过程（包括中间踩的雷的记录），真正实现上述目的的代码和结果附在最后面，可以直接运行。

【数据与目的】

使用的数据为GPM降水量数据，空间分辨率为0.1°，赤道上约合11132m（但实操时需要经纬度转换）。目的是将其提取到北京的范围，并根据降水栅格对应的每个像元生成完全重叠（大小一致，空间位置也完全相同）的格网。

【尝试1（踩的第一个坑）】

第一次尝试时大致的思路是：1.将降水栅格提取到北京的四至范围；2.提取每个栅格像元的中心点；3.对于每个中心点，向外生成宽度为0.1°（这里直接用11132m）的外接正方形格网。

gridsize=11132;

gridcell=point.buffer(gridsize/2).bounds();

再转换为Feature。

然后就得到了这个结果（错误的代码就不附了，最后附上正确的）：

可以发现基本思路是对的，但是有一个细节问题：格网南北方向的边界都与栅格像元完全一致，但东西方向有一定的重叠。这是因为地理坐标系的0.1°并不时刻等于111132米，这个关系仅在赤道上成立，而随着纬度升高，东西方向的距离实质上会缩减。比如我们在GEE里导出栅格时设置其分辨率为11132m，导出坐标系设置为WGS84，然后在ArcGIS中打开时，实际测量可以发现其南北方向分辨率确实为11km左右，但是东西方向分辨率仅为8km多。这就启示我们要注意度和米之间的换算，不能直接使用11132km作为边长。

【尝试2（还是重叠，但没想明白为什么不对）】

所以下一步尝试保持南北方向长度不变，将东西方向的宽度根据纬度进行距离换算：

使用 经度方向的缩放因子，即：cos(纬度)。

对应的公式为：实际经度上的距离 = 111320 米 * cos(纬度)。

这个方法从理论上应该是最直观和正确的，但不知为什么生成的结果还是和第一步一样，现在也没想明白。接下来只能稍微绕一下弯子，用另外的办法了。

【尝试3（踩的另外一个坑）】

那么有什么其他方法呢？这里可以用一个稍微复杂一点的方法：

任意选取两个相邻的点，先计算出可能重叠的距离，再用（两点间的东西距离-重叠距离/2）作为东西方向的宽度生成格网。

然后生成的结果是这样的......：

可以发现最北的一行的东西方向真的不重叠了，但是越往南重叠的区域越大。说明还是不能以一个固定的米数作为宽度，即使在北京这样大的一个区域这种由单位换算带来的误差仍不能忽略。所以问题是出在不能任意选择两个东西相邻的点计算重叠距离，而应该严格地从每一行中任意选择两个东西相邻的点计算重叠距离，再按照刚才的方法减去这个距离生成格网。

【尝试4（这回终于成功了）】

那接下来就按照正确的思路，先从每一行中任意选择两个东西相邻的点计算重叠距离，再按照刚才的方法减去这个距离的一半生成格网。

最后生成的中心点和格网如下：

可以看到这次的结果终于实现目的了，每个网格的边界都与对应的栅格像元完全重合。

导入到ArcGIS里看一下，结果也是一样的（蓝色为栅格像元，灰色为格网）：

至此终于完全解决了如何在GEE中生成与栅格像元完全重叠的格网的问题！

当然这个方法有些取巧，之后可能还会改进更好的方法。

【更新：最后的方法优化】

在之后修改代码时发现前面一直是卡在单位换算上，即先要把度换算为米再用buffer的方法从中心点向外生成格网。但如果不用这种向外缓冲的方法是否就能避免这样的单位换算呢？这里回归到另外一个最简单的方法：先确定每个像元的中心点坐标（以度为单位），接着获取每个点上下左右0.05度（边长的一半）的经纬度坐标（依然以度为单位），将这四个坐标作为对应网格的坐标四至直接生成矩形，中间没有涉及任何单位换算，即得到的格网矩形边长都为准确的0.1度。生成的结果和上面完全一样。

附上正确的示例代码（链接如下：https://code.earthengine.google.com/ac60284dd941d658eefb75ce9538c587?accept_repo=users%2Fsofiaermida%2Flandsat_smw_lst，可以直接运行，有需要自取）

（ps：B站现在发专栏插入不了代码块了）

// 加载北京市的矢量文件

var BJ = ee.FeatureCollection("users/pengzijie0611/Beijing_WGS84");

// 获取北京市的四至坐标

var BJ_bounds = BJ.geometry().bounds();

print('Beijing Bounds:', BJ_bounds);

// 加载 GPM 数据集（全球降水日数据）

var GPM = ee.ImageCollection("projects/climate-engine-pro/assets/ce-gpm-imerg-daily")

        .filterDate('2020-07-01', '2021-01-05') // 使用一段时间的数据，确保有数据

        .mean(); // 计算平均降水

// 获取 GPM 栅格的投影信息

var proj = GPM.projection();

// 获取栅格的分辨率（单位：米）

var scale = proj.nominalScale(); // 使用 nominalScale() 获取栅格分辨率

print('Scale (meters):', scale);

// 将北京市的四至坐标作为感兴趣区域（ROI）

var geometry = BJ_bounds;

// 获取每个像元的中心点和降水值信息

var points = GPM.sample({

 region: geometry, // 在北京市区域内采样

 scale: 11132, // 使用适当的分辨率

 projection: proj, // 保持与 GPM 数据的投影一致

 geometries: true // 生成包含几何（即中心点）的 FeatureCollection

});

// 打印中心点的坐标

print('Sampled points:', points);

// 四舍五入纬度值到小数点后三位

var roundedPoints = points.map(function(feature) {

 var lat = ee.Number(feature.geometry().coordinates().get(1)).multiply(1000).round().divide(1000); // 四舍五入到小数点后三位

 return feature.set('rounded_lat', lat); // 添加四舍五入后的纬度作为属性

});

// 将点按照纬度和经度分别排序

var sortedPoints = roundedPoints.sort('rounded_lat').sort('longitude');

// 遍历每个点，生成网格

var grid = sortedPoints.map(function(point) {

 var lon = ee.Number(point.geometry().coordinates().get(0)); // 获取当前点的经度

 var lat = ee.Number(point.geometry().coordinates().get(1)); // 获取当前点的纬度

 // 计算东西方向的缓冲距离，假设为 0.1 度

 var eastWestBuffer = ee.Number(0.1); // 0.1度的缓冲距离

 // 计算南北方向的缓冲距离，设定为 11132米

 var northSouthBuffer = ee.Number(11132).divide(111320); // 南北方向对应的度

 // 计算矩形格网的边界

 var west = lon.subtract(eastWestBuffer.divide(2));

 var east = lon.add(eastWestBuffer.divide(2));

 var south = lat.subtract(northSouthBuffer.divide(2));

 var north = lat.add(northSouthBuffer.divide(2));

 // 创建矩形

 var gridCell = ee.Geometry.Rectangle([west, south, east, north]);

 // 返回生成的格网

 return ee.Feature(gridCell);

});

// 将格网显示在地图上

Map.centerObject(BJ, 8); // 将地图居中到北京市

Map.addLayer(GPM.clipToCollection(BJ), {}, 'GPM Data (clipped)'); // 显示裁剪后的 GPM 数据

Map.addLayer(sortedPoints, {color: 'red'}, 'Pixel Center Points'); // 显示像素中心点

Map.addLayer(grid, {color: 'blue'}, 'Grid Cells'); // 显示生成的格网

// 打印生成的格网信息

print('Generated Grid:', grid);

// 导出矢量网格为 Shapefile

Export.table.toDrive({

 collection: grid, // 要导出的 FeatureCollection

 description: 'Beijing_Grid_Export', // 导出任务的名称

 fileFormat: 'SHP' // 文件格式，Shapefile

});

// 导出 GPM 栅格数据

Export.image.toDrive({

 image: GPM.clipToCollection(BJ), // 将影像裁剪到北京市范围

 description: 'Beijing_GPM_Export', // 导出任务的名称

 scale: 11132, // 分辨率

 region: BJ.geometry().bounds(), // 裁剪区域

 fileFormat: 'GeoTIFF', // 文件格式

 crs: 'EPSG:4326' // 投影为 WGS84

});