简介本资源为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, dtypenp.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(dtyperasterio.uint8, count1, compresslzw) with rasterio.open(Yunnan_2019_reclass.tif, w, **profile) as dst: dst.write(reclassed, 1)这段代码的核心是查找表映射比逐像元循环快几个数量级。参数上profile.copy()保留了原始的地理变换和投影信息确保输出栅格能和其它图层对齐compresslzw能显著减小文件体积对全省数据尤其有用。注意重分类前一定要确认nodata值否则背景值可能被误映射成某个地类导致面积统计虚高。常见做法是先用np.unique看一下原始编码的实际分布再决定映射关系不要照搬网上的对照表。3. 面积统计与变化检测从栅格到可汇报的数字3.1 按行政区统计各地类面积拿到重分类结果后最常被问到的就是“云南省2019年耕地有多少公顷”。如果手头有县级或乡镇级的行政边界矢量可以用rasterstats做分区统计一次性输出每个行政单元内各地类的像元数和面积。下面这个脚本假设你已经有了一个yunnan_counties.shpimport 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, categoricalTrue, # 按类别统计像元数 nodata0, affinerasterio.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, indexFalse, encodingutf-8-sig)这里的关键参数是categoricalTrue它让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_countsTrue) 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(dtyperasterio.int16, nodatanodata, compresslzw) with rasterio.open(Yunnan_lulc_change_2019_2022.tif, w, **profile) as dst: dst.write(change, 1)变化编码用“前一期类别×10后一期类别”的方式能在一张图里同时表达变化方向和类型比如21表示从农业用地变成生态用地。参数上astype(np.int16)是为了防止乘法溢出原始uint8乘以10后最大到60int16绰绰有余。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年把草地编码为32022年编码为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(fYunnan_{y}_reclass.tif, maskedTrue) da da.squeeze(band, dropTrue) # 去掉波段维 da da.assign_coords(yeary) arrays.append(da) # 沿时间维合并 cube xr.concat(arrays, dimyear) # 计算耕地编码2像元数逐年变化 cropland (cube 2).sum(dim[x, y]) print(cropland.values) # 输出每年的耕地像元数 # 计算每个像元上耕地存在与否的线性趋势 cropland_mask (cube 2).astype(np.float32) trend cropland_mask.polyfit(dimyear, deg1) slope trend.sel(degree1)[polyfit_coefficients] # slope为正表示该像元从非耕地转为耕地为负则相反 slope.rio.to_raster(cropland_trend_2019_2022.tif)这段代码的关键在于xr.concat沿年份维合并形成一个带时间坐标的数据立方体。polyfit对每个像元做一阶线性拟合斜率就是变化速率。参数上maskedTrue让nodata自动变成NaN避免参与运算deg1表示线性拟合如果想捕捉非线性趋势可以改成2但解释起来更复杂。输出的趋势栅格可以直接在QGIS里用发散色带渲染一眼就能看出哪些区域在退耕哪些在复耕。我一般会把这个趋势图和行政区边界叠加快速定位变化热点乡镇再针对性去查原始影像。从那以后我每次拿到新的土地覆盖年份包都会先跑一遍这个趋势脚本确认数据连续性没问题再往下做。希望帮到你。本文还有配套的精品资源点击获取
简介本资源为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, dtypenp.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(dtyperasterio.uint8, count1, compresslzw) with rasterio.open(Yunnan_2019_reclass.tif, w, **profile) as dst: dst.write(reclassed, 1)这段代码的核心是查找表映射比逐像元循环快几个数量级。参数上profile.copy()保留了原始的地理变换和投影信息确保输出栅格能和其它图层对齐compresslzw能显著减小文件体积对全省数据尤其有用。注意重分类前一定要确认nodata值否则背景值可能被误映射成某个地类导致面积统计虚高。常见做法是先用np.unique看一下原始编码的实际分布再决定映射关系不要照搬网上的对照表。3. 面积统计与变化检测从栅格到可汇报的数字3.1 按行政区统计各地类面积拿到重分类结果后最常被问到的就是“云南省2019年耕地有多少公顷”。如果手头有县级或乡镇级的行政边界矢量可以用rasterstats做分区统计一次性输出每个行政单元内各地类的像元数和面积。下面这个脚本假设你已经有了一个yunnan_counties.shpimport 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, categoricalTrue, # 按类别统计像元数 nodata0, affinerasterio.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, indexFalse, encodingutf-8-sig)这里的关键参数是categoricalTrue它让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_countsTrue) 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(dtyperasterio.int16, nodatanodata, compresslzw) with rasterio.open(Yunnan_lulc_change_2019_2022.tif, w, **profile) as dst: dst.write(change, 1)变化编码用“前一期类别×10后一期类别”的方式能在一张图里同时表达变化方向和类型比如21表示从农业用地变成生态用地。参数上astype(np.int16)是为了防止乘法溢出原始uint8乘以10后最大到60int16绰绰有余。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年把草地编码为32022年编码为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(fYunnan_{y}_reclass.tif, maskedTrue) da da.squeeze(band, dropTrue) # 去掉波段维 da da.assign_coords(yeary) arrays.append(da) # 沿时间维合并 cube xr.concat(arrays, dimyear) # 计算耕地编码2像元数逐年变化 cropland (cube 2).sum(dim[x, y]) print(cropland.values) # 输出每年的耕地像元数 # 计算每个像元上耕地存在与否的线性趋势 cropland_mask (cube 2).astype(np.float32) trend cropland_mask.polyfit(dimyear, deg1) slope trend.sel(degree1)[polyfit_coefficients] # slope为正表示该像元从非耕地转为耕地为负则相反 slope.rio.to_raster(cropland_trend_2019_2022.tif)这段代码的关键在于xr.concat沿年份维合并形成一个带时间坐标的数据立方体。polyfit对每个像元做一阶线性拟合斜率就是变化速率。参数上maskedTrue让nodata自动变成NaN避免参与运算deg1表示线性拟合如果想捕捉非线性趋势可以改成2但解释起来更复杂。输出的趋势栅格可以直接在QGIS里用发散色带渲染一眼就能看出哪些区域在退耕哪些在复耕。我一般会把这个趋势图和行政区边界叠加快速定位变化热点乡镇再针对性去查原始影像。从那以后我每次拿到新的土地覆盖年份包都会先跑一遍这个趋势脚本确认数据连续性没问题再往下做。希望帮到你。本文还有配套的精品资源点击获取
