★ FEATURED ARTICLE 简介广东Landcover数据10米分辨率是面向GIS研究与规划应用的2020年土地覆盖资源包含栅格与矢量两种格式。栅格部分以tif为主配有tfw坐标配准、ovr金字塔及xml元数据便于在ArcGIS、QGIS中直接显示与处理矢量部分提供广东省界shp及其dbf、shx、sbn等附属文件可用于缓冲区、叠加及分区统计。数据同时提供UTM投影与WGS地理坐标系两个版本满足不同制图和空间分析需求。资源包共16个文件压缩后约163.29MB文件类型覆盖tif/shp/dbf/xml等常用GIS格式。目前已有277人学习下载。借助该数据用户可快速获取广东10米尺度的森林、农田、水体、建设用地等地表覆盖信息结合边界矢量开展地类面积统计、区域对比与生态环境分析为学术研究、城市规划和自然资源管理提供基础底图与分析支持。1. 为什么广东的10米Landcover数据值得自己再处理一遍做广东的项目需要一张讲得清来路、对得上日期的10m分辨率Landcover底图这在2020年前很尴尬GlobeLand30和CLCD都只有30m粤北山区的零星果园、珠三角的基塘、粤西的养殖水面根本分不出来国土变更调查的覆盖成果精度高但属于管理数据拿不到也不能随便用。2020年前后ESA WorldCover和Esri先后开源了10m全球地表覆盖产品广东终于有了可以自由下载的10m分辨率数据。但全球产品在广东不能直接用——果园被归成林地、养殖塘被归成水体、红树林被切成一堆碎像元面积统计出来和统计年鉴差一大截。这篇笔记按「选源→裁剪重映射→避坑→验证」的路线把一份能用于业务分析的广东省10m Landcover数据怎么从全球瓦片里抠出来讲清楚。2. 10m分辨率广东Landcover的数据源选型ESA WorldCover 与 Esri LC 的取舍2.1 先看表主流产品在广东的适配性真正能下载到10m分辨率的公开产品就那几个先把底牌摆出来。这张表是我在广东项目里常用到的选型依据覆盖了多数从业者会遇到的主流选项产品分辨率基准年份类别口径广东使用的注意点ESA WorldCover10m2020 / 202111类树、灌丛、草、农田、建成区、裸地、水、湿地、红树林、冰雪、苔藓类别稳定、有全球精度报告果园和灌丛经常混淆Esri Land Cover10m2017 起逐年发布10类农田、森林、草地、水体等更新快适合做年际变化养殖水塘容易被并进水体FROM-GLC10m / 30m2020偏生态制图的类别体系国内团队生产训练样本公开但后续更新慢GLC_FCS3030m1985 起逐年细分类别含旱地/水田等亚类对农地评估友好但分辨率只有30m破碎地块吃亏CLCD30m1990 起逐年9类做长时间序列分析首选类别少城市化过程看得清GlobeLand3030m2000 / 2010 / 202010类官方发布做宏观对比没问题小地块基本丢光选型结论很直白如果你的目标是「全省一张底图面积统计说得清」ESA WorldCover 2021版是默认选择它的类别定义和全球一致性最好如果你的目标里有「逐年变化监测」那就用 Esri 的年度产品它每年出一版10m省去自己配准时间序列的工作。注意任何全球产品在广东都不能直接用原因见下一节。2.2 类别口径差异广东的果园、基塘和红树林为什么总被错分全球产品的类别体系是按「全球普遍地物」设计的落到广东就会出现三类典型错分。第一类是果园广东的柑橘、荔枝、龙眼成片种植冠层结构和天然林在Sentinel-2影像上几乎一个样ESA和Esri都会把大片果园归到乔木覆盖或森林导致林地面积虚高、耕地园地面积少一大块。第二类是基塘和养殖塘珠三角桑基鱼塘、粤西高位池养殖在影像上是「水面上有细碎堤坝和棚屋」分类器常整片判成水体或湿地实际按国土口径这是农业设施用地两者面积含义完全不同。第三类是红树林10m分辨率下红树林冠层窄通常只有1到3个像元多数产品能标出大块红树林但破碎边缘漏判严重。理解这些错分之后选型就不能只看分辨率。我的做法是把 ESA 和 Esri 两套产品同时下载做一次类别交叉对比凡是两边结果不一致的区域就是本地物种最集中的地方后续单独处理。千万不要二选一就匆匆进入裁剪不然面积表的解释成本会非常高。2.3 10m分辨率在广东丘陵地区的斑块破碎化问题10m相对30m最直观的代价是图斑变碎。广东地形以山地丘陵为主粤北的梯田、珠三角的河网水渠、粤东的台地果园在10m像元下会被田埂、阴影、窄路切成大量细碎斑块。30m产品里一个完整果园可能是一个面10m产品里同一个果园可能有几十个碎面。这直接影响到后面制图和统计如果按最小图斑逐个统计广东全省的斑块数量会达到百万级制图符号化时屏幕上一片麻点。应对办法在生产上一般是引入「最小成图单元」控制常见做法是把面积小于0.1公顷大约10个10m像元的斑块并入其邻域占多数的类别。这个步骤要放在类别重映射之后、面积统计之前否则碎斑会被累加成可观的面积误差。具体平滑参数在第6章里给命令先记住这个概念即可。3. 从全球分幅到广东边界裁剪、投影与类别重映射流程3.1 边界准备为什么用省级行政边界而不是经纬度外接框这一步最容易翻车。很多人在ESA下载页里直接按广东省的经纬度范围框选瓦片出来一个矩形区域然后交给分类器统计。问题在于矩形框会把湖南南部、广西东部、福建西部一部分连带海南岛北侧的部分岛礁全圈进来统计面积时这些非广东区域全被计入。正确做法是先准备一份省级行政边界矢量再按边界裁剪。我一般去天地图或自然资源主管部门的公开数据服务里下载省级行政区划拿到一个全国SHP或GeoPackage再用ogr2ogr把广东省筛出来。注意不同来源的字段名不一样常见的是NAME或PACogr2ogr -f GPKG guangdong.gpkg china_admin.gpkg \ -sql SELECT * FROM admin WHERE NAME 广东省这条命令把全国面状边界里NAME等于「广东省」的记录单独导出为geopackage。逻辑说明ogr2ogr的-sql参数是直接在源数据上执行OGR SQL比先加载到GIS再导出更快也更可控NAME字段在大多数行政区划公开数据里都存在但如果你拿到的数据字段是PAC行政区划代码就改成WHERE PAC LIKE 44%。参数说明-f指定输出格式源数据格式由GDAL自动识别无需额外指定。边界文件要有明确的坐标系如果拿到的是经纬度EPSG:4326后面裁剪时GDAL会自动做坐标转换如果拿到的是CGCS2000投影坐标最好先统一到与影像一致的系统减少边界几何套合误差。3.2 坐标系选择度坐标下统计面积一定会翻车ESA WorldCover 原始瓦片是按UTM分带组织的广东恰好横跨UTM 49N和50N两个分带。如果不做任何处理直接拼接输出影像会在约东经114度附近出现一条接缝两边的平面坐标基准不一样如果直接以经纬度当作投影来统计面积结果会偏离物理真相因为你实际上用「度」当「米」在计算。这个坑我见过太多次项目汇报时面积和统计年鉴对不上回溯发现是坐标系统一没做。生产上常用做法是统一到Albers等积圆锥投影。为什么选Albers而不是Web Mercator因为面积统计要求等积Web Mercator在高纬度会被拉长广东虽然纬度不高但统计全省1500多万亩耕地时哪怕百分之几的变形也经不起推敲。Albers的两个标准纬线要覆盖广东南北范围我一般取21度和25度中央经线取114度gdalwarp -cutline guangdong.gpkg -crop_to_cutline -dstnodata 0 \ -t_srs projaea lat_121 lat_225 lat_00 lon_0114 datumWGS84 unitsm no_defs \ guangdong_esa.vrt guangdong_landcover_10m_aea.tif参数说明-cutline指向3.1生成的边界文件-crop_to_cutline表示让输出栅格范围严格贴合边界外接矩形外的裁剪线边界外的像元统一写0-dstnodata 0把0值声明为无效值避免后续统计把边界外大面积的0当成裸地-t_srs后面的proj字符串就是Albers等积投影定义。逻辑说明gdalwarp会把原始UTM瓦片重投影到目标坐标系同时完成裁剪文件输出为GeoTIFF。如果源瓦片数量多建议先按3.4节的VRT拼接再执行这条命令否则gdalwarp会逐个处理上百个文件效率极低。分两步做也可以先重投影成全省阿拉斯投影影像再用cutline裁剪。一步做的风险是当边界线和影像边缘在几何上高度重合时个别版本的GDAL会出现边界像元抖动但实际影响很小我习惯一步完成省一半磁盘IO。3.3 类别重映射把全球11类合成业务用的6类拿到广东范围的10m原始分类栅格后下一步是类别重映射。ESA WorldCover的类别代码是固定的常用的代码对应关系如下10乔木、20灌木、30草地、40耕地、50建成区、60裸地、80水体、90草本湿地、95红树林、70冰雪、100苔藓地衣。广东没有冰雪和苔藓裸地集中在湛江一带的火山台地。如果不合并直接拿11类做统计汇报时解释成本极高所以我一般合并成6个一级类。用Python做重映射比用QGIS的栅格计算器更可控尤其是当你需要反复调试合并规则时。下面这段脚本我直接贴在项目库里复用import rasterio import numpy as np src rasterio.open(guangdong_landcover_10m_aea.tif) data src.read(1) reclass np.zeros_like(data, dtypenp.uint8) # ESA WorldCover v200 官方类别代码映射到业务大类 reclass[np.isin(data, [10, 20])] 1 # 乔木灌木 → 林地 reclass[data 30] 2 # 草地 reclass[data 40] 3 # 耕地 reclass[data 50] 4 # 建设用地 reclass[data 60] 5 # 裸地 reclass[np.isin(data, [80, 90, 95])] 6 # 水体湿地红树林 → 水域湿地 profile src.profile.copy() profile.update(dtypeuint8, nodata0, compresslzw) with rasterio.open(guangdong_landcover_6class.tif, w, **profile) as dst: dst.write(reclass, 1)逻辑说明np.isin把数组里所有匹配类别代码的像元统一赋值比逐个写data 10再相加更清晰也避免嵌套条件把代码写成一团。原影像的nodata在读取时会被投影为0而重映射数组初始化为0因此新文件里nodata仍是0不会和任何类别混淆。参数说明林地合并了乔木和灌木这是按广东林业口径处理的——广东的灌木林比如粤东的桃金娘灌丛面积不小单独划出来业务上意义不大水域湿地合并了水体、草本湿地和红树林如果你的项目关注红树林碳汇必须把95单独保留成一个类别别合并掉。重映射前建议先看一次直方图确认原始数据的类别代码分布符合预期再跑脚本。直接跳过这一步的人往往会在最后发现某个类别的像元数明显异常追溯半天才发现是代码理解错了。3.4 拼接顺序先VRT后裁剪别一次性解压上百个分幅从ESA官网按经纬度框选广东范围通常会下载几十到上百个分幅瓦片每个瓦片是一个独立GeoTIFF。最差的做法是把它们全部解压后合并成一个巨型栅格再裁剪广东。广东范围虽然不大但全球产品分幅切入一组瓦片叠加起来动辄几十GB先拼接再裁剪会占掉可观磁盘空间。正确做法是先建VRT虚拟目录再让gdalwarp只处理广东边界相关的那些瓦片gdalbuildvrt guangdong_esa.vrt ESA_WorldCover_10m/*.tif gdalwarp -cutline guangdong.gpkg -crop_to_cutline -dstnodata 0 \ -t_srs projaea lat_121 lat_225 lat_00 lon_0114 datumWGS84 unitsm no_defs \ -wo NUM_THREADSALL_CPUS -co COMPRESSDEFLATE -co TILEDYES \ guangdong_esa.vrt guangdong_landcover_10m_aea.tif逻辑说明gdalbuildvrt生成的.vrt文件本身不复制任何像元数据只是记录了每个瓦片的地理位置和路径gdalwarp读取VRT时只访问那些与目标范围相交的瓦片内存和磁盘占用都小一个量级。裁剪和重投影在同一个命令里完成避免中间产物占用空间。参数说明-wo NUM_THREADSALL_CPUS让重投影过程多线程并行处理百张瓦片时速度差距明显-co COMPRESSDEFLATE和-co TILEDYES分别是压缩和瓦片化对后续反复读取有利。如果机器内存吃紧可以把线程数固定成4而不是ALL_CPUS不然GDAL会为每个线程开辟缓存块内存不足时反而变慢。4. 避坑清单广东10m Landcover最容易翻车的四个问题4.1 时相不统一2020版和2021版在同一地块上的类别边界对不上现象同一片潮汕平原用ESA WorldCover 2020版看是连片农田换成2021版再看变成农田和草地交错边界犬牙交错完全谈不上一致性。原因两个版本的分类训练样本和输入影像时相不同分类器对同一地物的判定阈值有差异。广东省云雨多影像有效观测次数少两版产品落在关键物候期的影像可能差了好几个月水稻田和收割后的裸地就容易摇摆。解决选一个版本作为主数据不要混用。做年际对比时要么都用同一产品的连续年份要么至少保证对比文件和历史文件来自同一套分类框架。混用ESA和Esri做变化监测是常见的坑因为它们的类别定义本身就不一致统计出来的「变化」有一半是口径差异。4.2 云掩膜和阴影让沿海大片区域掉进空洞或错判现象湛江雷州半岛沿海一带分类结果里出现大量空白区域或成片的裸地和实际的红树林、养殖塘地貌完全对不上像被人用橡皮擦掉一块。原因广东春季回南天、夏秋台风季的云覆盖影响了Sentinel-2有效影像数量全球产品在生成时对这些区域做了云掩膜掩膜区没有分类值或只能依赖质量较差的时相。雷州半岛是重灾区处理时容易让人误以为是投影问题。解决拿到数据后先用直方图统计nodata占比云掩膜区域超过5%就考虑换年份版本。对于小范围缺口可以用前一年或后一年产品对应区域做补值但要记录补值范围最终报告里列出来否则精度验证时会发现样本点落在空洞边缘混淆矩阵一团糟。4.3 行政边界和影像几何套合差几十米边界出现黑线白边现象按广东省界裁剪后沿省界的影像边缘出现一圈黑色或白色像元带局部乡镇级边界附近村庄地块和影像偏移肉眼可见放大看像拼贴画。原因行政边界数据来自国土变更调查或地图服务坐标系参考是CGCS2000ESA产品基于WGS84两者在理论上差别极小但实际生产中公开下载的行政边界经过多次数据加工几何精度可能数十米偏差。省界附近山区没有明显地物参照时偏差最明显。解决不要在原边界EPSG:4326状态直接裁剪先把边界和影像都转到同一投影再做。我的顺序是先做3.2节的Albers投影再在QGIS里叠加影像和边界检查几处明显拐点误差超过30米就换数据源。另外裁剪命令里加上-dstnodata 0并配合-overwrite可减少黑线出现的概率但根本解法还是边界源问题。4.4 红树林和滨海湿地被切成椒盐状斑块完整性尽失现象湛江红树林国家级自然保护区周边的红树林在10m分类结果里呈现大量孤立像元原本应该连片的林带像是被砂纸打磨过斑块破碎度高到无法计算有效面积。原因红树林冠层在影像上呈窄带状宽度常只有1到3个10m像元分类器对这种线性地物的空间一致性约束弱加之潮位影响同一片林子涨潮和退潮时的光谱完全不同模型会把边缘像元判成水体或湿地。解决分类后做一次连通性清洗把小于等于3个像元且与红树林类别相邻的碎斑并入红树林。更高精度的做法是在分类阶段加入红树林专题指数比如用GEE里现成的 mangrove index 计算后作为额外特征。如果只是想快速改进面积统计最简单的形态学开运算就能把孤立点抹掉代价是边界会圆滑一些但对省级统计完全可接受。5. 精度验证用分层随机抽样跑混淆矩阵别只看官方总精度5.1 为什么官方总体精度在广东不适用ESA和Esri随产品发布的精度报告总体精度数字看起来都不低但那是全球抽样验证的结果。官方验证样本的分布策略兼顾全球均匀性落到广东这样一个面积占全国约1.8%的省份有效样本可能只有几十个而且集中在珠三角城市群粤北山区、粤西滨海基本没覆盖。你在广东做业务拿着全球总精度去汇报等于用一个黑匣子给自己背书经不起同行追问。做本地精度验证是必须的验证要点不是笼统看总体精度而是看每一类别的用户精度和生产者精度这两个指标能告诉你哪类被高估、哪类被低估。5.2 用QGIS生成每个类别至少50个的随机验证点实践上最小样本量是每类50个项目要求高就80个。类别多时可以按比例分配但保证每个类别有底线数量否则混淆矩阵里某一行全是0Kappa计算也会失真。我一般在QGIS里完成抽样步骤很固定先把重映射后的6类栅格矢量化成面用Raster菜单里的Conversion → Polygonize然后按照类别字段把面拆开每一类单独成层最后用Vector菜单下的Research Tools → Random Points inside Polygons对每个类别层生成固定数量随机点。生成的点表导出CSV包含坐标、类别代码两个字段。接下来还需要人工目视判读——拿高分影像去逐个确认点的真实地类这一步偷不得懒直接拿分类结果当真实值是自欺欺人。5.3 混淆矩阵与Kappa一份可直接跑的Python脚本验证点表长这样real列是人工判读的真实类别pred列是产品分类结果。用sklearn三行代码出矩阵和报告import pandas as pd from sklearn.metrics import confusion_matrix, classification_report, cohen_kappa_score df pd.read_csv(validation_points.csv) labels [林地, 草地, 耕地, 建设用地, 裸地, 水域湿地] cm confusion_matrix(df[real], df[pred], labelsrange(6)) print(pd.DataFrame(cm, indexlabels, columnslabels)) print(classification_report(df[real], df[pred], target_nameslabels)) print(Kappa:, cohen_kappa_score(df[real], df[pred]))逻辑说明confusion_matrix的行是真实类别、列是产品类别对角线上是判对的数量。classification_report里precision对应用户精度recall对应生产者精度。用户精度低说明该类在图上画多了实际是草地却被分到林地生产者精度低说明该类被漏画了实际是林地只画出一部分。参数说明labels数组顺序和类别代码0到5对应如果你的类别代码不同先按第3.3节重映射结果核对。如果某个类别样本数为0classification_report会输出nan建议回到抽样步骤补充样本不要硬读结果。Kappa小于0.6说明一致性不达标项目验收时基本过不去需要回到类别合并规则或数据源选择上重新考虑。5.4 验证完要做的事把误差反馈到类别合并规则里混淆矩阵不只是拿来写报告的更应该用来倒推哪里改。我见过最典型的案例耕地生产者精度只有0.5但用户精度有0.9说明大量真耕地被分到了其他类——查交叉表发现这些错分样本大多落在果园和灌丛里也就是第2.2节说的果园被归林现象。这时只改重映射脚本解决不了因为10m影像本身对果园和天然林的区分度不够。我的做法是反馈到两个地方一是如果错分集中在某个特定区域比如湛江的菠萝地、粤北的茶园就单独圈出该区域用更高分辨率的影像局部修正二是如果错分比例整体偏高说明当前数据源不适合这个业务口径回到第2章换候选产品重新处理一轮而不是硬用一套数据交差。6. 进阶用法多期10m数据做土地覆盖变化和最小图斑平滑多期变化监测是10m数据相对30m数据最明显的红利。用Esri年度产品或ESA的2020/2021两期数据逐像元比较类别变化输出「林地转耕地」「耕地转建设用地」等转移矩阵这对珠三角城市扩张监测、粤北撂荒地识别都有直接价值。不过逐像元比较的原始结果会惨不忍睹一个像元因为云影或季节差异被误判为变化整张变化图上全是椒盐噪声。常规做法是先做多数滤波再消除细碎斑块。多数滤波用scipy一行解决消除孤立像元gdal_sieve做连通性清洗把小于阈值的连通斑块并入邻居from scipy.ndimage import median_filter import rasterio with rasterio.open(change_2020_2021.tif) as src: data src.read(1) filtered median_filter(data, size3, modenearest)gdal_sieve -st 10 -8 change_filtered.tif change_sieved.tifmedian_filter的size3表示3×3窗口内取中值modenearest保证边界像元不会被填充成0gdal_sieve的-st 10表示小于10个像元的连通斑块并入邻域占多数的类别-8指定8邻域连通规则。执行完这两步变化图斑从几十万个骤降到几千个图面才具备制图输出和面积统计的条件。验证方法是把变化结果和城市建设用地审批数据、土地利用变更调查做交叉叠加抽查20个变化图斑高分辨率影像目视确认变化真实性。我个人的习惯是任何一版分类结果交付前先跑一遍第3.3节的直方图检查类别逻辑再做一次第5章的分层抽样验证这两步能挡掉一多半后续返工。十年遥感数据处理下来最贵的教训就是省掉验证环节直接出图。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
