做土地利用变化分析这些年我一直在等一套“既能看清细节、又不用自己从头训练模型”的全球土地覆盖数据。Esri在2020年放出的这套10米全球土地覆盖数据第一次把全球尺度的分类产品拉到了10米分辨率配合官方那套Land Cover Downloader下载入口从数据获取到出图统计整个流程比我想象中顺得多。这篇就把我从找数据、下载、投影处理、面积统计到踩坑排查的完整经历捋一遍给想用这套数据的人一个可以直接照着做的参考。这套数据适合谁做城市规划、生态评价、林业与农业宏观监测、气候变化研究的人基本都能用上。它帮你省掉下载海量Sentinel-2原始影像、自己跑深度学习模型分类的过程拿到手的已经是分好类的栅格结果。下面从数据本身的底细开始讲再到下载方式最后是实操中最容易翻车的坐标系和NoData问题。1. 这套数据是什么Esri 2020 10m全球土地覆盖的底细1.1 从原始影像到分类产品很多人第一次打开这套数据时会有个疑问这看起来不像卫星影像颜色一块一块的是不是把照片压扁了其实这正是它的核心价值所在。这套数据不是普通的遥感影像而是基于ESA Sentinel-2哨兵二号10米分辨率光学影像通过深度学习模型自动分类生成的“语义标签”栅格。具体来说Esri联合了Impact Observatory和微软的AI for Good团队用深度卷积神经网络对全球的Sentinel-2影像进行逐像元分类把每个10米×10米的像元归到10个地表类别之一。整个处理流程在云端完成最终以影像服务的形式发布在ArcGIS Living Atlas上。对普通用户来说你不需要自己下载任何原始影像不需要GPU不需要训练模型直接通过在线服务按范围导出即可。我看过官方披露的模型思路训练数据来自全球各地人工解译样本模型对光谱、纹理、植被指数等多维特征做融合分类。这套流程的厉害之处在于它把过去需要超级计算中心才能完成的任务变成了一个普通GIS用户可以随时调用的公开数据服务。也正是因为底层的自动化程度高它才能以10米这样的分辨率覆盖全球。1.2 10米分辨率到底意味着什么我先给几个对比数据传统的全球尺度土地覆盖产品比如ESA CCI-LC分辨率是300米GlobeLand30是30米而Esri这套是10米。10米这个档位意味着什么拿一个具体的场景说一条农村公路两侧的防护林带宽度大约15到20米在30米分辨率下基本就和道路糊在一起了但在10米分辨率下能清晰分出一条一条的树带。一个小村庄周边的零散菜地30米分辨率只能看到破碎的混合像元10米分辨率则可以数出地块边界。当然分辨率高不等于一切都好。10米数据对存储、传输、处理的要求也成倍增长全球一张10米栅格的体量是30米产品的9倍左右。这也是为什么官方并不推荐用户直接下载全球“一整张”tif而是建议按研究区范围、按需要分块导出。关于这个后面实操部分会详细说。从行业影响来看这套数据让很多中小团队第一次拥有了“免费、全球、高分辨率”的土地覆盖底图。过去做一个省级的生态评价光是准备土地覆盖数据就要花大量时间协调数据源现在打开ArcGIS Pro加个图层就能开始分析这是它最大的贡献。1.3 十个类别别把编号搞错这套数据的像元值就是类别编码不是连续的反射率数值。你每看到一个值必须对应到官方分类体系里去理解。官方共定义了10个类别编码英文名称中文含义简要说明1Water水体河流、湖泊、海洋、水库等常年有水区域2Trees树木郁闭度较高的林地包括天然林和人工林3Grass草地天然草地和稀疏草被区别于人工牧草4Flooded vegetation洪泛植被红树林、沼泽、季节性淹没植被5Crops农作物耕地、农田含一年生和多年生作物6Built Area建成区建筑、道路、硬化地面等不透水面7Bare ground裸地裸露土壤、沙地、岩石等8Snow/Ice冰雪常年积雪、冰川、冰盖9Clouds云云遮挡区域属于有效类别但不是地表类别10Rangeland灌丛草地以灌木和草本为主的自然植被区别于人工草地这里最容易搞混的是第3类Grass和第10类Rangeland。官方定义里Rangeland更偏向自然状态的灌木-草本混合植被而Grass更多指草为主的区域农业上的人工牧草通常会被分到Crops那一类。实际使用中如果只是看宏观植被格局把3和10合并成“草地灌丛”也很常见但严格按官方定义做统计时不要随便合并。另外特别注意第9类Clouds。很多人在统计面积时会把9当异常值删掉但实际上云是这套数据的“正常组成部分”全球不少区域尤其是热带雨林和季风区云遮挡比例相当高。处理时必须把它单独拎出来看占比它在分类结果里是白色的很容易和背景混淆后面第4章我会专门讲怎么处理。1.4 与同类全球产品的横向对比做科研或者写报告经常会被人问为什么不用GlobeLand30为什么不用FROM-GLC我把主流几套全球土地覆盖产品放在一起比较一下方便你选型。产品名称分辨率主要数据源类别数优势局限Esri 2020 Land Cover10米Sentinel-210类分辨率高、获取即时、全球一致性好类目较粗、有云污染、单期为主GlobeLand3030米Landsat等10类精度验证扎实、应用广泛更新周期长、下载流程较繁琐FROM-GLC10米/30米Landsat/Sentinel-230米版10类左右清华团队长期维护版本较多需仔细选择对应年份ESA CCI-LC300米MERIS/SPOT等22类时间序列长、适合气候模式分辨率太低区域细节不可用CGLS-LC100100米PROBA-V10类左右每年更新、有连续植被覆盖度层分辨率不如Esri我的建议是如果你做的是城市群、县域、省级尺度的分析Esri这套10米数据在分辨率上优势明显如果你要做长时间序列的变化分析CCI-LC的20多年连续产品更合适如果你对精度要求较严且只关心特定年份GlobeLand30的30米产品可以作为交叉验证。没有一套数据是万能的组合使用才是常态。2. 下载前的准备先想清楚范围与输出方式2.1 你需要的到底是“全球一张图”还是“研究区一张图”第一次用这套数据的人最容易犯的错就是试图把全球数据一次性下载回来。我见过群里有人问“哪里有全套10米全球土地覆盖数据直接下载”然后去找了几百GB的种子文件下了三天还没下完。其实官方并不建议这么做原因很简单全球10米分类栅格完整铺开的体量非常庞大普通电脑光打开一次就要几分钟更别说做分析了。正确的思路是先明确研究区。你在ArcGIS Pro里打开数据源浏览时看起来是“全世界都在”但真正需要的是某个流域、某个省、某个县的范围。先用矢量边界圈定范围再按范围导出这样导出的文件可能只有几十MB到几个GB处理起来非常舒服。还有一个细节范围框选时尽量把矢量边界的坐标系统一到WGS84经纬度。因为在线服务的导出接口通常使用地理坐标EPSG:4326来定义范围你直接用WGS84的bbox去请求最不容易出错。如果你手头是投影坐标的shp例如中国2000大地坐标系EPSG:4490/4547等建议先投影转WGS84再框定导出范围或者导出后再裁。2.2 两个下载入口Living Atlas与REST导出目前获取这套数据主要有两个入口都需要联网使用本质上是同一套云端影像服务只是访问方式不同。第一个入口是ArcGIS Living Atlas。打开ArcGIS Pro在目录窗格里的“门户”节点下选择Living Atlas搜索“Esri 2020 Land Cover”就能找到对应条目。双击加载到地图后可以直接预览也可以在图层面板右键选择“数据—导出栅格”来保存本地的tif文件。这个入口适合交互式操作边看边选范围比较直观。第二个入口是直接调用REST影像服务。在Living Atlas页面里找到图层服务URL把它复制出来可以在QGIS的“ArcGIS REST Servers”连接中加载也可以用脚本方式请求导出。这个入口适合批量处理你有十个县的范围要分别下载写个循环脚本跑一晚上就完事不用手动点二十次。我个人的习惯是首次接触先走Living Atlas确认数据范围和类别没问题后再转用脚本化导出。关于“Land Cover Downloader”这个名字其实指的就是这套数据配套的下载导出体系——无论是官方页面里的导出按钮还是第三方封装的下载小工具底层调用的都是影像服务的Export Image能力。所以只要理解了REST导出任何封装工具在你眼里都是透明的。2.3 把“下载器”理解成一个导出参数问题我用过几个不同版本的Land Cover Downloader下载工具体验参差不齐。有的工具只是把导出请求包装了一个界面有的则支持批量和断点续传。但你只要理解这背后的核心就不会被工具牵着走。一次导出请求本质上就是在问服务端这样几个问题你要哪块范围bbox即最小经度、最大经度、最小纬度、最大纬度你要什么分辨率官方源数据是10米但你可以指定导出30米以减小文件体量你要什么输出格式tif、png、jpg等你要用什么重采样方式分类数据必须用最近邻nearest你要什么坐标系默认3857也可直接导出为UTM等投影理解了这些参数你就可以用任何工具、任何语言去请求数据。这也是为什么我会在后面实操部分给一套兼顾ArcGIS Pro和GDAL的流程。3. 实操流程从拿到数据到出图统计3.1 ArcGIS Pro手动导出如果你是ArcGIS Pro用户手动导出整套流程是这样走的打开Catalog面板展开Portal进入Living Atlas搜索“Esri 2020 Land Cover”把影像图层拖进地图视图中。第一次加载可能会因为全球范围较大而显示比较慢这是正常的等绘图层刷新出来后再放大到你的研究区。确定范围后在图层面板上右键你加载的栅格图层选择“数据—导出栅格”。在导出窗口里关键是下面几个设置范围选“当前地图范围”或者手动输入研究区边界一般习惯选当前地图范围配合“按地图范围裁剪”使用像元大小保持10米不建议在这里强行放大比如改成5米因为源数据本身只有10米放大分辨率只是插值不会增加真实信息重采样必须要选“最近邻”Nearest Neighbor。分类数据是离散标签用双线性或三次卷积会算出“2.5类”这种无效值像素类型选8位无符号整型NoData设置为0或者跟随源数据的NoData设置。确认后点导出会生成一个本地tif文件。这里有一个官方文档不会告诉你的细节在线影像服务导出是有单次像素上限的。如果你的研究区特别大比如一整个省一次导出几万×几万像素很容易报错。解决办法是拆成多个小块分别导出后面用拼接工具合并。3.2 QGIS与GDAL脚本化下载如果你习惯开源工具链用QGIS加GDAL也能完成同样的事情而且脚本化之后效率会明显提升。先用QGIS加载服务数据源管理器选择ArcGIS REST Servers新建连接把Living Atlas条目里的服务URL贴进去注意URL应以/ImageServer结尾。连接成功后把图层加到画布右键“导出—另存为”照样设置tif格式和10米分辨率。这是交互式做法。想做批量导出直接用GDAL的命令行更省事。核心思路是通过/vsicurl虚拟文件系统直接读取HDF或者云服务端点GDAL请求服务端的exportImage子路径把结果当作一个本地文件来处理。示意命令如下gdal_translate -of GTiff -projwin 115 40 122 30 -tr 10 10 -r near \ /vsicurl/https://服务地址/ImageServer/exportImage?bbox115,30,122,40bboxSR4326size7000,7000formattifffimage \ lc_2020_tile.tif命令里的bbox是按WGS84经纬度写的四个角点格式是左、下、右、上size如果不给服务端会按范围自动计算-r near对应最近邻重采样。需要注意不同版本的GDAL对ArcGIS影像服务的处理能力有差异如果直接用/vsicurl读不了可以考虑先用QGIS导出一块测试数据确认流程通顺后再批量操作。我实际测试下来脚本化导出的优势不仅仅是省时间更重要的是参数一致性。手动一次一次点导出很容易某次忘记改重采样方法或者NoData设置而脚本每次请求都是相同参数这个在后期数据处理时能省很多麻烦。3.3 拼接、裁剪与投影转换分块导出之后第一步是拼接。这里有个经验不要一上来就调用“镶嵌至新栅格”这种重型工具先尝试用gdalbuildvrt建立虚拟栅格再转一次tif。好处是内存占用小速度快拼接时也不会因为边缘重叠产生条痕。gdalbuildvrt lc_all.vrt tile_0.tif tile_1.tif tile_2.tif gdal_translate -co COMPRESSLZW -co BIGTIFFYES lc_all.vrt lc_2020_study.tifLZW压缩对分类栅格非常友好因为类别值重复度高压缩比很大。BIGTIFFYES是为了防止输出超过4GB时报错。如果你研究区特别大文件超过几十GB可能需要考虑分带处理不要强行拼一张全球图。投影转换建议放在拼接完成之后做。因为原始服务的网格是Web Mercator3857而多数应用场景需要转成UTM或其他投影。用gdalwarp转投影时同样要注意重采样方法gdalwarp -t_srs EPSG:32650 -r near -co COMPRESSLZW lc_2020_study.tif lc_2020_utm50n.tif这里EPSG:32650是WGS84 UTM 50N适合中国中部大部分地区。大家根据自己的研究区所在分带选择合适的EPSG代码就好。3.4 类别面积统计有了分类栅格最常做的操作就是统计各类别面积。这里要强调如果栅格还在3857投影下直接统计的面积是错的这点非常重要下一章会专门讲。假设你已经转到了合适的等积投影统计面积就很简单了。在ArcGIS Pro里可以直接用“Tabulate Area”工具输入栅格表里会自动输出每个类别码对应的像元数量和面积。在QGIS里可以打开属性表用栅格唯一值统计的方式或者用Raster Layer Unique Values Report。如果用栅格计算器提取某个类别推荐写成这样# QGIS raster calculator expression lc_20201 5这会生成一个二值栅格真为1、假为0然后在属性里查看有效像元数乘以单个像元面积10米×10米100平方米就是该类别面积。更严谨的做法是先排除云和NoData再统计占比这部分我在4.2节继续讲解。4. 最容易翻车的投影与NoData细节4.1 3857投影下的面积陷阱标题里写了“全球10m”很多人以为数据是标准的等距等面积投影拿到手就直接统计面积结果一国面积算出来比理论上大了好几倍还以为是数据错了。其实问题出在Web Mercator投影。原始服务使用的是EPSG:3857也就是Web Mercator。这种投影的优点是在线地图显示时非常好用形状保持也还可以但面积变形非常夸张。它的变形规律是越靠近高纬度面积放大越严重。变形系数大约是1除以纬度的余弦值再平方公式为sec²(φ)。举个直观的例子在北纬60度附近面积会被放大约4倍中国东北地区大约在北纬45到53度之间面积也会被放大2到3倍。也就是说你在黑龙江用3857统计一片林地的面积得出的数字可能是真实面积的近3倍。所以做面积统计、密度计算、碳汇估算这些涉及数值结果的工作一定先把栅格投影转换成等积投影。中国的全国尺度项目可以用Albers等积投影区域项目用对应的UTM分带即可。制图展示时用3857没问题但数值统计必须换坐标系这是这套数据使用中最重要的一个坑。4.2 0和9的区别NoData与云分类栅格下载回来后看一眼属性表你会发现除了1到10之外还有一个0值。这个0通常代表NoData即背景或无效区域。它不是一个真实的地表类别统计面积时如果把它算进去会让总量偏大。第9类Clouds则完全不同它是有效像元但代表云遮挡。云下面的地表是什么分类算法无法判断所以它被单独标出来。如果你做植被覆盖分析云既不是植被也不能当无效数据删除——正确做法是先把云在分析中排除统计有效面积时用“非云且非NoData”的像元数作分母。我常用gdal_calc来做这个处理gdal_calc.py -A lc_2020_study.tif --outfilevalid_mask.tif \ --calc((A0)*(A!9))*1 --NoDataValue0生成的有效掩膜里1表示有效像元0表示云或NoData。后续所有统计都在这个掩膜约束下进行结果才真正可比。用QGIS栅格计算器也能写同样的表达式。还有一个经验下载完顺手算一下云的占比。如果某个区域的云占比超过10%说明这期数据的局部质量不太好做分析时要谨慎有机会的话用相邻年份的数据或者另一套产品插补。4.3 官方样式文件与配色建议分类栅格直接显示时通常是一堆灰阶颜色很难看。Esri官方为这套数据提供了配套的样式文件在Living Atlas条目页面或者ArcGIS Pro的符号系统里可以找到“Esri 2020 Land Cover”配色方案加载后就是大家常见的那套标准配色水体深蓝、森林深绿、草地浅绿、建成区红色、裸土棕色等。这套官方配色有几个好处第一各类别色差明显读图效率高第二社区分享图件时大家都用同一套颜色沟通成本低第三官方对第9类云用白色第8类冰雪用浅蓝白色两者在视觉上有区分。如果要自定义配色我建议注意两点。一是第9类云不要做成纯白否则和背景NoData撞色最好加一点灰色或浅黄。二是避免使用红绿这样的组合作为唯一区分考虑色盲用户可以用蓝橙棕等替代。制图输出时图例名称尽量写上中文简写加编码比如“2-树木”读者对照官方说明时不容易对不上号。5. 常见问题与排查技巧实录5.1 导出超时与任务排队在线影像服务不是实时出结果的超算它有自己的任务队列。我遇到过两三次这种情况提交了一个范围很大的导出请求等了好几分钟结果返回一个“Failed”或者“Timeout”。这不是数据坏了而是服务端限制了单次任务的像素数和计算时间。解决办法很简单缩小范围分块导出。我个人的经验是单次导出控制在5000×5000像素左右也就是2500万像元以内稳定性和成功率最高。如果你的研究区很大写一个网格切分脚本按每块约0.3到0.5个经纬度来切分逐块下载最后拼接。另外碰到服务繁忙时段比如工作日白天导出排队时间可能很长。如果条件允许把批量任务放到夜间跑成功率会高不少。这个跟任何大型在线渲染服务的道理是一样的。5.2 影像发黑、发白与类别值错乱下载完数据拖进软件里呈现出一片纯黑或者一片纯白这是新手最容易慌张的情况。先别急着删数据绝大多数是因为显示渲染方式不对。分类栅格是离散整数但软件默认会按“连续色带”去拉伸显示结果类别1和类别10之间的差值被硬拉出一条渐变看起来要么黑乎乎要么白花花。解决办法是在图层符号系统里把渲染方式从“拉伸”改成“唯一值”以Value字段作为分类依据再套上官方配色画面立刻正常。还有一个隐蔽问题如果你在导出时选了双线性或三次卷积重采样像元值会被插值成小数比如2.7类、5.2类。这种文件本身就已经坏了必须重新导出。记住一句话凡是分类数据任何步骤的重采样都只能选“最近邻”。5.3 面积统计偏大偏小怎么看如果你算出的某个类别面积和统计年鉴对不上先不要怀疑数据错按这个顺序排查第一检查当前栅格的坐标系。如果还是3857面积偏大是必然的尤其在高纬度地区第二检查是否包含NoData和云的像元分母有没有处理好第三检查面积单位换算。10米像元面积是100平方米很多人直接拿像元数当平方米用最后数值大了100倍第四如果以上都没问题再考虑分类误差。Esri这套数据是自动分类结果和实地调查相比肯定有误差建成区在低密度郊区尤其容易漏分这个只能靠结合其他数据源修正。我自己的习惯是先用这套数据做宏观格局分析再对关键区域用高分辨率影像抽样验证做到心里有数。5.4 常见问题速查表现象可能原因解决办法导出请求超时或失败范围过大、像素数超限分块导出单块控制在5000×5000以内影像显示全黑/全白默认拉伸渲染不适合分类数据改为“唯一值”渲染套用官方配色像元值出现小数如2.7重采样用了双线性或三次卷积重新导出全程使用最近邻重采样面积统计明显偏大还在Web Mercator投影下转UTM或其他等积投影后再统计水体面积与常识不符混淆了NoData(0)与水(1)明确排除NoData单独检查类别1的分布热带区域有云斑块正常现象云占比高统计时排除云或与邻近年份数据融合补齐拼接线明显或不连续各分块投影或NoData设置不一致统一导出参数使用gdalbuildvrt拼接6. 一些个人工作流与扩展玩法6.1 我的处理模板把前面所有内容串起来我现在处理一个省级范围的Esri 2020土地覆盖数据工作流基本固定为五步第一步确定研究区在WGS84下的bbox第二步按UTM分带或经纬度网格切块用脚本批量调用导出接口第三步用gdalbuildvrt拼接所有分块第四步转成目标等积投影同时用gdal_calc生成有效掩膜第五步做类别统计和制图。这套流程全程脚本化之后一个中等省份从数据下载到出图统计大约两个小时就能完成。如果全程在ArcGIS Pro里手动点可能大半天就耗进去了还容易在重采样或投影上出幺蛾子。所以我要再强调一次分类数据的影像处理规范和自动化比手快更重要。6.2 时间序列变化检测Esri除了发布2020年版还发布了对应的月度版本和相邻年份的版本。把这些版本放在同一个坐标系下做差值可以用来做大尺度变化检测。比如对比2017版和2020版差值为正的地方可能就是新增建成区或新造林地。要注意的是不同年份的分类结果存在噪声直接逐像元比较会产生大量“假变化”。我的建议是先做一个众数滤波把孤立的类别像元清理掉再去做变化检测。这里的滤波必须用众数滤波Majority Filter或者类似保留原始类别的算法绝不能用均值滤波——均值滤波会把类别算术平均产生一堆无效值。另外变化检测得到的是“可能变化区域”需要进一步用高分辨率影像验证不能直接当最终结论。6.3 与矢量边界结合的应用建议实际项目中这套数据和行政区划、流域边界、保护地边界等矢量数据结合是最常见的玩法。比如统计某个县域的林地面积或者评估某个保护区里农田占比的变化操作上都是先用矢量边界做栅格裁剪掩膜再用分区统计工具输出各类别面积表。这里有一个我踩过的坑掩膜提取之后矢量边界外的区域会变成NoData统计时如果不管会把边界外的像元也带进来。所以每次裁剪后先检查一下生成的栅格属性确认最小值、最大值、NoData情况是否符合预期再继续做统计。还可以顺手把“类别码、面积_m2、占比”三列整理成标准输出格式后续做报表和汇报时省很多事。最后再分享一个小技巧不管用什么工具导出导完第一时间用QGIS或ArcGIS Pro打开检查属性表里最小值是否为1、最大值是否为10或9如果出现0、255或者小数就不要进入下一步分析直接回去排查下载参数。这套数据本身质量很稳大部分问题都出在下载和处理环节的参数设置上。养成这个“拿到数据先体检”的习惯能帮你避开绝大多数坑。