GDAL/rasterio实战：1km土地利用重分类与转移矩阵

简介面向地理学、环境科学、城乡规划等领域的研究者与高年级学生这套中国1km土地利用栅格数据覆盖1980至2015年共7个时点的遥感监测结果可用于长时间序列的土地用途变迁分析、城市扩张识别与生态状况评估。压缩包共126个文件以ArcInfo栅格分块存储为主49个adf、21个nit、21个dat与7个dir、7个001共同组成每年一组的图层另有doc说明文档、log处理日志和xml元数据辅助判读整体约21.83MB按年份取用方便。像素为1km×1km耕地、林地、草地、建设用地、水域等类型均有编码记录逐年份叠加即可提取转移矩阵与变化热点。目前已有2904人学习下载适合作为城市化与生态保护评估的基础底图也可与其他社会经济数据融合搭建土地利用变化模型开展趋势推演。1. 一张 1km 栅格能顶半年的底图活这套 1980–2015 土地利用数据适合谁接手一个跨度三十多年的土地利用变化分析最费时间的往往不是建模而是找一套时间连续、分类标准一致的底图。这套 1980–2015 的 1km 土地利用数据就是用来做这件事的整型编码的栅格切片1 公里像元覆盖多期打开就能算面积、算转移、算建设用地扩张强度。它的价值不在“精度高”而在“分类标准稳”。同一套分类体系跨几十年做趋势分析时省掉大量重投影、重解译和对齐的活1 km 的粒度对省级、流域级、全国级的格局判断够用对地块级当然不够。适合做国土空间、生态评估、碳储量估算、城市扩张的人拿来当背景图层或先验分区也适合课程设计数据量不大普通笔记本能跑完全流程。2. 解压之后先别画图用 gdalinfo 和 rasterio 把数据底细摸清刚拿到压缩包的人九成会先拖进 QGIS 里看一眼配色一拉花花绿绿挺像那么回事然后就以为数据没问题了。真正容易出问题的地方全在元数据里投影是地理坐标还是等积投影、NoData 到底设没设、存的是二级码还是一级码、七期数据的栅格范围是不是严格一致。这几项只要有一项不对后面算出来的面积和转移矩阵全部作废而且错得很隐蔽——图看着是对的数字是错的。所以我的习惯是先做一次“体检”不改数据、不出图只把底细摸清楚。2.1 元数据里到底要看哪几行gdalinfo输出很长但真正决定后续能不能算的只有几行。Size is是行列数全国 1km 大致在 9000×8000 这个量级如果拿到的是 900×800说明分辨率或范围不对。Pixel Size最关键如果单位是度常见 0.00833333 左右说明是地理坐标这个“1km”只是在赤道附近成立往高纬度走东西向像元实际边长会缩短直接拿像元数乘 1 km² 算面积会偏如果单位是米且接近 1000说明已经是投影坐标面积统计会省事很多。NoData Value有没有决定了 0 到底是“无数据”还是某个合法类别。Type是 Byte、Int16 还是 UInt16决定你能不能直接做 256 长度的查找表。还有一个容易忽略的点这套数据有的年份是 GeoTIFF有的年份给的是 ArcInfo GRID 目录一堆hdr.adf、w001001.adf那种。GDAL 对两者一视同仁但如果你写脚本用os.listdir找.tif就会漏掉一半年份。2.2 命令行三分钟体检# 1) 基础信息尺寸、波段、数据类型、NoData、投影 gdalinfo landuse_2015.tif # 2) 顺带把统计量算出来最小值/最大值/均值/标准差 gdalinfo -stats landuse_2015.tif # 3) 七期批量比对只抽关键行省得一行行翻 for f in landuse_*.tif; do echo $f gdalinfo $f | grep -E Size is|Pixel Size|Origin|NoData|Type|PROJCRS|GEOGCRS|PROJCS|GEOGCS | head -n 10 done第一行不用解释任何栅格到手都该先跑一遍。第二行的-stats会真实读一遍数据算统计量全国 1km 的整型栅格也就几秒钟但它给出的最小值非常有价值——如果最小值是 0 而 NoData 显示为 None那基本可以判定 0 是被当成“背景/无数据”用了只是没写进头文件。第三行的批量循环是血泪经验七期数据里混进一期投影不一样的你不批量比对根本发现不了。GDAL 3.x 输出的是PROJCRS/GEOGCRSGDAL 2.x 是PROJCS/GEOGCS两个都 grep 上换机器不至于抓瞎。提示如果gdalinfo里完全没有 NoData 那一行别假设“没有 NoData”要按“0 可能是无效值”来处理后面统计时显式做掩膜。2.3 Python 读进来做一次值分布元数据看完下一步是把像元值本身数一遍。这一步的意义在于确认三件事码值是二级码还是已合并的一级类、有没有 0 或 255 这类无效值混在里面、每一类占的比例是否合理。import numpy as np import rasterio path landuse_2015.tif with rasterio.open(path) as src: arr src.read(1) # 只读第 1 波段多波段数据要确认用哪个 print(shape :, arr.shape) print(dtype :, arr.dtype) print(nodata:, src.nodata) # None 表示头文件里没写 print(crs :, src.crs) print(res :, src.res) # 地理坐标下单位是度投影坐标下是米 print(bounds:, src.bounds) vals, counts np.unique(arr, return_countsTrue) total arr.size order np.argsort(counts)[::-1] # 按出现次数从多到少 for i in order[:15]: print(fcode{vals[i]:4} pixels{counts[i]:10} share{counts[i]/total:.4%})src.read(1)只取第一波段全国 1km 大约 900 多万个像元Byte 类型占 10 MB 左右Int16 也就 20 MB普通机器完全扛得住不需要分块。src.res是最容易被误读的参数在地理坐标系下它的单位是度输出 0.008333 不代表 8.3 米而是 30 弧秒。np.unique(..., return_countsTrue)一次拿到码值和频次比循环快两个数量级。输出怎么读如果码值列表里出现了 11、12、21、23 这种两位数说明是二级码必须重分类后才能做一级类统计如果只有 1 到 6说明数据已经合并过可以直接用。如果最大的几个码值是 0 且占比不小那 0 就是无效值或者背景必须排除。如果出现 127、255 这种在分类体系里找不到对应含义的值八成是填充值没写进 NoData 字段得自己手动掩掉。注意不要把np.unique的结果直接当成面积占比对外发布。它的分母是总像元数包含了无效值和边界外的填充像元必须先做掩膜再算比例。3. 分类码不是随便编的6 个一级类与二级码的映射与重分类3.1 编码体系长什么样为什么不能拿原始码直接统计这类 1km 土地利用产品普遍采用两级分类6 个一级类下面挂若干二级类用两位数编码十位表示一级类、个位表示二级类。常见标准如下以压缩包内附的编码说明为准不同年份可能有个别调整一级类一级码常见二级码二级类含义耕地111、12水田、旱地林地221、22、23、24、25、26有林地、灌木林、疏林地等草地331、32、33高覆盖、中覆盖、低覆盖草地水域441、42、43、44、45、46河渠、湖泊、水库坑塘、滩涂等建设用地551、52、53城镇用地、农村居民点、其他建设用地未利用地661–67沙地、戈壁、盐碱地、沼泽、裸地等直接用二级码做统计最典型的错误是把“耕地”当成一个码去筛。实际上耕地被拆成 11 和 12你只筛 11得到的耕地面积会少一大截而且少的比例在不同区域还不一样——南方水田多、北方旱地多这种偏差会伪装成“区域差异”让你在分析里得出完全错误的结论。另一种情况更隐蔽某些年份的编码做了微调比如某类从 24 挪到 25跨年比较时如果不重分类到一级类就会出现“林地突然减少、另一类突然增加”的假信号。重分类的另一个理由是可比性。二级类的定义在不同年份之间本身就有细微差别做长时间序列时一级类才是相对稳定的比较口径。做趋势、做转移矩阵、做景观格局指数都建议先把二级码压到一级类需要细节时再单独对某一年做二级分析。3.2 用查找表做重分类比循环快几十倍重分类最慢的写法是双层 for 循环加np.where嵌套900 万像元能跑几分钟。正确做法是建一个 256 长度的查找表LUT直接用数组下标取值一次运算完成。import numpy as np import rasterio # 二级码 - 一级类。以压缩包里的编码说明为准这里给的是常见口径 LUT np.zeros(256, dtypenp.uint8) LUT[[11, 12]] 1 # 耕地 LUT[[21, 22, 23, 24, 25, 26]] 2 # 林地 LUT[[31, 32, 33]] 3 # 草地 LUT[[41, 42, 43, 44, 45, 46]] 4 # 水域 LUT[[51, 52, 53]] 5 # 建设用地 LUT[[61, 62, 63, 64, 65, 66, 67]] 6 # 未利用地 # 下标 0 保持 0作为无效值二位数以外的码值也自然落到 0 def reclass(src_path, dst_path): with rasterio.open(src_path) as src: arr src.read(1).astype(np.uint8) # LUT 索引必须是 0-255 的整型 profile src.profile.copy() out LUT[arr] # 核心就这一行向量化完成 out[arr 0] 0 # 原数据 0 视为无效避免 0 被映射 out[out 6] 0 # 兜底任何越界结果都归 0 profile.update(dtypeuint8, nodata0, compressdeflate, tiledTrue, blockxsize512, blockysize512) with rasterio.open(dst_path, w, **profile) as dst: dst.write(out, 1) vals, counts np.unique(out, return_countsTrue) print(dict(zip(vals.tolist(), counts.tolist()))) reclass(landuse_2015.tif, landuse_2015_l1.tif)核心逻辑是LUT[arr]arr里每个像元值是 11就取LUT[11]得到 1是 12 就得到 1是 0 就取LUT[0]还是 0。整个过程是一次内存寻址900 万像元毫秒级完成。astype(np.uint8)不能省如果原始 dtype 是 Int16 而码值里混进了负数直接索引会绕回数组尾部取到错误的值这种 bug 极难排查。out[arr 0] 0是显式保护防止某些数据把 0 定义成“未分类”以外的含义。out[out 6] 0是兜底保证输出只有 0–6 这七个值后面做交叉表时不会冒出意料之外的类别。如果手里的数据码值超过 255比如某些衍生产品的编码上万LUT 就不适用了改用np.searchsorted配合有序的码值数组做映射思路一样只是把“下标查表”换成“二分查找”。3.3 重分类结果怎么落盘不要覆盖原始文件原始数据是唯一的任何时候都不要原地改值。落盘时有两件事要做对一是profile必须从源文件复制保证投影、变换、尺寸全部一致二是显式设定nodata0和压缩参数。compressdeflate对分类栅格压缩比很高全国 1km 一般能从几十 MB 压到几 MB方便传阅。tiledTrue加上 512 的块大小是为了后续按窗口读取时不会整幅解码——做局部区域分析时这个设置能省掉大量 IO。nodata0是给下游工具的约定QGIS、GDAL、rasterio 看到这个字段就会自动把 0 当透明处理。还有个小细节重分类后的 dtype 一定要落到uint8。如果沿用源文件的 Int16文件体积翻倍不说某些库在做类别判断时会因为符号位出岔子。分类数据的动态范围就那么几个值没必要占两个字节。4. 面积和转移矩阵像元计数、等积投影与七期对齐4.1 面积到底该用像元数乘多少这是整套流程里出错率最高的一步。如果数据已经是 1km 的等积投影一个像元就是 1 km²直接计数即可如果数据是地理坐标src.res单位是度像元在地面上的实际面积随纬度变化越往高纬东西向越窄全国范围粗略估算的偏差能到两三成。两种处理方式。第一种是重投影到等积坐标系Albers 是最常用的选择投影参数按研究区的中纬度设置# 重投影到等积坐标系同时把像元统一成 1000m # -r near 是必须的分类数据只能用最近邻双线性会造出不存在的类别 gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 \ x_00 y_00 datumWGS84 unitsm no_defs \ -tr 1000 1000 -r near \ -srcnodata 0 -dstnodata 0 \ landuse_2015_l1.tif landuse_2015_aea.tif-t_srs指定目标投影lat_1/lat_2是两条标准纬线一般取研究区南北边界的纬度lon_0取中央经线。-tr 1000 1000强制输出像元为 1000 米避免 GDAL 自己算出一个 987.3 米之类的奇怪分辨率。-r near是最近邻重采样这一条不能改分类栅格用双线性或三次卷积会把 5建设用地和 1耕地插值成 3草地凭空造出几万平方公里的假草地。-srcnodata/-dstnodata保证无效值在重投影后还是无效值不会被当成 0 类参与统计。重投影之后面积就是像元数 × 1 km²简单直接。第二种方式是不重投影按纬度带逐行计算像元实际面积每行乘以111.32 × cos(lat) × 111.32公里代码稍复杂但不用生成新文件适合只想快速估个量级的场景。我一般选第一种文件多占几百 MB但后面所有统计口径统一省心。数据类型像元面积算法注意事项等积投影1000m像元数 × 1 km²确认src.res确实是 1000地理坐标0.008333°逐行111.32² × cos(lat)高纬误差大不建议直接乘常数其他投影先转等积再统计不要用投影坐标直接乘4.2 转移矩阵行列方向反了结论就反了土地利用转移矩阵本质是一张交叉表行是前期类别、列是后期类别对角线是未变化的面积非对角线是转移量。听起来简单实际动手有三个必踩的点无效值没排除、类别没对齐导致行列缺项、行列方向搞反。import numpy as np import pandas as pd import rasterio names {1: 耕地, 2: 林地, 3: 草地, 4: 水域, 5: 建设用地, 6: 未利用地} codes list(names.keys()) def read_l1(path): with rasterio.open(path) as src: return src.read(1) a read_l1(landuse_1990_l1.tif) # 前期 b read_l1(landuse_2015_l1.tif) # 后期 assert a.shape b.shape, 两期尺寸不一致先做网格对齐 mask (a 0) (b 0) # 只统计两期都有效的像元 mat pd.crosstab( pd.Categorical(a[mask], categoriescodes), pd.Categorical(b[mask], categoriescodes), dropnaFalse, ) mat.index [names[i] for i in mat.index] mat.columns [names[i] for i in mat.columns] mat.index.name 1990 mat.columns.name 2015 print(mat) # 单位km²1 像元 1 km² print(总有效面积:, mat.values.sum())pd.Categorical(..., categoriescodes)是关键。不用它的话如果某一期里恰好没有“未利用地”交叉表就会少一行一列矩阵变成 5×5后面算转移比例时维度对不上或者更糟——行列错位把“耕地转林地”读成“耕地转草地”。dropnaFalse保证缺失组合补 0 而不是被丢弃。mask同时排除两期的无效值因为转移矩阵只能建立在“两期都有明确类别”的像元上单期有效的像元无论怎么处理都会引入偏差。算完的矩阵可以直接mat.to_csv(transition_1990_2015.csv, encodingutf-8-sig)带 BOM 是为了 Excel 打开不乱码。想要转出比例用mat.div(mat.sum(axis1), axis0)得到每行归一化的结果读法是“某一类中有多大比例转去了别的类”。4.3 七期对齐为什么 1980 和 2015 的像元不能直接相减多期栅格的shape一样不代表网格一样。transform里原点偏半个像元两期数据叠在一起就会出现整幅的条带状伪变化——每条带宽度正好一个像元。这种错误在出图时几乎看不出来但转移矩阵会凭空多出大量“耕地转建设用地”。判断方法很简单把七期的transform和crs拿出来逐个比对。import rasterio files [flanduse_{y}_l1.tif for y in (1980, 1990, 1995, 2000, 2005, 2010, 2015)] ref None for f in files: with rasterio.open(f) as src: key (src.crs.to_string(), tuple(round(v, 6) for v in src.transform[:6]), src.shape) if ref is None: ref key print(基准:, f, key) elif key ! ref: print(不一致:, f, key) # 任何一行输出都意味着需要重采样对齐round(v, 6)是因为不同软件写出的transform浮点尾数会有微小差异直接比较会误报。只要有一期不一致就用gdalwarp -te把它对齐到基准网格# 以 2015 期为基准网格把其他期统一过来最近邻重采样 gdalwarp -te 7300000 0 12300000 6000000 -tr 1000 1000 -r near \ -srcnodata 0 -dstnodata 0 \ landuse_1980_l1.tif landuse_1980_aligned.tif-te后面四个数是基准期的范围从gdalinfo的Upper Left/Lower Right或src.bounds取-tr是分辨率。这套组合的意义是不改变基准期把所有其他期强行套到同一张网格上误差只发生在边缘的半行半列中间区域严格逐像元对应。核对完再跑转移矩阵结果才对得上。5. 常见问题排查五个让结果对不上的坑这一章是回头看踩过的坑攒出来的。五个问题的共同点是图看着正常、代码不报错、结果说不通。5.1 读取与投影阶段的三个坑坑一看到Pixel Size 0.008333就当成 1km。现象是脚本里直接count * 1.0算面积最后全国耕地总面积比公开数据大了近三成原因在于0.008333的单位是度30 弧秒的像元在中纬度东西向只有约 0.8 km南北向接近 1 km实际面积约 0.8 km²你按 1 km² 算就系统性高估解决办法是先看src.crs.is_geographic为真就重投影到等积坐标系或者按纬度逐行乘cos(lat)修正。坑二重采样用了默认参数地图上多出一片“草地”。现象是重投影后的分类图在耕地和建设用地的交界处出现细碎的 3 类像元转移矩阵里“耕地转草地”异常高原因是gdalwarp不带-r时对某些格式默认用双线性分类码被插值成中间值5 和 1 之间插出来的 2、3、4 全是假类别解决办法是所有分类栅格的重采样一律写死-r near并且在脚本里加一行注释说明原因免得下次接手的人改回去。坑三0 被当成“未利用地”参与统计。现象是未利用地面积比预期大出一大块且集中在研究区边界外的矩形区域原因是数据本身没有写NoData Value0 是背景填充值但统计脚本把它当成了合法码值解决办法是在gdalinfo里确认 NoData 为空时显式按arr 0做掩膜落盘时写入nodata0并在重分类 LUT 里把下标 0 固定映射为 0。5.2 统计与矩阵阶段的坑坑四转移矩阵对角线占比 99% 以上。现象是七期数据两两做矩阵变化面积小得离谱跟文献里看到的趋势完全相反原因通常有两个一是两期网格差半个像元逐像元比对时大量像元被错配但这种错配往往会让变化量偏大而不是偏小所以更可能是第二个原因——重分类时把某一期单独处理过比如某年只保留了一级类另一年保留二级码交叉表的 categories 没对齐大量像元落到了“未匹配”里被丢掉解决办法是先把七期全部走一遍统一的重分类函数再用固定的codes [1,2,3,4,5,6]作为 Categorical 的类别集合。坑五两期shape一样相减全 0 或者出现规则条带。现象是变化检测结果要么一片空白要么出现宽度恰好一个像元的斜向条纹原因是shape相同但transform的原点差了半像元这是不同批次数据导出时的常见现象解决办法是用transform前六位做 round 后比对不一致就用gdalwarp -te -tr对齐到基准网格绝不要靠“看起来对得上”来判断。注意这五条里最容易被忽略的是坑三和坑五。前三条出错时结果会很离谱你一定会去查后两条出错时结果“看起来合理”可能一路带到论文或报告里。6. 把七期串成一条流水线批量骨架与三条自证单期跑通之后合理的下一步是把它函数化七期一次跑完避免手工操作带来的不一致。6.1 流水线骨架import os import numpy as np import pandas as pd import rasterio YEARS [1980, 1990, 1995, 2000, 2005, 2010, 2015] RAW raw/landuse_{}.tif # 原始文件 L1 out/landuse_{}_l1.tif # 重分类后 LUT np.zeros(256, dtypenp.uint8) LUT[[11, 12]] 1 LUT[[21, 22, 23, 24, 25, 26]] 2 LUT[[31, 32, 33]] 3 LUT[[41, 42, 43, 44, 45, 46]] 4 LUT[[51, 52, 53]] 5 LUT[[61, 62, 63, 64, 65, 66, 67]] 6 def to_l1(year): 二级码压到一级类统一 dtype 与 nodata if os.path.exists(L1.format(year)): return L1.format(year) with rasterio.open(RAW.format(year)) as src: arr src.read(1).astype(np.uint8) prof src.profile.copy() out LUT[arr] out[arr 0] 0 out[out 6] 0 prof.update(dtypeuint8, nodata0, compressdeflate) with rasterio.open(L1.format(year), w, **prof) as dst: dst.write(out, 1) return L1.format(year) # 1) 全部重分类 paths {y: to_l1(y) for y in YEARS} # 2) 一致性门禁网格不一致直接抛错不要往下跑 ref None for y, p in paths.items(): with rasterio.open(p) as src: key (src.shape, tuple(round(v, 6) for v in src.transform[:6])) if ref is None: ref key elif key ! ref: raise ValueError(f{y} 期网格与基准不一致先做 gdalwarp -te 对齐) # 3) 逐对算转移矩阵只出相邻期避免组合爆炸 names {1: 耕地, 2: 林地, 3: 草地, 4: 水域, 5: 建设用地, 6: 未利用地} rows [] for a_y, b_y in zip(YEARS[:-1], YEARS[1:]): a rasterio.open(paths[a_y]).read(1) b rasterio.open(paths[b_y]).read(1) m (a 0) (b 0) tab pd.crosstab(a[m], b[m]) for i in tab.index: for j in tab.columns: if i ! j and tab.loc[i, j] 0: rows.append({from_year: a_y, to_year: b_y, from: names[i], to: names[j], area_km2: int(tab.loc[i, j])}) pd.DataFrame(rows).to_csv(out/changes.csv, indexFalse, encodingutf-8-sig)骨架里最重要的是第二段的“一致性门禁”。很多人把校验写在最后跑完才发现数据没对齐白等半小时。放在中间错就立刻停。第三段只算相邻期是因为七期两两组合有 21 对实际分析里绝大多数只需要看相邻期和首末期把组合压下来输出表也清爽。6.2 三条自证**第一条总面积守恒。**每一期的 1–6 类面积加起来应该等于该期有效像元数乘以像元面积。如果某一期加起来比别的期少 5% 以上说明有类别在重分类时被丢掉了。**第二条首末对比方向一致。**1980 到 2015 的建设用地增量应该约等于所有相邻期增量之和允许边缘像元的微小差异。如果两者差得多多半是中间某两期网格不一致。**第三条抽一个小窗口人工看图。**随机裁一块 100×100 像元的区域把两期分类图并排出图肉眼扫一遍。建设用地扩张、水体变化这类信号是能看出来的如果图上明显在变而矩阵里几乎为零问题一定在数据处理而不是在现实世界。我第一次批量跑这套数据时没做门禁七期里有一期是别人后补的投影偏了几百公里转移矩阵算出来的“耕地转水域”面积大到能上新闻。从那以后我每次处理多期栅格第一件事就是把transform和crs打印出来逐个比对宁可多花两分钟也不相信“看着对得上”。希望帮到你。本文还有配套的精品资源点击获取