广东10米分辨率土地覆盖数据：选型、处理、统计与避坑指南

简介广东10米分辨率土地覆盖数据由ESRI于2020年发布面向地理信息系统分析与城市规划、生态研究等场景可清晰区分森林、农田、水体、建筑等土地覆盖类型。压缩包共16个文件大小约163.29兆字节核心为两种坐标系版本的栅格影像——UTM投影适合大比例尺制图与跨经度区域分析WGS地理坐标便于全球定位与多源数据叠加另有广东省边界矢量、投影、坐标配准及属性表等配套文件可直接用于缓冲区分析、重分类与制图输出。10米分辨率意味着每个像素代表实地10米乘10米的区域能较精细地刻画地物边界也为后续空间统计提供可靠输入。目前已有277人学习下载。该数据集兼顾栅格与矢量两类格式既支持大范围地表分类研究也支持行政区划边界内的空间统计与城市扩张分析适合需要广东省高分辨率土地利用底图的科研与规划人员快速取用在环境分析、资源管理和决策支持中具有实用价值。1. 广东10m Landcover数据一张图看清全省地类的家底如果你在规划院或高校遥感组待过一定遇到过这种尴尬领导要全省土地利用一张图手头只有30m的GlobeLand30珠三角的城中村和河涌糊成一片粤北山地阴影连水体都分不清。广东Landcover数据10m分辨率就是为这类需求准备的它用Sentinel-2等10m级卫星影像做分类给全省每个像元贴上地类标签水体、林地、草地、耕地、灌丛、建设用地都在一张栅格里对得上号。10m分辨率下村塘、香蕉林、连片宅基地大致可辨做生态本底、区县对比、变化监测比30m方案信息量大一截又比亚米级数据省掉大量标注成本。这篇按我实际做广东项目的流程讲产品怎么选、数据怎么拼、统计怎么做、坑在哪。2. 选对10m地类数据源从Esri、WorldCover到FROM-GLC的分类体系取舍2.1 三款主流10m产品在广东的选型对比现在能稳定拿到手的10m级公开Landcover产品我常接触的是三款Esri的10m全球土地覆盖基于Sentinel-2影像生产、ESA的WorldCover系列以及清华大学发布的FROM-GLC系列。三者都宣称10m分辨率、全球覆盖也都包含广东但各自的生产思路、类别定义和获取方式差异不小直接决定你后面统计和出图要补多少工作。产品基准年类别数广东适用性获取方式类型Esri 10m Land Cover2020另有逐年版本9类建成区、水体轮廓比较稳类别偏粗ArcGIS Living Atlas 在线按范围导出ESA WorldCover2020 / 20211112类有红树林类别对广东海岸带更友好ESA官网按10°分幅下载GeoTIFFFROM-GLC2017年起多期9类中国区域训练样本相对密官网按分幅/瓦片下载我一般优先用Esri 2020版做全省概览因为它的建成区和水体在珠三角表现不差而且按AOI导出时不需要自己拼太多图幅。做沿海红树林、滨海湿地相关课题时我会换成WorldCover它把红树林单独列为Mangroves类在湛江、汕头、深圳湾这些地方明显更贴实际。FROM-GLC的优势是持续更新和多期对比但下载和预处理对新手不太友好如果你不是做长期时序分析第一版数据不建议从它入门。选型还有一个关键点这三款产品的“10m”都指原始栅格像元尺寸不是“10m能认清所有地物”。在广东这种地块破碎、地形起伏大的地方10m能保证100平方米以上的连续地类基本成片但3米宽的河涌、单棵大树、田埂间的菜地依然会混进邻类。把它当成“区域尺度可信、地块边界需谨慎”的底图而不是精确到地块的现状图后面能少很多翻车。2.2 分类体系拆开看广东九类地物如何映射Esri产品的9类编号是水、树木、草地、被淹没植被、耕地、灌丛、建成区、裸地、冰雪。广东除个别高山冬季短暂积雪外冰雪类基本用不上实际有效的是前8类。对业务口径来说这套体系有个好理解的口诀树木、草地、灌丛、被淹没植被都算生态空间耕地是农业空间建成区是城镇空间水是水域。但要注意被淹没植被这类常出现在鱼塘、滩涂、莲藕塘上直接并入水域或生态都要说明口径。WorldCover的分类更细一些常见类别包含林地、灌丛、草地、耕地、建成区、裸地、水体、草本湿地、红树林等。红树林单独成一类对广东沿海的生态评估很有价值但草本湿地和耕地有光谱重叠雨季的稻田在高分影像上看是水在分类图里可能落到湿地类。做面积统计前我会把这些易混类人工归并湿地并入生态空间稻田归耕地实在分不开的按“植被”合并处理避免台账里出现“既不是水也不是田”的尴尬类别。分类映射不能只看类名硬套要回到产品文档确认类别代码。我遇到过同事把Esri第4类被淹没植被当成湿地直接拿来算珠三角湿地面积结果把连片鱼塘都划进湿地面积翻了两倍。正确做法是先导出一小块已知地类区域把类别代码和影像对照一遍再定映射表。10m产品的类别边界本来就是概率输出映射口径不统一后续所有统计都会跟着偏。2.3 时相合成与季节为什么广东这种多雨省份要特别留意广东一年到头云雨多Sentinel-2单景影像超过一半被云覆盖产品生产方不会拿单期影像分类而是把一年内多次晴好观测合成为一个“合成影像”再分类。这意味着所谓“2020年”的类别实际可能由2019年冬天到2021年春天的晴好像元拼出来。同一个图斑里山体南坡可能是夏季像元北坡可能是冬季像元直接导致阴坡植被被分成草地、阳坡被分成林地。我看Esri 2020版“广东部分”时有个明显体感粤北连州、乳源一带的常绿阔叶林大片落在Trees类但同样山区的采伐迹地、火烧迹地因为合成时相偏冬被分进Grass甚至Bare ground。处理办法是做年度对比或变化检测时不要混用不同产品、不同年份的图层Esri 2020就配Esri 2021WorldCover 2021就配WorldCover 2021更稳妥的是先做类别重编码和投影统一再逐像元比较。另外多时相合成还会留下“时间残影”。我曾在珠江口一张10m分类图里看到孤立的水体像元沿着海岸线排成细线放大看是船舶、浪花、潮滩在不同月份来回切换分类标签。这类像元没有空间连续性在后续统计里应通过最小面积滤波或众数滤波处理。总之一句话分类图上的每一个像元都是一种“平均状态”不是某个具体日期的真实地表理解这一点才能正确处理广东这种多云多雨区域的数据。3. 下载、拼接、裁剪、投影把10m Landcover变成广东一张图3.1 获取数据按图幅还是按省边界先明确一个原则拿数据先按“产品分幅”拿再按“省界线”裁不要幻想一个文件覆盖全广东。WorldCover按10°×10°分幅组织广东大致落在N20E110、N20E120、N30E110这几幅交界附近沿海岛屿还要留意分幅边缘的缝隙Esri产品在Living Atlas里可以上传广东边界矢量作为AOI后端自动拼好返回但返回结果常是Web墨卡托坐标系面积会变形。FROM-GLC同样是分幅下载文件名里带经纬度信息。广东的大致范围是东经109.66°到117.3°、北纬20.19°到25.53°。下载时把外包矩形向外扩0.5°不然海岸线、南澳岛、担杆列岛这些边缘区域容易落在图幅外。下载完成后第一件事不是打开看颜色而是用gdalinfo确认原始坐标系、位深、NoData值。三个产品的位深和NoData不统一有的用0有的用255有的不写NoData这一步记错了后面的面积统计全错。gdalinfo ESA_WorldCover_10m_2021_N20E110.tif这条命令会输出栅格大小、投影、像元尺寸和NoData信息。拿到文件后我习惯顺手把原始CRS和NoData记到项目笔记里因为后面拼接、裁剪要反复用到。3.2 VRT拼接与省界裁剪gdal两步完成覆盖广东的多个分幅不能直接丢进GIS做分析先建VRT再裁剪。VRT是个虚拟栅格文件不复制像素数据只记录各分幅的位置和范围拼接速度极快。直接gdal_merge.py也可以但遇到大范围、多文件时容易吃满内存而且分幅之间的黑边、重叠区处理不如VRT灵活。gdalbuildvrt guangdong_raw.vrt ESA_WorldCover_10m_2021_N20E110.tif ESA_WorldCover_10m_2021_N20E120.tif gdalwarp -cutline guangdong_boundary.shp -crop_to_cutline -dstnodata 0 -ot Byte -of GTiff -co COMPRESSDEFLATE -co BIGTIFFYES guangdong_raw.vrt guangdong_landcover.tif逻辑说明第一条命令把两幅邻接的WorldCover分幅串成一个虚拟文件第二条命令用广东省界矢量裁出最终成果。参数里-crop_to_cutline表示沿省界的折线裁剪而不是简单的矩形裁切-dstnodata 0把裁切后省界外的区域统一写成0便于后续统计时屏蔽-ot Byte是因为分类图用无符号整型存储最省空间DEFLATE压缩能把10m栅格压到原大小的一半左右BIGTIFFYES防止超过4GB写失败。这里有个容易忽略的细节裁剪后栅格的边界不再是矩形NoData值会以0填充在省界外。第4章统计面积前必须显式把0排除否则“全省耕地面积”里会混进大片省界外背景值。另外如果省界shp和栅格投影不一致gdalwarp默认会先重投影矢量速度慢且边界可能会变形建议先统一成同一坐标系再裁剪。3.3 投影与NoData处理面积统计前的三个必改参数广东横跨UTM 49N和50N两个分带直接用原始UTM投影做全省面积统计东西两侧会有不小的长度变形。我一般统一转Albers等积投影中央经线取111°E两条标准纬线取18°N和28°N中央纬线取22°N。等积投影保证每个像元代表的实地面积一致是做面积统计、密度计算的前提。Web墨卡托EPSG:3857只适合做在线地图展示用它算面积纬度越高虚胖越严重珠三角的区县面积能比真实值多出10%以上。gdalwarp -t_srs projaea lat_118 lat_228 lat_022 lon_0111 datumWGS84 -r near -dstnodata 0 -overwrite guangdong_landcover.tif guangdong_aea.tif逻辑说明-t_srs直接传Albers投影参数不依赖Proj4库版本-r near保持分类值不插值这点和连续变量处理完全不同不能用bilinear或cubic否则会出现“1.5类”这种非法值-overwrite强制覆盖已有文件。NoData检查也要放在这步。gdalinfo guangdong_landcover.tif | grep -i nodata能快速看到NoData值但WorldCover的NoData有时是0Esri导出结果是255FROM-GLC某些版本干脆没写NoData。我处理这类栅格时统一加一步“NoData归一化”读进Python后用numpy把非法值全置0再写回GeoTIFF。这样后续脚本只认0一个口径省掉无穷无尽的边界判断。4. 用Python做类别统计、区县汇总与制图渲染4.1 类别面积统计先屏蔽NoData再看直方图拿到一张广东10m分类栅格大多数人的第一个需求是“各类别到底有多少面积”。最简单的做法是直接np.unique统计每个值出现次数但如果不先屏蔽NoData省界外的0值会被算成一类列面积时莫名其妙多出几万平方公里的“未知地类”。正确步骤是先读NoData生成有效掩膜再做直方图。import numpy as np import rasterio with rasterio.open(guangdong_aea.tif) as src: data src.read(1) transform src.transform nodata src.nodata if nodata is not None: valid data ! nodata else: valid data ! 0 px abs(transform[0]) py abs(transform[5]) pixel_area px * py classes, counts np.unique(data[valid], return_countsTrue) for cls, count in zip(classes, counts): area_km2 count * pixel_area / 1_000_000 print(f类别 {cls}: 面积 {area_km2:.2f} km2占比 {count / valid.sum() * 100:.2f}%)逻辑说明transform[0]和transform[5]分别是X方向、Y方向的像元分辨率Albers投影下单位是米valid掩膜把所有NoData剔除保证求和分母真实输出单位换算成平方千米便于和统计年鉴对照。参数上唯一要改的是nodata的取值不同产品已经归一化到0的话直接用data ! 0即可。4.2 多产品分类体系重编码统计前必须做的归一化如果你要同时处理Esri和WorldCover两套产品或者和第三次国土调查数据做对照第一件事是重编码否则两套栅格的类别编号鸡同鸭讲。重编码的原则是“业务类别少而稳”我把广东归并为水域、生态空间、耕地、建设用地、其他五个一级类再把原始类别按产品文档映射进来。reclass np.zeros_like(data) reclass[data 1] 10 # Water - 水域 reclass[data 2] 20 # Trees - 生态空间 reclass[data 3] 20 # Grass - 生态空间 reclass[data 4] 20 # Flooded vegetation - 生态空间也可单独保留 reclass[data 5] 30 # Crops - 耕地 reclass[data 6] 20 # Scrub - 生态空间 reclass[data 7] 40 # Built - 建设用地 reclass[data 8] 20 # Bare ground - 生态空间 reclass[data 9] 0 # Snow - 剔除 reclass[~valid] 0逻辑说明目标值10、20、30、40是自定义编号方便后续交叉表统计被淹没植被放生态空间还是水域要按研究目标改我建议单独保留一个编号别过早合并。参数上唯一要注意的是确保data是整型数组如果读进来是float先data.astype(np.uint8)再执行比较避免浮点比较出意外。重编码完成后把结果写成一个新tif命名里带上体系名比如guangdong_reclass_lv1.tif。多个项目共用同一套重编码脚本能少走很多弯路。我吃过一次亏先统计后重编码全套类别面积表都作废得回头重跑。4.3 区县Zonal统计rasterstats的categorical模式省级汇总做完业务上通常还要分到地市或区县。用ArcGIS的Zonal Histogram也行但脚本化处理更省事。rasterstats库的zonal_stats支持分类栅格只要把categorical参数打开它会返回每个几何对象里各类别像元数而不是默认的平均值。from rasterstats import zonal_stats stats zonal_stats( guangdong_counties.shp, guangdong_reclass_lv1.tif, categoricalTrue, nodata0, geojson_outTrue ) for feature in stats: props feature[properties] county props.get(name) counts props.get(categories, {}) print(county, counts)逻辑说明categoricalTrue让每个区县输出一个dict键是类别编号值是该类别像元数nodata0确保省界外背景不参与geojson_outTrue把统计结果直接附加到原矢量属性方便导出GeoJSON或转shp。注意先保证shp和tif都是Albers投影如果shp还留在WGS84经纬度rasterstats会逐几何重投影大范围统计时会明显变慢。拿到各区县各类别像元数后面积像元数×10010m分辨率下像元面积100m²再除以1e6换算成km²。我习惯把这一步封装成一个函数输入区县矢量和分类栅格输出DataFrame这样不同区县、不同年份一键复用。4.4 出图参数颜色表、图例与DPI统计完成后的最后一步是出图。10m分辨率栅格动辄几千万像元出图时切忌用ArcGIS默认随机配色。我给广东常用的一套配色是水域用深蓝44,123,182、林地用深绿34,139,34、草地用浅绿152,213,143、耕地用黄240,210,80、建设用地用红214,72,63、裸地用灰180,180,180。这套色表接近第三次国土调查的制图规范甲方看着不陌生。import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap, BoundaryNorm colors [#2c7bb6, #228b22, #98d58f, #f0d250, #d6483f, #b4b4b4] cmap ListedColormap(colors) norm BoundaryNorm([10, 20, 30, 40, 50, 60], cmap.N)逻辑说明BoundaryNorm把重编码后的10、20、30等类别映射到对应颜色边界值设成半开区间避免类别间颜色溢出。参数上要注意分类编号必须连续如果重编码后出现空缺编号BoundaryNorm会把相邻类别混色。出图时DPI设300以上比例尺按尺度来全省用1:150万区县用1:10万到1:5万图例里写清楚“数据来源、年份、分辨率10m、投影”这三项是再审不看图都认。5. 广东Landcover数据使用避坑5条踩出来的经验5.1 山地阴影区被误判为水体粤北最常见的坏点现象粤北乳源、连州、乐昌一带的山谷和背阴坡上分类图里出现大片连贯的Water类图斑沿着山脊线断断续续延伸和真正的河流、水库完全对不上。用高分影像放大看这些地方是茂密植被的阴影沟谷。原因10m合成影像的阴影像元整体辐射值低归一化水体指数在阴影里偏高分类器把“暗”当成了“湿”。解决用SRTM坡度数据叠加把坡度大于10°的像元强制排除出水体类再用最小图斑面积滤波把小于1500m²的孤立水体直接删除。我一般处理完会专门出两张图对比一张原始分类、一张坡度修正后确认山脊线不再串水。5.2 图幅拼接带条带差异不是数据坏是时相不一致现象全省拼接图里图幅接缝处出现明显的类别断层比如同一条山谷东边是林地、西边是草地而高分影像上看两者没有本质差异。原因不同分幅使用了不同时相的合成影像季节差异被分类器放大成了类别差异。解决拿到多幅数据后先不做全省拼接而是按分幅统计类别占比如果相邻分幅同一类别占比差超过3个百分点就先停一下。日常处理时我通常把分幅重叠区取众数或者只保留其中一幅覆盖研究区别等拼完再找断层。5.3 耕地与草地混分10m的极限与物候补救现象粤西湛江的旱地、粤东丘陵的坡耕地在分类图里经常落到Grass类而不是Crops类冬闲期的水田被分到Flooded vegetation或湿地。原因10m像元在破碎田块上混合像元多冬春季节耕地光谱和草地几乎没有区分度。解决面积统计或转移矩阵前叠加三调耕地边界或农业确权边界做掩膜把界线内外的类别重新投票如果只是做宏观趋势干脆把Crops和Grass合并成“绿色植被”不为一个分类细节卡住整个分析。和甲方汇报时务必说清“10m产品不能精确识别轮作地块”。5.4 投影混用导致偏移几十米先gdalinfo再动手现象把分类tif拖进QGIS和区县矢量边界叠在一起边界错位几十到上百米面积统计结果和统计年鉴对不上。原因tif原始是Web墨卡托矢量是CGCS2000经纬度软件动态投影只负责“显示正确”一旦做叠加计算就按各自原始坐标运算自然错开。解决所有数据动手前先gdalinfo查CRS统一用3.3节的方式转Albers或UTM不要依赖数据框的实时投影。这个坑最气人因为图上看不出来只有算面积时才暴露。5.5 不跑混淆矩阵就下结论精度是最后一道闸现象从分类图统计出的某区县耕地面积和统计年鉴差了30%以上项目评审时被专家当场质疑。原因全球产品的训练样本在广东分布稀疏粤北、粤西局部精度可能远低于官方公布的全球总体精度。解决分层随机抽样每个地级市至少布50个验证点全省500个点起用高分影像或实地照片判读类别生成混淆矩阵计算总体精度、各类别用户精度和生产者精度以及Kappa系数。我个人的红线是总体精度低于0.8、Kappa低于0.6时这张图只准做趋势分析不许出具体面积台账。把官方精度当自己项目精度的做法是拿职业信誉开玩笑。6. 进阶用两期10m Landcover做变化检测与精度自检6.1 两期数据变化检测的最小流程与滤波参数如果手里有Esri 2020和Esri 2021两期广东分类图最简单可靠的变化检测是逐像元类别比较。注意先把两期数据分别重编码到同一分类体系再执行比较否则类别编号错位会让“变化”多到不可信。import numpy as np import rasterio from scipy.ndimage import label with rasterio.open(guangdong_2020_lv1.tif) as a, rasterio.open(guangdong_2021_lv1.tif) as b: arr20 a.read(1) arr21 b.read(1) profile a.profile valid (arr20 ! 0) (arr21 ! 0) new_built (arr21 40) (arr20 ! 40) valid new_built new_built.astype(uint8) structure np.ones((3, 3), dtypeint) labeled, num_features label(new_built, structurestructure) sizes np.bincount(labeled.ravel()) min_pixels 100 # 10m分辨率下约1公顷 big_change np.isin(labeled, np.where(sizes min_pixels)[0])逻辑说明arr21是2021年类别40代表建设用地new_built标记“2020年不是建设用地、2021年变成建设用地”的像元label函数按8邻域连通性把这些像元聚成图斑sizes统计每个图斑的像元数min_pixels100对应10m分辨率下的1公顷最小成图单元。这一步能滤掉大量椒盐式伪变化。业务上如果只关心建设占用耕地可以把条件改成arr20是30且arr21是40再加一层条件逻辑完全一样。6.2 自检习惯验证点、混淆矩阵与决策红线变化检测做完我必做的一件事是抽检。两期分类图本身都有误差逐像元相减会把两期的误差叠加起来“草地变林地”“林地变草地”这类伪变化经常混在结果里。我的做法是每年固定生成500个随机验证点存成GeoJSON并逐年追加用高分影像逐点判读输出混淆矩阵。确认2020年、2021年两期数据在“建设用地”和“耕地”这两个类上的用户精度都超过80%才敢把变化图斑写进报告。这套流程看似多花半天实则省掉无数解释成本。以前有一次我只看了直方图就下结论把草地到林地的伪变化当成造林成果递上去被甲方拿着高分影像逐块追问场面相当下不来台。后来凡是出广东Landcover相关的图我都先跑混淆矩阵再出图这个习惯成了我做遥感底图类项目最后的闸门。希望帮到你。本文还有配套的精品资源点击获取