简介：本资源为2019年10米空间分辨率的云南省土地覆盖与土地利用数据，面向地理信息、遥感、生态评估及国土空间规划等方向的研究者与从业者。数据基于哨兵影像与深度学习方法制作，经重投影由墨卡托直角坐标系转为WGS84地理坐标系，并按最新省市级行政边界裁剪，涵盖耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪冰等十类地物，可直接用于云南省各地市的土地利用分析与制图。压缩包共112个文件，约115.71MB，包含16个tif栅格数据、16个tfw坐标文件、16个dbf属性表、16个xml元数据、16个png预览图、16个cpg编码文件及16个xlsx统计表，按州市分文件组织，覆盖德宏、曲靖、文山、玉溪、大理、红河、丽江、迪庆、怒江、昭通等地。目前已有327人学习下载，适合需要快速获取云南省精细土地利用底图、开展空间统计或作为模型输入的用户。
1. 2019年10m精度云南省土地覆盖土地利用数据：这份栅格包到底能干什么
如果你手头正好有一个云南省的行政区划矢量边界，又需要做2019年前后的耕地变化、建设用地扩张或者林地扰动分析，那这份「2019年10m精度云南省土地覆盖土地利用.rar」大概率能省掉你至少两周的预处理时间。它本质上是一套已经裁切到云南省范围的10米空间分辨率土地覆盖栅格产品，按年份打包，解压后是常见的GeoTIFF格式，每个像元对应一个地类编码。10米分辨率意味着什么？云南省面积约39.4万平方公里，在这个精度下，一个像元代表地面上10米×10米的方块，全省大约有近40亿个像元——这个量级既能看清坝区里的建设用地斑块，又不至于让普通工作站直接爆内存。它适合做省域尺度的土地利用结构统计、多期变化检测、生态红线内的人类活动扰动评估，也适合作为遥感分类、地表参数反演的底图或验证数据。不适合拿来做地块级确权或者单栋建筑提取，那个精度需求得用亚米级影像。这份资源解决的核心问题就一个：让你跳过从原始影像分类到全省拼接的漫长流程，直接进入分析和出图阶段。
2. 数据拆包与坐标系统确认：从压缩包到可分析栅格
2.1 解压后的目录结构与文件格式
拿到压缩包后，我一般先不急着往GIS软件里拖，而是用命令行把目录树看清楚。常见的打包习惯是年份加省份作为顶层文件夹，里面再按地类或分幅存放。下面这段bash命令能帮你快速摸清家底：
# 创建独立工作目录，避免污染原始压缩包 mkdir -p /work/yunnan_lulc_2019 && cd /work/yunnan_lulc_2019 # 解压rar，注意路径按实际存放位置调整 unrar x /data/downloads/2019年10m精度云南省土地覆盖土地利用.rar # 查看解压后的目录层级和文件大小 find . -maxdepth 3 -type f -name "*.tif" -exec ls -lh {} \; # 统计tif文件数量，确认是否分幅 find . -name "*.tif" | wc -l这段命令的逻辑很直接：先建一个干净的工作目录，把压缩包内容释放进去，然后用find列出所有tif文件及其大小。如果返回的是单个几GB的大文件，说明全省已经镶嵌好了；如果返回几十个甚至上百个小文件，那就是分幅存储，后续需要先做镶嵌。参数上，-maxdepth 3控制目录遍历深度，防止在深层嵌套里浪费时间；-name "*.tif"只匹配栅格文件，忽略可能附带的说明文档或样式表。常见做法是解压后先检查有没有配套的元数据文件，比如.prj、.tfw或者一个readme文本，这些能直接告诉你坐标系和地类编码含义，比事后猜要靠谱得多。
2.2 用GDAL确认投影与地类编码
10米精度的土地覆盖产品，国内常见的有几套不同的分类体系，地类编码从1到几十不等，投影也可能是UTM或者Albers。不确认这两项，后面做面积统计时数字会差出十万八千里。用GDAL的gdalinfo命令可以一次性把关键信息拉出来：
# 查看单个tif的元数据，重点关注坐标系和波段数 gdalinfo -stats Yunnan_2019_LULC.tif | head -60 # 如果分幅，批量输出每个文件的投影信息到文本，便于核对是否一致 for f in *.tif; do echo "=== $f ==="; gdalinfo "$f" | grep -E "PROJCRS|GEOGCRS|Pixel Size"; done > proj_check.txtgdalinfo输出里，PROJCRS或GEOGCRS字段就是坐标系定义，Pixel Size告诉你像元大小是不是10米×10米。如果发现有的分幅是地理坐标系（度为单位），有的又是投影坐标系（米为单位），那就得先统一投影再镶嵌，否则接边处会错位。-stats参数会顺便计算波段统计值，对于单波段的地类栅格，你能看到最小值和最大值，大致判断编码范围。我一般会把gdalinfo的输出存一份到项目目录里，后面写报告时直接引用，省得反复查。
2.3 地类编码对照与重分类准备
拿到栅格只是第一步，真正要分析时得知道每个数字代表什么。这类产品通常附带一个分类对照表，常见的一级类包括耕地、林地、草地、水域、建设用地、未利用地等。如果压缩包里没有现成的对照表，就需要根据数据来源推断——国内10米产品多基于Sentinel-2或国产高分影像分类，分类体系往往参照《土地利用现状分类》国标。下面是一个典型的重分类Python脚本框架，用rasterio和numpy把原始编码映射成你需要的类别：
import rasterio import numpy as np # 打开原始栅格 with rasterio.open('Yunnan_2019_LULC.tif') as src: lulc = src.read(1) # 读取第一波段 profile = src.profile.copy() nodata = src.nodata # 假设原始编码：1耕地 2林地 3草地 4水域 5建设用地 6未利用地 # 重分类为：1生态用地 2农业用地 3建设用地 4水域 reclass_map = { 1: 2, # 耕地 -> 农业用地 2: 1, # 林地 -> 生态用地 3: 1, # 草地 -> 生态用地 4: 4, # 水域 -> 水域 5: 3, # 建设用地 -> 建设用地 6: 1, # 未利用地 -> 生态用地（按需调整） } # 构建查找表，加速映射 max_val = int(np.max(lulc)) lut = np.zeros(max_val + 1, dtype=np.uint8) for k, v in reclass_map.items(): lut[k] = v reclassed = lut[lulc] # 保持nodata不变 if nodata is not None: reclassed[lulc == nodata] = nodata # 写出重分类结果 profile.update(dtype=rasterio.uint8, count=1, compress='lzw') with rasterio.open('Yunnan_2019_reclass.tif', 'w', **profile) as dst: dst.write(reclassed, 1)这段代码的核心是查找表映射，比逐像元循环快几个数量级。参数上，profile.copy()保留了原始的地理变换和投影信息，确保输出栅格能和其它图层对齐；compress='lzw'能显著减小文件体积，对全省数据尤其有用。注意重分类前一定要确认nodata值，否则背景值可能被误映射成某个地类，导致面积统计虚高。常见做法是先用np.unique看一下原始编码的实际分布，再决定映射关系，不要照搬网上的对照表。
3. 面积统计与变化检测：从栅格到可汇报的数字
3.1 按行政区统计各地类面积
拿到重分类结果后，最常被问到的就是“云南省2019年耕地有多少公顷”。如果手头有县级或乡镇级的行政边界矢量，可以用rasterstats做分区统计，一次性输出每个行政单元内各地类的像元数和面积。下面这个脚本假设你已经有了一个yunnan_counties.shp：
import geopandas as gpd import rasterio from rasterstats import zonal_stats import pandas as pd # 读取行政边界 counties = gpd.read_file('yunnan_counties.shp') # 确认边界和栅格投影一致，不一致则重投影 with rasterio.open('Yunnan_2019_reclass.tif') as src: raster_crs = src.crs if counties.crs != raster_crs: counties = counties.to_crs(raster_crs) # 定义地类名称映射 class_names = {1: '生态用地', 2: '农业用地', 3: '建设用地', 4: '水域'} # 对每个地类分别做分区统计 results = [] for code, name in class_names.items(): # 创建只含目标地类的二值栅格用于统计 stats = zonal_stats( counties, 'Yunnan_2019_reclass.tif', categorical=True, # 按类别统计像元数 nodata=0, affine=rasterio.open('Yunnan_2019_reclass.tif').transform ) # 提取每个县的像元数并换算面积（10m像元=100平方米=0.01公顷） for idx, stat in enumerate(stats): pixel_count = stat.get(code, 0) area_ha = pixel_count * 0.01 results.append({ 'county': counties.iloc[idx]['NAME'], 'class': name, 'area_ha': area_ha }) df = pd.DataFrame(results) df.to_csv('yunnan_2019_lulc_area_by_county.csv', index=False, encoding='utf-8-sig')这里的关键参数是categorical=True，它让zonal_stats返回每个类别在区域内的像元计数，而不是默认的均值或总和。像元面积换算要注意：10米分辨率下单个像元面积是100平方米，即0.01公顷，所以像元数乘以0.01就是公顷数。如果行政区边界跨了多个分幅，zonal_stats会自动处理镶嵌后的栅格，前提是你已经把所有分幅合并成了一个文件。输出CSV时用utf-8-sig编码，能避免Excel打开中文乱码，这是血泪经验。
3.2 两期数据做变化检测的栅格运算
如果你手里还有另一年份的同类数据，就可以做变化检测了。最直接的方式是栅格代数运算，把两期重分类结果相减，非零像元就是发生变化的位置。下面用rasterio读两个年份的数据，生成变化类型图：
import rasterio import numpy as np # 读取两期重分类栅格 with rasterio.open('Yunnan_2019_reclass.tif') as src: lulc_2019 = src.read(1) profile = src.profile.copy() nodata = src.nodata with rasterio.open('Yunnan_2022_reclass.tif') as src: lulc_2022 = src.read(1) # 确保两期数据形状一致 assert lulc_2019.shape == lulc_2022.shape, "两期栅格尺寸不一致，需先重采样对齐" # 变化编码：年份A类别 * 10 + 年份B类别，便于追溯变化方向 change = lulc_2019.astype(np.int16) * 10 + lulc_2022.astype(np.int16) # 无变化像元置为0 change[lulc_2019 == lulc_2022] = 0 # 保留nodata if nodata is not None: change[(lulc_2019 == nodata) | (lulc_2022 == nodata)] = nodata # 统计主要变化类型的面积 unique, counts = np.unique(change[change > 0], return_counts=True) for u, c in zip(unique, counts): from_class = u // 10 to_class = u % 10 print(f"从 {from_class} 变为 {to_class}：{c * 0.01:.2f} 公顷") # 写出变化栅格 profile.update(dtype=rasterio.int16, nodata=nodata, compress='lzw') with rasterio.open('Yunnan_lulc_change_2019_2022.tif', 'w', **profile) as dst: dst.write(change, 1)变化编码用“前一期类别×10+后一期类别”的方式，能在一张图里同时表达变化方向和类型，比如21表示从农业用地变成生态用地。参数上，astype(np.int16)是为了防止乘法溢出，原始uint8乘以10后最大到60，int16绰绰有余。assert那行是后悔药，两期数据如果行列数不一致，后面的运算会直接报错或者产生错位结果，提前拦住比事后排查强。统计部分只输出变化像元，无变化区域被置零，这样面积汇总时不会把稳定地类混进来。
3.3 用QGIS快速出图与结果核验
命令行跑完统计，最终汇报往往需要一张像样的图。QGIS加载这类栅格很直接，但有几个设置能让你少走弯路。打开QGIS后，把重分类后的tif拖进去，右键图层属性，在“符号化”里选择“唯一值”渲染，然后点击“分类”按钮，它会自动读取所有唯一像元值。如果地类编码是1到6，你会看到六个类别，手动给每个类别指定颜色和标签即可。核验环节我一般会做两件事：一是用“识别要素”工具随机点几个像元，看属性值是否和预期地类吻合；二是叠加行政边界，检查边界处有没有明显的裁切错位。如果发现边界外还有数据，说明原始栅格没有按省界精确裁切，需要用“按掩膜提取”工具重新裁一次。出图时注意设置打印布局的分辨率不低于300dpi，否则10米精度的细节在打印稿上会糊成一片。
4. 避坑与常见问题排查：那些让你白跑一晚上的细节
4.1 现象：面积统计结果比官方公布数据大出一截
原因通常有两个：一是nodata值没有正确处理，背景区域被当成了某个地类；二是坐标系是地理坐标系，像元大小以度为单位，直接按10米换算面积会严重失真。解决方法是先用gdalinfo确认投影类型，如果是地理坐标系，必须先用gdalwarp重投影到等面积投影（如Albers）再统计。nodata方面，用np.unique检查原始栅格的值分布，找到背景值并在统计前排除。
4.2 现象：两期数据相减后出现大量异常变化斑块
这多半是两期数据的分类体系不一致导致的。比如2019年把草地编码为3，2022年编码为4，相减后就会产生虚假变化。解决方法是先分别确认两期的地类编码对照表，统一重分类到同一套体系后再做变化检测。如果对照表缺失，可以抽样对比两期影像的原始光谱特征，辅助判断编码含义。
4.3 现象：QGIS加载栅格后显示全黑或全白
这是拉伸设置的问题。10米土地覆盖栅格的像元值通常是离散的整数，默认的线性拉伸会把它们压缩到极窄的显示范围。解决办法是在图层属性里把“最小/最大值”手动设为实际编码范围，或者直接用“唯一值”渲染。如果还是不对，检查波段是否被正确识别——有些产品把地类码放在第二个波段，第一个波段是质量评价码。
4.4 现象：分区统计时提示内存不足
全省10米栅格数据量不小，如果一次性读入内存做分区统计，普通16GB内存的机器很容易崩。解决办法是分块处理：用rasterio的窗口读取功能，按行政区边界的外接矩形逐块读取，统计完再合并结果。或者先用gdal_translate把栅格切分成多个小文件，再逐个统计。我一般会在统计前用gdaladdo建立金字塔，能显著加快QGIS里的浏览速度。
4.5 现象：变化检测结果里建设用地减少
在快速城镇化地区，建设用地理论上只增不减。如果出现减少，大概率是两期数据中有一期的建设用地被误分为其他地类，或者两期数据的获取季节不同导致光谱差异被误判。解决办法是叠加原始影像抽查几个“减少”的斑块，确认是分类错误还是真实拆除。如果是分类错误，考虑引入更高精度的参考数据做局部修正。
5. 进阶用法：把土地覆盖栅格接入时序分析流程
这份2019年数据如果只做静态统计，多少有点浪费。我习惯把它作为时序分析的起点，配合后续年份的数据做趋势检验。一个实用的技巧是用xarray把多期栅格组织成时间维数组，然后按像元计算变化斜率。下面这段代码演示如何把多个年份的栅格堆叠成一个数据集，并计算耕地面积的逐年变化：
import xarray as xr import rioxarray import numpy as np # 假设有2019到2022四期重分类栅格 years = [2019, 2020, 2021, 2022] arrays = [] for y in years: da = rioxarray.open_rasterio(f'Yunnan_{y}_reclass.tif', masked=True) da = da.squeeze('band', drop=True) # 去掉波段维 da = da.assign_coords(year=y) arrays.append(da) # 沿时间维合并 cube = xr.concat(arrays, dim='year') # 计算耕地（编码2）像元数逐年变化 cropland = (cube == 2).sum(dim=['x', 'y']) print(cropland.values) # 输出每年的耕地像元数 # 计算每个像元上耕地存在与否的线性趋势 cropland_mask = (cube == 2).astype(np.float32) trend = cropland_mask.polyfit(dim='year', deg=1) slope = trend.sel(degree=1)['polyfit_coefficients'] # slope为正表示该像元从非耕地转为耕地，为负则相反 slope.rio.to_raster('cropland_trend_2019_2022.tif')这段代码的关键在于xr.concat沿年份维合并，形成一个带时间坐标的数据立方体。polyfit对每个像元做一阶线性拟合，斜率就是变化速率。参数上，masked=True让nodata自动变成NaN，避免参与运算；deg=1表示线性拟合，如果想捕捉非线性趋势可以改成2，但解释起来更复杂。输出的趋势栅格可以直接在QGIS里用发散色带渲染，一眼就能看出哪些区域在退耕，哪些在复耕。我一般会把这个趋势图和行政区边界叠加，快速定位变化热点乡镇，再针对性去查原始影像。从那以后我每次拿到新的土地覆盖年份包，都会先跑一遍这个趋势脚本，确认数据连续性没问题再往下做。希望帮到你。
本文还有配套的精品资源，点击获取
