简介面向GIS与遥感应用人员提供2019年10米分辨率云南省土地覆盖/土地利用栅格数据适合土地利用分析、城市扩展监测和生态评估等场景。数据源自全球陆地覆盖制图成果基于哨兵影像与深度学习方法制作涵盖耕地、林地、草地、灌木、湿地、水体、不透水面、裸地、冰雪等类别已统一投影为WGS84坐标系并按照云南省各州市行政边界裁剪可直接在ArcGIS、QGIS中加载省去自行下载、投影转换与裁剪的繁琐流程。压缩包共112个文件含16个州市的TIF栅格数据配套tfw坐标配准文件、xml元数据、dbf属性表、png预览图及xlsx面积统计表格整体约115.71MBTIF与tfw构成主数据xml/dbf记录来源与属性png/xlsx便于快速浏览与汇总。目前已有326人学习下载文件命名包含州市信息便于检索使用。1. 2019年10m精度云南省土地覆盖土地利用这份数据到底能不能直接用拿到这种命名规范、后缀还是.rar的数据包你的第一反应可能和我一样先解压看一眼再说。真正打开之后才发现它不是一张普通影像而是一份覆盖云南省全域的10米分辨率栅格分类图每个像元代表地面上10米乘10米的方块类别包括耕地、林地、草地、水体、不透水面等。相比大家更熟悉的30米Landsat分辨率产品它能把乡村道路、小型水塘、零散宅基地都画出来做县域尺度的土地利用现状调查和变化检测非常合适。但“10m精度”这个说法很容易被误读。它不是指每个类别的分类准确率都达到10米级的可靠程度而是空间采样间隔是10米分类本身依然是模型估算结果山区阴影、云覆盖、混合像元都会让局部类别出错。适合它的用户是GIS分析师、遥感应用工程师和规划领域从业者拿它做底图、做统计、做趋势判断都没问题但别指望它能替代高精度外业调查。下面的内容围绕“怎么打开、怎么用、怎么验、坑在哪”展开我把这几年处理这类产品的经验拆开讲。2. 打开压缩包之前先弄清这份数据是什么格式、什么坐标系、什么图例土地覆盖数据看着只是一个tif实际操作中百分之八十的问题都出在“你还没搞懂它是什么”就急着叠加分析。这个章节把解压后必做的三件事讲清楚看文件清单、读栅格元数据、理解坐标系。2.1 压缩包不是障碍rar解出来先看文件清单和栅格元数据.rar后缀在Linux服务器上偶尔会让不熟悉的人卡一下因为系统默认不带解压工具。常见做法是用unrar或7z解压Windows下用Bandizip或WinRAR都行。解压后先别急着拖进ArcGIS先看目录结构这是每份遥感数据都应该养成的习惯。# 解压命令x 表示保留完整路径 unrar x 2019年10m精度云南省土地覆盖土地利用.rar -d /data/yunnan_2019/ # 查看解压后的文件列表确认有没有说明文档 ls -lh /data/yunnan_2019/我见过太多次“数据打不开”的案例最后发现只是没解压或者解压后还有一个嵌套子目录。解压完成后用gdalinfo直接读取栅格头部信息这是判断数据底细最快的方式不需要打开任何图形界面。gdalinfo /data/yunnan_2019/LC_Yunnan_2019_10m.tif重点看四类输出Size后面的行列数可以估算数据覆盖范围Pixel Size如果是0.000089度左右说明是WGS84地理坐标下的10米像素Coordinate Reference System会明确写是GCS_WGS_1984还是某个投影坐标系NoData Value通常是255或0后续统计面积时要先把它排除。如果数据包里有配套的图例PDF或TXT也要一并打开因为你无法凭文件名猜出类别码的含义。很多公开产品喜欢用10、20、30这类整数代表类别没有图例对照就会在重分类时翻车。2.2 10m精度的真实含义它与30m、250m产品在分析算法上的差异先明确一个概念空间分辨率10米指的是单个像元在地面上对应的边长。10米比30米精细在哪一个30米像元覆盖900平方米一条3米宽的乡村道路在30米影像上完全消失10米像元只覆盖100平方米道路、田埂、小型水体都能以独立图斑的形式出现。这就是10m数据最直接的吸引力地物的空间形态更接近真实。从产品来源看这类10m土地覆盖数据一般是基于Sentinel-2 10米多光谱影像结合Landsat长时间序列和地形辅助数据用随机森林、梯度提升树这类分类器生成的。部分公开产品会把Landsat与Sentinel-2的融合特征一起输入模型以提升时间连续性。相比250米MODIS产品它能看到村级地物相比30米产品它在小图斑识别上明显占优。但代价也很现实分类噪声更大、图斑更碎、存储体积更大。一个覆盖云南省的10米栅格未压缩tif动辄几个GB压缩后依然不小。这里必须强调“分辨率”和“精度”不是一回事。10m空间分辨率描述的是像元大小而分类精度描述的是“图上类别和地面真实类别一致的程度”。我见过不少项目把“10m精度”理解成“精确到10米的土地利用调查成果”直接拿去做权属争议判据结果可想而知。云南多云多雨光学影像有效覆盖少山区地形阴影重分类产品在山区的可靠度天然低于平原地区。所以第4章专门讲验证流程那是决定这份数据能不能用的关键一步。2.3 坐标基准WGS84与投影坐标的区别以及为什么切片会错位拿到元数据后坐标系判断是第一道分水岭。如果gdalinfo显示GCS_WGS_1984说明栅格是经纬度网格每个像元的宽高单位是度不是米。这种数据做显示没问题但做面积统计、长度量算、与投影坐标矢量叠加时误差会悄悄出现。云南地处低纬度经度跨度约从东经97.5度到106.2度横跨UTM 47N和48N两个分带。如果你的分析范围在云南西部用EPSG:32647中部和东部用EPSG:32648。如果做全省尺度的面积统计我更推荐Albers等积投影因为等积投影能保证每个像元代表的面积在不同纬度保持一致这对土地覆盖面积占比分析很重要。还有一个经常被忽视的问题在线影像服务默认是Web Mercator坐标系而tif是WGS84或UTM。ArcGIS Pro打开多个图层时会做动态投影肉眼看着对齐了但实际叠加位置可能存在一个到几个像元的偏移。做目视解译时这种偏移还能接受一旦做栅格与矢量的空间连接或像元统计就会莫名出现错位。解决办法只有一个先把所有数据重投影到同一个坐标系再开始分析别依赖软件的实时投影。3. 把数据纳入分析流程重投影与图例重编码、裁剪预处理看起来简单实际上每一步都可能改变分析结果。重投影选错重采样方法类别码会变得面目全非图例映射搞错统计口径全乱裁剪边界不合适面积核查会出问题。这一章给出可直接复用的命令和代码。3.1 用GDAL重投影到适合云南的投影坐标系重采样方法千万别选错重投影对连续变量如NDVI、气温用双线性或三次卷积没问题但土地覆盖是类别变量像元值只是类别编号不是有物理意义的连续量。如果用bilinear重采样边界处会产生“1.5类”这种不存在的结果后续统计直接报废。正确做法是最近邻。# 从WGS84重投影到UTM 48N使用near重采样输出压缩tif gdalwarp -t_srs EPSG:32648 -r near -ot Byte -co COMPRESSLZW \ /data/yunnan_2019/LC_Yunnan_2019_10m.tif \ /data/yunnan_2019/LC_Yunnan_2019_UTM48.tif参数说明-t_srs EPSG:32648指定目标坐标系-r near是全篇最重要的参数它保证重采样后像元值仍然是原始类别整数-ot Byte保持8位整型输出避免无意间转成float64导致文件体积膨胀-co COMPRESSLZW是无损压缩选项能显著减小tif体积。命令跑完后再用gdalinfo检查一次Pixel Size是否接近10米、Coordinate Reference System是否已变成UTM不要跳过这一步。如果你的项目需要把云南全省拼到一张全国底图上可以统一到CGCS2000的Albers投影。但要注意一旦经过两次重投影类别边界会被重采样两次噪声会增加一轮。所以尽量一次到位确定好目标坐标系后后续所有数据都统一到它。3.2 重分类把数据集自带的分类体系映射到你自己的分类体系不同土地覆盖产品的分类体系差异很大有的分9类有的分24类。你的分析需求可能只关注耕地、林地、草地、水体、不透水面五类就需要把原始类别码重新映射。在做这一步之前先建立一个类别码-类别名对照表这是重分类的基准也是后续写报告时别人能看懂你统计口径的依据。原始类别码示例原始类别含义映射后类别码映射后含义10水田1耕地20旱地1耕地30阔叶林2林地40针叶林2林地50灌木3草地/灌丛60草地3草地/灌丛70水体4水体80不透水面5建设用地90裸地0剔除注意上面是示例映射实际类别码要以你手上数据包的图例为准不要照抄。很多数据集的草地和灌木是分开的如果项目需要区分灌丛和草地就不要合并。我用Python的rasterio做重分类因为脚本可重复执行换个数据集改一下字典就能跑。import rasterio import numpy as np # 示例原始类别码到目标类别码的映射 reclass_map { 10: 1, # 水田 - 耕地 20: 1, # 旱地 - 耕地 30: 2, # 阔叶林 - 林地 40: 2, # 针叶林 - 林地 50: 3, # 灌木 - 草地/灌丛 60: 3, # 草地 - 草地/灌丛 70: 4, # 水体 - 水体 80: 5, # 不透水面 - 建设用地 90: 0, # 裸地 - 剔除 } with rasterio.open(/data/yunnan_2019/LC_Yunnan_2019_UTM48.tif) as src: profile src.profile data src.read(1) out np.zeros_like(data, dtypenp.uint8) for old_val, new_val in reclass_map.items(): out[data old_val] new_val with rasterio.open(/data/yunnan_2019/LC_Yunnan_2019_reclassed.tif, w, **profile) as dst: dst.write(out, 1)这段代码先把原始栅格读成numpy数组再遍历映射字典把对应像元值替换掉最后写回新tif。逻辑很简单但有两个细节一是dtypenp.uint8要显式指定避免输出变成int32二是剔除类设为0后后续统计要记得把0排除否则面积统计会把裸地和非研究区混在一起。3.3 云南省边界裁剪用它与面积核查裁剪是预处理最后一步也是最容易出边界争议的地方。我一般用省级行政区划矢量做裁剪裁剪前先确认矢量边界的坐标系与栅格一致如果不一致先统一再裁。# 用云南省级边界裁剪重分类后的栅格 gdalwarp -r near -cutline /data/boundary/yunnan_province.shp -crop_to_cutline \ /data/yunnan_2019/LC_Yunnan_2019_reclassed.tif \ /data/yunnan_2019/LC_Yunnan_2019_clip.tif参数说明-cutline指定矢量边界文件-crop_to_cutline表示把栅格裁剪到矢量范围之外不留黑边。这里有个常见现象裁完的栅格沿着省界会看到锯齿状边缘这不是数据损坏因为省界矢量本身的精度远高于10米像元而栅格像元是方形网格边界必然呈锯齿。千万不要为此去做平滑滤波那会改变边界处的类别。裁剪后立刻做一次面积核查。云南全省面积约39.4万平方公里10米分辨率下理论上有效像元数应该在39亿左右用GIS统计该tif的有效像元数乘以100再除以1,000,000得到的结果如果在35万到40万平方公里区间说明裁剪正常。如果差距过大通常是NoData值没设对或裁剪边界叠加错位。4. 验证数据的可信度在云南复杂地形下如何抽样比对土地覆盖数据集最大的黑匣子在于你没参与它的训练样本构建就不知道它在你关注的区域到底靠谱不靠谱。我见过太多项目直接拿产品数据做面积统计然后出了奇怪的结论。验证不是学术研究的附属品是工程落地前必须走的一道工序。4.1 验证样本怎么布分区域分层随机抽样别全省撒点云南的地形是“十里不同天”横断山区、高原湖泊、喀斯特地貌、干热河谷并存。如果在全国或全省范围内随机抽点结果一定是森林和裸岩占大多数城市、农田、水体样本严重不足算出来的混淆矩阵看似漂亮但对你想用的类别没有参考价值。正确做法是分层抽样先按州市分区再按类别分层保证每个类别至少有50个验证点山区和坝区分开布点。import geopandas as gpd import numpy as np from shapely.geometry import Point # 读取重分类后的栅格类别分布简化演示用类别栅格的唯一值作分层依据 # 这里假设已经从栅格中提取出类别点列表 class_list [1, 2, 3, 4, 5] # 耕地、林地、草地、水体、建设用地 points [] # 每个类别按分层数量生成随机候选点 for cls in class_list: # 实际项目中应基于类别栅格生成掩膜再在掩膜范围内随机采样 n 60 # 每个类别至少60个点 x np.random.uniform(97.5, 106.2, n) y np.random.uniform(21.0, 29.5, n) for xi, yi in zip(x, y): points.append({class: cls, geometry: Point(xi, yi)}) gdf gpd.GeoDataFrame(points, crsEPSG:4326) gdf.to_file(/data/validation/random_samples.shp)这里的随机采样是演示骨架实际工作中我的习惯是用重分类后的栅格生成每个类别的掩膜在掩膜范围内做空间均匀随机采样避免所有点挤在一起再叠加坡度数据在坡度大于25度的区域额外增加20%的样本。布好点之后每个点用当年的高分影像或Sentinel-2真彩色影像做人工判读得到参考类别。自动判读软件可以作为辅助但最终参考类别必须经过人工核对否则验证就变成了两个模型之间的互相比较。4.2 混淆矩阵、总体精度与Kappa的计算验证不是看一个数字人工判读完成后把每个点的“分类结果类别”和“参考类别”整理成两列CSV直接计算混淆矩阵和Kappa系数。这一步用Python最顺手。import pandas as pd from sklearn.metrics import confusion_matrix, cohen_kappa_score df pd.read_csv(/data/validation/validation_samples.csv) # 列说明 # class_map 是10m数据集给出的类别 # class_ref 是人工影像判读得到的参考类别 cm confusion_matrix(df[class_ref], df[class_map], labels[1, 2, 3, 4, 5]) oa cm.diagonal().sum() / cm.sum() kappa cohen_kappa_score(df[class_ref], df[class_map]) print(混淆矩阵\n, cm) print(总体精度{:.2%}.format(oa)) print(Kappa系数{:.3f}.format(kappa))输出结果中混淆矩阵对角线是分类正确的样本数非对角线是混分情况。Kappa系数通常大于0.75认为一致性较好但这个阈值不能一刀切如果你把原始20类合并成5类精度天然会上升Kappa也会好看这不代表原始数据好只代表合并后的类别更容易分得开。除了总体精度更要看每个类别的用户精度和生产者精度。用户精度衡量的是“图上画成草地的像元真实是草地的比例有多大”生产者精度衡量的是“真实草地有多少被识别出来了”。用户精度低的类别在实际应用中会被高估生产者精度低的类别会被漏估。把混淆矩阵按行归一化就能得到用户精度按列归一化则是生产者精度两条线都要看。4.3 三类最容易翻车的类别草地、裸地、灌木的混分现象在云南验证这类10m数据我的血泪经验是别把时间花在林地和水体上它们通常分得比较准真正容易混的是草地、裸地和灌木。原因很直接稀疏灌丛在10米像元内可能是“草土石头”的混合光谱分类器很容易把它划到草地裸露的石灰岩表面光谱与建设用地中的裸土、屋顶材料相似山坡阴影区里的草地又会被误判为林地。所以验证报告里要单独统计这三个类别的混淆情况。如果发现草地的用户精度低于70%说明该产品在你的研究区不适合直接做草地面积统计要么换数据源要么在分析时把草地和灌木合并成一个“草地灌丛”类别再使用。这个决定应该在拿到混淆矩阵之后做而不是在项目设计阶段拍脑袋。数据产品的边界在哪验证之后才真正清晰。5. 避坑清单拿到这份数据后的5个高频踩坑记录这类数据包的坑大部分不在数据本身而在使用流程。下面5个问题是读者问得最多、也是我在实际项目里翻过车的地方按“现象-原因-解决”写清楚。5.1 坑1用ArcGIS直接打开.rar一直报错打不开栅格现象把.rar文件拖进ArcGIS Pro或ArcMap系统提示“Failed to open raster dataset”或根本识别不了文件。原因ArcGIS的栅格读取接口不支持rar容器格式它只认压缩包内的tif/img格式不会自动解压。很多人以为是压缩包损坏白白浪费时间重新下载。解决先解压到本地目录再打开。注意解压后可能还有一层嵌套文件夹需要解压到看见tif文件为止。如果是Linux环境用unrar x解压时保留完整路径避免中文目录名导致编码问题。5.2 坑2直接用WGS84经纬度栅格统计面积结果和官方数据对不上现象用ArcGIS的栅格属性表直接统计类别像元数再乘以像元面积得到的总面积比云南省官方面积明显偏大或偏小。原因WGS84地理坐标下的像元尺寸用度表示1度×1度的地面面积在不同纬度不一样。虽然云南纬度低经度方向的形变没有高纬度严重但统计精度要求高时这种误差依然不可接受。解决必须先把栅格重投影到等积投影或UTM投影再做面积统计。用gdalwarp -t_srs EPSG:32648或Albers投影均可重投影后像元尺寸变成米面积计算结果才可信。5.3 坑3叠加在线影像时“对不齐”地物错位半个像元到几十米现象把10m土地覆盖tif与在线卫星影像叠加道路和房屋的轮廓总是错开一段距离放大后很明显。原因在线影像默认是Web Mercator投影而tif是WGS84或UTM投影图层之间的动态投影只是显示层面的临时对齐不是真正的坐标统一另外不同影像的获取时间不同地物本身也可能发生了变化。解决不要依赖动态投影做目视判读。先把所有数据重投影到同一个坐标系再叠加。如果对齐后仍有系统性偏移可以检查数据包说明中是否提到原始影像几何校正精度部分产品存在几米到十几米的几何偏差这种偏差无法通过投影设置消除只能作为已知误差在报告中说明。5.4 坑4打开tif后整幅图黑白一片完全看不出类别现象单波段分类栅格导入GIS后显示为灰色渐变或全黑没有分类颜色。原因分类栅格的值是类别编号不是连续光谱值。GIS默认用拉伸渲染显示单波段数据拉伸后相近类别码显示灰度差异极弱看起来就是一片黑。解决在图层符号化设置里改成唯一值渲染按类别码赋予颜色。如果是QGIS双击图层打开样式面板渲染类型选择“单波段伪彩色”颜色条选择离散色带再把类别码逐一设置颜色。也可以把tif内嵌的颜色表导出使用部分产品本身携带Color TableGIS没有自动读取需要手动切换渲染方式。5.5 坑5拿不同年份的土地覆盖数据做变化检测结果全是碎斑现象用2015年和2019年两期数据做像素级差值输出结果全是密密麻麻的小图斑完全看不出有意义的区域变化。原因两期数据的分类算法、影像季节、空间对齐都可能不同同一条河流或同一片森林在不同年份可能被分到不同类别。像素级比较会把这种“分类不稳定”全部当成变化结果就是碎斑遍地。解决先做变化检测前的平滑处理对分类栅格执行多数滤波或众数滤波消除孤立像元再做类别转移矩阵分析而不是逐像元比较。统计层面用“转出面积”“转入面积”来描述变化比像素级视觉对比可靠得多。6. 用这份10m数据做分析一个云南某县域土地利用变化的工作流数据经过预处理和验证后可以支撑一个完整的县域土地利用变化分析。这里以2015年和2019年两期10m数据做类别转移矩阵为例给出一个可以落地的技术路径。先准备两个预处理完成、坐标系一致、类别体系一致的tif。用Python直接对两个栅格逐像元比较统计不同类别之间的转换面积。import rasterio import numpy as np with rasterio.open(/data/yunnan_2019/LC_2015_clip.tif) as a, \ rasterio.open(/data/yunnan_2019/LC_2019_clip.tif) as b: old a.read(1) new b.read(1) profile a.profile mask (old 0) (new 0) # 排除无效像元 n 6 # 类别数按实际类别数量调整 matrix np.zeros((n, n), dtypenp.float64) for i, j in zip(old[mask], new[mask]): matrix[i-1, j-1] 1 # 每个像元100平方米转为平方公里 matrix matrix * 100 / 1_000_000 print(类别转移矩阵平方公里) print(matrix)这段代码输出的矩阵里对角线是2015年到2019年没有变化的面积非对角线是类别转换面积。例如矩阵第1行第3列就表示2015年耕地到2019年草地转换了多少平方公里。我的习惯是做完转移矩阵后再套一次最小图斑过滤把面积小于0.1平方公里的碎斑剔除只保留有实际意义的连续变化区域。过滤后叠加到当年的高分影像上目视抽检确认是真实变化还是分类噪声。这一步不能省否则报告里出现的“草原转耕地”“林地转水体”很可能是虚假变化。我现在的习惯是拿到任何公开土地覆盖数据集先花半天做验证和预处理再决定它能否进入正式分析流程。每一份数据的坑都不一样但验证流程是通用的。希望帮到你。本文还有配套的精品资源点击获取