Copernicus LC100全球土地覆盖数据：从下载、投影校正到面积统计实战

前阵子给一个区域生态评估项目挑基础数据客户张口就要“100 m分辨率、全球覆盖、最好免费可用”。市面上的全球土地覆盖产品其实不少但把哥白尼计划的 Copernicus 土地利用数据摊开之后我基本就锁死了用它做主底图。这篇文章不是官网文档的翻译而是我从注册、下载、打开栅格到最终算出统计表格的完整记录顺带把新手最容易踩的投影、文件组织和精度理解上的坑都摆到台面上。正在挑全球尺度土地底图的朋友可以直接把这篇当操盘笔记来用。1. 先搞清楚100 m分辨率的全球尺度Copernicus拿什么传感器和算法来做1.1 数据源头是PROBA-V卫星不是你想的Sentinel这套数据全称叫Copernicus Global Land Cover / Land Use行业里习惯简称CGLS-LC100或直接叫 LC100。它属于欧盟 Copernicus 计划陆地服务板块中的全球产品核心数据源是PROBA-V卫星一颗 2013 年发射、2020 年退役的植被监测小卫星。为什么要单拎出来强调传感器因为我见过不少新手拿到数据后误以为它是 Sentinel-2 直接生产的 10 m 产品做小尺度分析时对精度期待落空。PROBA-V 的原始观测分辨率其实是 300 m附近通过多时相合成和超采样重建成 100 m 的全球产品。换句话说100 m 是它的发布格网分辨率不等于原始影像分辨率这一点直接关系到后面精度评估的预期管理。算法上官方采用了一套基于 FAO LCCS 分类体系的决策树结合全天候时序特征、物候参数和辅助地形数据把地表覆盖状态映射成离散类型和连续覆盖度。每隔一段时间还会发布质量层quality layer这是其他很多免费数据集不提供的配置。1.2 为什么在众多全球土地数据里100 m这个尺度反而“最香”对比一下市面上能免费拿到手的全球土地覆盖产品MODIS MCD12Q1 是 500 mGlobCover 是 300 mESA WorldCover 是 10 mCGLS-LC100 是 100 m。很多人只看分辨率觉得 10 m 当然碾压 100 m但放到真正做区域级以上分析时10 m 的数据量和处理成本会让你的笔记本直接跪。我实测过ESA WorldCover 全球一幅栅格在 10 m 分辨率下是几十 GB 级别做一次全国范围的分区统计光 IO 就要跑很久。CGLS-LC100 的全球离散分类图层压缩后通常只有几个 GB 到十几个 GB普通台式机配合按 tile 处理完全可以跑下来。100 m 这个尺度刚好卡在“能分辨县级土地利用格局”和“处理成本可控”的平衡点上对宏观生态规划、碳汇估算、大流域分析来说性价比非常高。而且它时间上给出 2015 到 2020 年的年度产品可以做逐年变化分析ESA WorldCover 目前只提供 2020 和 2021 两个时相。两者各有适用场景但谈“逐年长时序全球变化”LC100 是最顺手的选择。1.3 版本演进别下到旧版还不自知目前主流通用版本是V3.0.1提供2016—2020 年年度合成更早还有 V2.0.2提供 2015—2019 年。PROBA-V 退役之后这一系列产品线就不再滚动更新后续的全球覆盖产品交接给了基于 Sentinel-1/2 的 Copernicus 产品线。所以如果你做 2021 年以后的逐年分析光指望 LC100 是不够的得考虑与其他产品衔接。下载时一定看清版本号。曾有人在群里问“为什么我跑的 2019 年和 2020 年结果完全一样”排查半天发现他从不同入口下到了同一版衍生产品。版本和年份两列要对齐这是最基础但最容易翻车的检查项。2. 打开压缩包别懵离散分类、连续覆盖层和质量层的完整清单2.1 离散分类23类的LCCS体系比IGBP细致在哪LC100 的离散分类图层按 FAO LCCS 体系组织主要类别有二十多个核心覆盖类型和代码大致如下代码含义代码含义0未知/无数据11灌木林1闭合常绿针叶林12灌木林稀疏2闭合落叶针叶林13草地3闭合常绿阔叶林14有稀疏树木的草地4闭合落叶阔叶林15有稀疏灌木的草地5闭合混交林16湿季淹没的草本植被6开放常绿针叶林17稀疏植被7开放落叶针叶林18裸地8开放常绿阔叶林19建成区9开放落叶阔叶林20水体10开放混交林21永久积雪/冰跟 MODIS 的 IGBP 分类相比LC100 把“森林”这一大类拆得很细常绿/落叶、针叶/阔叶、开放/闭合是六个维度两两组合后代码一直排到 10。做林业、生态模型的人会很喜欢这种结构因为针叶林和阔叶林的碳参数、蒸散参数差异巨大直接按 IGBP 合并成“树林”会损失太多信息。这里有个实操注意点开放/闭合的划分阈值是树冠覆盖度 15% 和 65%。闭合林指树冠覆盖度大于 65%开放林在 15%—65% 之间低于 15% 则归入有稀疏树木的草地或灌木地。如果你做森林定义标准不同例如国内一些标准按 20% 郁闭度划线就需要用连续树覆盖层自己重新评估不能直接拿离散层一刀切。2.2 连续覆盖层做森林覆盖度、植被退化分析要用它除了离散分类LC100 还有个经常被忽略的宝藏——连续覆盖层Continuous Coverage Layers。每个像元保存 0—100 的百分比数值包括树木覆盖度Tree Cover Fraction灌木覆盖度Shrub Cover Fraction林下草本覆盖度Understory Herbaceous Cover全部植被覆盖度Total Vegetated Cover水体覆盖度Water Cover为什么要提供连续层因为离散分类把复杂连续的地表强行切成几类在过渡带会出现明显的“椒盐噪声”和边界锯齿。连续层的价值在于你可以自定义阈值例如把“树覆盖 30%”且“小于 60%”的区域定义为退化林地这在离散分类里做不到。做森林恢复潜力评估、草原退化监测、城市绿量核算时我一般优先拉连续层而不是只看离散层。打个比方离散分类是给你一张点菜菜单看着是宫保鸡丁还是麻婆豆腐连续层则是告诉你这盘菜里鸡肉占多少、花生占多少、辣椒占多少。很多宏观分析真正需要的是后者。2.3 文件命名和质量层看清后缀再决定要不要重新下载官方发布的文件命名通常包含产品名、版本、图层类型和年份例如离散分类图层会带有 DiscreteClassification 字样连续层会带 TreeCover 等说明。下完压缩包后建议先列出文件清单再对照官方产品文档逐条核验。另外每个年份都配了一个质量层Quality Layer记录了每个像元在分类时的可靠性等级包括“有效”“边缘”“无效”“覆盖水体”等状态。做严谨出图或统计时我建议把质量层中非有效像元mask掉。尤其是沿海潮汐带和高纬度积雪区不注意质量层统计结果里会把海冰和泥滩算进林地面积误差能到几个百分点。3. 拿数据的完整路径注册、账号权限、数据集入口和常见下载链路3.1 官方入口和账号注册LC100 的下载入口主要有两个一个是Copernicus Global Land Service 官网land.copernicus.eu/global另一个是Climate Data Storecds.climate.copernicus.eu。第一个走产品专题页适合按年份、按瓦片手动下载第二个适合用 API 脚本批量拉取也是我实际用得最多的入口。注册账号是必须的而且需要注意官网和 CDS 的账号体系有时并不完全通用需要分别注册。注册过程就是邮箱验证加简单信息填写不涉及任何费用。Copernicus 数据采用开放许可允许自由使用包括商业用途但要求在成果里注明数据来源和版本号。这个在学术论文和商业报告里都要体现别偷懒。3.2 按瓦片还是按全国范围我的建议是分块下载第一次下载的人总想“直接下全球一整份”结果发现一个文件动不动好几个 GB网络稍慢就中断还要重新排队。LC100 官方提供了按图幅tile切分的下载方式瓦片命名基于经纬度网格例如东经 20 度、北纬 0 度附近的图幅会命名为类似 E020N00 的形式。我的建议是如果你只要某个区域先算出经纬度范围再把覆盖该范围的最小瓦片集合列出来逐个下载。例如一个中等尺度省份大约涉及 6—12 个瓦片单个瓦片在几十 MB 到几百 MB 之间浏览器直接下载也没压力。下完之后用 GDAL 的gdal_merge.py拼接几秒钟就能合成一张完整区域图。3.3 使用 CDS API 批量下载CDS 还提供 API 下载方式适合要跨多年、多图层批量取数的场景。配置好 API Key 后核心请求类似import cdsapi c cdsapi.Client() c.retrieve( satellite-land-cover, { variable: all, year: [2019, 2020], version: v3.0.1, format: zip, }, land_cover_2019_2020.zip, )注意 CDS 请求的参数名和取值会随平台改版变化运行前要去产品文档页复制最新的规范否则会 400 报错。我第一次就是因为没看更新日志用旧参数调了半个上午。4. 真正落地时的硬骨头文件格式、投影、重投影和面积换算4.1 打开前先看懂 GeoTIFF 元数据LC100 发布格式为 GeoTIFF坐标系默认是WGS84 经纬度EPSG:4326像素尺寸对应到赤道约 0.0009 度。用 GDAL 或 Python 打开前先跑一下栅格信息确认波段数、数据类型和文件组织方式gdalinfo PROBAV_LC100_Discrete_Classification_global_100m_2019.tif离散分类图层通常是单波段整型数值范围 0—23连续层是单波段浮点或整型数值 0—100。确认完之后再开始处理不要想当然。4.2 经纬度投影直接算面积错得离谱这是我在实际项目里碰到最多的问题有人拿着 WGS84 经纬度的栅格直接统计每个像元面积然后用“像元数量 × 100 m × 100 m”算总面积。问题在于经纬度格网在不同纬度上对应的地面面积差异巨大。同一个 0.0009 度 × 0.0009 度的像元在赤道附近面积接近 1 公顷到北纬 60 度附近只有约 0.5 公顷纬度越高面积越小。如果研究区横跨多个纬度直接相乘会产生系统性偏差。正确做法是先重投影到一个等积投影再做面积统计。全球尺度推荐Lambert Azimuthal Equal Area或Mollweide区域尺度用对应国家的投影参数。以欧洲区域常用的 ETRS-LAEA 为例栅格重投影到 EPSG:3035 后再做逐像元面积统计就靠谱得多import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src rasterio.open(discrete_2019.tif) dst_crs EPSG:3035 transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height, }) with rasterio.open(discrete_2019_LAEA.tif, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest, )重采样方法也要注意离散分类图层用nearest保证类别值不产生插值模糊连续覆盖图层用bilinear或cubic平滑过渡避免出现阶梯状边界。4.3 计算量大怎么办用 VRT、分块和压缩全球尺度栅格尽管已经比 10 m 产品小很多但直接全幅读进内存还是会卡。推荐做法是给多个瓦片建 VRT 虚拟拼接做分析时逐窗口读取gdalbuildvrt input.vrt E020N00.tif E020N01.tif E021N00.tif E021N01.tif然后利用rasterio的窗口读取或直接让专为地理栅格设计的库去处理分块。这里有个我从项目里总结的稳妥组合VRT 做索引GDAL 做重投影分块统计后写回表格。整个过程内存占用可以控制在 2 GB 以内。5. 它到底准不准官方验证、我的实测观察和绕不开的局限5.1 官方验证精度和我的交叉验证CGLS 官方对 V3.0.1 产品做了全球分层随机采样验证整体精度在 80% 到 90% 以上其中把 23 类合并成 9 大主类后部分区域的总体精度可以超过 90%。具体数字因年份和地区有波动但横向对比下来它属于当前免费全球土地覆盖数据里精度靠前的一档。我自己在东亚季风区、欧洲中部平原和非洲萨赫勒带三个样区做过交叉验证。东亚季风区水田、旱地交错离散分类在“湿季淹没草本”和“农田”之间混淆比较明显欧洲中部平原的农田和草地边界清晰精度最高萨赫勒带受稀疏植被影响裸地和草地经常互相判错。结论很简单全球产品的精度永远是分区域的用之前一定要做局地验证。5.2 已知的几个最容易踩到的坑第一时间分辨率是“年度合成”不是瞬时观测。某个年份的图层反映的是该年多时相统计优势结果短期的火灾、洪灾、砍伐事件可能被平滑掉。第二建成区在100 m格网下被严重低估。城市内的小街区、绿地和道路混合像元会让建成区面积偏小同时连续树木覆盖层在城市里又容易偏高。做城市相关研究建议配合更高分辨率数据做校正。第三高纬度湿地和冻土区可靠性较差。因为 PROBA-V 在高纬度地区观测次数少云和雪覆盖影响大质量层里很多像元被标记为低可靠度。处理环北极区域时务必先看质量层。5.3 和相邻产品怎么选如果你纠结 LC100 和 ESA WorldCover 10 m 怎么选我的判断标准很简单对比维度CGLS-LC100ESA WorldCover分辨率100 m10 m年份序列2015/2016—20202020、2021分类体系23类 LCCS11类处理成本低高连续覆盖层有无适用场景长时序宏观分析、模型模拟现状精细制图、局部热点做碳汇核算、水文模型、宏观生态区划LC100 是我首选做城市内部绿地识别或村级地块调查则必须上 10 m 级数据。两者不是替代关系而是互补关系。6. 从数据到结论一个可复现的土地覆盖统计工作流6.1 第一步裁剪到研究区抽样检查先把下载的瓦片按研究区范围裁剪裁剪后立刻做一件事叠加同一区域的高分影像或在线影像随机取 20—30 个点肉眼确认分类对不对。这个过程只需要十分钟却能把后面的方向性错误全部拦住。import geopandas as gpd import rasterio from rasterio.mask import mask shp gpd.read_file(study_area.shp) geom [shp.geometry.values[0]] with rasterio.open(discrete_2019_LAEA.tif) as src: out_image, out_transform mask(src, geom, cropTrue) out_meta src.meta.copy()输出后的out_image就是研究区的离散分类数组后续分析全部基于它。6.2 第二步重分类并统计面积把 23 类合并成业务需要的类别比如“森林”“草地/灌木”“农田”“建成区”“水体”“其他”逐类统计面积import numpy as np import pandas as pd data out_image[0].astype(np.uint8) reclass_map { **dict.fromkeys(range(1, 11), 1), # 各类森林 - 森林 **dict.fromkeys(range(11, 16), 2), # 灌丛、草地 **dict.fromkeys([16, 17, 18], 3), # 稀疏植被、裸地 19: 4, # 建成区 20: 5, # 水体 **dict.fromkeys([21, 0], 6), # 积雪、未知 } reclassified np.zeros_like(data) for old, new in reclass_map.items(): reclassified[data old] new pixel_area abs(out_transform[0] * out_transform[5]) # 等积投影下像元面积 area_per_class pd.Series( [(reclassified c).sum() * pixel_area / 1e6 for c in range(1, 7)], index[森林, 草地灌木, 稀疏裸地, 建成区, 水体, 其他] ) print(area_per_class)这段代码里的面积换算有个重要前提栅格已经重投影到等积投影坐标系此时每个像元面积才可以用分辨率直接相乘。6.3 第三步多时相变化分析和出图把 2016 和 2020 两年的离散数据分别重分类再做逐像元差分就可以得到变化矩阵。变化矩阵比单纯出变化图更有用它能告诉你哪一类转入哪一类比如“草地转森林”的面积有多少“林地转农田”又有多少这是土地利用转移矩阵的核心。出图时推荐用官方推荐的 23 色配色方案可以避免随意配色导致类别之间难以区分。QGIS 里可以直接加载官方 QML 样式导入后一键出图比自己去调色板省太多时间。6.4 我的一个习惯所有中间结果都写进处理日志最后分享一个我的个人习惯。这套数据处理链路里有版本号、投影、掩模范围、重分类映射表、面积换算系数五个关键参数任何一个变了统计结果都会有偏差。我会开一个简单的文本日志每次跑数前把参数记下来数据交付时连同日志一起发出去。这个动作保证了下个月重算时能精确复现同结果也是应对数据审核最有效的护身符。如果你手头正好在做区域土地利用变化、生态修复评估或者水文模型输入准备这套流程基本可以无缝移植过去。数据下载加处理一个下午能出首版结果后续精度验证再花一天——整体效率在免费全球土地覆盖数据里算非常划算了。