云南10m土地覆盖数据处理与面积统计：从解压到NDVI自检全攻略

简介2019年10m精度云南省土地覆盖土地利用数据包基于哨兵影像与深度学习制作分类涵盖耕地、林地、草地、灌木、湿地、水体、建筑用地、裸地及雪/冰等十类。资源已完成坐标系转换与行政边界裁剪统一为WGS84地理坐标系并按云南省各地市分别输出便于研究或应用中直接调用。压缩包共112个文件包含16个TIF栅格文件及配套的tfw定位文件、xml元数据、dbf/cpg属性表、png预览图与xlsx表格合计约115.71MB结构清晰适合城乡规划、生态环境及GIS学习者作为底图数据或训练样本。目前已有327人学习下载可满足对省级尺度高精度土地覆盖数据的快速获取需求。1. 拿到云南10m土地覆盖数据之后先别急着出图2019年10m精度云南省土地覆盖土地利用.rar这类压缩包近年在GIS从业者手里流转得很多。它本质上是全球10米分辨率土地覆盖产品按云南省边界裁剪出来的本地化数据常见来源是欧空局或ESRI基于Sentinel-2制作的2019年分类结果。很多人下载完第一件事就是拖进ArcGIS或QGIS拉伸显示结果要么颜色乱成调色盘要么把水体显示成林地然后开始怀疑数据坏了。这里想先立一个判断10m精度指的是像元大小为10米而不是每个地类都精确到10米。分类结果受训练样本、云覆盖和地形阴影影响局部误差可能超过一个像元。它适合做宏观概览、面积统计、变化趋势不适合当作地块级取证依据。我一般拿到这类数据会先花半小时做四件事查分类体系、验投影、算有效覆盖、做一张能稳定出图的样式。这套流程能避掉后续大部分坑。下面从一个rar包的落地路径展开从数据底细讲到面积统计再落到云南地形带来的各种翻车点最后给一个用NDVI做自检的土办法。新手可以照着命令走熟手可以直接跳去第5章看踩坑清单。2. 数据底细与选型2019年10m土地覆盖数据到底是什么2.1 分类体系与栅格文件组织先看压缩包内部的结构。常见的10米土地覆盖产品是一个GeoTIFF加一个样式文件.lyr或.qml有的还带一个CSV格式的分类描述。云南省的文件名里通常有Yunnan、2019、LC等字样栅格值从0到11或从1到12不同产品对应关系完全不同。这里要特别警惕ESRI产品中1是水体、2是树木、3是草地但ESA WorldCover中10是树木、20是灌木。拿ESRI的图例去渲染ESA的产品整张图都会是错的。打开文件后用gdalinfo看波段数和像素类型。多数10米产品是单波段8bit分类栅格一个像素只占一个字节文件大小可控。云南面积约39万平方公里10米分辨率下理论像元数接近40亿但压缩后tif通常在1到2GB因为分类图有大量游程压缩。如果看到的是一个三波段RGB预览图而不是单波段分类图说明压缩包里带的是可视化版本真正的分类图层在子目录里。过去我懒得翻子目录直接用RGB预览图做统计结果全部作废这个亏吃过一次就不会再吃。2.2 10m与30m、250m的真实权衡选分辨率不能只看越细越好。30米Landsat数据历史长能回溯到上世纪80年代适合做长时间序列250米MODIS数据时间频次高适合做物候而10米Sentinel-2数据从2015年后才稳定获取2019年的产品已经是比较成熟的批次。对云南这种山高谷深、地块破碎的地形来说30米像元跨过一条河谷时混入的类别可能超过一半10米像元虽然也混但至少能把梯田、林窗、村寨边界大致分开。但10米也带来噪声问题。山区阴影、单次观测云遮挡、同物异谱都会让分类结果出现大量椒盐噪点。做面积统计时直接用原始10米栅格统计会发现灌木面积忽大忽小原因在于训练样本里灌木与草地本身就难分。实战中我常做一步3x3众数滤波把孤立像元去掉后再统计结果更接近业务口径。这一步属于常见做法不算严谨科学但确实有效。顺便说一句千万别用均值滤波处理分类栅格否则会出现0.6之类的地类值后期查都查不出来。2.3 云覆盖与有效像元云南数据的隐形参数云南地处低纬高原干季晴天多雨季云量常年在60%以上。2019年的合成产品理论上用到全年多时相合成但部分月份数据源缺口不小。判断数据能不能用首先要看有效像元比例看文件大小没用。有效像元指分类值落在合法范围内的像元如果某个类别像元占比低于0.1%且恰好集中在高海拔山区那多半是云阴影或雪被误分。这类区域在成图时建议标注为“数据不确定区”而不是强行修图。我见过不少项目把云南西北部的雪冰类像元直接归零理由是当地冬季确实有雪但分类产品中的雪冰面积会随年份浮动。如果你做土地利用变化而基准年恰好选了多雪年份森林减少面积会被严重高估。处理办法是统计前把坡度大于35度且海拔高于4000米的像元单独拎出来看这些地方的分类误差自带Buff不要直接并进森林或裸地。云南这种地形任何遥感分类产品都需要结合DEM做二次判读这是项目开始前就应该写进技术方案里的。3. 从rar到可用TIF解压、校验与转坐标系实操3.1 解压与完整性校验用命令而不是双击拿到一个rar包第一步不是双击解压到桌面而是用命令行做一次完整校验。rar文件在网盘里传几轮后文件头损坏概率不低。Windows上可以用WinRAR自带的rar tLinux下用unrar t。下面以Linux环境为例因为这后面处理GDAL大多在Linux或WSL下跑。# 校验压缩包完整性不释放文件 unrar t 2019年10m精度云南省土地覆盖土地利用.rar # 查看包内文件列表避免解压出奇怪的路径 unrar lb 2019年10m精度云南省土地覆盖土地利用.rar # 解压到data目录保留原路径结构 unrar x 2019年10m精度云南省土地覆盖土地利用.rar ./data/t参数是test只校验CRC不写盘lb列出文件路径x是完整解压并保留目录结构。逻辑很简单先测后解能省掉后期缺文件的麻烦。解压后立刻用ls -lh查看每个tif的大小若某个tif只有几KB那多半是空文件或损坏文件趁早重新下载。另外注意RAR包里的文件名如果包含中文或空格后续脚本处理前统一重命名成y2019_yunnan_lc.tif这类格式能避免大量编码问题。这不是玄学是Python读取中文路径时在Windows控制台下经常乱码。3.2 检查投影与波段gdalinfo是照妖镜解压出的tif很可能落在Web墨卡托坐标系下EPSG:3857因为很多在线发布源直接吐瓦片坐标。云南东西跨度大Web墨卡托虽然低纬变形不大但面积统计需要等积投影。先用gdalinfo看元数据再决定是否投影。# 查看栅格基本结构 gdalinfo ./data/yunnan_lc2019.tif # 如果栅格太大只看关键信息 gdalinfo -stats ./data/yunnan_lc2019.tif | head -60看输出的关键三行Size is后面的行列数Coordinate System is后面的EPSG代码NoData Value后面的值可能是255也可能是0。最容易翻车的就是NoData。很多产品把海洋或境外区域设为0而0又恰好是合法分类值水体可能为0或1导致裁剪边界出现一圈奇怪的黑边或白边。如果发现NoData和合法值重叠必须在下游处理里重新定义NoData。处理办法是gdal_translate -a_nodata 255把255设为新的NoData因为255在大多数8bit分类体系里不是合法地类。3.3 重投影与裁剪两个坑一次填平云南常用的投影是Albers等积圆锥投影或UTM 47NEPSG:32647。如果只是做省级统计我一般直接用UTM 47N它覆盖云南大部分区域像元面积接近常数统计方便。下面的命令先把数据从3857重投影到32647再用省界矢量裁剪。注意这里用-r near而不是cubic因为分类栅格不允许插值三次卷积会把地类值变成小数。# 重投影到UTM 47N最近邻重采样 gdalwarp -overwrite \ -t_srs EPSG:32647 \ -r near \ -dstnodata 255 \ ./data/yunnan_lc2019.tif ./data/yunnan_lc2019_utm47.tif # 用云南省界矢量裁剪同时保持投影不变 gdalwarp -overwrite \ -cutline ./data/yunnan_boundary.shp \ -crop_to_cutline \ -dstnodata 255 \ -of GTiff \ ./data/yunnan_lc2019_utm47.tif ./data/yunnan_lc2019_clip.tif参数说明-r near使用最近邻重采样分类值不会被插值这是处理分类栅格的基本纪律。-dstnodata 255把掩膜区设为255统计时直接剔除。-cutline指定边界矢量-crop_to_cutline让输出范围与边界完全一致不会留下矩形白底。如果数据本身已经是按云南裁剪好的重投影时不要再次裁剪否则会在边界处产生第二次重采样导致边界像元值变化。判断是否已经裁剪看gdalinfo里的角点坐标是否落在云南边界附近即可。重投影后影像的像元大小可能变成10.2米或9.8米这是地图投影变形造成的不是错了。如果后面要做变化检测所有年份的数据都必须重投影到同一坐标系否则像元错位会让你后悔没有做好这一步。3.4 为输出建金字塔与颜色表处理完成后建立金字塔能让你在QGIS里缩放不卡。分类栅格用平均采样会得到奇怪颜色所以要用最邻近采样。同时写一个简单的颜色表把地类颜色固定下来后续所有图都用这一套颜色省得每次调样式。# 建立金字塔最近邻采样 gdaladdo -r nearest ./data/yunnan_lc2019_clip.tif 2 4 8 16 # 用文本颜色表直接生成带颜色的渲染tif gdaldem color-relief -of GTiff \ -nearest_color_entry \ ./data/yunnan_lc2019_clip.tif \ ./data/lc_color.txt ./data/yunnan_lc2019_color.tif颜色表lc_color.txt的格式是每行一个值加对应R G B例如1 31 79 255表示水体蓝色。-nearest_color_entry保证像素值落在某区间时取最近的颜色不会插出中间色。这个带颜色的tif可以拖进软件当底图但项目交付时一定要保留一份未做渲染的原始分类tif。我见过有人只存了color tif结果业务方一问“你这里的类4是什么”他只能对着颜色猜非常被动。4. 分类图转矢量与面积统计让土地覆盖数据变成业务口径4.1 像元面积统计不要直接数像元数统计面积时很多人直接用像元数乘以100平方米这在局部小范围内成立但全省范围会累积投影变形误差。更稳的做法是读取栅格后统计每类像元数再乘以像元实际面积。UTM 47N下像元面积随纬度和位置变化很小但严谨起见我们可以用栅格的地理变换参数计算面积。from osgeo import gdal import numpy as np ds gdal.Open(./data/yunnan_lc2019_clip.tif) band ds.GetRasterBand(1) # 读取分类栅格把NoData之外的像元作为有效区域 clc band.ReadAsArray().astype(np.uint8) nodata band.GetNoDataValue() valid_mask clc ! nodata # 统计各类像元数 classes, counts np.unique(clc[valid_mask], return_countsTrue) # 获取像元尺寸UTM下单位为米 gt ds.GetGeoTransform() pixel_area abs(gt[1] * gt[5]) # 输出每个类别的面积平方公里 for cls, cnt in zip(classes, counts): area_km2 cnt * pixel_area / 1e6 print(f类 {cls}: {cnt} 像元, {area_km2:.2f} km²)参数说明gt[1]是东西方向像元宽gt[5]是南北方向像元高投影坐标系下通常为负乘积绝对值就是单像元面积。这个代码的隐含假设是投影后像元是规则矩形UTM下近似成立。统计前先剔除NoData否则边界外的黑色区域会被计入类0或类255。如果你发现某类面积大得离谱先回头检查NoData值而不是怀疑算法。4.2 栅格转矢量设置聚合参数避免碎面业务方经常要shp不要tif。但直接把分类栅格转矢量会产生密密麻麻的碎面云南这种陡峭地带更严重一个10米像元的变化就是一个多边形。常见做法是先做众数滤波再转矢量最后按面积过滤碎面。# 先做3x3众数滤波去掉孤立像元 gdal_fillnodata.py -md 3 -si 1 ./data/yunnan_lc2019_clip.tif ./data/yunnan_lc2019_fill.tif # 使用gdal_polygonize.py转矢量 gdal_polygonize.py \ ./data/yunnan_lc2019_fill.tif \ -f ESRI Shapefile \ ./data/yunnan_lc2019_poly.shp \ ./data/yunnan_lc2019_poly layer1这里有个隐藏坑gdal_polygonize.py输出的shp字段只有DN随后用ogr2ogr按面积过滤时如果shp没有定义投影坐标系面积字段可能是经纬度平方不是平方米。所以更稳的做法是转出后先定义投影再计算面积并过滤。# 定义投影并过滤碎面保留面积大于1公顷的图斑 ogr2ogr -overwrite \ -t_srs EPSG:32647 \ -dialect sqlite \ -sql SELECT *, ST_Area(geometry) AS area_m2 FROM yunnan_lc2019_poly WHERE ST_Area(geometry) 10000 \ ./data/yunnan_lc2019_poly_clean.shp \ ./data/yunnan_lc2019_poly.shpST_Area在投影坐标系下返回平方米-t_srs确保输出shp自带投影信息。过滤阈值10000平方米换来的是图面整洁和更快的渲染速度。代价是丢掉了小于1公顷的独立地类图斑如果你的项目关注小规模零散地块就不要过滤这么狠改到1000平方米或保留全部。这个权衡没有标准答案取决于业务口径。4.3 按行政区汇总统计表长什么样业务上最常要的统计是“云南省各州市土地利用面积”或“某流域地类构成”。用rasterstats库可以直接对矢量分区做分类统计输出一个CSV表。import geopandas as gpd import pandas as pd from rasterstats import zonal_stats # 读取州市界线矢量并统一投影到UTM 47N zones gpd.read_file(./data/yunnan_cities.shp).to_crs(EPSG:32647) # 对分类栅格做分区统计categoricalTrue会统计每个类别的像元数 stats zonal_stats( zones, ./data/yunnan_lc2019_clip.tif, categoricalTrue, nodata255, geojson_outTrue, ) rows [] for st in stats: props st[properties] row {市州: props[name]} for k, v in props.items(): # 过滤掉非统计字段和NoData类 if k.isdigit() and int(k) ! 255: row[f类{k}_km2] round(v * 100 / 1e6, 2) rows.append(row) df pd.DataFrame(rows) df.to_csv(./data/yunnan_lc2019_by_city.csv, indexFalse)categoricalTrue会让zonal_stats直接返回每个分类值的像元数。乘以100平方米再除以1e6得到平方公里。这段代码输出的表会非常整齐每行一个市州每列一个地类面积。报数据的时候最好同时输出一个面积占比列例如“类2占该州市总面积的比例”这样领导才不用自己拿计算器去按。占比可以直接计算row[f类{k}_km2] / total_area * 100。另外提醒一句行政边界矢量版本不同统计结果会有几个百分点的差异。如果项目跨年度对比所有年份必须使用同一版本的边界文件否则变化量会掺入边界修订的噪声。5. 避坑云南地形导致的五个常见问题与排查5.1 全黑或全白先看NoData再看直方图现象把tif拖进软件全黑拉伸无效。原因有两个一是NoData被设成00又是合法地类值渲染时整图被当成空值二是颜色表没有随文件加载软件找不到渲染映射。排查顺序是先看直方图。gdalinfo -hist ./data/yunnan_lc2019_clip.tif | grep -A 20 Histogram如果直方图显示0像元占比99%那文件本身可能下载错了。如果其他值正常只是渲染黑就用gdal_translate -a_nodata 255把NoData改掉再重新加载。注意gdal_translate会重写整个文件执行前先备份否则改了后悔没药可吃。5.2 高山积雪被分进水体现象德钦、香格里拉一带统计出的水体面积远高于常年水面。原因2019年产品中冰川、雪地与水体在部分合成场景中存在同物异谱阴坡积雪被分成了水体类。处理办法不是直接改分类栅格而是叠加DEM做掩膜。# 用DEM计算坡度和海拔把海拔4000米以上且坡度大于20度的水体像元重分类为高山不确定 gdaldem slope ./data/dem_utm47.tif ./data/dem_slope.tif \ -p -s 111120 -a 10 # 用gdal_calc.py做条件赋值 gdal_calc.py -A ./data/yunnan_lc2019_clip.tif \ -B ./data/dem_utm47.tif -C ./data/dem_slope.tif \ --outfile./data/yunnan_lc2019_masked.tif \ --calc((A1) * (B4000) * (C20)) * 200 (A!1) * A解释一下-A是分类图-p生成百分度坡度-s指定水平因子。gdal_calc.py中(A1)表示原分类为水体(B4000)和(C20)表示高海拔陡坡满足条件的像元重分类为200你可以定义成“高山不确定”。其他像元保留原值。加这样一层掩膜之后水体统计就不会被积雪地带污染了。5.3 边界出现连续不自然条带现象云南省界线外侧有一圈与界线平行的异常区块颜色介于两个地类之间。原因gdalwarp裁剪时边界外像元被赋予默认的srcnodata 0而0恰好是合法分类值导致边界带被填充成错误类别。解决方法是重投影时显式指定源NoData并在裁剪后重新检查。# 如果已经出现条带用缓冲矢量收边 ogr2ogr -dialect sqlite \ -sql SELECT ST_Buffer(geometry, -10) FROM yunnan_boundary \ ./data/yunnan_boundary_inner.shp ./data/yunnan_boundary.shp # 用内缩边界重新裁剪 gdalwarp -overwrite -cutline ./data/yunnan_boundary_inner.shp \ -crop_to_cutline -dstnodata 255 \ ./data/yunnan_lc2019_utm47.tif ./data/yunnan_lc2019_fixed.tifST_Buffer(geometry, -10)把边界向内收缩10米也就是一个像元把受污染的边缘像元裁掉。这样做会让统计面积比实际略小但换来的干净边界对出图更有价值。5.4 统计面积和官方年鉴对不上现象用栅格统计的耕地面积与省自然资源厅年鉴数字差距超过15%。原因分类器把大量撂荒地、梯田、园地归并到了草地或灌木同时年鉴口径来自国土调查分类体系完全不同。不要试图在这份数据上去“对齐”国土调查。正确的定位是10米覆盖产品适合看空间分布、相对占比和变化趋势不适合做绝对面积法定统计。如果你必须给一个估算值那就按自己的重分类表调整例如把“草地”中的25%估计为撂荒耕地然后单独写一个说明字段。这种做法属于模型估算不是数据修正千万别把栅格里的值直接改了。5.5 图例颜色和别人发布的风格不一致现象加载官方样式文件后水体是红色、森林是黑色。原因不同发布版本的颜色板顺序不一样文件名里都是2019但实际对应ESA或ESRI不同批次。最稳的做法是自己生成颜色表然后存成QGIS样式。颜色表的格式很简单每个分类值一行写RGB。我常用的配色是水体31-79-255树木10-107-53草地168-212-143耕地255-211-0湿地170-170-170。保存为lc_style.qml后整个项目组的出图效果就统一了不再有人“突然交出一张反色图”。5.6 一条命令快速体检写一个几十行的脚本太累直接用一个Python一行统计方式来做整体体检python3 -c from osgeo import gdal import numpy as np dsgdal.Open(./data/yunnan_lc2019_clip.tif) bds.GetRasterBand(1) ab.ReadAsArray() vals,cntsnp.unique(a,return_countsTrue) totala.size for v,c in zip(vals,cnts): print(v, round(c/total*100,2), %) 如果输出结果中NoData值255占比超过5%说明裁剪边界留白太多或者镶嵌时源数据就存在大量空洞。如果某个类别占比接近0%结合地理位置判断是合法还是异常。这个命令花不了几秒但能把分类分布、NoData比例、异常值一次看清。我习惯在每次重分类后都跑一遍形成一条快速自检的肌肉记忆。6. 进阶用NDVI做2019年土地覆盖的精度自检拿到这份10m数据如果没有实地样点怎么判断它到底靠不靠谱一个低成本的办法是用同一年的哨兵2号影像合成NDVI对分类结果做交叉验证。原理很简单不同地类的NDVI分布应该有明显区分。水体NDVI接近0或为负森林通常高于0.6草地和耕地集中在0.3到0.6裸地和建设用地低于0.2。如果某个“森林”像元对应位置的NDVI只有0.15那这个分类十有八九是错的。操作上我不建议下载整年影像做逐像元验证那样计算量太大。更实用的方案是分层随机抽样每个地类别随机抽200个点提取这些点位上的NDVI中位数和标准差画出每类的箱线图。如果某类别的NDVI分布与经验值完全偏离就该怀疑这类被混淆了。import geopandas as gpd import numpy as np from osgeo import gdal from rasterio.sample import sample_gen import rasterio # 读取分类栅格生成每个类别的随机采样点 src rasterio.open(./data/yunnan_lc2019_clip.tif) clc src.read(1) valid_mask (clc ! 255) # 用numpy随机采样200个位置 rows, cols np.where(valid_mask) np.random.seed(42) idx np.random.choice(len(rows), size2000, replaceFalse) samples [src.transform * (cols[i], rows[i]) for i in idx] class_vals clc[rows[idx], cols[idx]] # 读取对应的NDVI影像8月晴天的单景或合成 ndvi_ds rasterio.open(./data/yunnan_aug_ndvi.tif) ndvi_vals np.array([x[0] for x in sample_gen(ndvi_ds, samples)]) # 按类别统计NDVI中位数 for cls in np.unique(class_vals): mask class_vals cls if mask.sum() 10: print(类, cls, 中位数NDVI, np.median(ndvi_vals[mask]))这段代码用了rasterio的sample_gen直接读取采样点处的NDVI值。注意两个栅格必须保持地理坐标一致如果分类图已经投影到UTM 47NNDVI影像也要用同样投影否则采样点会落到错位位置。这也是为什么我在第3章坚持先把投影统一后面所有验证步骤都依赖坐标系一致。做完这步你会对这份数据的“坑位”有直观认识。比如我跑到云南南部发现常绿阔叶林的NDVI中位数只有0.55而北部针叶林能到0.75这其实是物候差异不一定分类错但如果某块“水体”的NDVI中位数是0.45那绝对是分错了。我个人的习惯是无论数据来自哪里永远保留一份未改动的原始tif然后把所有中间产物都命名为带处理后缀的版本。这样做是吃过亏之后的教训——有次我直接改原图结果后面想重新对比不同滤波参数时发现原始版本已经救不回来了只能重新解压。数据备份是最后的后悔药而这套NDVI自检流程则能在你向项目汇报前提前把明显错分类的区域识别出来。希望帮到你。本文还有配套的精品资源点击获取