去年帮一位做城市扩张研究的朋友处理人口栅格他拿到一份 2000—2023 年的 1km 精度人口分布栅格数据第一件事就是把研究区内的像元值全部加总想直接得到总人口曲线。结果算出来的数比公开统计口径少了将近三成他一度以为是数据本身有问题。后来一起排查问题出在两个地方栅格值的单位是“人/像元”而不是“人/平方公里”而他在裁剪时又用了all_touchedTrue把边界外侧的像元也带了进来正负抵消之后偏差就变得很随机。这件事很典型——1km 精度人口分布栅格数据这类东西看起来只是“下载一个 TIFF 然后加总”真正用起来坐标系、单位、nodata、边界版本、产品版本这几件事随便哪一件没对齐结论就会跑偏。这篇想聊的就是这类数据从下载到出图之间的那段路。它面向的是拿到全球、全国、分省、分市四种粒度人口栅格、准备做时间序列分析或空间统计的人不管你是刚接触栅格的新手还是已经用惯了 GDAL 的老手里面这些口径和坑大概率都会碰上。我会把 1km 这个分辨率到底意味着什么、数据是怎么被生产出来的、解压后该先看哪几个字段、按行政边界裁剪时误差从哪来、24 年叠成一条曲线要注意什么一项一项拆开讲。1. “1km 分辨率”不是精度只是格子大小1.1 一个像元里装了多少人决定了你能问什么问题很多人第一次看到“1km 精度”这个词下意识会理解成“每个格子里的人口数误差在 1 公里以内”这个理解是错的。1km 描述的是空间采样间隔也就是栅格的行列间距它跟数值准确度没有直接关系。真正决定一个像元代表多少人的是产品采用的分配模型在城市核心区一个 1km 像元可能对应几万人到了荒漠或高山同一个像元里可能是 0。这个差别会直接影响分析口径。如果你的研究对象是“某地级市人口总量”1km 栅格完全够用因为大数相加时随机误差会互相抵消但如果你要做“某街道内部人口密度梯度”1km 格子在建成区里只有一两个像元样本量太小任何一点分配误差都会被放大成显著的假梯度。我一般会跟人强调一句1km 产品适合回答“区域级”和“趋势级”的问题不适合回答“街区级”和“单一像元级”的问题。还有一层容易忽略经纬度坐标下的 1km 像元并不是正方形。在赤道附近0.008333333°也就是 30 弧秒的经纬跨度大约对应 927 米到了北纬 60 度经度方向的实地宽度收缩到大约 460 米纬度方向仍然是 927 米左右。所以同一个“1km 产品”在不同纬度上的实际覆盖面积能差一倍以上。如果你的栅格值是密度人/平方公里跨纬度求和必须先算每个纬度带的真实像元面积如果栅格值本身就是人数人/像元那直接求和没问题但换算成密度时同样要除以真实面积。1.2 人数型和密度型一行代码就能分出来拿到一份新数据我习惯先做一次“量纲体检”读一个像元值看它的数量级。全球人口约八十亿全球陆地面积约 1.5 亿平方公里平均密度约 50 人/平方公里。如果你读到的典型陆地像元值是几十、几百那大概率是密度型如果是零点几、几那大概率是人数型。这个判断不能靠猜最稳的办法是看产品附带的说明文档但很多二次分发的数据包会把文档丢掉。这时候可以用总量反推把全球或全国的像元值加总如果结果落在几十亿量级那就是人数型如果落在几千万到几亿量级那基本是密度型需要乘面积。我在实操中遇到过一份数据全球加总得到 78 亿直接就是人数型用起来最省事也遇到过加总得到 1300 多万的那个 1300 万其实是“人/平方公里”的总和没有任何物理意义必须乘面积。类型单位全球加总的量级换算方式人数型人/像元十亿量级如 70~85 亿直接求和密度型人/平方公里千万到亿量级每像元乘真实面积后求和归一化型无量纲或 0~1与像元数同量级需按区域总量缩放提示不要用“最小值是 0”来判断 nodata。很多产品把海洋和无人区都填成 0nodata 另有标记值常见的是 -99999、-3.4e38 或 NaN。把 nodata 当 0 参与求和结果不会报错但会悄悄拉低总量。1.3 分辨率数字背后的产品代际差异市面上能拿到 2000—2023 逐年序列的 1km 人口栅格通常不是一套模型从头跑到尾的产物。早年份2000—2010往往依赖当时的普查底账和较低分辨率的辅助数据中期2010—2015开始引入建筑轮廓和更细的土地利用近些年2016 以后则大量使用遥感建成区提取和机器学习分配。这意味着同一个产品内部不同年份的分配逻辑可能存在系统性差异。我在做长序列分析时踩过一次坑某产品 2012 年和 2013 年之间全国总量只涨了不到 1%但几个大城市的中心城区像元值却出现了明显跳变。查了说明才知道2013 年换了一版建成区底图城市内部的分配权重整体重算过。这种跳变不是人口真实变化是产品版本变化。所以做年际对比之前先做一遍“总量平滑性检查”很有必要把逐年总和画成折线正常的人口曲线应该是平缓的如果某一年出现尖峰或台阶先怀疑数据版本再怀疑现实。2. 从普查表到栅格这套数据是怎么“造”出来的2.1 底账是普查但普查最小单元通常比 1km 粗所有人口栅格产品的源头都是人口普查或户籍登记问题在于普查公布的最小空间单元往往比 1km 粗得多。有的地方能到社区级有的地方只能到乡镇级还有一些地区只公布县级汇总。1km 栅格要做的事就是把这些“面状底账”拆到“网格”上。这个拆分过程在业内叫分区密度制图dasymetric mapping。思路不复杂假设人口在普查单元内部不是均匀分布的而是跟着“能住人的地方”走。能住人的地方怎么识别靠辅助图层——建筑轮廓、建成区范围、夜间灯光、土地利用分类、道路密度、坡度、水体掩膜。把普查单元的总人口按这些辅助图层给出的权重分配到单元内的各个像元上就得到了一张看起来连续的栅格。理解这一点很关键因为它解释了这类数据的误差结构总量在普查单元级别是准的越往像元级别拆误差越大。换句话说把一整个省的所有像元加起来结果会非常接近该省普查总量但单独看某一个 1km 像元它可能只是“看起来合理”的一个估计值不代表真实观测。2.2 辅助图层选什么决定了城乡边界的锐利程度辅助图层的选择直接决定产品的“质感”。只用土地利用分类的城乡过渡带会显得很糊因为分类本身是 30 米到 300 米不等的粗糙多边形引入建筑轮廓的城市内部会出现非常锐利的高低值结构因为建筑物确实有明确的占地面积引入夜间灯光的城市外围的稀疏聚落会被适度提亮但也会把一些工业区和主干道误判成居住区。我在对比两套产品时做过一个简单的检验在同一片城乡结合部取一条 20 公里的剖面看像元值的衰减曲线。只用土地利用的产品曲线是阶梯状下降中间有好几段平台用建筑轮廓的产品曲线是连续衰减能看出村落一个个散开的节奏。前者适合做总量统计后者适合做形态分析。选哪套取决于你要回答的问题不存在绝对更好的产品。这里还有一个常被忽略的细节水体掩膜。有些产品会把大型水库、宽阔河道直接从可居住区域里剔除有些则保留。如果你研究的区域有大面积水域两套产品的总量差异可能就来自这里。做跨产品对比前先把水体掩膜对齐不然比出来的差异全是口径差。2.3 时间序列里的“插值年份”和“实测年份”要分开看2000—2023 是 24 个年份但普查不是每年都做。很多产品的做法是在有普查的年份用真实底账在没有普查的年份用插值或回归外推。有些产品说明书里会标出哪些年份是“基准年”哪些是“估计年”有些则完全不标。我一般会自己做一个粗略判断看逐年总量曲线的二阶差分。基准年之间的插值段曲线会异常平滑几乎是一条直线真实年份则会有小幅波动。如果发现某几段过于平滑就要在分析里降低这些年份的权重或者至少在做趋势拟合时加一个“年份来源”的虚拟变量。另外要留意行政区划调整。2000 年和 2023 年的县级边界不是同一套如果产品是按“当年行政区划”组织的分省分市文件那么同一个文件名在两个年份里可能对应不同的地理范围。做长序列时我通常会用一套固定的边界比如最近年份的边界去裁所有年份虽然会有少量错位但至少口径统一。3. 解压后的第一小时元数据体检清单3.1 用五行代码把关键字段全打出来拿到 GeoTIFF 之后别急着在 QGIS 里拖进去看颜色。先跑一段元数据检查把这几个字段确认清楚坐标系、像元大小、行列数、数据类型、nodata 值、是否分块、是否带金字塔。这六项里任何一项不对后面的统计都可能白做。import rasterio path pop_2020_1km.tif with rasterio.open(path) as src: print(CRS :, src.crs) print(SIZE :, src.width, x, src.height) print(PIXEL :, src.transform.a, src.transform.e) print(BOUNDS :, src.bounds) print(DTYPE :, src.dtypes) print(NODATA :, src.nodata) print(BLOCKS :, src.block_shapes[0]) print(OVERVIEWS:, src.overviews(1))重点看PIXEL那一行。如果是0.008333333说明是 30 弧秒的经纬度栅格需要按纬度带算真实面积如果是1000.0附近的小数说明用的是米制投影像元接近正方形换算简单得多。这两类产品的处理代码不能混用。再看NODATA。有些产品这里返回None但实际数据里存在-99999这样的哨兵值。这种情况必须手工指定 nodata否则求和时会把 -99999 当成真实人口加进去总量直接变成负数。我建议在所有读取操作里统一用maskedTrue然后显式检查掩膜比例with rasterio.open(path) as src: data src.read(1, maskedTrue) valid_ratio 1 - data.mask.mean() print(有效像元占比:, round(float(valid_ratio), 4)) print(总量:, float(data.sum()))全球陆地占比大约在 29% 左右如果你算出来的有效像元比例远高于这个数说明海洋没被掩掉远低于这个数说明掩膜过度可能有陆地被误判成 nodata。3.2 全球单文件有多大内存怎么安排全球 1km 经纬度栅格的行列数大致是 21600 行 × 43200 列约 9.3 亿个像元。Float32 存下来是 3.7 GB 左右Int32 也是这个量级。单文件读进内存对普通工作站是有压力的尤其是你要同时开好几年的数据做对比。更稳妥的做法是分块流式读取用block_windows按窗口遍历边读边累加内存占用只跟窗口大小有关def block_sum(path, band1): total 0.0 count 0 with rasterio.open(path) as src: for _, window in src.block_windows(band): arr src.read(band, windowwindow, maskedTrue) total float(arr.sum()) count int((~arr.mask).sum()) return total, count total, n block_sum(pop_global_2020_1km.tif) print(全球总量:, total, 有效像元:, n)如果文件是按 strip 存的block_shapes返回(1, 43200)这种说明没做分块优化读取效率会很差。这种情况我一般先用gdal_translate重写成 tiled 格式加上压缩体积能降一半以上读取速度也能提上来gdal_translate -of GTiff \ -co TILEDYES -co COMPRESSDEFLATE -co PREDICTOR3 \ -co BIGTIFFIF_SAFER \ pop_global_2020_1km.tif pop_global_2020_1km_opt.tifPREDICTOR3对浮点数据的压缩效果明显这是个小技巧很多人只用默认的PREDICTOR1压缩率差不少。3.3 色带断点没设对会把正常数据看成异常数据这一条属于“看着不像坑的坑”。把 1km 人口栅格拖进 GIS 软件默认色带通常按最小最大值线性拉伸。人口分布的极值差异巨大——城市核心区像元可能是几万人乡村像元可能是个位数线性拉伸的结果就是全国一片深色只有几个点在发亮看起来像是数据坏了。正确的做法是用分位数断点或者对数断点。我常用的是 2%、10%、25%、50%、75%、90%、98% 这七档能同时看清稀疏区和密集区。也可以直接对值取对数再渲染。这个操作本身不影响数据但会影响你对数据的判断。我见过有人因为默认渲染太难看误以为下载的文件损坏重新下载了三遍。提示分层设色之前先算一遍分位数值域里如果出现明显不合理的极大值比如单个像元十几万人先查是不是单位错了或者混入了非人口图层。4. 按行政边界裁剪那一刀切在哪里误差就从哪里来4.1all_touched这个参数能改变百分之几的总量用rasterio.mask或gdalwarp -cutline裁剪时边界像元的处理有两种策略只保留中心点落在多边形内的像元或者只要像元与多边形有任何接触就保留。前者是all_touchedFalse后者是all_touchedTrue。这个参数对小区域的影响非常显著。一个形状不规则的地级市边界周长可能有好几百公里涉及的边界像元数以千计。用all_touchedTrue会把边界外侧的像元也算进来总量偏高用False会把边界内侧只覆盖一部分的像元丢掉总量偏低。对于面积只有几千平方公里的市这个偏差能到 3%~8%。我的建议是做总量统计用False做可视化或面积覆盖分析用True并且在整个项目里保持统一。如果追求更精确可以对边界像元做面积加权也就是按像元与多边形的重叠比例分配人口但这需要额外的几何运算一般项目里用不上。import geopandas as gpd from rasterio.mask import mask import rasterio bound gpd.read_file(admin_boundary.shp) # 边界坐标系必须与栅格一致 with rasterio.open(pop_2020_1km.tif) as src: if bound.crs ! src.crs: bound bound.to_crs(src.crs) geom bound.loc[bound[name] A市, geometry].values out_img, out_transform mask( src, geom, cropTrue, nodatasrc.nodata, all_touchedFalse ) meta src.meta.copy() meta.update({ height: out_img.shape[1], width: out_img.shape[2], transform: out_transform }) with rasterio.open(pop_2020_A市.tif, w, **meta) as dst: dst.write(out_img)命令行版本更省事批量处理时我基本都用这个gdalwarp -cutline admin_boundary.shp \ -cwhere nameA市 \ -crop_to_cutline \ -dstnodata -99999 \ -co COMPRESSDEFLATE -co TILEDYES \ pop_2020_1km.tif pop_2020_A市.tif注意-crop_to_cutline会顺带把输出范围收紧到边界外接矩形文件小很多但输出栅格的原点会变后续做多图层叠加时要重新对齐。4.2 边界版本不一致是两个数据源对不上的头号原因分省分市数据最常见的用途就是把栅格统计值跟统计年鉴对比。对不上的时候先别怀疑栅格先怀疑边界。行政区划每年都在调整撤县设区、合并乡镇、跨市划转。2000 年的“A市”和 2023 年的“A市”地理范围可能差了百分之几。如果你用 2023 年的边界去裁 2000 年的栅格再用 2000 年的年鉴数字去校验得到的差异里混进了区划调整的成分根本说不清是数据误差还是口径差。我通常的做法是准备两套边界一套是“最新边界”用于所有年份的裁剪保证空间范围一致另一套是“逐年边界”只在需要严格对应官方统计时使用。两套边界的统计结果差值本身就是一个有价值的指标能反映出区划调整带来的影响量级。边界方案适用场景优点代价固定最新边界长序列趋势分析空间口径统一可比性强与早年官方统计有系统偏差逐年边界单年校验、年鉴对比与官方口径一致跨年不可比需额外对齐自定义研究区流域、都市圈分析贴合研究问题需自行构建工作量大4.3 分幅文件之间的缝隙和重叠全球或全国产品有时会按经纬度分幅提供比如每 5 度 × 5 度一个文件。批量处理时有两个细节要盯住一是分幅之间可能有 1 个像元的重叠直接拼接会导致边界处重复计数二是分幅的 nodata 填充方式可能不一致有的用 -99999有的用 0。拼接前先检查相邻文件的bounds是否严格相接。如果有重叠用gdalwarp的重采样会自然处理但rasterio.merge默认取后一个文件的值不会做加权。做总量统计时我倾向于先把所有分幅读成掩膜数组用全局坐标判断每个像元归属避免重复。5. 24 层栅格叠成时间序列年际波动里哪些是真信号5.1 先把“总量曲线”画出来这是最便宜的质量检查24 个年份的文件都到位之后第一件事不是做空间分析是把逐年总量算出来画一条线。这条线能暴露很多问题某年突然翻倍可能是单位错了某年突然腰斩可能是 nodata 处理不一致曲线整体过于平滑可能是插值年份占了大头。import glob import os import re files sorted(glob.glob(pop_*_1km.tif)) years, totals [], [] for f in files: year int(re.search(r(\d{4}), os.path.basename(f)).group(1)) total, _ block_sum(f) years.append(year) totals.append(total) for y, t in zip(years, totals): print(y, round(t / 1e8, 3), 亿)正常的人口总量曲线年际变化率一般在千分之几到百分之一之间。如果某一年相对上一年的变化超过 3%就需要单独查原因。常见的三类原因产品版本切换、边界文件更换、单位或 nodata 处理不一致。5.2 把“总量”和“结构”分开看总量曲线平稳不代表数据没问题。人口在空间上一直在重组中心城区可能流失外围新区在增长两者相互抵消总量看起来没变。所以我一般会同时算两个指标总量以及“高密度像元占比”比如每平方公里超过 5000 人的像元比例。这两个指标的变化方向经常不一致。我看到过一个案例某都市圈总量十年只涨了 4%但高密度像元占比从 12% 降到了 7%说明人口在向外围摊开。如果只看总量这个结构性变化完全被掩盖了。做结构分析时阈值的选择要跟研究问题匹配。做城乡划分可以用 1500 人/平方公里作为分界做城市蔓延可以用 5000 人/平方公里。阈值一旦定了整个序列都要用同一个不能因为某年数据分布不同就调整。5.3 年际变化率的空间可视化比单个年份的图有用得多单看某一年的分布图信息量其实有限因为人口分布形态十几年变化不大。真正有信息量的是“变化图”把 2023 减去 2000或者算年均变化率用发散色带渲染红色表示流失蓝色表示增长。做差值之前有两个预处理必须做一是把两年数据对齐到同一个网格分辨率和原点完全一致二是把 nodata 区域统一掩掉。如果两年数据的投影或范围有细微差异直接相减会产生大量假的“边缘变化”。import numpy as np import rasterio with rasterio.open(pop_2000_1km.tif) as s0, \ rasterio.open(pop_2023_1km.tif) as s1: a0 s0.read(1, maskedTrue) a1 s1.read(1, maskedTrue) # 两年都有效的像元才参与计算 valid (~a0.mask) (~a1.mask) diff np.where(valid, a1.filled(0) - a0.filled(0), np.nan) growth np.where(valid (a0.filled(0) 0), diff / a0.filled(0), np.nan)计算增长率时要特别小心分母为 0 的像元。荒漠区从 0 人变成 3 人增长率是无穷大这种像元会污染整个色带。我的做法是设一个最小基数阈值比如基数低于 10 人的像元不参与增长率计算只参与绝对差值计算。6. 和官方数字对不上时的排查链路6.1 从大到小逐级定位别一上来就抠像元统计结果和官方口径有差异是常态1km 产品的省级总量误差通常在 ±5% 以内市级在 ±10% 左右区县级可能到 ±20%。差异超出这个范围才需要认真排查。我习惯按“全球 → 大区 → 省 → 市 → 县”的顺序逐级缩小范围每一步都记录总量和差异率。如果差异在省级就出现了说明是全局性问题单位、nodata、边界版本如果省级对得上、市级对不上说明是局部性问题某个市的边界文件错了或者该市有大面积水体被误判。def compare(unit_name, raster_total, official_total): diff raster_total - official_total rate diff / official_total print(f{unit_name}: 栅格{raster_total:,.0f} 官方{official_total:,.0f} f差异{diff:,.0f} ({rate:.2%}))这个函数看着简单但逐级跑一遍问题基本就定位了。6.2 五个高频原因按出现概率排序排序原因典型表现定位方法1边界版本不一致全省对得上个别市偏差大换用同年边界重算2单位混淆差异在 100 倍或 1e6 量级检查是人数型还是密度型3nodata 被当 0总量系统性偏低统计掩膜比例检查哨兵值4裁剪参数不当小区域偏差明显大区域不明显对比all_touched两种结果5水体/冰川掩膜差异特定区域固定偏差叠加水体图层核对第一条我遇到过太多次。有个研究要做某市 2005—2020 的人口变化用的边界是 2020 年的结果 2005 年那几年的数字总是偏低。后来发现该市在 2012 年并入了一个邻近的县用 2020 年边界去裁 2005 年数据等于把当年还不属于该市的区域也算了进去而那个区域当时人口密度很低拉低了整体均值。第三条也很隐蔽。有些产品的海洋区域填的是 0nodata 值是 -99999但只在部分分幅文件里设置了nodata元数据另一些分幅没设。批量处理时如果不显式指定-99999就会被当成真实值加进去。全球数据里只要有几百万个这样的像元总量就会偏差几千万。6.3 校准什么时候该做什么时候不该做发现系统偏差之后很多人的第一反应是按比例缩放让总量对上官方数字。这个操作要谨慎。如果偏差来自明确的口径差异比如边界版本不同、水体掩膜不同按比例缩放在数学上能让数字好看但会把误差均匀摊到每个像元上掩盖真实的空间误差结构。比如某个市的实际误差集中在城区按比例缩放会让乡村像元也一起被调整反而更失真。我的原则是只在做跨区域对比时才校准且必须记录校准系数。如果研究目的是看单个区域的绝对量宁可保留原始值并在报告里说明误差范围也不要做无依据的缩放。校准系数本身也是有信息量的——如果连续多年系数稳定在 1.05 附近说明是系统性口径差如果系数逐年漂移说明产品版本在变这时候校准反而会引入假趋势。6.4 一个容易被忽略的校验维度小区域的表现大部分校验都在省级和市级做因为官方数字好找。但 1km 产品真正的软肋在小区域。我建议至少抽 5 到 10 个面积较小、地形有特点的区县做抽查比如山区县、水域占比高的县、快速城镇化的新区。抽查时会发现一些规律平原农业县的总量通常对得比较准因为人口分布均匀分配模型的误差互相抵消山区县容易偏高因为模型倾向于把人口分配到坡度较缓的像元而实际山谷里的聚落可能比模型判断的更集中水域占比高的县容易偏低因为水体掩膜把一些临水聚落也误删了。这些规律没法写进通用文档只能靠一个个区域试出来。我在做跨省对比时会先给每个省标记一个“地形类型”然后在同类型省份之间比较避免把地形导致的系统偏差误读成区域差异。7. 我自己的几条使用习惯用这类数据几年下来有几条习惯是一直保留的。第一任何一次统计之前先把单位、nodata、边界版本这三项写成一行注释贴在脚本开头强迫自己确认一遍。第二所有中间结果都保留原始未校准版本校准版本另存文件名带_cal后缀避免几个月后自己都分不清哪个是哪个。第三做长序列时一定会画两张图总量折线和高密度像元占比折线两条线放在一起看比任何单项指标都能说明问题。还有一个不算技巧但很省事的做法把 24 个年份的数据先统一重采样到一个固定的分析网格上而不是每次分析都从原始文件读。重采样会有一次精度损失但换来的是后续所有空间运算都可以直接做数组运算速度提升非常明显。重采样时用最近邻还是双线性取决于你的栅格值类型——人数型用最近邻更稳妥避免插值造出小数人口密度型用双线性更平滑。这个选择做错了后续所有统计都会带上一层看不见的偏差。