尧图网站建设 
简介2019年10m精度云南省土地覆盖土地利用RAR包面向GIS、遥感、测绘及城乡规划从业者提供按云南省及下辖16个州市行政边界裁剪完成的栅格土地覆盖数据。数据基于10米哨兵影像与深度学习方法制作覆盖耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪/冰等类别原始墨卡托坐标已统一转为WGS84坐标系可直接叠加其他地理图层服务于资源调查、生态评估与国土空间规划。压缩包共112个文件以16个州市为单位分幅存放每个行政区均包含tif栅格、tfw坐标定位、xml元数据、dbf属性统计、png快视图、cpg字符编码及xlsx面积表格总大小115.71MB便于按区域检索与分州处理。使用者无需自行拼接或投影转换用ArcGIS/QGIS打开即可开展面积量算、地类统计与专题制图。已有327人学习下载是开展高分辨率土地覆盖研究与制图的实用基础数据。1. 刚拿到这个2019年10m精度云南省土地覆盖土地利用.rar先别急着解压刚拿到这份2019年10m精度云南省土地覆盖土地利用.rar的人第一反应大多是赶紧解压、拖进ArcGIS看一眼云南省的森林和水体。这个动作本身没错但翻车率不低有人叠加到矢量上发现地物边界偏移了上百米有人统计出的面积比官方年鉴差出一倍还有人没读随包类别表把建设用地认成耕地。这份数据本质上是把Sentinel-2这类10m多光谱影像分类成离散的土地覆盖/土地利用栅格每个像元保存一个类别整数编码同时配套水体、森林、草地、耕地、建设用地、裸地等映射表。如果你只想要一张大体正确的示意它够用但如果你想拿“10m精度”去算地块面积、做建设用地变化监测或当野外调查底图那必须先把这个压缩包里的黑匣子打开做一轮完整的体检和预处理。这篇文章把我平常处理这类数据包的流程走一遍从rar包一路做到能出图、能统计、能验证的成果并把最容易踩的四个坑直接写在明处。2. 解读rar包里的数据结构类别体系、坐标系统和文件格式拿到rar包别急着双击解压先打开压缩包浏览文件清单。用7-Zip或WinRAR临时打开看是否有分卷、是否有注释。很多从网盘分享出来的数据包会带一个密码说明或者把数据放在子目录里。如果压缩包损坏后续全是白忙。正常情况下一个土地覆盖数据包会包含栅格主体、元数据、类别说明和配色样式这些文件的完整程度直接决定你后面能不能把事情做好。2.1 解压rar后的典型文件清单解压后目录里通常至少有三种东西GeoTIFF栅格文件.tif、描述数据的XML元数据、类别对照表.csv或.xlsx。运气好的话还有配色样式文件.sld / .qml和面积统计样本。如果只给你一个孤零零的tif而没有元数据那类别编码就只能靠猜后面统计面积时很容易张冠李戴。在Linux下我一般用unar解压它对中文文件名支持比unrar更稳能避免把“云南省”解压成乱码。Windows用户用7-Zip右键解压即可但要注意目标路径不要太深否则后续GDAL处理会报路径过长。# 推荐用 unar处理中文文件名更稳 unar -f -o ./yunnan_lc ./2019年10m精度云南省土地覆盖土地利用.rar # -f 强制覆盖同名文件-o 指定输出目录解压完成后先查看文件列表和目录结构find . -maxdepth 2 -type f | sort这一步能看到是否存在多个分幅tif、是否有类别CSV和XML。常见问题是不看清单直接拖tif进GIS最后发现类表没载入渲染全是单色。下面这个表是标准数据包应该具备的文件缺少时建议回到数据源去要文件用途缺失时的影响.tif 栅格存放每个像元的类别编码无数据主体其他都白搭.xml / .txt 元数据记录投影、传感器、生产日期、版本无法确认坐标系和类别来源.csv / .xlsx 类别表定义类别编码与中文名称统计结果无法解释容易算错地类.qml / .sld配色和渲染样式出图时颜色乱套仍可自己配.rar 校验文件校验下载完整性解压出错时找不到原因2.2 类别体系怎么读从整型像元值到中文地类名的映射表土地覆盖/土地利用数据是单波段整型栅格每个像元的值对应一个类别。不同生产机构用的编码完全不同有的用0到9有的用10、20、30这样的间隔还有的把NoData设为0、背景也设为0。所以第一步永远是读类别表而不是直接用GIS里的“唯一值”去猜。类别表通常是CSV因为生产方多为国内团队编码常见GBK。用Python读时先试gbk报错再试utf-8import pandas as pd # 类别表常用gbk编码读失败再换utf-8 df pd.read_csv(class_legend.csv, encodinggbk) print(df.head(10))代码逻辑很简单把类别表读进来确认Value列与你想用的地类一一对应。重点看两处一是0是不是背景或NoData二是“不透水层”这类编码是否被单独拆出。拿到类别表后下一步用栅格的直方图核对表里的值是否真实出现。如果类别表有5个类但直方图出现第6个峰说明要么NoData没设置要么这个产品实际类别比你拿到的表多。也可以用gdalinfo快速看直方图gdalinfo -hist -mm landcover.tif | tail -40这里-hist会输出像元值分布直方图-mm会输出最小最大值。观察直方图有没有异常的单峰尖刺特别是0值数量巨大时基本可以断定NoData没有正确识别后续面积计算必须先清洗。2.3 坐标系与投影的判定方法10m精度的土地覆盖产品很多直接给WGS84经纬度坐标还有一些按UTM分带存储。云南省横跨东经97度到106度横跨UTM 47和48两个带如果你拿到的是全云南省单幅数据且用UTM 47N西侧边缘的变形和误差就会明显偏大。所以我拿到tif后的第一件事是用gdalinfo确认坐标系而不是直接叠加做判断。gdalinfo landcover.tif | egrep ^(PROJCRS|GEOGCRS|AUTHORITY|Size is|Pixel Size|Origin)这条命令会输出投影名称、椭球体、权威EPSG编码、行列数和像元尺寸。如果输出里带有AUTHORITY[EPSG,4326]那就是WGS84经纬度如果是PROJCRS[WGS 84 / UTM zone 47N]就是投影坐标。还有一种情况是CGCS2000地理坐标EPSG编码会是4500、4526之类不要看到“GEOGCRS”就默认是WGS84。如果压缩包里包含多个分幅tif建议循环检查所有文件的坐标系是否一致for f in *.tif; do echo $f; gdalinfo $f | grep -E AUTHORITY | head -1; done这个脚本会把每个tif的EPSG编码打出来。若分幅之间坐标系统一说明是生产方已经预处理过若不一致就标记下来等到第4章统一重投影。3. 用GDAL完成解压后第一轮体检元数据读取、空值统计、边界核对解压完成只是第一步我习惯把“拿到tif”到“敢放心用”之间至少做三轮检查元数据、类别统计、边界范围。很多人在这一步为了省时间直接跳过去最后在制图或面积统计时发现问题再回来重做反而更花时间。3.1 用gdalinfo确认行列数、分辨率、投影与NoData副本栅格信息是后续所有处理的基础。重点看四件事像元总数、像元尺寸、坐标系、NoData值。用gdalinfo一条命令就能拿全。gdalinfo -mm -stats yunnan_lc.tif输出里的Size is 71230, 81230代表列数和行数Pixel Size (9.98, -9.98)代表像元尺寸是10m左右最下面NoData Value0代表0是背景。如果像元尺寸出现明显的小数比如9.95或10.03说明原始数据重采样时做过几何校正后续做面积统计之前最好重新对齐网格。用Python rasterio读取更直观尤其适合写脚本批量检查import rasterio with rasterio.open(yunnan_lc.tif) as src: print(CRS:, src.crs) print(尺寸:, src.width, src.height) print(分辨率:, src.res) print(NoData:, src.nodata) data src.read(1) print(有效值范围:, data[data ! src.nodata].min(), data.max())注意rasterio的src.nodata可能返回None但数据里实际存在0值背景。遇到这种情况建议手动把NoData指定为0否则后续gdalwarp会保留0值而不会被当背景裁掉。3.2 统计土地利用/覆盖类别面积分布有了类别表就可以统计各类别面积。最稳的方法是把栅格读成numpy数组剔除NoData后用np.unique统计像元数再乘以单像元面积。这里有一个容易忽略的问题如果数据还是WGS84经纬度坐标像元尺寸是度而不是米直接相乘出来的结果是“平方度”不是面积。更合理的做法是重投影成等积投影后再统计。下面这段代码假设数据已经投影到米制坐标系比如UTM或Albers。import numpy as np import rasterio with rasterio.open(yunnan_lc.tif) as src: data src.read(1) nodata src.nodata # 剔除NoData统计各类别像元数 if nodata is not None: valid data[data ! nodata] else: valid data vals, counts np.unique(valid, return_countsTrue) pixel_area abs(src.res[0] * src.res[1]) # 面积单位为平方米 for val, cnt in zip(vals, counts): print(f类别{int(val)}: {cnt} 个像元, 面积 {cnt * pixel_area:.2f} 平方米)这里用abs是因为rasterio的y方向分辨率通常是负值代表像元从北往南排列src.res[0] * src.res[1]直接乘会得到负数用abs取绝对值。如果数据是经纬度坐标先不要用这段代码到第4章重投影后再统计。此外还要注意像元面积分辨率是10m但实际统计面积会受边界裁切影响。如果数据边界没有贴合行政区边界面积统计是“矩形范围内的面积”不是真正的云南省面积。所以要配合行政区矢量做裁切。3.3 边界与影像范围核对tif的四至范围和行政区边界是否吻合直接决定你能不能把它作为基础底图使用。我一般用gdaltindex生成一个范围轮廓再叠加到QGIS里看。gdaltindex -t_srs EPSG:4326 -f GeoJSON bounds.geojson yunnan_lc.tifgdaltindex会为每个输入的栅格生成一个包围框多边形这里指定输出WGS84坐标方便和县级、省级区划叠加。叠加后发现边界与云南边界相差几十公里说明数据范围是标准的矩形还需要按省界裁切。如果想看实际有数据的区域而不是矩形边界可以先把有效像元变成掩膜再转成面gdal_calc.py -A yunnan_lc.tif --calcA0 --outfilevalid_mask.tif --NoDataValue255 gdal_polygonize.py valid_mask.tif -mask valid_mask.tif valid_area.shp第一条gdal_calc.py生成二值掩膜有效像元为10值背景被排除第二条gdal_polygonize.py把值为1的区域转成面要素。这样就能在GIS里清晰看到这个产品真正覆盖了多少地方是否存在空洞。4. 裁剪、重投影与格式转换把10m分幅数据做成能用的本地底图体检完以后通常需要对原始数据做三件事按行政区裁剪、重投影到目标坐标系、把多幅数据拼接成统一栅格。这三步的顺序一般建议先裁剪再重投影再拼接如果有多幅且跨带可以先拼接再重投影但这样会放大内存消耗。我习惯先按省界/市界裁掉无关区域数据量小了重投影和拼接都更快。4.1 按行政区边界裁剪避免跨带孤岛云南边界曲折数据四周的0值背景如果不裁掉后续面积统计会自动把0当背景处理但出图时黑色边框很难看。裁剪使用gdalwarp的cutline参数配合crop_to_cutline让输出范围严格贴合矢量边界。gdalwarp -overwrite -cutline yunnan_province.shp -crop_to_cutline \ -dstnodata 0 -of GTiff -co TILEDYES -co COMPRESSDEFLATE \ yunnan_lc.tif yunnan_cut.tif说明-cutline指定省界矢量这里建议先用ogr2ogr把矢量投影到和栅格一致的坐标系避免自动转换带来额外误差-crop_to_cutline让其输出范围正好是矢量边界而不是矩形裁剪-dstnodata 0将背景统一写成0。如果你把NoData设为0但类表里的0也是真实类别后面就容易出错因此要提前确认类表遇到0是真实类别的产品把背景改成255。裁剪后检查一下边界是否圆滑如果出现大量锯齿或空洞多半是矢量与栅格坐标系没对齐。比如省界用的是CGCS2000栅格用的是WGS84两者虽然在多数地区差别不大但在高山区域会造成几十米偏移。稳妥做法是先把省界转成与栅格一致ogr2ogr -t_srs EPSG:4326 yunnan_wgs84.shp yunnan_province.shp再用转换后的矢量做cutline。这样裁剪边界的坐标参考和栅格完全一致不会因为投影转换在边界产生毛刺。4.2 重投影到CGCS2000/UTM并控制像素尺寸云南地形起伏大如果只做目视分析用WGS84地理坐标也可以但如果要计算面积必须使用等积投影。我推荐用Albers等积投影中央经线105°E双标准纬线25°N和47°N基本覆盖云南全境并保持面积守恒。制图输出时再按目标片区转成UTM 47N或48N。分类栅格重投影有一个关键参数重采样方式必须用near最近邻不能用bilinear或cubic否则会在类别之间插出本不存在的类别比如在林地和水体之间凭空生成草地。gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs \ -r near -tr 10 10 -dstnodata 0 -overwrite \ yunnan_cut.tif yunnan_cut_aea10.tif这里-t_srs可以接受一个PROJ字符串也可以写EPSG:xxxx。如果你需要输出到CGCS2000建议把datumWGS84换成ellpsGRS80并加上中央经线参数具体EPSG编号则要按所在经度选带。-tr 10 10强制输出像元为10m的正方形如果省略输出分辨率会根据源数据计算可能出现9.98这种小数影响后续统计的口径。制图用UTM投影的命令略有不同gdalwarp -t_srs EPSG:32647 -r near -tr 10 10 -dstnodata 0 \ -overwrite yunnan_cut.tif yunnan_cut_utm47.tifEPSG:32647是WGS84 UTM 47N。注意UTM 47N只覆盖东经99度到105度如果数据涉及云南西北部东经97-99度该区域应使用48NEPSG:32648。如果整个云南省统一出图不如直接用Albers避免一个图幅里出现两套投影。4.3 分幅拼接与格式转换把多块tif合并成一个vrt如果数据是按标准分幅发布的比如每个经纬度一格一个云南省有几十个分幅。直接逐个处理很痛苦一般先用gdalbuildvrt生成一个虚拟目录再物化成单一GeoTIFF。gdalbuildvrt -srcnodata 0 -vrtnodata 0 yunnan_lc.vrt *.tif gdal_translate -co COMPRESSDEFLATE -co BIGTIFFYES -co TILEDYES \ yunnan_lc.vrt yunnan_lc_merge.tifVRT本身不复制数据只是记录所有分幅的路径、位置和投影信息打开速度特别快。生成VRT后如果分幅之间有重叠默认会显示先添加到命令的那个文件如果重叠区域的类别不一致就需要回到源数据检查接边问题。最后用gdal_translate把它真正输出成一个tif。BIGTIFFYES是必须的云南省全幅10m数据很容易超过4GB不加这个参数gdal会报“创建失败”或自动转成分卷文件。拼接完成后用gdalinfo -stats再看一次统计值确认拼接过程中NoData是否被正确保留。如果发现接缝处出现黑色斑块或异常类别多半是分幅之间存在NoData间隙需要回到第5章的避坑部分处理。5. 避开四类常见坑类别错位、精度虚标、坐标系漂移与文件损坏排查处理这种网盘流出的数据包大多数人翻车不在技术难度而在默认数据没问题。下面四个坑是我反复遇到的每条按“现象—原因—解决”写清楚你可以直接对照排查。5.1 类别错位统计结果和官方面积对不上现象把栅格统计的林地面积和年鉴一对比差出三四成甚至翻倍或者统计表里出现一个莫名的类数量巨大。原因最常见的是NoData没有正确识别0值背景被当成真实类别统计比如类表里0正好代表裸地背景也写成0裸地瞬间膨胀。另一个原因是类别表编码错位比如类表把10定义为耕地但这个产品其实用1表示耕地、2表示林地你拿着另外一份相似产品的类表硬套。解决不要相信文件名或README里的“标准类表”直接用gdalinfo -hist看直方图把直方图里的像元值和随包类表逐一比对。统计前先打印所有唯一值gdalinfo yunnan_lc.tif | grep -A 100 Histogram如果输出里的直方图显示0值有上亿像元而其他类别只有几百万基本可以判断0是背景。这时在裁剪和重投影命令里把-dstnodata改成255、0留给真实类别或用一个条件计算把0转为NoDatagdal_calc.py -A yunnan_lc.tif --outfileclean.tif --calc(A0)*255(A0)*A --NoDataValue255这条命令把等于0的像元变成255把大于0的像元原值保留然后指定255为NoData。之后统计面积前再看一眼直方图确认异常峰消失。5.2 精度虚标10m分辨率不等于90%分类精度现象拿着10m分类结果和无人机影像对比发现山区的细碎地块大量错分林地里面夹杂草地建设用地漏掉很多新小区。原因10m是空间分辨率不代表分类精度。卫星影像的云遮挡、时相差异、山体阴影都会让分类算法出错在云南这种高差大、地块破碎的地方精度打八折都不奇怪。生产方给的“总体精度90%”往往是抽样区域的结果不代表你家山沟里也那么准。解决使用前先做局部验证。拿几个已经知道底细的小区域和高分辨率影像叠在一起数一下错分像元比例。如果发现山区阴影区错分率高建议把那部分像元在分析中标记为“低置信”不要直接参与面积统计。更严谨的做法是第6章的随机分层抽样用小样本算一个实际混淆矩阵再决定能不能用这个数据做业务。不要因为标题写“10m精度”就默认它是准的。5.3 坐标系漂移WGS84与CGCS2000混用导致几十米偏移现象数据和其他测绘成果叠加边界和高精度影像错位50到100米在山谷区域尤其明显。原因很多国际产品用WGS84国内测绘数据用CGCS2000两者在大部分地区点位差异只有几米到十几米但在边界或局部高程起伏大的地区肉眼可见。还有一种是产品虽标着WGS84但实际生产时做了局部几何校正和标准WGS84网格对不上。解决先确认产品元数据和最终使用参考框架。如果产品是WGS84而你要叠加县域国土数据就把产品重投影到CGCS2000不要直接在GIS里改“坐标系显示”那只是换了帽子像元位置没动。用第4章的gdalwarp命令把源数据从EPSG:4326转到CGCS2000相应投影例如gdalwarp -t_srs EPSG:4535 -r near -tr 10 10 yunnan_cut.tif yunnan_cgcs2000.tif但EPSG:4535具体对应哪个CGCS2000带要看你的区域不确定就不要乱猜编号。稳妥做法是先用gdalsrsinfo查一个已知CGCS2000带确认范围覆盖云南后再转。对于历史数据如果偏移是整体性的还可以用控制点做一次仿射校正这里不展开。5.4 文件损坏与分幅接边缝隙现象解压到一半报“CRC错误”或“文件头损坏”拼完的VRT出现接边黑条、细长空洞有的分幅纹路不对。原因网盘下载不完整或上传时文件损坏尤其是rar分卷少了一个卷。接边黑条则往往是分幅数据边缘有NoData拼接后形成了缝隙。解决rar包解压前一定先测试完整性rar t 2019年10m精度云南省土地覆盖土地利用.rar这条命令会遍历所有分卷报错就重新下载不要心存侥幸。接边缝隙可以用gdalwarp先给每个分幅加少量缓冲再重新拼接如果缝隙已经形成可以用邻域填充不过分出图比如用gdal_fillnodata.pygdal_fillnodata.py -md 10 -si 0 -o yunnan_filled.tif yunnan_merge.tif-md 10表示最大填补距离为10个像元-si 0关闭平滑迭代只填补空洞不改变其他像元值。填补后还需要把边界重新裁剪一次防止填充到图幅外的背景区域。6. 数据落地后的进阶用法专题制图、分区统计与精度验证从rar包到干净的栅格只是开始。真正让这份10m数据产生价值的是后面的专题制图和统计。这一章讲我用的三个进阶手段都能直接抄。6.1 用样式文件快速出图sld/qml与配色土地覆盖图用单一色带渲染很难看因为类别是离散值必须用唯一值配色。如果随包没有提供样式可以根据类别表生成QGIS可用的QML文件。下面这段脚本会遍历类别CSV生成一个最简单的唯一值渲染QMLimport pandas as pd df pd.read_csv(class_legend.csv, encodinggbk) with open(landcover.qml, w, encodingutf-8) as f: f.write(qgispipe-datarenderer-v2 typesinglebandpseudocolor) for _, row in df.iterrows(): value int(row[Value]) r, g, b int(row[R]), int(row[G]), int(row[B]) f.write(fitem value{value} color{r},{g},{b} label{row[Name]}/) f.write(/renderer-v2/pipe-data/qgis)这个脚本依赖CSV里有R、G、B列如果没有你可以自己在QGIS里手动双击每个类别调色。调好后右键图层→导出→保存QML下次直接加载省得每次重复配颜色。QML只改显示不改原始栅格值出图时非常方便。6.2 分区统计按州市矢量统计各类面积如果要做“云南省各州市森林覆盖面积排行榜”用rasterstats一行就能搞定。它会对每个矢量面计算栅格各类别像元数再乘以像元面积from rasterstats import zonal_stats stats zonal_stats( cities.shp, yunnan_cut_aea10.tif, categoricalTrue, nodata0 ) for i, st in enumerate(stats): print(i, st)categoricalTrue返回一个字典键是类别值值是像元数。注意nodata0要和你的栅格一致否则0值会被当成类别统计进去。得到像元数后因为栅格已经是10m等积投影面积就是像元数乘以100平方米再换算成平方千米。这个统计速度快适合批量输出表但不能替代空间分析因为它是按多边形边界切分的。6.3 精度验证的抽样策略和我的收尾最后说验证。我一直坚持任何土地覆盖数据在用于报告之前都必须做一次小样本验证否则就是拿运气赌结果。做法是分层随机抽样先按类别把栅格多边形化在每个类别内部随机布设10到30个点然后独立判读这些点的高分辨率影像记录“一致/不一致”最后算混淆矩阵。一个简化版验证脚本可以这样import numpy as np from sklearn.metrics import confusion_matrix, classification_report # predicted 是栅格中抽样点的类别reference 是人工判读类别 predicted np.array([1, 2, 2, 3, 1, 2, 3, 3, 1, 2]) reference np.array([1, 2, 1, 3, 1, 2, 3, 3, 1, 2]) print(confusion_matrix(reference, predicted)) print(classification_report(reference, predicted))这里的预测类别从土地覆盖栅格提取人工判读类别从现地调查或高分辨率影像得到。两列对齐后输出混淆矩阵会直接告诉你哪些类别容易混。我自己做云南山地项目时最常见的混类是草地和低矮灌木如果验证显示这两类分不清我宁可把它们合并成“灌草混合”而不是硬分。收尾说一句我自己的习惯每次从这类数据包里取数我都会把验证点结果、数据时相和原始rar的校验信息写在一个CSV里放在处理目录下。没有时相记录的验证点过了三个月再看就废了文件来源不清晰的栅格换了电脑就可能被覆盖。数据能用多久往往取决于你把“元信息”留得多完整。希望帮到你。本文还有配套的精品资源点击获取 10m精度云南省土地覆盖数据处理全流程：解压、投影、裁剪到验证
简介2019年10m精度云南省土地覆盖土地利用RAR包面向GIS、遥感、测绘及城乡规划从业者提供按云南省及下辖16个州市行政边界裁剪完成的栅格土地覆盖数据。数据基于10米哨兵影像与深度学习方法制作覆盖耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪/冰等类别原始墨卡托坐标已统一转为WGS84坐标系可直接叠加其他地理图层服务于资源调查、生态评估与国土空间规划。压缩包共112个文件以16个州市为单位分幅存放每个行政区均包含tif栅格、tfw坐标定位、xml元数据、dbf属性统计、png快视图、cpg字符编码及xlsx面积表格总大小115.71MB便于按区域检索与分州处理。使用者无需自行拼接或投影转换用ArcGIS/QGIS打开即可开展面积量算、地类统计与专题制图。已有327人学习下载是开展高分辨率土地覆盖研究与制图的实用基础数据。1. 刚拿到这个2019年10m精度云南省土地覆盖土地利用.rar先别急着解压刚拿到这份2019年10m精度云南省土地覆盖土地利用.rar的人第一反应大多是赶紧解压、拖进ArcGIS看一眼云南省的森林和水体。这个动作本身没错但翻车率不低有人叠加到矢量上发现地物边界偏移了上百米有人统计出的面积比官方年鉴差出一倍还有人没读随包类别表把建设用地认成耕地。这份数据本质上是把Sentinel-2这类10m多光谱影像分类成离散的土地覆盖/土地利用栅格每个像元保存一个类别整数编码同时配套水体、森林、草地、耕地、建设用地、裸地等映射表。如果你只想要一张大体正确的示意它够用但如果你想拿“10m精度”去算地块面积、做建设用地变化监测或当野外调查底图那必须先把这个压缩包里的黑匣子打开做一轮完整的体检和预处理。这篇文章把我平常处理这类数据包的流程走一遍从rar包一路做到能出图、能统计、能验证的成果并把最容易踩的四个坑直接写在明处。2. 解读rar包里的数据结构类别体系、坐标系统和文件格式拿到rar包别急着双击解压先打开压缩包浏览文件清单。用7-Zip或WinRAR临时打开看是否有分卷、是否有注释。很多从网盘分享出来的数据包会带一个密码说明或者把数据放在子目录里。如果压缩包损坏后续全是白忙。正常情况下一个土地覆盖数据包会包含栅格主体、元数据、类别说明和配色样式这些文件的完整程度直接决定你后面能不能把事情做好。2.1 解压rar后的典型文件清单解压后目录里通常至少有三种东西GeoTIFF栅格文件.tif、描述数据的XML元数据、类别对照表.csv或.xlsx。运气好的话还有配色样式文件.sld / .qml和面积统计样本。如果只给你一个孤零零的tif而没有元数据那类别编码就只能靠猜后面统计面积时很容易张冠李戴。在Linux下我一般用unar解压它对中文文件名支持比unrar更稳能避免把“云南省”解压成乱码。Windows用户用7-Zip右键解压即可但要注意目标路径不要太深否则后续GDAL处理会报路径过长。# 推荐用 unar处理中文文件名更稳 unar -f -o ./yunnan_lc ./2019年10m精度云南省土地覆盖土地利用.rar # -f 强制覆盖同名文件-o 指定输出目录解压完成后先查看文件列表和目录结构find . -maxdepth 2 -type f | sort这一步能看到是否存在多个分幅tif、是否有类别CSV和XML。常见问题是不看清单直接拖tif进GIS最后发现类表没载入渲染全是单色。下面这个表是标准数据包应该具备的文件缺少时建议回到数据源去要文件用途缺失时的影响.tif 栅格存放每个像元的类别编码无数据主体其他都白搭.xml / .txt 元数据记录投影、传感器、生产日期、版本无法确认坐标系和类别来源.csv / .xlsx 类别表定义类别编码与中文名称统计结果无法解释容易算错地类.qml / .sld配色和渲染样式出图时颜色乱套仍可自己配.rar 校验文件校验下载完整性解压出错时找不到原因2.2 类别体系怎么读从整型像元值到中文地类名的映射表土地覆盖/土地利用数据是单波段整型栅格每个像元的值对应一个类别。不同生产机构用的编码完全不同有的用0到9有的用10、20、30这样的间隔还有的把NoData设为0、背景也设为0。所以第一步永远是读类别表而不是直接用GIS里的“唯一值”去猜。类别表通常是CSV因为生产方多为国内团队编码常见GBK。用Python读时先试gbk报错再试utf-8import pandas as pd # 类别表常用gbk编码读失败再换utf-8 df pd.read_csv(class_legend.csv, encodinggbk) print(df.head(10))代码逻辑很简单把类别表读进来确认Value列与你想用的地类一一对应。重点看两处一是0是不是背景或NoData二是“不透水层”这类编码是否被单独拆出。拿到类别表后下一步用栅格的直方图核对表里的值是否真实出现。如果类别表有5个类但直方图出现第6个峰说明要么NoData没设置要么这个产品实际类别比你拿到的表多。也可以用gdalinfo快速看直方图gdalinfo -hist -mm landcover.tif | tail -40这里-hist会输出像元值分布直方图-mm会输出最小最大值。观察直方图有没有异常的单峰尖刺特别是0值数量巨大时基本可以断定NoData没有正确识别后续面积计算必须先清洗。2.3 坐标系与投影的判定方法10m精度的土地覆盖产品很多直接给WGS84经纬度坐标还有一些按UTM分带存储。云南省横跨东经97度到106度横跨UTM 47和48两个带如果你拿到的是全云南省单幅数据且用UTM 47N西侧边缘的变形和误差就会明显偏大。所以我拿到tif后的第一件事是用gdalinfo确认坐标系而不是直接叠加做判断。gdalinfo landcover.tif | egrep ^(PROJCRS|GEOGCRS|AUTHORITY|Size is|Pixel Size|Origin)这条命令会输出投影名称、椭球体、权威EPSG编码、行列数和像元尺寸。如果输出里带有AUTHORITY[EPSG,4326]那就是WGS84经纬度如果是PROJCRS[WGS 84 / UTM zone 47N]就是投影坐标。还有一种情况是CGCS2000地理坐标EPSG编码会是4500、4526之类不要看到“GEOGCRS”就默认是WGS84。如果压缩包里包含多个分幅tif建议循环检查所有文件的坐标系是否一致for f in *.tif; do echo $f; gdalinfo $f | grep -E AUTHORITY | head -1; done这个脚本会把每个tif的EPSG编码打出来。若分幅之间坐标系统一说明是生产方已经预处理过若不一致就标记下来等到第4章统一重投影。3. 用GDAL完成解压后第一轮体检元数据读取、空值统计、边界核对解压完成只是第一步我习惯把“拿到tif”到“敢放心用”之间至少做三轮检查元数据、类别统计、边界范围。很多人在这一步为了省时间直接跳过去最后在制图或面积统计时发现问题再回来重做反而更花时间。3.1 用gdalinfo确认行列数、分辨率、投影与NoData副本栅格信息是后续所有处理的基础。重点看四件事像元总数、像元尺寸、坐标系、NoData值。用gdalinfo一条命令就能拿全。gdalinfo -mm -stats yunnan_lc.tif输出里的Size is 71230, 81230代表列数和行数Pixel Size (9.98, -9.98)代表像元尺寸是10m左右最下面NoData Value0代表0是背景。如果像元尺寸出现明显的小数比如9.95或10.03说明原始数据重采样时做过几何校正后续做面积统计之前最好重新对齐网格。用Python rasterio读取更直观尤其适合写脚本批量检查import rasterio with rasterio.open(yunnan_lc.tif) as src: print(CRS:, src.crs) print(尺寸:, src.width, src.height) print(分辨率:, src.res) print(NoData:, src.nodata) data src.read(1) print(有效值范围:, data[data ! src.nodata].min(), data.max())注意rasterio的src.nodata可能返回None但数据里实际存在0值背景。遇到这种情况建议手动把NoData指定为0否则后续gdalwarp会保留0值而不会被当背景裁掉。3.2 统计土地利用/覆盖类别面积分布有了类别表就可以统计各类别面积。最稳的方法是把栅格读成numpy数组剔除NoData后用np.unique统计像元数再乘以单像元面积。这里有一个容易忽略的问题如果数据还是WGS84经纬度坐标像元尺寸是度而不是米直接相乘出来的结果是“平方度”不是面积。更合理的做法是重投影成等积投影后再统计。下面这段代码假设数据已经投影到米制坐标系比如UTM或Albers。import numpy as np import rasterio with rasterio.open(yunnan_lc.tif) as src: data src.read(1) nodata src.nodata # 剔除NoData统计各类别像元数 if nodata is not None: valid data[data ! nodata] else: valid data vals, counts np.unique(valid, return_countsTrue) pixel_area abs(src.res[0] * src.res[1]) # 面积单位为平方米 for val, cnt in zip(vals, counts): print(f类别{int(val)}: {cnt} 个像元, 面积 {cnt * pixel_area:.2f} 平方米)这里用abs是因为rasterio的y方向分辨率通常是负值代表像元从北往南排列src.res[0] * src.res[1]直接乘会得到负数用abs取绝对值。如果数据是经纬度坐标先不要用这段代码到第4章重投影后再统计。此外还要注意像元面积分辨率是10m但实际统计面积会受边界裁切影响。如果数据边界没有贴合行政区边界面积统计是“矩形范围内的面积”不是真正的云南省面积。所以要配合行政区矢量做裁切。3.3 边界与影像范围核对tif的四至范围和行政区边界是否吻合直接决定你能不能把它作为基础底图使用。我一般用gdaltindex生成一个范围轮廓再叠加到QGIS里看。gdaltindex -t_srs EPSG:4326 -f GeoJSON bounds.geojson yunnan_lc.tifgdaltindex会为每个输入的栅格生成一个包围框多边形这里指定输出WGS84坐标方便和县级、省级区划叠加。叠加后发现边界与云南边界相差几十公里说明数据范围是标准的矩形还需要按省界裁切。如果想看实际有数据的区域而不是矩形边界可以先把有效像元变成掩膜再转成面gdal_calc.py -A yunnan_lc.tif --calcA0 --outfilevalid_mask.tif --NoDataValue255 gdal_polygonize.py valid_mask.tif -mask valid_mask.tif valid_area.shp第一条gdal_calc.py生成二值掩膜有效像元为10值背景被排除第二条gdal_polygonize.py把值为1的区域转成面要素。这样就能在GIS里清晰看到这个产品真正覆盖了多少地方是否存在空洞。4. 裁剪、重投影与格式转换把10m分幅数据做成能用的本地底图体检完以后通常需要对原始数据做三件事按行政区裁剪、重投影到目标坐标系、把多幅数据拼接成统一栅格。这三步的顺序一般建议先裁剪再重投影再拼接如果有多幅且跨带可以先拼接再重投影但这样会放大内存消耗。我习惯先按省界/市界裁掉无关区域数据量小了重投影和拼接都更快。4.1 按行政区边界裁剪避免跨带孤岛云南边界曲折数据四周的0值背景如果不裁掉后续面积统计会自动把0当背景处理但出图时黑色边框很难看。裁剪使用gdalwarp的cutline参数配合crop_to_cutline让输出范围严格贴合矢量边界。gdalwarp -overwrite -cutline yunnan_province.shp -crop_to_cutline \ -dstnodata 0 -of GTiff -co TILEDYES -co COMPRESSDEFLATE \ yunnan_lc.tif yunnan_cut.tif说明-cutline指定省界矢量这里建议先用ogr2ogr把矢量投影到和栅格一致的坐标系避免自动转换带来额外误差-crop_to_cutline让其输出范围正好是矢量边界而不是矩形裁剪-dstnodata 0将背景统一写成0。如果你把NoData设为0但类表里的0也是真实类别后面就容易出错因此要提前确认类表遇到0是真实类别的产品把背景改成255。裁剪后检查一下边界是否圆滑如果出现大量锯齿或空洞多半是矢量与栅格坐标系没对齐。比如省界用的是CGCS2000栅格用的是WGS84两者虽然在多数地区差别不大但在高山区域会造成几十米偏移。稳妥做法是先把省界转成与栅格一致ogr2ogr -t_srs EPSG:4326 yunnan_wgs84.shp yunnan_province.shp再用转换后的矢量做cutline。这样裁剪边界的坐标参考和栅格完全一致不会因为投影转换在边界产生毛刺。4.2 重投影到CGCS2000/UTM并控制像素尺寸云南地形起伏大如果只做目视分析用WGS84地理坐标也可以但如果要计算面积必须使用等积投影。我推荐用Albers等积投影中央经线105°E双标准纬线25°N和47°N基本覆盖云南全境并保持面积守恒。制图输出时再按目标片区转成UTM 47N或48N。分类栅格重投影有一个关键参数重采样方式必须用near最近邻不能用bilinear或cubic否则会在类别之间插出本不存在的类别比如在林地和水体之间凭空生成草地。gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs \ -r near -tr 10 10 -dstnodata 0 -overwrite \ yunnan_cut.tif yunnan_cut_aea10.tif这里-t_srs可以接受一个PROJ字符串也可以写EPSG:xxxx。如果你需要输出到CGCS2000建议把datumWGS84换成ellpsGRS80并加上中央经线参数具体EPSG编号则要按所在经度选带。-tr 10 10强制输出像元为10m的正方形如果省略输出分辨率会根据源数据计算可能出现9.98这种小数影响后续统计的口径。制图用UTM投影的命令略有不同gdalwarp -t_srs EPSG:32647 -r near -tr 10 10 -dstnodata 0 \ -overwrite yunnan_cut.tif yunnan_cut_utm47.tifEPSG:32647是WGS84 UTM 47N。注意UTM 47N只覆盖东经99度到105度如果数据涉及云南西北部东经97-99度该区域应使用48NEPSG:32648。如果整个云南省统一出图不如直接用Albers避免一个图幅里出现两套投影。4.3 分幅拼接与格式转换把多块tif合并成一个vrt如果数据是按标准分幅发布的比如每个经纬度一格一个云南省有几十个分幅。直接逐个处理很痛苦一般先用gdalbuildvrt生成一个虚拟目录再物化成单一GeoTIFF。gdalbuildvrt -srcnodata 0 -vrtnodata 0 yunnan_lc.vrt *.tif gdal_translate -co COMPRESSDEFLATE -co BIGTIFFYES -co TILEDYES \ yunnan_lc.vrt yunnan_lc_merge.tifVRT本身不复制数据只是记录所有分幅的路径、位置和投影信息打开速度特别快。生成VRT后如果分幅之间有重叠默认会显示先添加到命令的那个文件如果重叠区域的类别不一致就需要回到源数据检查接边问题。最后用gdal_translate把它真正输出成一个tif。BIGTIFFYES是必须的云南省全幅10m数据很容易超过4GB不加这个参数gdal会报“创建失败”或自动转成分卷文件。拼接完成后用gdalinfo -stats再看一次统计值确认拼接过程中NoData是否被正确保留。如果发现接缝处出现黑色斑块或异常类别多半是分幅之间存在NoData间隙需要回到第5章的避坑部分处理。5. 避开四类常见坑类别错位、精度虚标、坐标系漂移与文件损坏排查处理这种网盘流出的数据包大多数人翻车不在技术难度而在默认数据没问题。下面四个坑是我反复遇到的每条按“现象—原因—解决”写清楚你可以直接对照排查。5.1 类别错位统计结果和官方面积对不上现象把栅格统计的林地面积和年鉴一对比差出三四成甚至翻倍或者统计表里出现一个莫名的类数量巨大。原因最常见的是NoData没有正确识别0值背景被当成真实类别统计比如类表里0正好代表裸地背景也写成0裸地瞬间膨胀。另一个原因是类别表编码错位比如类表把10定义为耕地但这个产品其实用1表示耕地、2表示林地你拿着另外一份相似产品的类表硬套。解决不要相信文件名或README里的“标准类表”直接用gdalinfo -hist看直方图把直方图里的像元值和随包类表逐一比对。统计前先打印所有唯一值gdalinfo yunnan_lc.tif | grep -A 100 Histogram如果输出里的直方图显示0值有上亿像元而其他类别只有几百万基本可以判断0是背景。这时在裁剪和重投影命令里把-dstnodata改成255、0留给真实类别或用一个条件计算把0转为NoDatagdal_calc.py -A yunnan_lc.tif --outfileclean.tif --calc(A0)*255(A0)*A --NoDataValue255这条命令把等于0的像元变成255把大于0的像元原值保留然后指定255为NoData。之后统计面积前再看一眼直方图确认异常峰消失。5.2 精度虚标10m分辨率不等于90%分类精度现象拿着10m分类结果和无人机影像对比发现山区的细碎地块大量错分林地里面夹杂草地建设用地漏掉很多新小区。原因10m是空间分辨率不代表分类精度。卫星影像的云遮挡、时相差异、山体阴影都会让分类算法出错在云南这种高差大、地块破碎的地方精度打八折都不奇怪。生产方给的“总体精度90%”往往是抽样区域的结果不代表你家山沟里也那么准。解决使用前先做局部验证。拿几个已经知道底细的小区域和高分辨率影像叠在一起数一下错分像元比例。如果发现山区阴影区错分率高建议把那部分像元在分析中标记为“低置信”不要直接参与面积统计。更严谨的做法是第6章的随机分层抽样用小样本算一个实际混淆矩阵再决定能不能用这个数据做业务。不要因为标题写“10m精度”就默认它是准的。5.3 坐标系漂移WGS84与CGCS2000混用导致几十米偏移现象数据和其他测绘成果叠加边界和高精度影像错位50到100米在山谷区域尤其明显。原因很多国际产品用WGS84国内测绘数据用CGCS2000两者在大部分地区点位差异只有几米到十几米但在边界或局部高程起伏大的地区肉眼可见。还有一种是产品虽标着WGS84但实际生产时做了局部几何校正和标准WGS84网格对不上。解决先确认产品元数据和最终使用参考框架。如果产品是WGS84而你要叠加县域国土数据就把产品重投影到CGCS2000不要直接在GIS里改“坐标系显示”那只是换了帽子像元位置没动。用第4章的gdalwarp命令把源数据从EPSG:4326转到CGCS2000相应投影例如gdalwarp -t_srs EPSG:4535 -r near -tr 10 10 yunnan_cut.tif yunnan_cgcs2000.tif但EPSG:4535具体对应哪个CGCS2000带要看你的区域不确定就不要乱猜编号。稳妥做法是先用gdalsrsinfo查一个已知CGCS2000带确认范围覆盖云南后再转。对于历史数据如果偏移是整体性的还可以用控制点做一次仿射校正这里不展开。5.4 文件损坏与分幅接边缝隙现象解压到一半报“CRC错误”或“文件头损坏”拼完的VRT出现接边黑条、细长空洞有的分幅纹路不对。原因网盘下载不完整或上传时文件损坏尤其是rar分卷少了一个卷。接边黑条则往往是分幅数据边缘有NoData拼接后形成了缝隙。解决rar包解压前一定先测试完整性rar t 2019年10m精度云南省土地覆盖土地利用.rar这条命令会遍历所有分卷报错就重新下载不要心存侥幸。接边缝隙可以用gdalwarp先给每个分幅加少量缓冲再重新拼接如果缝隙已经形成可以用邻域填充不过分出图比如用gdal_fillnodata.pygdal_fillnodata.py -md 10 -si 0 -o yunnan_filled.tif yunnan_merge.tif-md 10表示最大填补距离为10个像元-si 0关闭平滑迭代只填补空洞不改变其他像元值。填补后还需要把边界重新裁剪一次防止填充到图幅外的背景区域。6. 数据落地后的进阶用法专题制图、分区统计与精度验证从rar包到干净的栅格只是开始。真正让这份10m数据产生价值的是后面的专题制图和统计。这一章讲我用的三个进阶手段都能直接抄。6.1 用样式文件快速出图sld/qml与配色土地覆盖图用单一色带渲染很难看因为类别是离散值必须用唯一值配色。如果随包没有提供样式可以根据类别表生成QGIS可用的QML文件。下面这段脚本会遍历类别CSV生成一个最简单的唯一值渲染QMLimport pandas as pd df pd.read_csv(class_legend.csv, encodinggbk) with open(landcover.qml, w, encodingutf-8) as f: f.write(qgispipe-datarenderer-v2 typesinglebandpseudocolor) for _, row in df.iterrows(): value int(row[Value]) r, g, b int(row[R]), int(row[G]), int(row[B]) f.write(fitem value{value} color{r},{g},{b} label{row[Name]}/) f.write(/renderer-v2/pipe-data/qgis)这个脚本依赖CSV里有R、G、B列如果没有你可以自己在QGIS里手动双击每个类别调色。调好后右键图层→导出→保存QML下次直接加载省得每次重复配颜色。QML只改显示不改原始栅格值出图时非常方便。6.2 分区统计按州市矢量统计各类面积如果要做“云南省各州市森林覆盖面积排行榜”用rasterstats一行就能搞定。它会对每个矢量面计算栅格各类别像元数再乘以像元面积from rasterstats import zonal_stats stats zonal_stats( cities.shp, yunnan_cut_aea10.tif, categoricalTrue, nodata0 ) for i, st in enumerate(stats): print(i, st)categoricalTrue返回一个字典键是类别值值是像元数。注意nodata0要和你的栅格一致否则0值会被当成类别统计进去。得到像元数后因为栅格已经是10m等积投影面积就是像元数乘以100平方米再换算成平方千米。这个统计速度快适合批量输出表但不能替代空间分析因为它是按多边形边界切分的。6.3 精度验证的抽样策略和我的收尾最后说验证。我一直坚持任何土地覆盖数据在用于报告之前都必须做一次小样本验证否则就是拿运气赌结果。做法是分层随机抽样先按类别把栅格多边形化在每个类别内部随机布设10到30个点然后独立判读这些点的高分辨率影像记录“一致/不一致”最后算混淆矩阵。一个简化版验证脚本可以这样import numpy as np from sklearn.metrics import confusion_matrix, classification_report # predicted 是栅格中抽样点的类别reference 是人工判读类别 predicted np.array([1, 2, 2, 3, 1, 2, 3, 3, 1, 2]) reference np.array([1, 2, 1, 3, 1, 2, 3, 3, 1, 2]) print(confusion_matrix(reference, predicted)) print(classification_report(reference, predicted))这里的预测类别从土地覆盖栅格提取人工判读类别从现地调查或高分辨率影像得到。两列对齐后输出混淆矩阵会直接告诉你哪些类别容易混。我自己做云南山地项目时最常见的混类是草地和低矮灌木如果验证显示这两类分不清我宁可把它们合并成“灌草混合”而不是硬分。收尾说一句我自己的习惯每次从这类数据包里取数我都会把验证点结果、数据时相和原始rar的校验信息写在一个CSV里放在处理目录下。没有时相记录的验证点过了三个月再看就废了文件来源不清晰的栅格换了电脑就可能被覆盖。数据能用多久往往取决于你把“元信息”留得多完整。希望帮到你。本文还有配套的精品资源点击获取