简介2019年10m精度云南省土地覆盖土地利用数据源于基于哨兵影像与深度学习生成的全球陆地覆盖产品经重新投影为WGS84地理坐标系并按照最新省市级行政边界裁剪形成可直接使用的云南省栅格数据集。面向GIS与遥感分析人员、国土空间规划研究者提供耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、冰雪等十类地表覆盖信息免去自行下载、拼接、投影转换与分幅裁剪的繁琐环节。资源包共112个文件压缩后约115.71MB主要包含tif栅格数据、tfw地理配准文件、xml元数据、dbf属性表、cpg编码文件、png快速预览图及xlsx分类说明结构清晰便于直接加载、属性查询与面积统计。全省各州市均单独成图数据坐标信息完整可在ArcGIS、QGIS等主流平台中无缝使用。已有327人学习适合用于土地利用变化监测、生态环境评估、国土空间规划及科研教学等场景。1. 2019年云南10米土地覆盖数据包先搞清这是什么精度和值不值得下2019年10m精度的云南省土地覆盖土地利用数据这个压缩包我拆过不止一次里面主体是一张十米分辨率的分类栅格把云南的林地、灌丛、水体、农田、房屋、裸地都按统一编码切好。很多人看到“10m”第一反应是“比30m清楚一点”但真正用来算面积、做转移矩阵时才体会到这个精度能把山地破碎地块的边界画得准不少前提是你先弄懂它的坐标系和分类体系。这份数据适合做县域国土空间底图、生态保护评估、耕地非粮化监测或者作为训练样本的标签。它不适合只想知道“云南森林覆盖率多少”这种宏观问题因为那样的数字用30m数据加统计年鉴就够了。真正值得下载的人是要把栅格裁剪成研究区、按类别统计出可信面积甚至和另一年数据做变化分析的从业者。2. 数据包里到底有什么分类体系、坐标系与文件清单的拆解2.1 分类体系从全球产品到云南区域裁剪的取舍要判断这份数据能不能用于你的场景第一步不是找软件打开而是弄清分类体系。市面上流通的2019年10米云南土地覆盖栅格绝大多数源自基于Sentinel-2影像的全球级公开土地覆盖产品比如ESA的WorldCover系列、Esri的全球土地覆盖产品生产时用随机森林或深度模型对逐像元分类再按云南行政区边界裁剪。省级裁剪的好处是省去你自己下载全球数据和做掩膜的时间坏处是裁剪过程可能丢失原始元数据导致你拿到手的分类代码没有对应的中文图例。云南的地理特征决定了分类难点集中在“林地”和“灌丛”的区分上。高山峡谷区大量存在稀疏林地、灌丛、草甸的过渡带全球训练样本在平坦区域表现好到了云南就容易把一个山坡分成上半截林地、下半截灌丛。所以你拿到数据后首先要做的是把栅格的所有唯一值读出来列一张清单。举个例子常见编码可能如下表所示但不同来源的数字可能完全不同以下只是示意分类代码中文名称说明1乔木林地郁闭度较高2灌木林地包括天然灌丛3草地含高山草甸4耕地含水田和旱地5建设用地房屋、道路、工矿6水体河流、湖泊、水库7裸地裸土、裸岩、沙地我一般建议在QGIS里打开后先用工具栏的“识别”功能点几个典型位置比如滇池、昆明城区、哀牢山森林确认编码和地物是否对得上。如果包里的README缺失那就要自己通过高分辨率影像比对来推断每个数字的含义。这一步不做后面的面积统计都是黑匣子出了错你都不知道往哪找。2.2 坐标系与投影为什么打开后形状有点“歪”原始分类栅格一般默认采用WGS84经纬度坐标系EPSG编码是4326。像元尺寸在参数中通常写作0.00009度左右但这并不意味每个像元都是10米乘10米。在赤道附近0.00009度经度对应的地面距离约10米但在云南北纬21度到29度区域经度方向的真实距离要乘以纬度的余弦值。以北纬25度计算cos25°约等于0.906也就是说经度方向一个像元实际只有约9.06米而纬度方向仍是约11米左右导致单像元面积从“标称的100平方米”变成了实际约99.7平方米看起来差别不大但全省数亿个像元累计后面积误差能达到几千平方公里。打开后云南的轮廓看起来有点“歪”也是因为经纬度网格在低纬地区并不等距屏幕显示时默认做了等距投影所以弧线和直线都会变形。计算面积和做统计之前必须先重投影到等积投影或UTM投影。常用的方案如下投影EPSG特点适用场景WGS84 经纬度4326原始格式像元不等面积查看、拼接UTM 47N32647六度分带局部变形小单县或局部精确量算Albers 等积投影视参数而定面积守恒适合全区域汇总全省统计、制图需要特别提醒的是如果你用Albers投影参数里的中央经线和标准纬线不同输出像元大小也会不一样。最稳妥的做法是直接在GDAL的Warp里指定目标投影然后设置输出分辨率这样后续每一个像元面积都是你设定的值例如100平方米。2.3 文件清单从tif到辅助文件哪些必须保留解压这个rar包后里面通常不会只有一个孤零零的tif。一个合格的成果包至少包含这几个文件主栅格文件、坐标配准文件、元数据文件和样式文件。主栅格往往是单波段UInt8类型因为分类数最多不超过几十个一个字节足够文件体积一般在几百MB到几个GB之间云南全省10米分辨率大约有十亿量级像元压缩成LZW或DEFLATE后能小很多。坐标配准文件.tfw非常关键它记录了栅格左上角坐标、像元尺寸以及旋转参数。如果没有它很多软件会把tif当成无坐标文件无法套到正确位置。元数据文件.xml或.json包含数据源、投影、生产日期、分类体系这是解释数据最权威的凭证。样式文件.qml或.lyr是给QGIS或ArcGIS用的配色方案能让你打开时直接看到彩色分类图而不是灰度或黑屏。我习惯的做法是解压后先执行一次完整性检查用GDAL的gdalinfo命令看看能否正常读取元数据和波段信息。命令示例如下gdalinfo LC_2019_Yunnan.tif输出里重点看Coordinate System、Metadata、Band 1的NoData Value和Minimum/Maximum。如果Coordinate System显示为空说明坐标文件丢失需要从同源数据中找回如果NoData Value没有定义后面统计面积时就会把背景值当成真实类别这是一个极其隐蔽的坑。还有一点容易被忽略不要只保留tif文件.tfw文件只有几KB一旦丢了栅格就变成一张“无定位图片”等于白下载。3. 把数据用起来QGIS和Python读取、裁剪、统计面积的完整流程3.1 QGIS中快速查看和符号化拿到数据第一件事不要急着做统计先在QGIS里加载看一眼睛。很多人在这一步就开始翻车加载后整幅图像要么全黑要么灰蒙蒙。根源在于QGIS默认用连续渐变色带渲染而分类栅格的值是离散类别代码比如1、2、3它们之间没有物理上的连续关系默认的拉伸渲染会把它们当成灰度值显示出来自然不是你要的样子。正确操作是在“图层”面板右键当前图层选择“属性”进入“符号化”选项卡把渲染类型从“单波段灰度”切换为“单值”。然后点击“分类”按钮QGIS会自动列出所有唯一值。接着你需要根据元数据里的分类表手动修改每个类的颜色和标签。如果包内含.qml样式文件直接点击样式下拉框里的“加载样式”选择该qml文件QGIS会自动完成全部渲染配置。这里有一个关键点建议把背景值设为透明。在单值渲染的列表中找到代表NoData或背景的类别常见的是0或255把填充色设置为“透明”。如果你不透明地显示背景整幅图的边缘会出现一个巨大的黑色或白色方块干扰你对有效范围的判断也影响后续的目视检查。3.2 PythonGDAL读取和重投影图形界面适合快速查看批量处理还是得靠脚本。Python环境下最常用的是GDAL的osgeo模块或Rasterio。下面这段代码演示读取一张10米土地覆盖tif输出基本信息和唯一值列表from osgeo import gdal import numpy as np ds gdal.Open(LC_2019_Yunnan.tif) band ds.GetRasterBand(1) print(投影:, ds.GetProjection()) print(仿射变换:, ds.GetGeoTransform()) print(数据类型:, gdal.GetDataTypeName(band.DataType)) print(栅格尺寸:, ds.RasterXSize, x, ds.RasterYSize) data band.ReadAsArray() unique_vals, counts np.unique(data, return_countsTrue) for val, cnt in zip(unique_vals, counts): print(f类别 {val}: {cnt} 个像元)逻辑说明GetGeoTransform()返回六个参数顺序是左上角x坐标、像元宽度、旋转项、左上角y坐标、旋转项、像元高度。这里旋转项通常为0像元高度为负值代表栅格从左上角开始逐行向下。ReadAsArray()默认读入全图如果文件有好几个GB内存会非常吃紧建议改成band.ReadAsArray(col_offset, row_offset, col_count, row_count)分块读取或者用gdal.Warp先生成一块缩略图。np.unique能快速暴露数据中是否存在意外值比如某类的代码是128而不是0~10之间那说明数据可能被错误转换过。重投影是面积统计前必须做的事。推荐用gdal.Warpout_tif LC_2019_Yunnan_UTM47.tif gdal.Warp(out_tif, ds, dstSRSEPSG:32647, resampleAlgnear, xRes10, yRes10, formatGTiff) print(重投影完成)参数说明dstSRSEPSG:32647把坐标系直接换成UTM 47NresampleAlgnear表示最近邻重采样这一点必须写死。分类栅格的值是标号不是连续值如果使用bilinear或cubic重采样类别数值会被插值成小数比如3.7这既不是草地也不是水体整个栅格就作废了。xRes10, yRes10强制输出像元尺寸为10米这样后续面积计算可以直接用100平方米乘像素数不必再动态读取仿射参数。3.3 按行政区裁剪从省到县的矢量掩膜提取如果研究范围是某个县或多边形区域需要从全省tif中裁剪出子集。常见做法是使用gdal.Warp的cutline参数配合一个边界矢量文件。下面这段代码把云南省A县边界裁剪到土地覆盖图层上gdal.Warp( A县_LC_2019.tif, LC_2019_Yunnan_UTM47.tif, cutlineDSNameA县边界.shp, cropToCutlineTrue, dstNodata255, formatGTiff )逻辑说明cutlineDSName指定矢量边界的路径cropToCutlineTrue表示输出范围裁剪到矢量边界的外接矩形并对边界外的部分填充NoData。但请注意gdal.Warp默认裁剪的是外接矩形而不是严格按矢量边界形状。如果矢量边界是凹多边形外部但仍然在外接矩形内的区域会被填充NoData算面积时你需要排除255。如果你需要精确到矢量边界内部更合适的是用Rasterio的mask函数它支持像素级掩膜import rasterio from rasterio.mask import mask import json with rasterio.open(LC_2019_Yunnan_UTM47.tif) as src: with open(A县边界.geojson) as f: geojson json.load(f) out_image, out_transform mask(src, geojson[features], cropTrue, nodata255) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: 255 }) with rasterio.open(A县_LC_2019_precise.tif, w, **out_meta) as dst: dst.write(out_image)参数说明mask函数读入栅格和GeoJSON要素cropTrue让输出范围收缩到要素边界内。注意掩膜操作会把矢量边界外的像元设置为nodata255但边界内的NoData也会被保留为255后续统计时统跳过即可。这种方式比gdal.Warp更精确代价是内存占用更高大区域操作时尽量分块或先做好矢量简化。3.4 像元面积统计不要直接用像素数乘100统计各类别面积是使用这份数据最核心也是最容易出错的环节。先重投影到UTM或Albers再用以下代码统计from osgeo import gdal import numpy as np ds gdal.Open(LC_2019_Yunnan_UTM47.tif) band ds.GetRasterBand(1) data band.ReadAsArray() # 读取像元实际尺寸 gt ds.GetGeoTransform() pixel_width abs(gt[1]) pixel_height abs(gt[5]) pixel_area pixel_width * pixel_height # 统计唯一值及其出现次数 classes, counts np.unique(data, return_countsTrue) results [] for cls, cnt in zip(classes, counts): if cls 255: # 跳过NoData continue area_m2 cnt * pixel_area area_km2 area_m2 / 1e6 results.append((cls, area_km2)) print(f类别 {cls}: {area_km2:.2f} km²)逻辑说明gt[1]是经向像元尺寸gt[5]是纬向像元尺寸在UTM投影下两者通常都是10米但你不能假设一定是10因为重投影时四舍五入可能导致9.999或10.001。最稳妥的方式就是这里动态读取。np.unique返回两个数组顺序一致用zip配对即可。注意我把NoData值定为255如果某个有效类别恰好也是255就得改成其他值这是为什么在上一步重投影时设置dstNodata255之前先检查原始分类最大值的道理。如果想输出成CSV加一个np.savetxt或pandas即可。这里建议用pandasimport pandas as pd from osgeo import gdal import numpy as np ds gdal.Open(LC_2019_Yunnan_UTM47.tif) data ds.GetRasterBand(1).ReadAsArray() gt ds.GetGeoTransform() pixel_area abs(gt[1]) * abs(gt[5]) classes, counts np.unique(data, return_countsTrue) df pd.DataFrame({ class: classes, pixel_count: counts, area_km2: counts * pixel_area / 1e6 }) df df[df[class] ! 255] df.to_csv(Yunnan_LC2019_area.csv, indexFalse) print(df)这样处理之后输出的面积表格可以直接和统计年鉴对照。需要说明的是栅格分类天然存在混合像元10米分辨率下梯田边界、道路两侧的细碎地物仍然会被错分所以面积统计不是物理实测而是一个基于最大似然分类的估算值使用在趋势分析上完全没问题。4. 土地利用转移矩阵从静态分类到动态变化分析4.1 转移矩阵的计算逻辑这份数据本身是2019年单时期的但做土地利用变化的同行通常会把2019年作为基准期再找一期相邻年份的数据构建转移矩阵。转移矩阵的基本思想是统计每个像元在初期和末期的类别组合并把结果填入二维交叉表。矩阵的行表示初期的类别列表示末期的类别对角线上的数值表示“未变化面积”非对角线上的数值表示“从某种类型转变为另一种类型”的面积。举个例子假设初期只有草地3和耕地4末期也只有这两类。矩阵中第一行第一列是“草地转为草地”的面积第一行第二列是“草地转为耕地”的面积。矩阵中最关键的一步是保证两期栅格的像元一一对应。要求坐标系、分辨率、裁剪范围完全一致否则一个像元对不上另一个像元的状况矩阵里会出现大量虚假变化。由于原始数据已经是2019年你需要准备第二期。常见做法是下载同源产品的2020年或2021年图幅再按同一个裁剪范围处理。如果直接拿不同来源的两期数据做矩阵分类代码不一致会让矩阵无法解释必须先重分类统一编码。这一步比矩阵计算本身更耗时但也更重要。4.2 用Python构建2019年分类的转移矩阵当两期数据准备好后下面的代码构建一个最直观的转移矩阵。为方便理解这里使用一个循环实现虽然效率不是最高但逻辑清晰适合小面积试点。from osgeo import gdal import numpy as np def read_band(path): ds gdal.Open(path) band ds.GetRasterBand(1) return band.ReadAsArray() data_2019 read_band(LC_2019_Yunnan_UTM47.tif) data_2020 read_band(LC_2020_Yunnan_UTM47.tif) # 保证两期数据形状一致且在相同位置 assert data_2019.shape data_2020.shape, 两期栅格尺寸不一致 # 构造有效像元掩膜 valid (data_2019 ! 255) (data_2020 ! 255) a data_2019[valid] b data_2020[valid] # 收集所有类别 classes np.unique(np.concatenate([a, b])) index {cls: i for i, cls in enumerate(classes)} matrix np.zeros((len(classes), len(classes)), dtypenp.int64) for a_val, b_val in zip(a, b): matrix[index[a_val], index[b_val]] 1 print(矩阵行2019年矩阵列2020年) for i, cls in enumerate(classes): print(f{cls}: {matrix[i, :]})逻辑说明valid掩膜同时排除了两期中的NoData。np.concatenate是为了从两期里收集全部类别编号保证矩阵行列覆盖所有可能值。循环遍历每个有效像元把组合计数写入矩阵对应位置。dtypenp.int64务必保留如果云南全省像元大约10亿累计计数很容易超过int32上限。对于大型影像这个循环会非常慢。改进办法是使用np.add.attrans (np.searchsorted(classes, a) * len(classes) np.searchsorted(classes, b)) matrix_opt np.zeros((len(classes), len(classes)), dtypenp.int64) np.add.at(matrix_opt, trans, 1)searchsorted把类别值映射成从0开始的索引然后通过一维索引展开二维矩阵的位置np.add.at在扩展位置进行累加。这种方法比Python循环快两个量级适合全云南处理。4.3 结果解读与常见误区转移矩阵得出后常见的解读错误有两个。第一是把行和列的方向搞反。建议在输出表格的注解里直接写“行是2019年列是2020年”并且导出到CSV时保留类别码不用阿拉伯数字代替类别名称。第二是直接把两期分类图相减例如data2020 - data2019这种做法的结果完全无法解释因为5变成3和3变成5都会得到2而且类别编号没有数值含义减法得到的数字根本不是物理量。如果要让矩阵变成面积把每个计数乘以像元面积后放在矩阵里。计算面积时用3.4节的动态像元面积而不是手动假定100。做完后可以输出一个面积转移表再用热力图展示哪种转换最频繁。我在实际项目中见过一个常见陷阱周边省份的边界像元因为云掩膜存在导致NoData不重合矩阵里出现大量“有到无”的变化这时候必须先把两期的云掩膜并集作为有效区再计算矩阵否则数十万个虚假变化会让结果失真。5. 避坑记录10米土地覆盖数据的五个典型坑和一个后悔药这一章写我和同行踩过的五个具体坑按“现象→原因→解决”的顺序列出来后面还有一个后悔药希望你能绕过这些地方。5.1 黑屏和花屏渲染方式不对现象加载tif后画面全黑或者只有少数几个颜色完全看不出地形和地物的边界。原因分类栅格的值是离散代码不是连续灰度QGIS默认使用彩色拉伸渲染把分类数字当成灰度值或连续色带显示导致没有过渡色的区域表现为黑色。解决在图层属性“符号化”中选择“单值”渲染点击“分类”让软件自动枚举值。如果有.qml样式文件直接加载样式。如果还是不行检查数据是否损坏用gdalinfo的Minimum... Maximum...看是否能读取有效范围如果读取失败大概率栅格文件本身出了问题。5.2 面积统计结果偏小5%~10%现象用原始tif统计云南林地面积和《云南省统计年鉴》或三调数据对比总是少一块差距有时接近十分之一。原因原始WGS84坐标下像元在经度方向的实际长度随纬度升高而变短。云南大部分处于北纬21°~29°按北纬25°计算经度方向一个像元大约只能对应9.05米而不是10米纬度方向约11.04米实际像元面积约100平方米但分布不均衡。如果你直接用0.00009度换算成10米误差就出来了。解决面积统计前强制重投影至UTM 47N或Albers等积投影并设置输出像元尺寸为10米。用gdal.Warp加xRes10, yRes10再用动态读取仿射参数计算面积。从那以后我看到任何一张土地覆盖tif第一件事就是看坐标系如果显示EPSG:4326我绝不会直接用像元数乘以100。5.3 NoData被当成有效类别现象类别统计时出现一个面积几万平方公里的“类别0”占掉总量一半整个统计表失去意义。原因数据生产时把无效区域或云掩膜设置为0但很多软件不会自动识别0为NoData。有些数据集则把NoData设为255你在原始WGS84坐标下可能没注意。解决在QGIS里用“标识”工具点击影像外框看该区域的值是什么然后统一重编码。最好在重投影时用dstNodata255并在后续统计代码里显式跳过255。如果原始NoData就是0则跳过0。千万别盲目跟网上教程默认“0是背景”一定要自己确认。5.4 不同来源数据图例混用现象把两个不同年份或不同来源的tif放在一起发现同一地点类别变化频繁很多区域在“林地”和“灌木”之间来回跳变。原因不同产品分类体系不同。比如某个产品把“疏林”归为林地另一个产品把它归为灌木编码更是天差地别。直接用第二套图例去解释第一套数据必然混乱。解决使用前先分别读取唯一值做成两张“原始编码→统一分类”的重分类映射表。在Python中用np.copy做重编码例如import numpy as np def recode(data, mapping): out data.copy() for old, new in mapping.items(): out[data old] new return out mapping_2019 {1: 10, 2: 20, 3: 30, 4: 40, 5: 50} mapping_2020 {2: 10, 3: 20, 5: 30, 7: 40, 9: 50} data_2019_recoded recode(data_2019, mapping_2019) data_2020_recoded recode(data_2020, mapping_2020)这里把两期类别统一到同一套编号体系这样才能保证转移矩阵有实际意义。这个步骤最费时间但它是变化分析可信度的根基。5.5 裁剪后边界出现锯齿和杂斑现象用矢量边界裁剪后边界地带出现很多细碎的类别突变仔细观察会发现某些像元在边界外但被填入了类别值。原因裁剪本质是按像元对齐矢量边界穿过像元时该像元会被判定为“部分落入”。gdal.Warp的cropToCutline默认外接矩形填充如果使用-cutline但未-crop_to_cutline输出不裁剪如果用cropToCutlineTrue边界外的像元处理取决于重采样方法有时会产生边缘噪点。解决对于边界问题接受栅格本身的锯齿。不要试图把栅格矢量化成光滑多边形那会丢失像元精度。如果要做精确边界统计先用Rasteriomask进行像素级掩膜再统计不要在边缘像元上反复纠结。对于杂斑如果某类面积占比小于0.1%且形态上不符合地物分布可以用多数滤波或形态学开闭运算清理但必须评估是否影响真实细小地物。5.6 后悔药及时备份原始rar包现象裁剪或重投影过程中误覆盖原始tif想重置某个区域的分类时只能重新下载整个压缩包耗时几天。原因很多人解压后为了节省空间删掉了rar包或者直接在工作目录里把原始tif作为输出文件名覆盖。解决下载后先创建一个“原始备份”目录把rar包复制进去并设置文件夹只读。工作目录放一份解压副本所有脚本处理都指向工作副本。我习惯在脚本开头强制检查原始文件是否存在import os raw_path 原始备份/LC_2019_Yunnan.tif assert os.path.exists(raw_path), 请先检查原始数据是否存在这个后悔药听起来太基础但我在实际中已经看到不止一个同事因为误删原始文件加班重新下载数据。别让最基础的事成为最大的坑。6. 进阶验证用随机样本点检验数据精度把“10米”变成可信数字无论这份分类数据来自哪个公开产品拿到手后最好自己做一轮精度验证尤其是你要把它作为成果交付时。最常用的方法是生成随机检验点用高分辨率影像逐点判读真实类别再与栅格提取值比较计算整体精度和Kappa系数。十米精度听起来可信但实际分类结果在云南山区往往比平原地区虚高验证一下才能给你兜底。在QGIS中生成本地随机点按类别分层抽样。每个类别至少抽80个点样本太少混淆矩阵会空行或空列Kappa值也会失真。导出CSV后用sklearn计算混淆矩阵import pandas as pd from sklearn.metrics import confusion_matrix, accuracy_score, cohen_kappa_score df pd.read_csv(validation_samples.csv) cm confusion_matrix(df[reference], df[predicted], labelssorted(df[reference].unique())) oa accuracy_score(df[reference], df[predicted]) kappa cohen_kappa_score(df[reference], df[predicted]) print(混淆矩阵:\n, cm) print(f总体精度: {oa:.3f}, Kappa: {kappa:.3f})参数说明reference列是人工判读的类别predicted列是栅格上提取的类别。labels参数让矩阵行列覆盖所有可能类别避免某类缺失导致列数不足。cohen_kappa_score会计算Kappa系数一般Kappa大于0.75说明可用于业务制图0.6到0.75之间可用于宏观趋势小于0.6就要谨慎了。验证结果出来后查看哪一类最容易混淆。云南最常见的混淆是“草地”和“灌木”以及高海拔裸岩和建筑。如果你发现某一类的用户精度特别低可以在后续统计中把它并进更粗略的一级类别例如把“灌丛”和“草地”合并成“草灌地”。这样损失一点细节但换来更可信的数字。那一次我在云南某项目里帮客户统计不同坡度带的林地面积原始数据Kappa只有0.58在我合并了草灌类后提升到0.74客户认账了。从此以后我每次拿到新的10米土地覆盖数据都会先抽出两三百个点验证一轮再决定是否用于业务。这个习惯帮我避开了很多“看起来精度高然并卵”的公开数据。希望这份拆解能帮你把10米数据用到实处也希望你别跳过验证这一步。本文还有配套的精品资源点击获取