数据这行干久了最容易出问题的不是模型而是喂进去的那张底图。前阵子帮一个做灾害风险评估的团队看结果全国尺度的暴露人口算出来比普查数少了近两成查了两天最后发现是栅格单位搞错了——人家给的是人/平方公里的密度值他们直接按像元求和还漏了纬度方向上的面积变形。这类坑在1km精度人口分布栅格数据上特别常见因为这套数据看起来太标准了GeoTIFF、WGS84、一年一个文件谁拿到都觉得自己会用。实际上从2000到2023这一长串年份、从全球到分市这四级尺度里面藏着单位口径、总量校准、掩膜规则、像元对齐一堆细节任何一个没对齐后面的统计和建模全是错的。这篇就把我这些年用1km人口栅格做分析的整套流程摊开讲数据是怎么生产出来的、四种尺度版本的本质差别、下载后第一件事该做什么检查、裁剪重采样和分区统计怎么写代码不翻车以及那些只有真正跑过几十个年份数据才会遇到的坑。不管你是刚接触空间数据的学生还是要把人口栅格接进业务系统的开发者都能从中找到能直接抄的部分。1. 先把1km人口栅格这件事讲清楚1.1 它到底是什么从普查表格到连续表面的那一步转换人口普查给出的原始结果是一张以行政单元为行的表格某个县多少人、某个街道多少人。这种数据最大的问题是空间粒度被行政边界锁死了一个面积几千平方公里的县所有人口被压缩成一个数字你没法知道人是聚在县城还是散在山沟里。1km人口分布栅格数据要解决的正是这个问题——把行政单元上的总量拆解重分配到规则网格上让每个1km×1km的像元都有一个估计人口值从而得到一个空间上连续、边界上无缝的人口表面。这个拆解过程在学术上有个专门的名字叫分区密度制图dasymetric mapping。它的核心逻辑不复杂人口不是均匀铺在一个县的土地上的而是高度集中在建成区、居民点、道路沿线这些人真正待着的地方。所以只要找到能指示人待在哪的辅助信息就能把总量按权重摊下去。举个生活化的例子一个县有50万人县域总面积5000平方公里如果简单按面积平均每平方公里100人这显然是错的——城区每平方公里可能上万人山区可能每平方公里不到10人。分区密度制图做的事就是先识别出哪些格子是建成区、哪些是林地水域给前者高权重、后者零权重然后按权重把50万人分下去。这里有个特别容易被忽略的概念区分栅格值到底是人数还是密度。绝大多数全球尺度的产品输出的是每个像元内的估计人数count单位是人/像元但也有部分产品、部分衍生版本输出的是人口密度人/平方公里。这两种口径在做求和统计时完全不同——前者直接相加就是总人口后者必须先乘以像元实际面积再相加。而1km网格如果用的是等经纬度划分通常是30弧秒约0.008333度像元的实际面积是随纬度变化的赤道附近约0.93平方公里纬度40度左右降到约0.71平方公里纬度60度只剩约0.46平方公里。这意味着同一份密度数据在不同纬度求和误差可以差出一倍。所以拿到任何一份人口栅格第一件事永远是打开元数据确认单位和投影而不是急着算。1.2 四种空间尺度的版本本质差别比你想的大标题里列的全球/全国/分省/分市四级很多人以为是同一套东西切了四刀实际上它们的生产逻辑往往完全不同用之前必须搞清楚手里拿的是哪一类。全球尺度的产品通常是把各国普查数据、多源遥感特征、机器学习模型统一跑一遍覆盖完整、时间序列长优点是能跨国比较、口径相对统一缺点是国别之间的校准精度参差有些地区因为普查数据陈旧或缺失误差会明显偏大。国别尺度的产品一般是用本国更细的普查单元比如乡镇街道级作为控制总量辅助数据也更贴近本地实际精度通常优于全球产品在这一区域的切片。至于省级和市级版本情况就比较微妙了——相当一部分所谓分省分市数据并不是针对该省该市重新建模而是拿国别或全球产品按行政边界裁剪出来的再按该行政单元的普查总量做一次比例缩放。这种处理方式本身没问题但它意味着市级版本的空间分布格局其实还是原始产品给的只是总量被重新锚定了你在做市内精细分析时不能指望它比上游产品更准。我在实际项目里的判断原则是这样的做跨国或大区域对比用全球版做单一国家的深入分析优先找国别版只关心某个省或市的汇总数字用分省分市版省事但要研究城市内部的人口分布格局这四级都不够得去找更高分辨率的产品或者本地建筑体量、手机信令之类的数据补充。这个判断不是拍脑袋而是由数据生产的控制总量和辅助特征决定的——分辨率是表象控制单元和特征质量才是精度的真正来源。2. 数据是怎么造出来的生产链路拆解2.1 基准总量普查数据的年份错位与总量校准任何人口栅格的第一层地基都是普查或人口登记数据。问题在于普查不是每年都做多数国家是五年或十年一次而产品要输出2000到2023的连续年度序列中间年份怎么办常见做法是在两个普查年之间做线性或非线性插值再叠加人口登记、抽样调查等辅助序列做修正。这就带来一个必须知道的特性所谓2023年数据在普查间隔期内其实是插值加模型外推的产物它的可靠性明显低于普查当年。更关键的是总量校准环节。模型输出的是相对权重权重归一化后乘以控制总量才能保证格子里的人口加起来等于普查总数。这个校准一般在哪个层级做直接决定了数据的可信边界。如果校准在省级做那么省级汇总数一定是准的市级就可能偏如果校准在市级做市级准区县就可能偏。有些产品为了兼顾采用分层校准先在国家层面对齐再逐级下推工程量大但一致性更好。你在用数据做分析时如果结论依赖某个层级的汇总最好先验证一下该层级的汇总值和官方统计差多少差得离谱就说明校准层级不匹配。还有个隐蔽的坑是行政边界变更。2000到2023这二十多年区县一级的撤并、改区、新设非常频繁一个地级市在2010年和2020年的管辖范围可能都不一样。如果你拿2023年的行政边界去统计2005年的人口遇到边界调整过的区域数字必然对不上而且这种偏差不会报错只会安静地错在那里。我的处理习惯是涉及跨年比较时统一用一套固定年份的边界并且对发生过调整的单元单独标注避免把边界变化误读成人口变化。2.2 空间权重辅助变量怎么选、怎么处理权重层决定了人口在空间上摊得准不准这是整个生产链路里最见功力的部分。目前主流做法会堆叠几类特征。第一类是建成区与不透水面信息来自遥感影像分类直接告诉模型哪里有建筑、哪里是裸地农田这是最核心的负向/正向判别特征。第二类是夜间灯光灯光强度和人口密度有较强相关性尤其在城市化程度高的区域但它对农村人口的指示能力弱容易出现农村被严重低估、城市被高估的系统性偏差所以现在多数产品只把它当辅助而非主力。第三类是土地利用分类把林地、水体、耕地、建设用地的像元做差异化处理水域和保护区一般直接给零权重。第四类是道路网络和POI密度用来刻画可达性和人类活动强度POI密集的地方通常人口也密集。第五类是建筑体量与高度信息来自开放建筑数据集或雷达数据这是近些年提升精度最明显的一类特征因为它能区分一片低矮平房和一片高层住宅这两者的人口容量差好几倍。这些特征进模型前都要做处理灯光和POI密度这类长尾分布变量通常取对数不同量纲的特征要归一化分类变量要做独热编码或者直接作为分层条件。还有一个容易被忽略的步骤是掩膜——水系、冰川、自然保护区、军事管制区这些理论上不该有常住人口的区域要在权重层就置零否则模型会把人口漏到湖面上。我见过最典型的问题就是某个湖泊中心出现了一个非零人口像元一查是水体掩膜用的年份和人口数据年份不一致那年湖区水位低露出了一个岛。2.3 模型与重分配从随机森林到分区密度制图特征准备好之后就是建模建立特征—人口的关系。目前主流是随机森林一类的集成树模型原因很实际人口与辅助变量的关系高度非线性且在不同区域差异很大树模型不需要预设函数形式对共线性也不敏感训练和调参成本低。训练样本一般来自有细粒度普查数据的区域模型学到映射关系后再推广到只有粗粒度数据的区域。模型输出的是每个像元的相对权重接下来的重分配环节才是真正决定总量守恒的地方。标准流程是按掩膜筛掉不可居住像元剩余像元权重归一化乘以该行政单元的控制总人口得到每个像元的估计人数。这个过程在数学上保证了格子之和等于控制总量但也带来一个副作用——如果权重估计有系统偏差误差会在单元内部被强行摊平表现为某些像元被高估、某些被低估但总量看起来完美。这就是为什么不能只看总量校验通过就放心空间格局的合理性得单独看。还有一个技术选择值得说明是先做高分辨率再聚合到1km还是直接建模1km。前者通常先做出100m甚至更细的权重层再聚合求和好处是保留了城市内部的异质性聚合后的1km值更接近真实的空间集中程度后者直接以1km为单元训练计算量小但对城市核心区的刻画会偏平滑。从数据表现上能看出差别前者在城市中心往往有更尖锐的峰值后者更平缓。做城市内部研究时这个差别会影响结论值得在元数据里确认。2.4 2000到2023的时序一致性为什么相邻两年会跳时间序列数据最怕的不是单年误差而是年际之间的不一致。如果每一年都独立训练一次模型、独立做一次校准那么即使真实人口几乎没变相邻两年也会因为特征更新、样本变动、模型随机性产生肉眼可见的跳变。这种伪变化在做趋势分析时会变成噪声甚至被误读为人口迁移信号。主流产品对这个问题的处理有几套思路。一是固定掩膜和固定特征集只让控制总量随时间变化这样年际差异主要反映总量变化空间格局相对稳定。二是在基准年锚定其他年份做增量调整保证序列平滑。三是做时序平滑后处理用移动窗口抑制异常波动。每种都有代价固定掩膜会让新增建成区无法体现锚定基准年会让后期年份的空间格局偏保守平滑处理会削弱真实的突变信号比如新区突然建成入住。我在做长序列分析时的经验是先算一遍全国或全省的总量年序列看曲线是否平滑合理如果出现某一年突然掉一截又弹回来基本可以判定是数据版本或方法切换而不是真实人口变化。然后再挑几个已知发生过快速城镇化的区域看栅格的空间格局是否体现了建成区扩张如果二十多年过去格局几乎没动说明这份数据的时间动态信息很弱只能当静态底图用不能做变化检测。这个前置检查花不了半小时但能避免后面几个月的分析白做。3. 实操把数据用起来的完整流程3.1 下载后的第一件事文件体检与目录组织数据下载下来通常是一堆GeoTIFF按年份分文件全球版单个文件可能好几个GB。我建议在写任何分析代码之前先做一轮体检用几十行代码把关键元数据全部打印出来形成一张清单。import rasterio import glob import os files sorted(glob.glob(./pop/*.tif)) for f in files[:3]: with rasterio.open(f) as src: print(os.path.basename(f)) print( CRS :, src.crs) print( Size :, src.width, x, src.height) print( Res :, src.res) print( Bounds :, src.bounds) print( NoData :, src.nodata) print( Dtype :, src.dtypes[0]) print( Min/Max :, src.read(1, maskedTrue).min(), src.read(1, maskedTrue).max())这段代码要重点看四样东西。第一是CRS确认是地理坐标经纬度还是投影坐标后者要知道具体投影才能算面积。第二是Res确认像元大小是不是你以为的1km——有些产品名义上1km实际是30弧秒换算到地面约0.93km边界框对不上会导致后续裁剪错位。第三是NoData这个值五花八门有-99999、-32768、255、0如果NoData被当成真实值参与求和结果会离谱到离谱必须先确认再统一处理。第四是Dtype和Min/Max如果最大值出现几万甚至几十万那基本是密度值而不是人数或者NoData没被识别。目录组织上我习惯按产品名/版本/尺度/年份四级建目录文件名里带上分辨率、单位和生成日期。听起来啰嗦但当你手上同时有五六个来源的人口数据、每个来源好几个版本时这套命名能救命。另外强烈建议对原始文件做只读备份所有裁剪、重采样、格式转换的操作都输出到新目录永远不要覆盖原始数据——空间数据的坑一大半来自我以为我改的是副本。3.2 坐标系、NoData与像元对齐三个必须处理的底层问题坐标系问题在面积统计里最致命。经纬度坐标下的度不是等面积单位直接算面积会随纬度严重失真。正确做法是先把栅格重投影到一个等积投影比如按研究区中央经线定制的Albers等积投影或者全球范围的等积圆柱投影再做面积相关的计算。如果你的分析只是按行政单元汇总人数且数据单位是人/像元那不重投影也能做因为求和不需要面积但一旦涉及密度、人均面积、建筑密度这类指标重投影就是必须的。NoData处理有个容易翻车的细节很多栅格用0表示无数据但0也可能是合法的这里确实没人。这两种情况在语义上完全不同混淆会导致大片区被误判为无数据。我的做法是先看数据的值分布直方图如果0的像元数量异常大且集中在海洋、荒漠区域那多半是掩膜值如果0零星分布在内陆那可能是真实的零人口像元。确认之后统一用masked array处理让NoData在运算中被正确忽略而不是当0参与。像元对齐是三个问题里最隐蔽的。当你把人口栅格和行政边界、土地利用图、洪水淹没范围叠加时如果两张栅格的网格原点或像元大小差一点点就会产生半个像元的偏移导致边界上的统计出现系统性偏差。判断方法很简单打印两张栅格的transform参数逐项对比不一致就先做对齐。对齐时重采样方法的选择也有讲究——人口数据重采样到更粗分辨率用求和sum能保持总量重采样到更细分辨率用最近邻nearest能保持总量但会出现块状用双线性bilinear会平滑但破坏总量守恒。如果需要平滑又需要总量准确那就平滑之后再按行政单元做一次比例校准把总量拉回去。from rasterio.warp import calculate_default_transform, reproject, Resampling dst_crs EPSG:6933 # 等积投影示例 with rasterio.open(pop_2020_1km.tif) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs src.meta.copy() kwargs.update({crs: dst_crs, transform: transform, width: width, height: height}) with rasterio.open(pop_2020_1km_eqarea.tif, w, **kwargs) as dst: reproject(sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest)注意重投影到等积坐标系后像元不再严格是1km×1km如果后续要按每平方公里报数务必用重投影后的实际像元面积不要沿用1.0这个常数。3.3 裁剪与分区统计把栅格变成能进表格的数字最常见的需求是给我每个行政区的人口数。这件事看起来简单实际有几个决策点。第一是像元归属规则一个像元跨在边界上算给谁常规做法是按像元中心点落在哪个单元就归给谁这是rasterstats等工具的默认行为速度快但边界单元会有偏差。更严格的做法是按像元与多边形的重叠面积做加权分配计算量大但边界精度高适合面积小的城市单元。第二是球形 vs 平面的问题地理坐标下做面积加权会引入误差最好在投影坐标下做。下面是一个不依赖rasterstats、自己写的分区统计函数好处是能清楚控制NoData和像元中心规则也方便加面积加权import numpy as np import geopandas as gpd import rasterio from rasterio.features import rasterize def zonal_sum(raster_path, zones, zone_idadcode, all_touchedFalse): with rasterio.open(raster_path) as src: arr src.read(1, maskedTrue).astype(float64) transform src.transform nodata src.nodata # 用行政区 id 栅格化得到每个像元归属 shapes ((geom, int(v)) for geom, v in zip(zones.geometry, zones[zone_id])) zone_raster rasterize(shapes, out_shapearr.shape, transformtransform, fill0, all_touchedall_touched, dtypeint32) result {} for zid in zones[zone_id].unique(): sel (zone_raster int(zid)) vals arr[sel] result[zid] float(vals.sum()) if vals.size else 0.0 return result这段代码有两个地方要提醒。一是all_touched参数设为True会把所有接触边界的像元都算进来相邻单元会重复计数导致总量膨胀设为False只算中心点落入的像元总量守恒但小面积单元可能一个像元都没覆盖到结果变成0。对于面积小于1平方公里的城市街道这个问题很现实解决办法是提高栅格分辨率或者改用面积加权。二是arr.sum()在masked array上会自动忽略NoData但如果你提前用了.filled(0)NoData就变成了0混进去了这两种写法结果不同务必确认自己用的是哪一种。如果统计单元数量很多比如全国几千个乡镇上面这种逐单元循环会慢得难以接受。优化思路是把zone_raster拉平用pandas的groupby做聚合速度能提升一到两个数量级import pandas as pd flat_zone zone_raster.ravel() flat_val arr.filled(np.nan).ravel() df pd.DataFrame({zone: flat_zone, val: flat_val}) df df[df[zone] 0] summary df.groupby(zone)[val].sum()3.4 时间序列拼接把23年数据变成一条可分析的曲线单年数据处理完接下来是把2000到2023拼起来。这里推荐用xarray配合rioxarray它能把时间维度当成一个坐标轴做时序运算非常顺手import rioxarray import xarray as xr import glob files sorted(glob.glob(./pop_clip/pop_*_clip.tif)) das [] for f in files: year int(f.split(_)[-2]) da rioxarray.open_rasterio(f, maskedTrue).squeeze() da da.assign_coords(yearyear) das.append(da) cube xr.concat(das, dimyear) # 全区总人口年序列 total cube.sum(dim[x, y], skipnaTrue) # 逐年增长率 growth total.diff(year) / total.shift(year1)拼接时最常遇到的问题是对齐。如果每年的栅格范围或网格有细微差别xr.concat会报错或者按并集扩展出一堆NaN。解决办法是拼接前统一裁到一个固定的边界框并用同一个模板栅格做reproject_match保证所有年份的transform完全一致。还有一点如果某些年份缺文件xr.concat不会自动补得手动构造完整的年份索引再reindex否则算增长率时会把不相邻的两年当成相邻年来算得出的增长率完全没有意义。做趋势分析时我还习惯额外算两个东西一是分位数曲线看人口分布是越来越集中还是越来越分散这比总量更能反映城镇化过程二是高值像元的重心迁移轨迹用像元的加权平均坐标逐年计算能直观看到人口重心往哪个方向移动。这两个指标对判断数据质量也有帮助——如果重心一年跳几百公里那基本可以确定是数据问题而不是真实迁移。4. 几种典型应用场景的落地打法4.1 灾害暴露评估把人口叠到风险范围上这是人口栅格最经典也最刚需的用法。洪水、台风、地震、山火任何灾害风险评估都绕不开一个问题有多少人处在危险范围内。流程本身不复杂拿到灾害范围淹没深度栅格、影响半径缓冲区、历史灾害多边形和人口栅格做空间叠加统计范围内的总人口。但实操中有几个点决定了结果可不可信。第一是时间匹配。灾害发生在2020年你用的是2023年的人口数据那暴露人口就偏了。正确做法是用灾害发生当年或最近年份的数据如果只有普查年数据做一次线性内插到灾害年份。第二是空间范围匹配洪水淹没范围通常只覆盖河道及两岸而人口栅格是连续表面叠加时要注意淹没范围的边缘效应——范围边界上的像元人口可能只有一部分在淹没区内直接整像元计入会高估。第三是暴露的定义是范围内常住人口还是范围内夜间人口还是范围内需要疏散人口这三个数差别很大得在报告里说清楚口径。实际项目中我一般会输出三档结果范围全覆盖的乐观估计、按像元中心点落入的中间估计、按重叠面积加权的保守估计三个数一起给让决策方知道不确定性区间在哪。这比给一个精确到个位的数字实际得多也更诚实。4.2 公共服务设施可达性人口加权才是关键医疗、教育、养老设施的布局分析如果只用设施数量和行政区面积结论往往和真实感受差很远因为人口分布极度不均。把1km人口栅格接进来就能算出每个设施实际服务了多少人多少人能在30分钟内到达这类分析在城市规划和公共政策评估里用得非常多。具体做法上两步半法2SFCA是最常用的框架第一步对每个设施算供需比用设施容量除以服务范围内的人口第二步对每个需求点通常是栅格像元或小区质心算可达性把它能触达的所有设施的供需比按距离衰减加权求和。这里面人口栅格的作用是提供分母也就是服务范围内的人口基数。操作上的关键是把人口栅格聚合成需求点同时保留位置坐标聚合单元的大小要和出行距离尺度匹配——分析15分钟步行可达性用1km像元就太粗了最好细到250m或直接用建筑单元分析县域医疗资源覆盖1km就够用。还有一个容易出错的地方是距离的计算方式。直线距离欧氏距离会严重高估可达性尤其在有大河、山体、高速阻隔的区域。如果没有路网数据至少要用成本距离或者加入地形阻力如果有路网就用网络分析算实际通行时间。这个选择对结果的量级影响很大我在一次县域医疗可达性分析里对比过直线距离下90%的人口能在15分钟内到达换成路网实际通行时间后降到62%差距足够改变政策建议的方向。4.3 商业选址与市场容量估算商业场景里人口栅格的用法更直接以候选店址为圆心画1公里、3公里、5公里三个圈分别统计圈内常住人口、人口密度、家庭数估算值作为市场容量的第一层筛选。这套方法的好处是快能在几小时内把几百个候选点位筛到几十个坏处是它只看了有多少人住在这没看这些人会不会来。我一般会在人口栅格之外再叠三层修正。第一层是竞争密度统计圈内同类门店数量把市场容量按竞争强度打折。第二层是通行阻隔用路网判断圈内的实际可达范围一条高架或者一条河能把3公里圈切成两个世界。第三层是人流方向人口栅格是居住地分布但消费行为发生在通勤路上和工作地如果有白天人口数据或者手机信令的驻留数据应该和居住人口做交叉把住在这里和白天在这里区分开。很多选址失败案例的根因就是把居住人口当成了到店客流。还有一个技术细节值得提商业分析里常用1公里内人口这样的表述但到底按圆形缓冲区还是按步行等时圈结果差别巨大。步行等时圈更贴近实际消费行为现在有开源的路网数据和等时圈计算工具成本不高建议直接用等时圈替代圆形缓冲区说服力和准确度都更好。4.4 环境暴露与健康研究人口加权浓度才是暴露量空气污染、高温、绿地覆盖这类环境暴露研究如果只算区域平均浓度会完全忽略人到底在哪。一个城市平均PM2.5是40微克/立方米但人口密集的城区可能是55人口稀疏的郊区可能是25用40去做健康效应评估会系统性低估实际暴露水平。正确做法是用人口栅格对环境浓度栅格做加权平均得到人口加权浓度这才是暴露量。计算方法很简单把两张栅格对齐到同一网格逐像元相乘再求和除以总人口。但前提是两张栅格的分辨率和网格必须严格一致否则必须先重采样对齐。环境数据往往是3km、5km甚至更粗的网格重采样到1km时用双线性会引入平滑但这是可接受的因为大气污染本身空间连续。反过来把人口重采样到粗网格时一定要注意用求和而不是平均否则人口总量会丢失。另外提醒一个伦理和方法上的问题健康研究往往需要分年龄、分性别的人口结构而人口栅格只给总量。这时候需要拿普查的人口结构比例按行政单元乘到栅格上得到一个近似的分组人口表面。这种近似会在结构差异大的区域比如大学城、工业园区产生明显误差做精细的健康效应估算时要在局限性里写清楚不能当成精确值用。5. 常见问题与排查实录5.1 数据类问题速查表下面这张表是我这些年积累下来的高频问题清单基本覆盖了八成以上的数据看起来不对的情况。遇到异常先对着表排查一遍能省下大量时间。现象最可能的原因处理办法求和结果和普查总数差很多单位是密度而非人数NoData被当0计入查元数据确认单位用masked array排除NoData高纬度区域人口明显偏少等经纬度网格未做面积校正转等积投影后重算或按cos(纬度)校正像元面积湖面、海面出现非零人口水体掩膜年份与数据年份不一致用同期水体数据重新掩膜或直接置零相邻两年某区域数值突变基准年切换、模型版本更替、边界调整查版本说明做时序平滑固定边界重算分幅数据拼接后有接缝网格原点或像元大小不一致用统一模板做reproject_match后拼接分区统计总量膨胀用了all_touchedTrue导致跨边界重复计数改为中心点规则或改用面积加权分配小面积单元统计结果为0单元面积小于一个像元中心点未落入提高分辨率或改用面积加权最大值异常巨大NoData值未被识别参与求和确认nodata定义重新mask重采样后总量变了用了双线性或平均而非求和人数数据重采样到粗网格必须用sum和其他栅格叠加错位半个像元transform不完全一致对齐前逐项比对transform统一后运算表格里最需要强调的是第一行和最后一行。单位问题出现的频率远超想象因为很多数据分享页面只写1km人口数据不写单位下载者默认按人数用。像元错位则是最难发现的因为软件不会报错图看起来也正常只有当你把两张图放大到像元级别对比才会发现偏移。5.2 几个我实际踩过的坑第一个坑是关于分市数据的。有次做某个地级市的人口空间分析用了标称分市的数据结果发现市辖区内部的人口分布极其平滑完全没有体现出老城区和新区的差异。后来对比了国别版才发现这份分市数据就是把国别产品按市界裁了一刀再按市级普查总量缩放空间格局一点没变。而国别产品本身在那一带的建模分辨率就偏粗导致市区内部细节丢失。这个教训让我养成了一个习惯拿到任何分省分市数据先和一个上游产品做像元级对比如果除了比例系数之外完全一致那就知道它是裁剪版用它做精细分析要谨慎。第二个坑是关于零值的。某次统计一个西部县域的人口结果比普查数少了三成查了很久发现该县有大片区域在原始栅格里的值是0而我的代码把0当成了无数据给排除了。实际上那些0是真实的零人口像元荒漠和保护区排除它们不影响总和但问题在于我的代码同时还把NoData-99999也当成了有效值参与运算两个错误互相抵消了一部分导致结果既不对也看不出规律。后来我改成显式声明nodata并全程用masked array问题才彻底解决。这件事让我明白处理空间数据时隐式约定是最危险的所有假设都必须写进代码里。第三个坑是关于时间序列的。做某区域二十年人口变化分析曲线画出来发现2008年和2015年各有一个明显的台阶。一开始以为是行政区划调整查了边界没变又怀疑是普查数据更新查了普查年份也不对。最后翻了数据发布说明才发现这两个年份正好是产品大版本更新的年份模型和特征集都换了所以产生了方法性断点。这个发现直接改变了分析结论——原本准备写的人口阶段性变化特征其实是数据生产方法的痕迹。从那以后我做任何长序列分析都会先去查产品的版本更新日志把版本切换年份在图上标出来避免把方法变化误读成现实变化。5.3 精度这件事怎么说才不算误导最后想认真说一下精度表述。1km精度这个说法本身是有歧义的它可以指空间分辨率是1km也可以指精度误差是1km。实际上前者是确定的后者完全不是——1km人口栅格的像元值误差可能达到百分之几十在城市核心区和人口稀疏区还不一样。把1km分辨率说成1km精度在技术文档里是常见但不太严谨的表述。更合适的表达方式是分开说空间分辨率1km单个像元的相对误差在什么量级汇总到多大范围后误差能收敛到什么水平。经验上单个1km像元的人口估计误差可以很大但汇总到一个县、一个市之后正负误差会相互抵消总量误差通常能控制在个位数百分比。所以这份数据适合做区域级、城市级的分析不适合对单个像元做精确断言更不适合拿去做某小区有多少人这种精细判断。我在给别人交付结果时一定会附上这句话因为一旦使用者误以为1km栅格能精确到小区后面所有的决策都会建立在一个错误的精度假设上。还有一个常被忽略的不确定性来源是时间。前面提过普查间隔期内的年份是插值外推的越远离普查年不确定性越大。如果做趋势分析最好在图上把普查年标出来让读者知道哪些点是实测锚定的、哪些是模型推出来的。这个细节不增加任何计算量但能显著提升结论的可信度。我个人在实际操作中的体会是人口栅格这类数据真正的门槛不在会不会打开文件而在知不知道它的边界在哪——它由什么控制总量、用什么特征分配、在哪个层级校准、哪一年换了方法。把这四件事搞清楚一份普通的数据就能用得很稳搞不清楚再高分辨率的数据也只是好看的热力图。最后分享一个提高效率的小习惯给每一份用过的数据写一个不超过十行的元数据卡片记下来源、版本、单位、坐标系、NoData、校准层级、已知问题存在同目录下。一年之后再回来用这张卡片能帮你省掉整整一天的重新摸索。