GEE中为研究区添加经纬网格的完整教程与代码实现

做遥感或地理空间分析的朋友应该都有过这种经历论文或项目报告里的研究区图必须带经纬网格GEE里处理完数据想展示某个区域的分布情况也需要一个清晰的坐标参考。之前有不少人问我Google Earth Engine里到底怎么给选定的研究区加上经纬网格划分正好最近在做一个北京市的案例就把整个流程完整梳理出来。这篇内容以北京市为例但代码逻辑完全通用改成你手头的行政区、流域、样地范围都行。这里的经纬网格划分本质上是用GEE生成一组线要素经线和纬线叠加在底图和研究区边界上既能当制图底网也能为后续的网格采样、分块统计做准备。1. 为什么要在GEE里给研究区画经纬网格以及整体方案选型1.1 经纬网格的实际应用场景很多人第一反应是“经纬网格不就是画几条线嘛”但在实际项目里这个操作解决的往往是比较具体的问题。第一个场景是论文出图。大多数学术期刊对研究区位置图都有坐标要求配一张带经纬度注记的底图是基本操作。GEE里处理完遥感影像或矢量结果之后直接在图上叠一层经纬网格导出图片时就有清晰的坐标参考省去在PS里手工描线的麻烦。第二个场景是空间抽样设计。比如我想在全市范围内均匀布设50个采样点最常见的方法就是先做一个经纬网格然后按网格中心点或者按每个网格内的随机点来采样。这时候网格不只是“画线”它是一套可计算的空间框架。第三个场景是分块统计。举个例子你想统计北京市不同区域的NDVI均值、土地利用面积或者降水量分布完全可以按经纬网把研究区切成若干块然后逐块做reduceRegion或reduceRegions。网格本身可以带编号、带面积属性后续分析直接按属性分组就行。第四个场景比较容易被忽略检查数据和结果的位置偏差。在GEE里切换不同数据源时我习惯把经纬网格叠上去快速核对比如某个矢量文件是不是偏了某景影像的覆盖范围到底在哪有网格一眼就能看出来。1.2 方案对比不是只有一种“网格”可以选在GEE里做经纬网格我见过几种常见做法这里先把方案选型讲清楚免得一上来就写代码做完才发现选错方向。方案优点缺点适合场景线状经纬网LineString集合要素少、加载快、样式干净出图效果好没有面属性不能直接做分区统计制图底网、位置展示、坐标参考面状渔网Rectangle集合每个网格是完整Polygon能做面积统计、采样分区要素数量多绘制相对慢复杂边界下需要再裁剪网格化统计、抽样框构建外部GIS生成后导入QGIS/ArcGIS生成网格可批量添加经纬度注记排版精细需要来回导数据不能实时联动出版级专题图、论文正式插图用内置样式或在线WMS图层省事直接加载定制能力弱网格间距和范围不可控快速看个大概不适合正式分析我在这篇案例里做的是“线状经纬网”也就是只生成经线和纬线。原因是北京市的面积不大不小出图展示坐标参考线状网格足够而且加载速度非常快。如果是要做格网统计那就要用第二种面状网格我在第3章里也会给出扩展方案。1.3 GEE实现经纬网格的核心思路整个实现逻辑其实不复杂核心就四步第一步确定研究区边界并拿到它的外接矩形范围也就是经纬度南北东西的极值第二步根据设定的步长在经度方向上生成一组垂直线在纬度方向上生成一组水平线第三步把生成的所有线合并成一个FeatureCollection第四步用Map.addLayer加载到地图上。但这里有一个初学者最容易翻车的地方GEE分“客户端”和“服务端”两套对象。你在浏览器里正常写JavaScript的for循环循环变量必须是普通的JavaScript数字不能是ee.Number或ee.List否则循环根本跑不起来。换句话说需要先把研究区边界的经纬度极值“取出来”变成JavaScript数字才能进入循环生成线。这一步用getInfo()或者evaluate()来实现。虽然getInfo()在大量数据下不推荐用但网格生成的边界坐标只有四个数字用一次完全没问题。2. 核心代码拆解网格生成逻辑与参数细节2.1 获取研究区边界先解决一个问题北京市的边界从哪来。GEE里现成的全球行政区划数据集有很多我用的是FAO/GAUL这个经典数据集。它是分级的level 1相当于我国的省级所以过滤出北京的FeatureCollection再取geometry就是北京市的边界了。var beijing ee.FeatureCollection(FAO/GAUL/2015/level1) .filter(ee.Filter.eq(ADM1_NAME, Beijing)) .geometry(); Map.centerObject(beijing, 8);如果你手头不是省级边界而是自己做的样地、流域或者特定区域那就把矢量上传到GEE Assets里然后通过ee.FeatureCollection(users/你的用户名/你的文件名)加载再取geometry即可。这里有个建议上传矢量时尽量使用WGS84坐标系经纬网格本身是基于经纬度的如果原始边界是投影坐标系叠加时容易出问题。2.2 从外接矩形到经纬线的循环生成先拿北京市边界的外接矩形范围。这一步我用bounds()得到最小外接矩形然后从它的坐标列表里抽出所有点的横纵坐标再分别求最小值和最大值就得到了西、东、南、北四个边界值。var coords beijing.bounds().coordinates().get(0); var xs coords.map(function(pt) { return ee.List(pt).get(0); }); var ys coords.map(function(pt) { return ee.List(pt).get(1); }); var west ee.List(xs).reduce(ee.Reducer.min()).getInfo(); var east ee.List(xs).reduce(ee.Reducer.max()).getInfo(); var south ee.List(ys).reduce(ee.Reducer.min()).getInfo(); var north ee.List(ys).reduce(ee.Reducer.max()).getInfo();注意这里reduce返回的是ee.Number后面加了一个getInfo()作用就是把它变成JavaScript数字。到这一步你可以看到类似west115.4, east117.5这样的结果后面循环里才可以用。生成经线和纬线的核心逻辑很简单就是for循环。我习惯把起点的经纬度校准到步长的整数倍比如步长是0.5度那就从Math.floor(west/step)*step开始这样网格线会落在整数或半整数刻度上出图好看也方便读坐标。var step 0.5; var gridLines ee.FeatureCollection([]); for (var lon Math.floor(west / step) * step; lon east; lon step) { var line ee.Feature( ee.Geometry.LineString([[lon, south], [lon, north]]), {type: lon, value: lon} ); gridLines gridLines.merge(ee.FeatureCollection([line])); } for (var lat Math.floor(south / step) * step; lat north; lat step) { var line ee.Feature( ee.Geometry.LineString([[west, lat], [east, lat]]), {type: lat, value: lat} ); gridLines gridLines.merge(ee.FeatureCollection([line])); }这里我把每条线的属性和value写入Feature属性里type区分经线和纬线value记录这条线是东经多少度或者北纬多少度。这个设计后面会有用比如你想单独给经线设置不同颜色或者按经纬度值筛选某些线直接filter属性就行。最后加载到地图Map.addLayer(beijing, {color: black, width: 2}, 北京边界); Map.addLayer(gridLines, {color: red, width: 1}, 经纬网格 step °);2.3 网格间距怎么选不是越小越好我把间距问题单独拎出来说是因为见过太多人一上来就设0.01度结果地图上密密麻麻全是线完全没法看。网格间距的选择主要看两个因素研究区大小和出图比例尺。北京市大约跨经度2.1度、纬度1.9度用0.5度步长会生成大概5条经线和四五条纬线图上看起来干净清爽坐标参考也够清楚。如果做全国尺度的图我通常用1度或者2度如果是街道、园区这种小范围0.01度甚至0.001度才合适。另一个经验是网格只是辅助参考不要盖过数据本身。配色上建议用半透明的浅色灰或浅红线宽保持1左右。行政边界用黑色粗线突出研究区轮廓经纬网用细线弱化存在感。这样出图后主次分明审稿人看着也舒服。2.4 经纬度标注GEE目前的限制与替代方案这里必须说实话GEE Code Editor里给矢量线要素添加文字注记并不像桌面GIS那样方便它没有一个直接让线的旁边显示“116°E”这种标签的参数。很多教程默认不标经纬度只画线但论文里往往需要标注。我的做法有两种。第一种是导出到桌面GIS标注。把网格线导出成Shapefile或GeoJSON然后在QGIS或ArcGIS里批量加注记再配图。这个流程虽然多了几步但可控性最强线旁标注的字体、位置、旋转都可以精细调。第二种是在GEE里直接叠一个ui.Label把经纬度和网格信息作为图面信息显示在地图角落。如果你只需要告诉看图的人“这张图的网格间距是多少”这个办法最简单。Map.add(ui.Label(网格间距0.5°, {position: bottom-left}));不过如果确实想在地图本身上标注“116°E”“40°N”我的建议是导出网格后到QGIS处理这个我会在第3章导出部分再详细说。3. 完整可直接运行的案例北京市经纬网格加载3.1 完整代码可直接复制运行把上面散落的代码整合成一个完整脚本并封装成函数方便以后换研究区直接调用。这段代码在GEE Code Editor里可以一次跑通。// 1. 生成经纬网格的函数 function makeGrid(region, step) { var coords region.bounds().coordinates().get(0); var xs coords.map(function(pt) { return ee.List(pt).get(0); }); var ys coords.map(function(pt) { return ee.List(pt).get(1); }); var west ee.List(xs).reduce(ee.Reducer.min()).getInfo(); var east ee.List(xs).reduce(ee.Reducer.max()).getInfo(); var south ee.List(ys).reduce(ee.Reducer.min()).getInfo(); var north ee.List(ys).reduce(ee.Reducer.max()).getInfo(); var gridLines ee.FeatureCollection([]); for (var lon Math.floor(west / step) * step; lon east; lon step) { var vLine ee.Feature( ee.Geometry.LineString([[lon, south], [lon, north]]), {type: lon, value: lon} ); gridLines gridLines.merge(ee.FeatureCollection([vLine])); } for (var lat Math.floor(south / step) * step; lat north; lat step) { var hLine ee.Feature( ee.Geometry.LineString([[west, lat], [east, lat]]), {type: lat, value: lat} ); gridLines gridLines.merge(ee.FeatureCollection([hLine])); } return gridLines; } // 2. 导入北京市边界 var beijing ee.FeatureCollection(FAO/GAUL/2015/level1) .filter(ee.Filter.eq(ADM1_NAME, Beijing)) .geometry(); Map.centerObject(beijing, 8); // 3. 生成并显示网格 var grid makeGrid(beijing, 0.5); Map.addLayer(beijing, {color: black, width: 2}, 北京边界); Map.addLayer(grid, {color: red, width: 1}, 经纬网格 0.5°); print(网格线数量:, grid.size()); print(网格类型及数值示例:, grid.first());运行后地图上会显示北京市边界和一组红色经纬网。控制台会输出网格线数量和第一条线的属性方便你检查网格是否正常生成。3.2 换成任意研究区两处改动即可想换研究区只需要改函数调用部分。第一种方式是自己上传矢量var roi ee.FeatureCollection(users/你的用户名/你的shp文件名).geometry(); Map.centerObject(roi, 8); var grid makeGrid(roi, 0.1); Map.addLayer(roi, {color: black, width: 2}, 研究区边界); Map.addLayer(grid, {color: red, width: 1}, 经纬网格);第二种方式是用GEE公共数据集里的现成区域比如流域、生态区、气象站点范围。只需要把filter条件换成你需要的字段就行。网格步长也要跟着研究区大小调整我一般先print一下外接矩形的面积心里有个数再定步长。这里有个容易犯的错当你用非矩形的roi生成网格时网格覆盖的是roi的外接矩形所以边界外的区域也会有网格线。如果只想保留研究区内部的线需要加一步裁剪我在下一小节讲。3.3 进阶玩法把网格裁剪进研究区给每个格子算面积如果把线状网格升级成面状渔网就能做空间统计。封装一个生成面状网格的函数原理和线状网格几乎一样只是把LineString换成Polygon。function makeGridPolygon(region, step) { var coords region.bounds().coordinates().get(0); var xs coords.map(function(pt) { return ee.List(pt).get(0); }); var ys coords.map(function(pt) { return ee.List(pt).get(1); }); var west ee.List(xs).reduce(ee.Reducer.min()).getInfo(); var east ee.List(xs).reduce(ee.Reducer.max()).getInfo(); var south ee.List(ys).reduce(ee.Reducer.min()).getInfo(); var north ee.List(ys).reduce(ee.Reducer.max()).getInfo(); var features []; for (var x west; x east; x step) { for (var y south; y north; y step) { features.push(ee.Feature(ee.Geometry.Rectangle([x, y, x step, y step]))); } } return ee.FeatureCollection(features); } var gridPoly makeGridPolygon(beijing, 0.5); var gridClipped gridPoly.map(function(f) { return f.intersection(beijing); });用intersection把每个网单元裁剪到北京市边界内之后每条Feature的面积就是真实的市内面积。这里注意intersection在边界复杂时计算量不小如果网格数量上百个可能会比较慢。一个优化办法是先用一个简化版本的研究区边界beijing.simplify(100)能明显加快速度精度损失在网格统计场景下通常可以接受。给每个网格加上面积属性var gridWithArea gridClipped.map(function(f) { var area f.area(); return f.set(area_m2, area); }); print(带面积属性的网格, gridWithArea);有了面积属性后你可以筛选面积大于某阈值、等于最大面积的网格也可以和影像数据做reduceRegions统计每个网格内的NDVI均值、降水量等。这一步是把“画网格”真正变成“用网格”的关键。3.4 导出网格用于正式制图网格画好之后如果想在别的软件里精细排版加注记可以导出。GEE导出矢量最常用的方式是Export.table.toDrive。Export.table.toDrive({ collection: grid, description: beijing_grid_0_5deg, fileFormat: SHP });如果要GeoJSON把fileFormat改成GeoJSON。导出后去Google Drive下载然后拖进QGIS。在QGIS里选中线图层右键“属性”-“标签”选择“按value字段显示标签”value字段就是每条线对应的经纬度值。这样出来的图就是带完整经纬度注记的出版级底图。我自己的习惯是GEE里快速出网格导出SHP到QGIS做注记和排版最后矢量结果用Print Layout出图。这样两边优势都发挥出来了。4. 常见问题与避坑指南4.1 网格太密导致地图卡顿用0.01度步长给北京市画网格会生成几百条线Map.addLayer加载起来明显卡顿缩放时更明显。这不是GEE坏了是要素数量太多了。解决办法很简单预览时用大间距网格比如0.5度或1度真正做小尺寸网格统计时不要把所有网格都addLayer而是先算好结果再导出数据。如果必须在图上展示小网格可以考虑把所有线合并成一个MultiLineString要素大幅减少要素数量。4.2 网格线和研究区边界不贴合网格覆盖了外接矩形这个问题很常见。记住一个原则网格线的生成范围是研究区的外接矩形bounds不是研究区本身。所以如果研究区是不规则形状网格自然会延伸到边界之外。解决方法上面提过了用裁剪把网格限制在研究区范围内。对于线状经纬网可以用line.intersection(region)来裁剪但裁剪线要素会产生很多小线段不太好看。实际制图时我通常不会裁剪经线而是直接用外接矩形的网格做底网再用行政边界压盖在上面这样反而干净。4.3 经纬网格在高纬度地区看着不对劲GEE默认的地图投影是Web墨卡托越靠近两极面积变形越严重。经纬网格是按固定经纬度间隔生成的所以它在地图上是均匀的网格但每个格子的实际面积随着纬度升高在缩小。这时候要注意如果做格网统计用等经纬度网格不等于等面积网格。想在北京这样范围不大的区域做统计误差很小问题不大。但如果研究区跨度很大比如从海南到黑龙江那就要考虑用投影坐标系下的等面积网格比如按公里网生成格网而不是简单的经纬度网格。4.4 GEE登录不上、连接失败还有一些“下载收费”的问题这里也顺便回应一下网上经常搜到的问题。Google Earth Engine本身不是一个需要下载安装的桌面软件它是一个在线云平台直接在浏览器里访问就可以。个人用于科研、教育、非商业用途是免费的前提是注册账号并通过申请。如果遇到连接不上、登录超时这些问题我一般先检查本地网络是否稳定然后清理浏览器缓存、更换浏览器或换一个网络环境再试再不行就等官方服务恢复后再操作。有一点想提醒大家网上有些所谓的“下载器”“离线版”“破解工具”基本都不靠谱风险很大完全没有必要去碰。正规使用路径只有一条就是官方网站注册使用。5. 我个人在实操中的一点体会整套流程跑下来最花时间的往往不是写代码而是调试网格间距和样式。我前后试过很多次总结出一个通用习惯先跑通最小示例确认边界和网格数量合理再逐步调整到目标步长。比如给北京市出图我一般先用1度步长快速看位置然后切到0.5度看细节最后决定用0.25度还是0.5度出图完全取决于图上需要标注哪些经纬度节点。另外一个让我比较受益的小技巧给网格线的Feature设置type和value属性虽然只是多写两行但后续按属性筛选会非常方便。比如想只显示东经116.5度的线直接grid.filter(ee.Filter.and(ee.Filter.eq(type,lon), ee.Filter.eq(value,116.5)))不用重新生成网格。最后再说一句经纬网格看似基础但在GEE里把它做通透很多高级应用就顺了。无论是出图、采样还是分块统计这套代码都能直接当底座用。以后换研究区、换步长改两行参数就行。如果你也在做区域尺度的空间分析建议把这套函数存成GEE代码库真的能省不少时间。