简介：本资源为2019年10米空间分辨率的云南省土地覆盖土地利用数据，面向地理信息、遥感、生态评估及国土空间规划等方向的研究者与从业者。数据基于10米哨兵影像，采用深度学习方法制作，共分耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、雪/冰等十类，并已由原始墨卡托投影重投影为WGS84地理坐标系，按最新省市级行政边界裁剪，覆盖云南省各地级市与自治州。压缩包共112个文件，约115.71MB，包含16个tif栅格主数据、16个tfw坐标文件、16个dbf属性表、16个xml元数据、16个cpg编码文件、16个png预览图及16个xlsx统计表，便于直接加载、符号化与面积统计。目前已有327人学习下载，适合用于土地利用变化分析、生态遥感建模与区域规划制图等场景。
1. 10m精度云南省土地覆盖数据：从压缩包到可分析图层的完整路径
拿到「2019年10m精度云南省土地覆盖土地利用.rar」这个包的人，十有八九是冲着两件事来的：一是云南这个地形气候都极端的省份，10米分辨率到底能看清什么；二是手里这份数据能不能直接喂进自己的模型或统计流程。云南从热带雨林到高山草甸垂直分布，30米产品在怒江峡谷、滇中高原坝子这些地方经常把耕地和林地混成一团，10米精度意味着地块级边界终于能分开了。这个包本质上是一套栅格分类成果，核心工作不是「打开看看」，而是把它变成坐标系正确、分类体系明确、能参与面积统计和变化检测的可用图层。适合做国土空间规划、生态评估、农业遥感以及需要省级高分辨率底图的研究者。下面按「先搞懂数据长什么样，再动手处理，最后避开那些让人返工的坑」这条线走一遍。
2. 先搞清楚10m土地覆盖栅格到底装了什么
2.1 分类体系与编码：别把DN值当面积用
10米土地覆盖产品绝大多数走的是二级分类体系，一级类通常包括耕地、林地、草地、水域、建设用地、未利用地这几大项，二级类会细到水田、旱地、有林地、灌木林、疏林地、高覆盖草地、中覆盖草地、低覆盖草地、河渠、湖泊、水库坑塘、城镇用地、农村居民点、裸地等。栅格里的每个像元存的是一个整数编码，不是反射率也不是面积。很多人第一次打开看到值域是1到十几，就以为是某种连续指数，直接拿去做归一化，结果全错。
云南的特殊性在于，同一套全国分类体系落到这里，会出现大量「混合像元」被强行归到某一类的情况。比如橡胶林在编码上可能归到有林地，但它的物候和天然林完全不同；滇池周边的设施农业可能被归到耕地，实际是温室大棚。所以拿到数据第一件事是找到配套的分类对照表，确认每个编码对应的地类名称，而不是凭经验猜。
常见做法是：如果压缩包里带了说明文档或.txt编码表，直接读；如果没有，就按该产品公开的全国分类体系对照，但要在报告里注明「编码解释依据XX体系，未获官方随包说明」。这一步不做，后面面积统计出来的数字没法解释。
2.2 坐标系与投影：面积算不准的根源
10米栅格产品在国内通常有两种投影：一种是地理坐标系（如CGCS2000或WGS84），单位是度；另一种是投影坐标系（如Albers等面积投影），单位是米。如果拿到的是地理坐标系的栅格，直接算面积会得到一个以「平方度」为单位的荒谬数字。云南跨了多个纬度带，地理坐标系下高纬度像元实际面积比低纬度小，不做投影转换就统计面积，误差能到百分之几甚至更多。
判断方法很简单：用gdalinfo看一眼栅格的地理变换参数。如果像素大小是0.0001度左右，那就是地理坐标系；如果是10米左右，那就是投影坐标系。这一步决定了后面要不要做重投影。
gdalinfo 2019_yunnan_landcover.tif重点看输出里的Coordinate System和Pixel Size。如果Pixel Size显示的是类似(0.0001, 0.0001)，说明是地理坐标，需要转到等面积投影再做面积统计。如果显示(10, 10)，单位是米，可以直接用。
提示：不要用「度」直接乘111000来估算面积，云南东西跨度和南北跨度都大，这种粗略换算在省级尺度上误差不可接受。
2.3 数据完整性检查：云南边境和云层遮挡区
云南西北部与西藏、四川交界处，以及南部边境一带，部分产品会存在数据缺失或填充值。常见填充值是0、-9999或255。如果不做掩膜处理，这些填充值会被当成某一类参与统计，导致未利用地面积虚高。
检查方法是统计栅格值分布：
import rasterio import numpy as np with rasterio.open('2019_yunnan_landcover.tif') as src: data = src.read(1) unique, counts = np.unique(data, return_counts=True) for u, c in zip(unique, counts): print(f"值 {u}: {c} 像元, 占比 {c/data.size*100:.2f}%")如果看到某个异常值占比很小且集中在边境线附近，基本可以判定是填充值。处理方式是在后续分析前用np.where把它设为np.nan或单独掩膜掉。这一步不做，后面做变化检测时会出现「2019年有值、2020年无值」的假变化。
3. 从压缩包到可分析图层：解压、裁剪与重投影
3.1 解压与文件结构确认
.rar包在Linux下用unrar或7z解压，Windows下用7-Zip或WinRAR。解压后常见文件包括：一个或多个.tif栅格文件、一个.txt或.xlsx分类说明、有时还有一个.prj投影文件。如果只有.tif没有说明文件，需要根据产品名称去查对应分类体系。
# Linux下解压 mkdir yunnan_landcover_2019 unrar x "2019年10m精度云南省土地覆盖土地利用.rar" yunnan_landcover_2019/ ls -lh yunnan_landcover_2019/解压后先确认栅格文件是否分幅。省级10米产品有时会按图幅切分，比如按1:100万标准分幅。如果是多幅，需要先镶嵌成一整幅再裁剪。判断方法：看每个.tif的范围是否只覆盖云南的一部分。
3.2 按云南省行政边界裁剪
如果拿到的是全国或大区域数据，第一步是按云南边界裁剪。用gdalwarp配合一个云南省的矢量边界（.shp或.geojson）即可。注意裁剪时要用-crop_to_cutline，否则输出范围还是原始范围，只是边界外被设为NoData。
gdalwarp -cutline yunnan_boundary.shp \ -crop_to_cutline \ -dstnodata 0 \ -of GTiff \ input_landcover.tif \ yunnan_landcover_2019.tif参数说明：-cutline指定裁剪矢量；-crop_to_cutline让输出范围紧贴矢量边界；-dstnodata 0把边界外设为0，方便后续掩膜；-of GTiff指定输出格式。如果矢量边界和栅格坐标系不一致，gdalwarp会自动做重投影，但建议先把矢量转到和栅格一致的坐标系，避免意外。
3.3 重投影到等面积坐标系
如果裁剪后的栅格还是地理坐标系，需要转到等面积投影。国内常用Albers投影，参数为：中央经线105°E，双标准纬线25°N和47°N，椭球体CGCS2000或WGS84。用gdalwarp一行搞定：
gdalwarp -t_srs "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs" \ -tr 10 10 \ -r near \ -of GTiff \ yunnan_landcover_2019.tif \ yunnan_landcover_2019_albers.tif-t_srs指定目标投影；-tr 10 10指定输出分辨率为10米；-r near表示用最近邻重采样，分类栅格必须用最近邻，不能用双线性或立方卷积，否则会出现不存在的类别值。这一步做完，栅格每个像元就是10米×10米，面积统计可以直接用像元数乘以100平方米。
注意：重投影后一定要再跑一次值分布检查，确认没有因为重采样产生新的异常值。最近邻虽然不会创造新值，但如果原始数据有填充值，重投影后填充值可能扩散到边界外。
4. 面积统计与分类提取的实操细节
4.1 按地类统计面积：从像元数到平方公里
面积统计的核心逻辑是：每个像元面积=10m×10m=100平方米，统计每个编码的像元数，乘以100再除以10^6得到平方公里。用rasterio和numpy可以几行搞定：
import rasterio import numpy as np import pandas as pd with rasterio.open('yunnan_landcover_2019_albers.tif') as src: data = src.read(1) nodata = src.nodata # 掩膜掉NoData和填充值 valid = data != nodata data_valid = data[valid] unique, counts = np.unique(data_valid, return_counts=True) # 像元面积100平方米，转平方公里 area_km2 = counts * 100 / 1e6 df = pd.DataFrame({'code': unique, 'area_km2': area_km2}) df['percent'] = df['area_km2'] / df['area_km2'].sum() * 100 print(df.sort_values('area_km2', ascending=False))这段代码的关键点：先掩膜NoData，再统计；面积换算用100平方米/像元；最后算百分比。如果分类说明表里有编码对应的地类名称，用pd.merge关联上去，输出就是可读的面积表。
参数说明：src.nodata可能为None，如果为None需要手动指定填充值，比如data != 0。云南边境填充值常见为0，但0也可能对应某个真实地类（如未利用地），所以最好结合值分布判断。
4.2 提取单一地类做专题图
做生态评估或规划时，经常需要单独提取某一类，比如林地或水域。用numpy布尔索引即可：
import rasterio import numpy as np with rasterio.open('yunnan_landcover_2019_albers.tif') as src: data = src.read(1) profile = src.profile # 假设林地编码为2（有林地）、3（灌木林）、4（疏林地） forest_mask = np.isin(data, [2, 3, 4]).astype(np.uint8) profile.update(dtype=rasterio.uint8, nodata=255) with rasterio.open('yunnan_forest_2019.tif', 'w', **profile) as dst: dst.write(forest_mask, 1)np.isin用来匹配多个编码；输出为0/1二值栅格，1表示林地，0表示非林地。profile.update里把数据类型改为uint8，NoData设为255，避免和0/1混淆。如果要做更细的二级类提取，把编码列表换成对应值即可。
4.3 与行政区划叠加统计
省级数据经常需要落到州市甚至县级。做法是用rasterstats库的zonal_stats，或者用geopandas加rasterio做分区统计。rasterstats更直接：
from rasterstats import zonal_stats import geopandas as gpd gdf = gpd.read_file('yunnan_prefectures.shp') stats = zonal_stats(gdf, 'yunnan_landcover_2019_albers.tif', categorical=True, nodata=0) gdf['landcover_stats'] = stats print(gdf[['name', 'landcover_stats']].head())categorical=True表示按类别统计像元数，返回的是字典，键是编码，值是像元数。再乘以100/1e6就是每个州市各地类的面积。注意矢量边界和栅格坐标系必须一致，不一致先用gdf.to_crs(src.crs)转一下。
提示：
zonal_stats默认统计所有像元，如果矢量边界和栅格边界不完全重合，边缘像元会被算入。对精度要求高的场景，可以加all_touched=False，只统计中心点落在多边形内的像元。
5. 避坑与排查：10米土地覆盖处理的5个血泪教训
5.1 现象：面积统计出来全省总面积对不上
原因：最常见的是没做重投影，直接在地理坐标系下用像元数乘100平方米。地理坐标系下像元大小是度，不是米，乘100得到的数字没有物理意义。另一个原因是没掩膜NoData，填充值被当成某一类参与统计。
解决：先gdalinfo确认像素大小单位，如果是度，先转Albers等面积投影；统计前用src.nodata或手动指定填充值做掩膜。做完这两步，全省总面积应该和云南省官方公布面积（约39.4万平方公里）在一个合理范围内，差太多说明还有问题。
5.2 现象：重投影后出现原始数据里没有的类别编码
原因：用了双线性或立方卷积重采样。这两种方法会在地类边界处插值出中间值，比如编码1和编码3之间插出2，凭空多出一个类别。
解决：分类栅格重投影必须用最近邻-r near。如果已经用了其他方法，重新跑一遍。检查方法是重投影前后各跑一次np.unique，对比编码集合是否一致。
5.3 现象：裁剪后边界外像元值为0，但0也是某个真实地类
原因：gdalwarp的-dstnodata 0把边界外设为0，但分类体系里0可能对应「无数据」或某个真实类。如果后续统计不区分，边界外像元会被算入。
解决：裁剪时用一个原始数据里不存在的值作为NoData，比如255或-9999。如果原始数据值域是1-20，用255做NoData最安全。然后在统计时用data != 255掩膜。
5.4 现象：不同来源的10米产品分类编码不一致
原因：国内10米土地覆盖产品有多套，分类体系和编码各不相同。A产品的「林地」编码是2，B产品可能是21。直接混用会导致类别错乱。
解决：拿到任何一套产品，先找它的分类说明。如果没有说明，不要凭经验猜编码含义。可以在小范围内用高分辨率影像或实地知识验证几个典型地类的编码，确认后再大面积使用。混用多套产品做变化检测时，必须先做编码统一。
5.5 现象：云南边境和西北部出现大片异常值
原因：原始产品在边境线附近可能有填充值或缺失，西北部高山峡谷区可能因为云层遮挡或地形阴影导致分类失败。
解决：先做值分布统计，找出占比异常小且空间上集中的编码。用rasterio掩膜掉这些值，或者在报告里注明「XX区域数据缺失，面积统计不含该部分」。不要直接把这些值归到最近的类别，那会引入无法量化的误差。
6. 进阶：用10米数据做变化检测与精度验证
6.1 两期数据做变化检测的最小流程
如果手里有2019年和另一年的同源10米产品，可以做变化检测。核心逻辑是：两期栅格逐像元对比，编码不同的地方就是变化。但前提是两期数据的分类体系、投影、分辨率完全一致。不一致就先统一。
import rasterio import numpy as np with rasterio.open('yunnan_landcover_2019_albers.tif') as src1: data2019 = src1.read(1) profile = src1.profile with rasterio.open('yunnan_landcover_2020_albers.tif') as src2: data2020 = src2.read(1) # 掩膜NoData mask = (data2019 != 255) & (data2020 != 255) change = np.zeros_like(data2019, dtype=np.uint8) change[mask & (data2019 != data2020)] = 1 profile.update(dtype=rasterio.uint8, nodata=255) with rasterio.open('change_2019_2020.tif', 'w', **profile) as dst: dst.write(change, 1)输出是0/1二值栅格，1表示发生变化。进一步可以统计「林地转耕地」「耕地转建设用地」等具体转移矩阵，用pandas的crosstab即可。
参数说明：mask同时掩膜两期NoData，避免「一期有值一期无值」被误判为变化。如果两期NoData值不同，分别掩膜后再取交集。
6.2 用高分辨率影像做抽样精度验证
10米产品在云南复杂地形下，精度不是均匀的。做验证时，用Google Earth或天地图的高分辨率影像，随机抽100-200个点，目视判读真实地类，和产品编码对比，算混淆矩阵和Kappa系数。抽样时要分层：坝区、山区、边境各抽一定比例，否则坝区精度高会拉高整体精度，掩盖山区的问题。
验证结果如果Kappa低于0.7，说明这套数据在你的研究区需要谨慎使用，或者需要做局部修正。修正方法包括：用更高分辨率数据替换错误区域、用DEM和坡度辅助区分耕地和林地、用时间序列NDVI辅助区分草地和灌木。
6.3 一个我常用的习惯
每次拿到新的土地覆盖产品，我第一件事不是打开ArcGIS看，而是跑三行命令：gdalinfo看投影和像素大小，np.unique看值分布，gdalwarp裁到研究区。这三步做完，数据能不能用、要怎么处理，心里就有数了。云南这种地形复杂的地方，10米数据比30米强很多，但强的前提是你把投影、编码、NoData这三件事都处理对了。希望帮到你。
本文还有配套的精品资源，点击获取
