简介中国1km土地利用数据1980—2015压缩包面向地理学、环境科学、城乡规划与资源管理等领域的研究者和学生提供覆盖1980至2015年多个时点的1km格网土地利用类型信息可用于分析近四十年土地变迁、城市扩张、生态保护与农业动态。包内共126个文件以adf、nit、dat、dir等ArcGIS栅格单元文件为主并含log、doc、xml等说明与辅助文件压缩包约21.83MB文件按年份与数据集组织。已有2904人学习下载。数据以栅格形式记录耕地、林地、草地、建设用地、水域等类型便于对比不同年份识别显著变化也可与社会经济数据结合构建土地利用变化模型。对开展区域案例研究、政策评估和可持续分析具有较高参考价值。1. 1km 分辨率能干什么先把这套土地利用栅格的适用边界说清楚一个做县域生态评估的 A 同学拿到 1980–2015 的中国 1km 土地利用数据之后直接按像元数乘以 1 平方公里去算耕地转出的面积。第一次跑出县城建成区扩张 8%换一个投影再跑变成 3%。他以为是数据坏了其实这份土地利用数据里的问题几乎都集中在“单位”和“编码”这两件事上。这套资源是 1980 到 2015 年逐期的全国 1 公里栅格按年份分包解压之后就是可以直接读进 GDAL 的类别栅格省掉了自己去拼接、配准、对齐的整段流程。它适合做长时序、大范围的地类结构分析比如耕地流向、城镇化扩张、林地草地的年际波动也适合给生态模型当驱动数据。但它不适合地块级权属判断也分不清城市内部的用地功能。1 公里就是一个像元一平方公里边界模糊、细小地物会被吞掉这是分辨率本身决定的跟数据质量无关。2. 解压之后先别画图把目录、编码表和元数据核一遍拿到一个压缩包最省时间的做法不是先把它塞进 QGIS 看一眼出图效果而是花二十分钟把三件事确定下来文件是怎么组织的、编码是什么意思、栅格自己的元数据靠不靠谱。这三件事没定下来就动手算面积后面所有的数字都是悬空的。我把这一步叫体检体检做完后面写脚本才有依据。2.1 摊开压缩包年份目录与文件命名的识别规律先看清单再解压不要一口气全铺开。很多全国尺度的栅格包体量大解压完才发现里面是 ArcInfo Grid 那种一堆.adf的目录结构单独拷一个文件出去就坏了得整个目录一起搬。用清单先把结构摸清能省掉一次返工。# 先看清单确认里面是单文件栅格tif/img还是 ArcInfo Grid 目录 unzip -l 中国1km土地利用数据1980-2015.zip | head -40 # 解压到独立目录别和别的项目数据混在一起 unzip -q -d ./lu_raw 中国1km土地利用数据1980-2015.zip # 列出所有栅格候选按路径排序后年份顺序一眼可见 find ./lu_raw -type f \( -iname *.tif -o -iname *.img \) | sortWindows 上没有unzip就用 7-Zip 的图形界面做同样的事重点是先看清单、再解压到独立目录。清单里如果出现.prj、.tfw、.aux.xml这些小文件说明栅格带坐标和统计附属信息搬家的时候要跟主文件一起走漏了.prj有些软件会默认按像素坐标处理位置直接跑偏。解压完把年份抽出来建一张“年份 → 文件路径”的映射表。后面所有脚本都靠这张表循环不要在每个脚本里重新拼文件名。import re from pathlib import Path files sorted(Path(./lu_raw).rglob(*.tif)) year_map {} for f in files: m re.search(r(19|20)\d{2}, f.stem) # 从文件名里抓四位年份 if m: year_map[int(m.group(0))] f print(sorted(year_map)) # 期望看到 1980、1990、1995、2000、2005、2010、2015 之类的序列这段代码的逻辑很直白遍历解压目录下所有.tif用正则从文件名主干里抓年份建成字典。如果打印出来的年份少于预期说明命名规则不是纯年份可能带了版本后缀或者分幅编号这时候先把files的前十个路径打出来看规律别硬套正则。年份数量对不上是最容易被忽略的第一个坑少一期数据做时序分析结论能差出好几年。2.2 分类编码表一级类看趋势二级类看结构这一类全国土地利用栅格主流用的是两级分类一级类 6 个二级类二十多个用一个整数编码表示。一级类够做趋势和总量二级类才用来看内部结构比如耕地里水田和旱地的比例、林地里是有林地还是灌木林。一级类编码常见二级类与编码耕地111 水田、12 旱地林地221 有林地、22 灌木林、23 疏林地、24 其他林地草地331 高覆盖度草地、32 中覆盖度草地、33 低覆盖度草地水域441 河渠、42 湖泊、43 水库坑塘、44 永久性冰川雪地、45 滩涂、46 滩地城乡工矿居民用地551 城镇用地、52 农村居民点、53 其他建设用地未利用土地661 沙地、62 戈壁、63 盐碱地、64 沼泽地、65 裸土地、66 裸岩石质地、67 其他这张表要当成“常见口径”来用不能当成“这份数据一定如此”来用。不同分发批次里二级类可能做过合并也可能某些年份只给了一级类编码。我一般会先跑一遍值域统计看实际出现了哪些整数再和这张表对照出现的编码是不是都在表里、有没有超出 67 的值。有超出就说明这一期用了别的口径得单独处理。提示先跑值域统计再决定用哪套编码解释数据比对着文档猜快得多也准得多。把二级编码折算回一级最省事的写法是整除def to_level1(code: int) - int: 两位二级编码折算成一级编码11-143-467-6 return code // 10 if code 10 else code注意这个函数只对两位编码成立。如果某期数据里混进了三位编码整除会算出 10 以上的值这时候要改成查表。这也是为什么我坚持先看值域函数写得再漂亮口径不对一样算错。2.3 用 gdalinfo 和 rasterio 给数据做一次体检体检的核心是三件事坐标系是什么、像元多大、NoData 是怎么声明的。这三项决定了后面面积能不能直接乘。gdalinfo -stats LU_2015.tif-stats会顺带算出最小最大值和直方图输出里的Coordinate System is、Pixel Size、NoData Value三个字段要重点看。如果Pixel Size显示的是类似(0.00833, -0.00833)这种小数说明这份数据是经纬度网格不是米制投影——这一点非常关键后面第 3 章会展开。再用 rasterio 批量体检一遍把关键属性整理成能直接比较的结构import rasterio import numpy as np with rasterio.open(LU_2015.tif) as src: print(CRS :, src.crs) print(宽高 :, src.width, src.height) print(像元尺寸 :, src.res) # (x, y)单位跟随 CRS print(NoData :, src.nodata) print(数据类型 :, src.dtypes[0]) arr src.read(1) vals, counts np.unique(arr, return_countsTrue) for v, c in zip(vals, counts): print(f编码 {int(v):3} 像元 {int(c):12,})src.res返回的是像元在x、y方向上的尺寸单位跟 CRS 走米制投影就是米经纬度就是度。src.nodata如果打印出None说明栅格没有声明无效值这是最常见的隐患之一海洋、境外区域可能是 0 或者 255 这种占位数字会被当成真实地类统计进去。np.unique那一段是看值域分布顺便能发现数据里到底有没有二级编码如果只出现 1 到 6 这六个数那这份就是一级类版本别按二级类去拆。体检做完应该能回答三个问题这是等积投影还是经纬度网格、像元实际物理面积是多少、无效值怎么处理。三个问题都有答案了再往下算面积才不会翻车。3. 面积算不对的根源投影、像元面积与 categorical 重采样面积统计是这套数据被用得最多的场景也是错得最多的场景。错法通常有三种拿经纬度网格直接当等面积网格用、按像元计数乘一个固定系数、重采样的时候用了插值算法。这三种错误各有各的隐蔽性跑出来的数字看着都很像那么回事。3.1 Albers 等积投影与经纬度网格同一块地的两个答案全国尺度的土地利用栅格常见做法是采用 Albers 等积圆锥投影一套典型参数是双标准纬线 25°N 和 47°N、中央经线 105°E。选等积投影的理由很实际土地利用统计只关心面积不关心形状和距离是否保真等积投影能保证任意区域的面积比例不变做汇总时才不会被投影变形带偏。如果拿到的是发布方给的 Albers 版本像元是 1000 米 × 1000 米一个像元就是一平方公里面积直接乘就行。麻烦的是经纬度版本。经纬度网格里一个像元覆盖的东西向距离随纬度收缩同样“1 度 × 1 度”的方格在 30°N 和 45°N 覆盖的实际面积能差两成左右。粗算像元物理面积的公式是这样import math def cell_area_km2(px_x_deg: float, px_y_deg: float, lat_center: float) - float: 经纬度网格单个像元的近似面积km² dx abs(px_x_deg) * 111.320 * math.cos(math.radians(lat_center)) # 经向实际距离 dy abs(px_y_deg) * 111.132 # 纬向实际距离 return dx * dypx_x_deg、px_y_deg是像元的经纬度尺寸lat_center是像元中心纬度。这个近似在单个像元上误差不大但做全国汇总时不同纬度带要用不同系数逐像元算一遍才准写起来不难但容易漏。我一般会直接绕开这个问题先在等积投影下统计实在要用经纬度版本就先重投影过去再算面积。注意判断依据是src.crs.is_geographic。返回True就是经纬度网格这时候把像元数乘以一个固定的平方公里数一定会偏。3.2 用 rasterio 批量统计分类面积确定了投影和像元面积统计就是纯计数。这段代码我用了很多次改的主要是输入输出路径import numpy as np import rasterio def class_area_km2(tif_path: str, pixel_area: float None) - dict: 返回 {编码: 面积km²}。pixel_area 为 None 时按米制像元自动推算。 with rasterio.open(tif_path) as src: arr src.read(1) nodata src.nodata # 米制投影下像元面积 res_x * res_y除以 1e6 换成 km² if pixel_area is None: pixel_area src.res[0] * src.res[1] / 1e6 if nodata is not None: arr arr[arr ! nodata] # 剔除无效值别让占位数字混进统计 vals, counts np.unique(arr, return_countsTrue) return {int(v): float(c * pixel_area) for v, c in zip(vals, counts)}逐行说明一下。src.res[0] * src.res[1]得到像元面积单位是平方米除以1e6换成平方公里如果是 1000 米分辨率结果正好是 1.0。arr ! nodata这一步看起来简单但前提是nodata有值所以前面体检时发现nodata为None的要在这之前手动指定一个占位值。最后返回字典键是编码、值是面积后面聚合一级类、跨年份比较都基于这个结构。拿到二级类面积之后聚合到一级类def aggregate_level1(area_dict: dict) - dict: out {} for code, area in area_dict.items(): k code // 10 if code 10 else code out[k] out.get(k, 0.0) area return out这套写法在图幅边界上会有个细节跨图幅的像元被重复计数或者漏计。全国数据如果是一整幅不用管如果是分幅拼接出来的统计前要先做一次镶嵌别对着单幅挨个加。3.3 重采样到别的分辨率categorical 只用最近邻有人为了“让图更好看”或者“和别的 500 米数据对齐”会把这份 1 公里数据重采样到 500 米。技术上没问题但必须用最近邻。双线性会在两个编码之间插值比如 1耕地和 5建设用地中间插出 3而 3 在编码表里是草地——凭空多出一个不存在的草地像元面积统计立刻失真。import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path, dst_path LU_2015.tif, LU_2015_500m.tif with rasterio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, src.crs, src.width, src.height, *src.bounds, resolution500) meta src.meta.copy() meta.update(transformtransform, widthwidth, heightheight) with rasterio.open(dst_path, w, **meta) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crssrc.crs, resamplingResampling.nearest, # 类别数据只能最近邻 )如果要顺便把投影换成 Albers把上面的dst_crs从src.crs改成目标 CRS 就行calculate_default_transform会自动算出新的尺寸和变换参数。这里有个语义上的坑要说清楚重采样到 500 米不代表精度变成 500 米。原始 1 公里像元里本来就没有的边界信息插值插不出来只是把一张粗网格画成了细网格图好看了信息量没增加。做面积统计的时候重采样还会引入边界像元的类别替换误差全国范围能差出千分之几跨年份比较时这个误差如果方向一致趋势会被轻微带偏。我的习惯是统计永远在原始分辨率上做重采样只用于出图和叠置显示。4. 避坑排查五个把结果带偏的现场这套数据用久了翻车的场景会重复出现。下面五个是我自己和身边人踩过最多次的现象看起来各不相同根子上大多是“没对齐”或者“没声明清楚”。每一条按现象、原因、解决写遇到类似的可以先照着对一遍。4.1 两期数据地图上几乎没差统计出来差了一大截现象某省 1980 和 2015 两期叠在一起看地类边界几乎没有位移算出来的耕地转出却是两位数百分比。原因两期的网格原点或者尺寸对不上。有些年份的裁剪窗口差了半个像元地图上肉眼看不出来数组下标一对齐就成了错位比较每个像元都在和自己的邻居比。解决做比较之前先断言两期网格一致不一致就用重投影对齐到同一套 transform 上。import rasterio def assert_same_grid(a_path: str, b_path: str, tol: float 1e-6) - None: with rasterio.open(a_path) as a, rasterio.open(b_path) as b: assert a.crs b.crs, 两期坐标系不一致 for i in range(6): # transform 六个分量逐个比 assert abs(a.transform[i] - b.transform[i]) tol, ftransform[{i}] 不一致 assert (a.width, a.height) (b.width, b.height), 行列数不一致tol是米制下的容差1e-6 米对栅格来说可以认为完全相等如果两期本来就是不同投影这个断言会先报 CRS 不一致那就先做重投影再比。4.2 某县面积逐年跳一下别的县却平滑现象全省趋势平滑单独看某个县某一类面积每年跳 5% 到 10%像有人在改数据。原因裁切时没有统一 NoData。按行政区裁剪之后边界外的区域可能被填成 0也可能被填成 255还有可能保留原值。填 0 的话会被当成一个真实存在的“编码 0 地类”算进去面积自然跳。解决裁切时显式指定输出 NoData统计时用掩膜读取。with rasterio.open(county_2015.tif) as src: data src.read(1, maskedTrue) # maskedTrue 会把 nodata 自动掩掉 valid_ratio 1 - data.mask.mean() # 有效像元占比低于 0.8 就要查裁切参数 print(有效像元占比:, round(float(valid_ratio), 3))maskedTrue依赖栅格自己的 nodata 声明声明缺失时这个掩膜是空的所以裁切那一步必须把 NoData 写进输出文件不能只靠统计端补救。4.3 草地面积是别人软件算的两倍现象和同行对数同一范围同一年的草地面积自己算出来差不多是人家的两倍。原因一级类和二级类混着算。数据里是二级编码用arr 3去筛草地一个像元都筛不到改用(arr 30) (arr 33)才是对的。反过来如果数据是一级类版本用arr // 10聚合会把 1 到 6 全部映射成 0统计直接归零。解决先确认数据是哪个层级然后用统一口径聚合别在同一份统计里混用两种写法。判断方法在 2.3 的值域统计那一步已经做了值域里最大到 67 就是二级最大到 6 就是一级。4.4 转格式或者建金字塔之后冒出 13、27 这种不存在的编码现象把栅格转成别的格式、或者在 GIS 软件里建完金字塔再看直方图出现了一堆编码表里没有的值。原因概览图和格式转换默认用了平均或者双线性重采样。类别栅格的概览层必须是最近邻否则每一层金字塔都是一次插值污染。解决命令行建金字塔时显式指定算法转格式时也一样。# 显式用最近邻建金字塔层级 2/4/8/16 gdaladdo -r nearest LU_2015.tif 2 4 8 16 # 转 GeoPackage 时同样指定别用默认 gdal_translate -of GPKG -r nearest LU_2015.tif LU_2015.gpkg已经建错的先删掉.ovr附属文件再重建不删的话旧的污染概览层还在。4.5 NoData 没声明全国汇总面积明显偏大现象把全国所有地类面积加起来比国土面积大出一截或者某一年莫名其妙多出几万平方公里。原因栅格没有 nodata 属性海洋、境外、无数据区被填成了 0 或者某个占位数字统计时当成真实像元计入了。解决先看值域找出那个不该出现的编码再在脚本里把它当无效值剔除。如果占位值是 0而编码表里本来就没有 0那arr arr[arr 0]就够了。如果有多个占位值就显式列出白名单只保留编码表里存在的值VALID_CODES {11, 12, 21, 22, 23, 24, 31, 32, 33, 41, 42, 43, 44, 45, 46, 51, 52, 53, 61, 62, 63, 64, 65, 66, 67} arr arr[np.isin(arr, list(VALID_CODES))]白名单写法比黑名单稳因为占位值可能不止一个。代价是每次数据口径变化都要改白名单所以我把它放在配置区不散落在业务逻辑里。5. 从两期栅格到转移矩阵变化检测的完整链路面积统计回答的是“有多少”转移矩阵回答的是“从哪来到哪去”。后者才是土地利用变化分析真正有价值的部分耕地减少了多少不重要重要的是减少的这部分去了建设用地还是转成了林地两者的政策含义完全不同。这一章的链路是先对齐网格再用交集掩膜最后做交叉表。5.1 掩膜先对齐再谈变化两期数据直接做逐像元比较前提是它们在同一套网格上。第 4 章那个断言只能告诉你一不一致不能帮你对齐。实际做法是把后一期重采样到前一期的网格上最近邻然后取两期都有效的交集。import numpy as np import rasterio from rasterio.warp import reproject, Resampling def load_on_grid(path: str, ref_path: str): 把 path 对齐到 ref_path 的网格返回 (参考期数组, 对齐后数组) with rasterio.open(ref_path) as ref: h, w ref.height, ref.width transform, crs ref.transform, ref.crs ref_arr ref.read(1) with rasterio.open(path) as src: dst np.zeros((h, w), dtypesrc.dtypes[0]) reproject( sourcerasterio.band(src, 1), destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crscrs, resamplingResampling.nearest, src_nodatasrc.nodata, dst_nodata0, # 无效区统一填 0方便后面掩掉 ) return ref_arr, dst函数的第二个返回值是已经对齐到参考网格的数组dst_nodata0把重采样产生的空白统一成 0后面用 0就能一次性掩掉两期的无效区。如果两期本来就是同一套网格这一步等于做了一次恒等重采样开销可以忽略但保留了代码路径的一致性我不喜欢为“应该一致”的情况写两条分支。5.2 用 crosstab 生成转移矩阵有了对齐的两期数组交叉表一行代码就能出。但行列标签必须统一不然某些编码只在第一期出现、不在第二期出现时对角线会错位算出来的“未变化”面积是假的。import numpy as np import pandas as pd def transfer_matrix(a: np.ndarray, b: np.ndarray, valid: np.ndarray) - pd.DataFrame: a、b 已对齐的两期数组valid 为两期都有效的布尔掩膜 va, vb a[valid].ravel(), b[valid].ravel() codes np.union1d(np.unique(va), np.unique(vb)) # 两期编码的并集 tab pd.crosstab(pd.Series(va, namet1), pd.Series(vb, namet2)) # 用并集重排行列缺失的组合补 0保证行列标签顺序一致 return tab.reindex(indexcodes, columnscodes, fill_value0)np.union1d拿到的是两期所有编码的并集reindex把交叉表强行对齐到这个并集上缺的组合填 0。这一步看起来多余但没有它某些年份某类地类完全消失时行列长度不同np.trace取到的对角线就不是真正的“未变化”面积。调用的时候掩膜这样构造y1980, y2015 load_on_grid(LU_2015.tif, LU_1980.tif) valid (y1980 0) (y2015 0) # 两期都有效的像元才参与比较 tab transfer_matrix(y1980, y2015, valid)5.3 把矩阵读成人话净变化和总动量要分开看矩阵出来之后最容易犯的错是只看净变化。净变化等于行和减列和正数表示这类地增加了看起来直观但它会把方向相反的变化互相抵消。total int(tab.values.sum()) unchanged int(np.trace(tab.values)) print(f有效像元 {total:,}未变化 {unchanged:,}变化占比 {1 - unchanged/total:.2%}) net tab.sum(axis1) - tab.sum(axis0) # 每类净变化正数为增加 print(net.sort_values())假设某地一年内 100 平方公里耕地转成了林地同时 90 平方公里林地转成了耕地净变化只剩“耕地减少 10 平方公里”。报上去的数字看着很温和实际上发生了 190 平方公里的地类互换。所以做区域评估时我一般同时报三个数净变化、总转出、总转入。总转出等于行和减对角线总转入等于列和减对角线这三个数放一起读者才知道这块地到底“动”了多少。另外转移矩阵最好在一级类和二级类两个层级各做一次。一级类矩阵适合对外汇报结论干净二级类矩阵能看出耕地内部水田旱地的互换这种变化在一级类里是看不见的。6. 收尾一个脚本把全部年份体检完顺便做三个守恒检查一份数据七个年份、三十来个分幅手工一个个打开看元数据看完第一个小时人就麻了而且一定会漏。我现在拿到新一期的土地利用数据会先跑一个批量体检把每年每幅的关键属性摊成一张表再用几条断言卡住明显的异常。表里的字段就是前面几章反复提到的那些年份、坐标系、行列数、像元面积、有效像元数、总面积。import pandas as pd def audit(year_map: dict) - pd.DataFrame: rows [] for year, path in sorted(year_map.items()): with rasterio.open(path) as src: arr src.read(1) nodata src.nodata if src.nodata is not None else 0 valid arr ! nodata px_km2 src.res[0] * src.res[1] / 1e6 rows.append({ year: year, crs: src.crs.to_string() if src.crs else None, w: src.width, h: src.height, pixel_km2: round(px_km2, 4), valid_px: int(valid.sum()), total_km2: round(valid.sum() * px_km2, 1), }) return pd.DataFrame(rows) df audit(year_map) print(df.to_string(indexFalse))三个检查我每次都跑成本几乎为零但救过我好几次。第一所有年份的crs必须只有一个取值出现两个就说明混进了不同投影的版本后面统计全废。第二total_km2的年际波动要在一个合理区间内全国范围用最大最小比小于 1.05 做粗筛超过就去看是不是某一年边界外被算进来了。第三pixel_km2在所有年份必须一致出现 0.98 或者 1.02 这种说明有年份被重采样过了原始分辨率信息已经不可逆地丢了。assert df[crs].nunique() 1, 存在多套坐标系 assert df[total_km2].max() / df[total_km2].min() 1.05, 总面积年际波动异常 assert df[pixel_km2].nunique() 1, 像元尺寸不一致可能被重采样过断言失败的时候不要急着改阈值先去看失败那一年的原始文件绝大多数情况是分发批次不同、参数被打乱重新取一份对应版本的就行。阈值本身是经验值全国尺度上海洋和境外区域基本不变1.05 的宽松度足够覆盖真实的地类变化。血泪经验是最开始我图省事把面积统计和变化检测写在一个脚本里结果中间某一期投影不对一路算到最后才报错排查花了整整一个下午。从那以后我每次拿到新一期的土地利用数据都强制先跑一遍这个体检表和三条断言通过了再进分析流程。看着啰嗦其实是给自己留了后悔药。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/11/5 5:37:53 · 阅读 1.2k 1km土地利用栅格：编码、投影与转移矩阵避坑
简介中国1km土地利用数据1980—2015压缩包面向地理学、环境科学、城乡规划与资源管理等领域的研究者和学生提供覆盖1980至2015年多个时点的1km格网土地利用类型信息可用于分析近四十年土地变迁、城市扩张、生态保护与农业动态。包内共126个文件以adf、nit、dat、dir等ArcGIS栅格单元文件为主并含log、doc、xml等说明与辅助文件压缩包约21.83MB文件按年份与数据集组织。已有2904人学习下载。数据以栅格形式记录耕地、林地、草地、建设用地、水域等类型便于对比不同年份识别显著变化也可与社会经济数据结合构建土地利用变化模型。对开展区域案例研究、政策评估和可持续分析具有较高参考价值。1. 1km 分辨率能干什么先把这套土地利用栅格的适用边界说清楚一个做县域生态评估的 A 同学拿到 1980–2015 的中国 1km 土地利用数据之后直接按像元数乘以 1 平方公里去算耕地转出的面积。第一次跑出县城建成区扩张 8%换一个投影再跑变成 3%。他以为是数据坏了其实这份土地利用数据里的问题几乎都集中在“单位”和“编码”这两件事上。这套资源是 1980 到 2015 年逐期的全国 1 公里栅格按年份分包解压之后就是可以直接读进 GDAL 的类别栅格省掉了自己去拼接、配准、对齐的整段流程。它适合做长时序、大范围的地类结构分析比如耕地流向、城镇化扩张、林地草地的年际波动也适合给生态模型当驱动数据。但它不适合地块级权属判断也分不清城市内部的用地功能。1 公里就是一个像元一平方公里边界模糊、细小地物会被吞掉这是分辨率本身决定的跟数据质量无关。2. 解压之后先别画图把目录、编码表和元数据核一遍拿到一个压缩包最省时间的做法不是先把它塞进 QGIS 看一眼出图效果而是花二十分钟把三件事确定下来文件是怎么组织的、编码是什么意思、栅格自己的元数据靠不靠谱。这三件事没定下来就动手算面积后面所有的数字都是悬空的。我把这一步叫体检体检做完后面写脚本才有依据。2.1 摊开压缩包年份目录与文件命名的识别规律先看清单再解压不要一口气全铺开。很多全国尺度的栅格包体量大解压完才发现里面是 ArcInfo Grid 那种一堆.adf的目录结构单独拷一个文件出去就坏了得整个目录一起搬。用清单先把结构摸清能省掉一次返工。# 先看清单确认里面是单文件栅格tif/img还是 ArcInfo Grid 目录 unzip -l 中国1km土地利用数据1980-2015.zip | head -40 # 解压到独立目录别和别的项目数据混在一起 unzip -q -d ./lu_raw 中国1km土地利用数据1980-2015.zip # 列出所有栅格候选按路径排序后年份顺序一眼可见 find ./lu_raw -type f \( -iname *.tif -o -iname *.img \) | sortWindows 上没有unzip就用 7-Zip 的图形界面做同样的事重点是先看清单、再解压到独立目录。清单里如果出现.prj、.tfw、.aux.xml这些小文件说明栅格带坐标和统计附属信息搬家的时候要跟主文件一起走漏了.prj有些软件会默认按像素坐标处理位置直接跑偏。解压完把年份抽出来建一张“年份 → 文件路径”的映射表。后面所有脚本都靠这张表循环不要在每个脚本里重新拼文件名。import re from pathlib import Path files sorted(Path(./lu_raw).rglob(*.tif)) year_map {} for f in files: m re.search(r(19|20)\d{2}, f.stem) # 从文件名里抓四位年份 if m: year_map[int(m.group(0))] f print(sorted(year_map)) # 期望看到 1980、1990、1995、2000、2005、2010、2015 之类的序列这段代码的逻辑很直白遍历解压目录下所有.tif用正则从文件名主干里抓年份建成字典。如果打印出来的年份少于预期说明命名规则不是纯年份可能带了版本后缀或者分幅编号这时候先把files的前十个路径打出来看规律别硬套正则。年份数量对不上是最容易被忽略的第一个坑少一期数据做时序分析结论能差出好几年。2.2 分类编码表一级类看趋势二级类看结构这一类全国土地利用栅格主流用的是两级分类一级类 6 个二级类二十多个用一个整数编码表示。一级类够做趋势和总量二级类才用来看内部结构比如耕地里水田和旱地的比例、林地里是有林地还是灌木林。一级类编码常见二级类与编码耕地111 水田、12 旱地林地221 有林地、22 灌木林、23 疏林地、24 其他林地草地331 高覆盖度草地、32 中覆盖度草地、33 低覆盖度草地水域441 河渠、42 湖泊、43 水库坑塘、44 永久性冰川雪地、45 滩涂、46 滩地城乡工矿居民用地551 城镇用地、52 农村居民点、53 其他建设用地未利用土地661 沙地、62 戈壁、63 盐碱地、64 沼泽地、65 裸土地、66 裸岩石质地、67 其他这张表要当成“常见口径”来用不能当成“这份数据一定如此”来用。不同分发批次里二级类可能做过合并也可能某些年份只给了一级类编码。我一般会先跑一遍值域统计看实际出现了哪些整数再和这张表对照出现的编码是不是都在表里、有没有超出 67 的值。有超出就说明这一期用了别的口径得单独处理。提示先跑值域统计再决定用哪套编码解释数据比对着文档猜快得多也准得多。把二级编码折算回一级最省事的写法是整除def to_level1(code: int) - int: 两位二级编码折算成一级编码11-143-467-6 return code // 10 if code 10 else code注意这个函数只对两位编码成立。如果某期数据里混进了三位编码整除会算出 10 以上的值这时候要改成查表。这也是为什么我坚持先看值域函数写得再漂亮口径不对一样算错。2.3 用 gdalinfo 和 rasterio 给数据做一次体检体检的核心是三件事坐标系是什么、像元多大、NoData 是怎么声明的。这三项决定了后面面积能不能直接乘。gdalinfo -stats LU_2015.tif-stats会顺带算出最小最大值和直方图输出里的Coordinate System is、Pixel Size、NoData Value三个字段要重点看。如果Pixel Size显示的是类似(0.00833, -0.00833)这种小数说明这份数据是经纬度网格不是米制投影——这一点非常关键后面第 3 章会展开。再用 rasterio 批量体检一遍把关键属性整理成能直接比较的结构import rasterio import numpy as np with rasterio.open(LU_2015.tif) as src: print(CRS :, src.crs) print(宽高 :, src.width, src.height) print(像元尺寸 :, src.res) # (x, y)单位跟随 CRS print(NoData :, src.nodata) print(数据类型 :, src.dtypes[0]) arr src.read(1) vals, counts np.unique(arr, return_countsTrue) for v, c in zip(vals, counts): print(f编码 {int(v):3} 像元 {int(c):12,})src.res返回的是像元在x、y方向上的尺寸单位跟 CRS 走米制投影就是米经纬度就是度。src.nodata如果打印出None说明栅格没有声明无效值这是最常见的隐患之一海洋、境外区域可能是 0 或者 255 这种占位数字会被当成真实地类统计进去。np.unique那一段是看值域分布顺便能发现数据里到底有没有二级编码如果只出现 1 到 6 这六个数那这份就是一级类版本别按二级类去拆。体检做完应该能回答三个问题这是等积投影还是经纬度网格、像元实际物理面积是多少、无效值怎么处理。三个问题都有答案了再往下算面积才不会翻车。3. 面积算不对的根源投影、像元面积与 categorical 重采样面积统计是这套数据被用得最多的场景也是错得最多的场景。错法通常有三种拿经纬度网格直接当等面积网格用、按像元计数乘一个固定系数、重采样的时候用了插值算法。这三种错误各有各的隐蔽性跑出来的数字看着都很像那么回事。3.1 Albers 等积投影与经纬度网格同一块地的两个答案全国尺度的土地利用栅格常见做法是采用 Albers 等积圆锥投影一套典型参数是双标准纬线 25°N 和 47°N、中央经线 105°E。选等积投影的理由很实际土地利用统计只关心面积不关心形状和距离是否保真等积投影能保证任意区域的面积比例不变做汇总时才不会被投影变形带偏。如果拿到的是发布方给的 Albers 版本像元是 1000 米 × 1000 米一个像元就是一平方公里面积直接乘就行。麻烦的是经纬度版本。经纬度网格里一个像元覆盖的东西向距离随纬度收缩同样“1 度 × 1 度”的方格在 30°N 和 45°N 覆盖的实际面积能差两成左右。粗算像元物理面积的公式是这样import math def cell_area_km2(px_x_deg: float, px_y_deg: float, lat_center: float) - float: 经纬度网格单个像元的近似面积km² dx abs(px_x_deg) * 111.320 * math.cos(math.radians(lat_center)) # 经向实际距离 dy abs(px_y_deg) * 111.132 # 纬向实际距离 return dx * dypx_x_deg、px_y_deg是像元的经纬度尺寸lat_center是像元中心纬度。这个近似在单个像元上误差不大但做全国汇总时不同纬度带要用不同系数逐像元算一遍才准写起来不难但容易漏。我一般会直接绕开这个问题先在等积投影下统计实在要用经纬度版本就先重投影过去再算面积。注意判断依据是src.crs.is_geographic。返回True就是经纬度网格这时候把像元数乘以一个固定的平方公里数一定会偏。3.2 用 rasterio 批量统计分类面积确定了投影和像元面积统计就是纯计数。这段代码我用了很多次改的主要是输入输出路径import numpy as np import rasterio def class_area_km2(tif_path: str, pixel_area: float None) - dict: 返回 {编码: 面积km²}。pixel_area 为 None 时按米制像元自动推算。 with rasterio.open(tif_path) as src: arr src.read(1) nodata src.nodata # 米制投影下像元面积 res_x * res_y除以 1e6 换成 km² if pixel_area is None: pixel_area src.res[0] * src.res[1] / 1e6 if nodata is not None: arr arr[arr ! nodata] # 剔除无效值别让占位数字混进统计 vals, counts np.unique(arr, return_countsTrue) return {int(v): float(c * pixel_area) for v, c in zip(vals, counts)}逐行说明一下。src.res[0] * src.res[1]得到像元面积单位是平方米除以1e6换成平方公里如果是 1000 米分辨率结果正好是 1.0。arr ! nodata这一步看起来简单但前提是nodata有值所以前面体检时发现nodata为None的要在这之前手动指定一个占位值。最后返回字典键是编码、值是面积后面聚合一级类、跨年份比较都基于这个结构。拿到二级类面积之后聚合到一级类def aggregate_level1(area_dict: dict) - dict: out {} for code, area in area_dict.items(): k code // 10 if code 10 else code out[k] out.get(k, 0.0) area return out这套写法在图幅边界上会有个细节跨图幅的像元被重复计数或者漏计。全国数据如果是一整幅不用管如果是分幅拼接出来的统计前要先做一次镶嵌别对着单幅挨个加。3.3 重采样到别的分辨率categorical 只用最近邻有人为了“让图更好看”或者“和别的 500 米数据对齐”会把这份 1 公里数据重采样到 500 米。技术上没问题但必须用最近邻。双线性会在两个编码之间插值比如 1耕地和 5建设用地中间插出 3而 3 在编码表里是草地——凭空多出一个不存在的草地像元面积统计立刻失真。import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path, dst_path LU_2015.tif, LU_2015_500m.tif with rasterio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, src.crs, src.width, src.height, *src.bounds, resolution500) meta src.meta.copy() meta.update(transformtransform, widthwidth, heightheight) with rasterio.open(dst_path, w, **meta) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crssrc.crs, resamplingResampling.nearest, # 类别数据只能最近邻 )如果要顺便把投影换成 Albers把上面的dst_crs从src.crs改成目标 CRS 就行calculate_default_transform会自动算出新的尺寸和变换参数。这里有个语义上的坑要说清楚重采样到 500 米不代表精度变成 500 米。原始 1 公里像元里本来就没有的边界信息插值插不出来只是把一张粗网格画成了细网格图好看了信息量没增加。做面积统计的时候重采样还会引入边界像元的类别替换误差全国范围能差出千分之几跨年份比较时这个误差如果方向一致趋势会被轻微带偏。我的习惯是统计永远在原始分辨率上做重采样只用于出图和叠置显示。4. 避坑排查五个把结果带偏的现场这套数据用久了翻车的场景会重复出现。下面五个是我自己和身边人踩过最多次的现象看起来各不相同根子上大多是“没对齐”或者“没声明清楚”。每一条按现象、原因、解决写遇到类似的可以先照着对一遍。4.1 两期数据地图上几乎没差统计出来差了一大截现象某省 1980 和 2015 两期叠在一起看地类边界几乎没有位移算出来的耕地转出却是两位数百分比。原因两期的网格原点或者尺寸对不上。有些年份的裁剪窗口差了半个像元地图上肉眼看不出来数组下标一对齐就成了错位比较每个像元都在和自己的邻居比。解决做比较之前先断言两期网格一致不一致就用重投影对齐到同一套 transform 上。import rasterio def assert_same_grid(a_path: str, b_path: str, tol: float 1e-6) - None: with rasterio.open(a_path) as a, rasterio.open(b_path) as b: assert a.crs b.crs, 两期坐标系不一致 for i in range(6): # transform 六个分量逐个比 assert abs(a.transform[i] - b.transform[i]) tol, ftransform[{i}] 不一致 assert (a.width, a.height) (b.width, b.height), 行列数不一致tol是米制下的容差1e-6 米对栅格来说可以认为完全相等如果两期本来就是不同投影这个断言会先报 CRS 不一致那就先做重投影再比。4.2 某县面积逐年跳一下别的县却平滑现象全省趋势平滑单独看某个县某一类面积每年跳 5% 到 10%像有人在改数据。原因裁切时没有统一 NoData。按行政区裁剪之后边界外的区域可能被填成 0也可能被填成 255还有可能保留原值。填 0 的话会被当成一个真实存在的“编码 0 地类”算进去面积自然跳。解决裁切时显式指定输出 NoData统计时用掩膜读取。with rasterio.open(county_2015.tif) as src: data src.read(1, maskedTrue) # maskedTrue 会把 nodata 自动掩掉 valid_ratio 1 - data.mask.mean() # 有效像元占比低于 0.8 就要查裁切参数 print(有效像元占比:, round(float(valid_ratio), 3))maskedTrue依赖栅格自己的 nodata 声明声明缺失时这个掩膜是空的所以裁切那一步必须把 NoData 写进输出文件不能只靠统计端补救。4.3 草地面积是别人软件算的两倍现象和同行对数同一范围同一年的草地面积自己算出来差不多是人家的两倍。原因一级类和二级类混着算。数据里是二级编码用arr 3去筛草地一个像元都筛不到改用(arr 30) (arr 33)才是对的。反过来如果数据是一级类版本用arr // 10聚合会把 1 到 6 全部映射成 0统计直接归零。解决先确认数据是哪个层级然后用统一口径聚合别在同一份统计里混用两种写法。判断方法在 2.3 的值域统计那一步已经做了值域里最大到 67 就是二级最大到 6 就是一级。4.4 转格式或者建金字塔之后冒出 13、27 这种不存在的编码现象把栅格转成别的格式、或者在 GIS 软件里建完金字塔再看直方图出现了一堆编码表里没有的值。原因概览图和格式转换默认用了平均或者双线性重采样。类别栅格的概览层必须是最近邻否则每一层金字塔都是一次插值污染。解决命令行建金字塔时显式指定算法转格式时也一样。# 显式用最近邻建金字塔层级 2/4/8/16 gdaladdo -r nearest LU_2015.tif 2 4 8 16 # 转 GeoPackage 时同样指定别用默认 gdal_translate -of GPKG -r nearest LU_2015.tif LU_2015.gpkg已经建错的先删掉.ovr附属文件再重建不删的话旧的污染概览层还在。4.5 NoData 没声明全国汇总面积明显偏大现象把全国所有地类面积加起来比国土面积大出一截或者某一年莫名其妙多出几万平方公里。原因栅格没有 nodata 属性海洋、境外、无数据区被填成了 0 或者某个占位数字统计时当成真实像元计入了。解决先看值域找出那个不该出现的编码再在脚本里把它当无效值剔除。如果占位值是 0而编码表里本来就没有 0那arr arr[arr 0]就够了。如果有多个占位值就显式列出白名单只保留编码表里存在的值VALID_CODES {11, 12, 21, 22, 23, 24, 31, 32, 33, 41, 42, 43, 44, 45, 46, 51, 52, 53, 61, 62, 63, 64, 65, 66, 67} arr arr[np.isin(arr, list(VALID_CODES))]白名单写法比黑名单稳因为占位值可能不止一个。代价是每次数据口径变化都要改白名单所以我把它放在配置区不散落在业务逻辑里。5. 从两期栅格到转移矩阵变化检测的完整链路面积统计回答的是“有多少”转移矩阵回答的是“从哪来到哪去”。后者才是土地利用变化分析真正有价值的部分耕地减少了多少不重要重要的是减少的这部分去了建设用地还是转成了林地两者的政策含义完全不同。这一章的链路是先对齐网格再用交集掩膜最后做交叉表。5.1 掩膜先对齐再谈变化两期数据直接做逐像元比较前提是它们在同一套网格上。第 4 章那个断言只能告诉你一不一致不能帮你对齐。实际做法是把后一期重采样到前一期的网格上最近邻然后取两期都有效的交集。import numpy as np import rasterio from rasterio.warp import reproject, Resampling def load_on_grid(path: str, ref_path: str): 把 path 对齐到 ref_path 的网格返回 (参考期数组, 对齐后数组) with rasterio.open(ref_path) as ref: h, w ref.height, ref.width transform, crs ref.transform, ref.crs ref_arr ref.read(1) with rasterio.open(path) as src: dst np.zeros((h, w), dtypesrc.dtypes[0]) reproject( sourcerasterio.band(src, 1), destinationdst, src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crscrs, resamplingResampling.nearest, src_nodatasrc.nodata, dst_nodata0, # 无效区统一填 0方便后面掩掉 ) return ref_arr, dst函数的第二个返回值是已经对齐到参考网格的数组dst_nodata0把重采样产生的空白统一成 0后面用 0就能一次性掩掉两期的无效区。如果两期本来就是同一套网格这一步等于做了一次恒等重采样开销可以忽略但保留了代码路径的一致性我不喜欢为“应该一致”的情况写两条分支。5.2 用 crosstab 生成转移矩阵有了对齐的两期数组交叉表一行代码就能出。但行列标签必须统一不然某些编码只在第一期出现、不在第二期出现时对角线会错位算出来的“未变化”面积是假的。import numpy as np import pandas as pd def transfer_matrix(a: np.ndarray, b: np.ndarray, valid: np.ndarray) - pd.DataFrame: a、b 已对齐的两期数组valid 为两期都有效的布尔掩膜 va, vb a[valid].ravel(), b[valid].ravel() codes np.union1d(np.unique(va), np.unique(vb)) # 两期编码的并集 tab pd.crosstab(pd.Series(va, namet1), pd.Series(vb, namet2)) # 用并集重排行列缺失的组合补 0保证行列标签顺序一致 return tab.reindex(indexcodes, columnscodes, fill_value0)np.union1d拿到的是两期所有编码的并集reindex把交叉表强行对齐到这个并集上缺的组合填 0。这一步看起来多余但没有它某些年份某类地类完全消失时行列长度不同np.trace取到的对角线就不是真正的“未变化”面积。调用的时候掩膜这样构造y1980, y2015 load_on_grid(LU_2015.tif, LU_1980.tif) valid (y1980 0) (y2015 0) # 两期都有效的像元才参与比较 tab transfer_matrix(y1980, y2015, valid)5.3 把矩阵读成人话净变化和总动量要分开看矩阵出来之后最容易犯的错是只看净变化。净变化等于行和减列和正数表示这类地增加了看起来直观但它会把方向相反的变化互相抵消。total int(tab.values.sum()) unchanged int(np.trace(tab.values)) print(f有效像元 {total:,}未变化 {unchanged:,}变化占比 {1 - unchanged/total:.2%}) net tab.sum(axis1) - tab.sum(axis0) # 每类净变化正数为增加 print(net.sort_values())假设某地一年内 100 平方公里耕地转成了林地同时 90 平方公里林地转成了耕地净变化只剩“耕地减少 10 平方公里”。报上去的数字看着很温和实际上发生了 190 平方公里的地类互换。所以做区域评估时我一般同时报三个数净变化、总转出、总转入。总转出等于行和减对角线总转入等于列和减对角线这三个数放一起读者才知道这块地到底“动”了多少。另外转移矩阵最好在一级类和二级类两个层级各做一次。一级类矩阵适合对外汇报结论干净二级类矩阵能看出耕地内部水田旱地的互换这种变化在一级类里是看不见的。6. 收尾一个脚本把全部年份体检完顺便做三个守恒检查一份数据七个年份、三十来个分幅手工一个个打开看元数据看完第一个小时人就麻了而且一定会漏。我现在拿到新一期的土地利用数据会先跑一个批量体检把每年每幅的关键属性摊成一张表再用几条断言卡住明显的异常。表里的字段就是前面几章反复提到的那些年份、坐标系、行列数、像元面积、有效像元数、总面积。import pandas as pd def audit(year_map: dict) - pd.DataFrame: rows [] for year, path in sorted(year_map.items()): with rasterio.open(path) as src: arr src.read(1) nodata src.nodata if src.nodata is not None else 0 valid arr ! nodata px_km2 src.res[0] * src.res[1] / 1e6 rows.append({ year: year, crs: src.crs.to_string() if src.crs else None, w: src.width, h: src.height, pixel_km2: round(px_km2, 4), valid_px: int(valid.sum()), total_km2: round(valid.sum() * px_km2, 1), }) return pd.DataFrame(rows) df audit(year_map) print(df.to_string(indexFalse))三个检查我每次都跑成本几乎为零但救过我好几次。第一所有年份的crs必须只有一个取值出现两个就说明混进了不同投影的版本后面统计全废。第二total_km2的年际波动要在一个合理区间内全国范围用最大最小比小于 1.05 做粗筛超过就去看是不是某一年边界外被算进来了。第三pixel_km2在所有年份必须一致出现 0.98 或者 1.02 这种说明有年份被重采样过了原始分辨率信息已经不可逆地丢了。assert df[crs].nunique() 1, 存在多套坐标系 assert df[total_km2].max() / df[total_km2].min() 1.05, 总面积年际波动异常 assert df[pixel_km2].nunique() 1, 像元尺寸不一致可能被重采样过断言失败的时候不要急着改阈值先去看失败那一年的原始文件绝大多数情况是分发批次不同、参数被打乱重新取一份对应版本的就行。阈值本身是经验值全国尺度上海洋和境外区域基本不变1.05 的宽松度足够覆盖真实的地类变化。血泪经验是最开始我图省事把面积统计和变化检测写在一个脚本里结果中间某一期投影不对一路算到最后才报错排查花了整整一个下午。从那以后我每次拿到新一期的土地利用数据都强制先跑一遍这个体检表和三条断言通过了再进分析流程。看着啰嗦其实是给自己留了后悔药。希望帮到你。本文还有配套的精品资源点击获取