简介一份覆盖1990—2015年的中国人口空间分布GeoTIFF栅格数据集空间分辨率1km每5年一期共6期数据。面向地理信息系统、城乡规划、人口地理与数据可视化领域的研究人员、工程师及高校师生可满足省级到全国尺度的人口制图与时空分析需求。压缩包共24个文件每个时期均含tif主影像、ovr金字塔、xml元数据与tfw坐标配准文件主影像用于栅格计算与出图辅助文件用于快速显示和地理定位整体约179.57MB。已有259人学习下载数据覆盖1990、1995、2000、2005、2010、2015年六个时间断面可直接用于人口密度制图、城市化进程研究与区域发展评估。配合ArcGIS、QGIS等主流软件可完成时空演变分析和三维可视化省去自行收集、裁剪与配准的步骤是论文研究、教学演示和专题图制作的便捷基础数据。1. 一份 tif 格式的 1km 人口数据不是图是能算的栅格中国人口空间分布数据 tif 格式 1km名字像一张图实际是一层栅格每个像元有一个数值代表这个公里格网里住了多少人。第一次拿它做项目时我也当普通图片打开结果一片黑、叠不上、对不齐最后才意识到问题出在拉伸、坐标系和 NoData 三个环节。这份资源的价值是把普查汇总的人口数字拆到 1km 网格上让你能按任意多边形切出人口、做覆盖分析、做选址评分。适合刚接触栅格数据的学生、要出业务图的分析师以及被要求把人口叠到地块上的开发。这篇笔记按实操顺序走看懂元数据裁剪重投影可视化出图分区统计落表踩过的坑都写在对应章节里。2. 读懂这份 1km 人口栅格格式、坐标系与像元含义拿到 tif 格式的人口文件第一件事不是打开看图而是弄清三件事这个 GeoTIFF 是怎么组织的、像元值是什么口径、坐标系是哪一套。三件事里任何一件弄错后面统计和可视化都会白做这一章把三个问题一次讲透。2.1 GeoTIFF 的组织方式TIF 不只是一张图普通 TIF 只有像素值和压缩方式GeoTIFF 在文件头里多写了一套地理标签记录坐标原点、像元尺寸、投影定义和 NoData 值。所以同样一张 tif 格式的图用看图软件打开可能正常用 GIS 软件打开却要求你指定坐标系——这就是地理标签缺失或读不出来的表现。人口栅格一般用 Float32 单波段存储因为人口是一个连续数值一个格网算出 3.7 人这种值很正常用整数存就会丢总量。我拿到文件后的第一个动作永远是跑一遍 gdalinfo不猜gdalinfo CHN_ppp_v2b_2020.tif输出大致长这样示例实际以你拿到的文件为准Driver: GTiff/GeoTIFF Size is 7440, 4320 Coordinate System is: GEOGCRS[WGS 84, DATUM[World Geodetic System 1984]] Origin (73.000000000000000,54.000000000000000) Pixel Size (0.008333333333333,-0.008333333333333) Data Type: Float32 NoData Value-99999逐行解释Size 是宽和高中国范围的 1km 网格切片大约在 7000×4000 量级Pixel Size 是 0.0083333° 的经纬度步长即 1/120 度这就是1km的来历NoData Value 是无效值标记常见的有 -99999、0、65535不同来源不一样必须先确认。如果 NoData 没设对后面所有统计都会把无效区当真实数值算进去。在 Python 里对应读元数据的写法是这样import rasterio with rasterio.open(CHN_ppp_v2b_2020.tif) as src: profile src.profile # 波段数、类型、压缩、NoData 等完整参数 print(CRS:, src.crs) print(分辨率(度):, src.res) print(尺寸:, src.width, x, src.height) print(NoData:, src.nodata) print(边界:, src.bounds)profile返回一个字典包含驱动、数据类型、NoData、压缩方式等全部写入参数src.crs是投影对象src.res是像元尺寸src.bounds是地理范围。这一步花不了一分钟但能省掉后面数小时的返工。2.2 像元值含义人口数还是人口密度一字之差差出几十倍这一步最容易翻车。人口格网数据有两类口径先看对比表口径像元值含义求和结果常见来源每像元人数这个格网内约住多少人区域加总即人口总数WorldPop ppp、LandScan人口密度每平方公里多少人需乘像元面积才是人口国内机构发布的 1km 密度格网判断方法有两种。第一看文档里写的是 persons per pixel 还是 persons per km²。第二实测把全省范围的像元加起来和七普常住人口对比。加总结果是千万到亿级那就是每像元人数结果只有几百万且像元均值明显偏大那很可能是密度。如果是密度还要看投影投影后的 1000m 格网每个像元面积接近 1km²密度和人数数值上几乎相等但经纬度坐标下的 0.0083333° 格网像元实际面积随纬度变化东北地区一个格网约 0.7km²直接拿密度当人数会明显偏低。还有一个容易被忽略的点很多 1km 人口数据不是整数。人口分布在网格上是平滑的一个格网里有 3.7 人很正常。读取时不要自作主张四舍五入用 Float32 原值直接计算四舍五入到整数全省加总可能损失 2%3%对比年鉴时会造成错觉。2.3 坐标系先确认WGS84、CGCS2000 还是 Albers人口栅格常见坐标系有三类WGS84 经纬度EPSG:4326、CGCS2000 经纬度EPSG:4490、投影坐标如 Krasovsky 1940 Albers 或 CGCS2000 Albers 等积投影。世界公开数据多用 WGS84国内科研机构发布的 1km 人口格网很多用 Albers 等积投影中央经线 105°E、双标准纬线 25°N 和 47°N。用等积投影是因为人口统计要算面积等积投影下每个像元面积恒定统计才准确。为什么必须在裁剪前确认坐标系如果你的省界矢量是 CGCS2000栅格是 WGS84两者经纬度数值相差很小直接叠加肉眼几乎看不出来但裁剪时切出来的范围会整体偏移人口统计也随之偏差。更麻烦的是自定义 Albers 投影没有标准 EPSG 编码src.crs可能显示一长串 PROJCS 字符串这时最稳妥的做法是把矢量边界转换到栅格坐标系而不是反过来重投影整幅大栅格import geopandas as gpd prov gpd.read_file(province_boundary.geojson) print(矢量 CRS:, prov.crs) with rasterio.open(CHN_ppp_v2b_2020.tif) as src: print(栅格 CRS:, src.crs)两行代码就能对比。不一致时用prov.to_crs(src.crs)把边界转到栅格坐标系再做裁剪避免了对整个中国范围的大文件重投影速度和安全性都更好。3. 把人口 TIF 变成可用成果读取、裁剪、重投影一气呵成先读出数组、再裁剪、再重投影这三步是任何人口格网分析的地基。每一步都会遇到看着对、细算错的问题所以每一步的校验方法我一起写出来。3.1 用 rasterio 读取把第一波段读成二维数组读取是整个流程里最没技术含量、但最容易出错的一步错在 NoData 处理和数据类型上import rasterio import numpy as np with rasterio.open(CHN_ppp_v2b_2020.tif) as src: band src.read(1) # 读第 1 波段得到 (height, width) 的二维 ndarray nodata src.nodata # 人口数据不可能为负把负值统一处理掉 band_clean np.where(band 0, 0, band) # 粗校验全国加总应接近当次普查人口 total band_clean.sum() print(全国格网加总:, int(total))src.read(1)返回 Float32 数组内存可以估算宽 7440 × 高 4320 × 4 字节约 128MB普通机器完全放得下如果是更大范围的分块文件要分窗口读。把负值改成 0 而不是删掉是为了后面按多边形统计时数组形状不变量无人区像元值为 0不参与求和即可。粗校验加总如果和普查人口差一个数量级停在这一步回去查 2.2 节的口径问题不要继续往下做。对大文件还有一个省内存的写法用src.read(1, out...)复用缓冲区或者用src.windows()分窗口读取。对 1km 全国数据通常用不上但你要叠加多份变量一起算的时候这个方法能救急。3.2 按行政区裁剪先对齐坐标系再下刀裁剪最常见的翻车点是坐标系没对齐。我一般把裁剪写成固定流程每步都校验import geopandas as gpd import rasterio from rasterio.mask import mask city gpd.read_file(city_boundary.geojson) with rasterio.open(CHN_ppp_v2b_2020.tif) as src: # 第一步坐标系不一致就转换矢量 if city.crs ! src.crs: city city.to_crs(src.crs) # 第二步裁剪多边形外部填 nodata0 out_img, out_transform mask( src, city.geometry, cropTrue, # 输出贴合多边形范围 nodata0, # 外部区域标记为 0 filledTrue # 把外部像素填成 nodata避免混入真实值 ) # 第三步更新元数据并写盘 meta src.meta.copy() meta.update(driverGTiff, heightout_img.shape[1], widthout_img.shape[2], transformout_transform, nodata0) with rasterio.open(city_pop_1km.tif, w, **meta) as dst: dst.write(out_img)cropTrue决定输出范围是否贴合多边形filledTrue保证多边形外部的像素全部填成 0不让周边真实人口值混进裁剪文件这是统计正确的关键。nodata0给输出文件打上外部无数据标记QGIS 里打开外部区域不会被当成真实 0 值上色。裁剪完之后做一次验证把裁剪结果加总理论上等于整幅数据落在这个市范围内的那部分。如果总和明显偏少但边界又没错多半是边缘像元归属问题回到第 5 章 5.4 的解法处理。3.3 重投影与重采样人口总量要守恒别用 nearest 硬转当你要把 WGS84 的人口栅格转成 CGCS2000或把 1km 转成 5km 时重采样方法选错人口总量会直接崩掉。nearest 适合分类数据bilinear 适合连续表面而人口这种总量型数据降到更粗分辨率时必须用 sum 重采样保证新网格上叠加的源像元值加起来等于原来的总量import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path CHN_ppp_v2b_2020.tif dst_path pop_wgs84_5km.tif dst_crs EPSG:4326 with rasterio.open(src_path) as src: # 按源范围计算目标网格的变换参数 transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) profile src.profile.copy() profile.update(crsdst_crs, transformtransform, widthwidth, heightheight) with rasterio.open(dst_path, w, **profile) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.sum # 人口总量守恒 )calculate_default_transform根据源范围和目标坐标系算出合适的宽高与仿射变换避免目标网格过大造成冗余Resampling.sum是 rasterio 1.3 才有的选项适合人口、降水这类总量数据。版本不支持 sum 时退而求其次用 average再把结果乘上源像元面积 / 目标像元面积做换算但这只是近似不如 sum 干净。重采样方法适用数据类型对人口总量的影响nearest分类、土地利用重复或丢失格网总量失真大bilinear连续表面抹平峰值总量略变average降尺度通用总量须按面积比例换算sum人口等总量型守恒推荐升尺度转更细分辨率时没有完美方案任何插值都是想象出的分布总量保持但空间格局被平滑输出时要在说明里标注重采样方法和日期否则后手会把你平滑过的东西当原始数据用。4. 人口数据可视化分级设色、百分位拉伸与边界叠加数据处理完到出图环节。1km 人口数据可视化的难点不在工具而在人口分布极端不均匀绝大多数格网接近 0少数超大城市一个格网就好几万。默认的最小值-最大值拉伸会把色带全压在城市上全国图几乎全白所以出图的核心是拉伸策略。4.1 QGIS 流程三分钟出一张全国人口分布图我常用 QGIS 做快速出图操作路径固定拖入 TIF右键图层属性Symbology 选 Singleband pseudocolor色带选 YlOrRd 反转深色表示高值Mode 选 Quantile 分成 8 类Min/Max 不勾自动改用手动范围并按 2%98% 百分位裁剪。参数推荐值理由渲染类型Singleband pseudocolor连续数值按色带上色比灰度可读色带YlOrRd 反转低值浅黄、高值深红符合直觉ModeQuantile8 类等频分级避免大量格网挤在一个色级拉伸范围2%98% 百分位排除极小/极大值干扰Clip out of range勾选超范围像元显示为透明不被截断人口数据高度右偏等间距分级会让 80% 的格网落在第一级里图没法看分位数分级保证每级格网数量大致相同读图时看到的是相对稠密区而不是绝对数值。想进一步压缩动态范围可以对数据做 log1p 变换后再分级但出图说明里要写清楚变换方式。4.2 Python 出图百分位拉伸 色带的完整脚本不依赖 QGIS 的自动化出图用 matplotlib 三件事就能说清读数组、定拉伸范围、画图。完整脚本如下import matplotlib.pyplot as plt import numpy as np import rasterio with rasterio.open(CHN_ppp_v2b_2020.tif) as src: band src.read(1) # imshow 的 extent 顺序是 left, right, bottom, top extent [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top] crs src.crs # 无效值和非人口区域统一置为 nan让它们不上色 band np.where(band 0, np.nan, band) # 百分位拉伸避开少数大城市把色阶拉平 p2, p98 np.nanpercentile(band, [2, 98]) fig, ax plt.subplots(figsize(12, 9)) im ax.imshow( band, extentextent, cmapYlOrRd, vminp2, vmaxp98, interpolationnearest ) ax.set_title(中国人口空间分布1km 格网, fontsize14) cbar plt.colorbar(im, axax, fraction0.03, pad0.02) cbar.set_label(人 / 格网) # 若是密度数据请改成 人 / 平方公里 plt.savefig(china_pop_1km.png, dpi300, bbox_inchestight)extent必须按 left/right/bottom/top 顺序组装直接用src.bounds会把边界顺序弄混图像坐标就错了。interpolationnearest是关键matplotlib 默认插值会把人口栅格平滑成热力图误导人以为人口连续渐变。nanpercentile计算百分位时自动跳过 nanNoData 和负值不参与拉伸。vmin/vmax取 2% 和 98% 分位是为了出图对比度而不是为了数值准确——这跟统计是两回事。4.3 叠加边界与色标让图和业务对得上光有栅格没有行政界线业务没法用。接着上一段代码叠加省界import geopandas as gpd prov gpd.read_file(province_boundary.geojson) if prov.crs is not None and prov.crs ! crs: prov prov.to_crs(crs) prov.boundary.plot(axax, color#444444, linewidth0.4) plt.savefig(china_pop_1km_boundary.png, dpi300, bbox_inchestight)叠加前必须确认两侧坐标系一致否则在 1km 尺度上偏移几百米一般看不出但统计对不上。输出印刷级图建议直接在 QGIS 打印布局里做矢量边界和栅格一起导出图例、比例尺、指北针都齐备Python 脚本更适合批量出图或做多期对比。5. 常见问题与避坑五个人口栅格高频翻车点这几条是我和同事在实际项目里真实踩过、并且每条都赔过时间的坑。按现象 → 原因 → 解决的顺序写方便对号入座。5.1 现象TIF 拖进软件全黑或全白什么都看不见排第一的坑。1km 人口数据像元值从 0 到几万90% 的像元集中在低值区只有少数城市格网到几万。默认渲染用全幅最小最大值做线性拉伸色带全被极端高值占据剩下区域在色带里分不到一级看上去就是全黑或全白。原因不是数据坏是渲染策略不适合右偏分布。解决QGIS 里 Symbology 选 Singleband pseudocolorMode 选 QuantileMin/Max 手动改成 2%98% 百分位Python 出图用np.nanpercentile算vmin/vmax。想进一步拉开低值区差异可以做 log1p 变换后渲染但图例要标注变换方式。5.2 现象区域加总人口只有年鉴的一成甚至出现负数这个坑我印象最深。第一次用这份数据算一个地级市结果只有年鉴的 8%我一度以为是数据坏了。逐格排查才发现 NoData 是 -99999统计时没排除无效值几十万无效像元全被当成真实人口求和直接把结果拉崩。另一种常见原因是口径混用把密度数据当每像元人数直接加总或者把户籍口径数据拿去和常住人口年鉴对比差出 20% 是常态。解决分三步先 gdalinfo 确认 NoData 值和数据类型再全省加总和七普常住人口对比确定这份文件属于哪类口径最后统计时显式传nodata参数不要依赖函数默认值。5.3 现象裁剪结果和行政边界错位像元好像漂了半格矢量边界是 CGCS2000栅格是 WGS84两者经纬度数值只差零点几弧秒叠图在省级尺度几乎看不出放大到城市、区县尺度几百米偏移就非常明显。还有一类更隐蔽国内机构发布的人口栅格用自定义 Albers 投影没有标准 EPSG 编码读出来是一长串 PROJCS 字符串如果偷懒不处理偏移量会到公里级。解决每次处理前打印两侧 CRS 对比不一致就to_crs把矢量转到栅格坐标系不要凭文件名猜坐标系一切以元数据输出为准。5.4 现象山区或建成区边界出现一排空白像元裁剪用cropTrue时多边形边界是弧线像元是方块弧线切过像元时边缘像元归属哪边存在取舍如果边界精度高于栅格边界上只有小半落在多边形内的像元会被丢弃视觉上就是边界锯齿加一圈空白。解决统计环节用rasterstats时传all_touchedTrue把与多边形相交的所有像元纳入计算裁剪环节先给边界做缓冲0.01° 或一个像元尺寸裁完再用原始边界做可视化遮罩。对全省层面影响可忽略但对市区级别的小范围统计一排像元可能代表上万人。5.5 现象重投影后全国总人口变了 20% 以上WGS84 经纬度网格转成 CGCS2000 或 Albers 投影时目标网格和源网格像元位置不可能完全对齐重采样必须决定每个新像元从哪些源像元取值。nearest 取最近的一个会重复取某些像元、丢掉另一些人口总量瞬间失真bilinear 对四邻加权平均城市峰值被抹平形状对了但总量变。解决降尺度用Resampling.sum保证新网格加总等于原网格加总升尺度没有守恒方案必须标注插值来源涉及面积或密度计算先把栅格转入 Albers 等积投影再用真实像元面积算不要在经纬度坐标下用固定系数糊弄。注意重采样参数要在成果目录里留一份说明标明源分辨率、目标分辨率、重采样方法和日期防止后手把插值结果当实测数据。6. 进阶分区统计、总量校验与一张可用的表处理完栅格、出完图业务最终要的是这个省多少人、那个市多少人。分区统计和总量校验是收尾的关键两步用 rasterstats 一次跑完import geopandas as gpd from rasterstats import zonal_stats provinces gpd.read_file(province_boundary.geojson) stats zonal_stats( provinces, CHN_ppp_v2b_2020.tif, stats[sum, mean, max], nodata-99999, all_touchedTrue, geojson_outTrue ) df gpd.GeoDataFrame.from_features(stats) result df[[NAME, sum, mean, max]].sort_values(sum, ascendingFalse) result.to_csv(province_pop_1km.csv, indexFalse) print(result.head(10))stats[sum, mean, max]中 sum 才是人口总数mean 是平均每格网人数max 是最高密度格网nodata-99999必须和文件实际 NoData 一致否则负数参与求和all_touchedTrue把边界相交像元算进去避免区县级边界切漏人口。输出 GeoDataFrame 直接写 CSV这张表就能进 PPT 或对接业务系统。统计完必须做总量校验把各省 sum 加总和全国普查常住人口对比。世界公开数据有 UN adjusted 和 unadjusted 两个版本后者没按联合国分年龄性别数据校准低龄和老年人口偏差明显国内机构的格网还分常住、户籍、夜间灯光修正版。误差在 ±5% 以内可以放心用超过 20% 先回查 NoData 和坐标系再查口径。1km 格网在人口稀疏区本来就有空间平滑拿来比较城市间相对密度没问题但拿某个偏远县的绝对人口去做预算审批风险自负。从那以后我每次拿到人口 TIF 都强制走一遍固定三步gdalinfo 看元数据、全省加总对口径、边界 CRS 对齐截图确认然后才允许进入可视化和统计。这三分钟换来的是一次做对而不是反复返工。希望帮到你。本文还有配套的精品资源点击获取
1km人口栅格数据处理全指南：读取、裁剪、重投影与统计
简介一份覆盖1990—2015年的中国人口空间分布GeoTIFF栅格数据集空间分辨率1km每5年一期共6期数据。面向地理信息系统、城乡规划、人口地理与数据可视化领域的研究人员、工程师及高校师生可满足省级到全国尺度的人口制图与时空分析需求。压缩包共24个文件每个时期均含tif主影像、ovr金字塔、xml元数据与tfw坐标配准文件主影像用于栅格计算与出图辅助文件用于快速显示和地理定位整体约179.57MB。已有259人学习下载数据覆盖1990、1995、2000、2005、2010、2015年六个时间断面可直接用于人口密度制图、城市化进程研究与区域发展评估。配合ArcGIS、QGIS等主流软件可完成时空演变分析和三维可视化省去自行收集、裁剪与配准的步骤是论文研究、教学演示和专题图制作的便捷基础数据。1. 一份 tif 格式的 1km 人口数据不是图是能算的栅格中国人口空间分布数据 tif 格式 1km名字像一张图实际是一层栅格每个像元有一个数值代表这个公里格网里住了多少人。第一次拿它做项目时我也当普通图片打开结果一片黑、叠不上、对不齐最后才意识到问题出在拉伸、坐标系和 NoData 三个环节。这份资源的价值是把普查汇总的人口数字拆到 1km 网格上让你能按任意多边形切出人口、做覆盖分析、做选址评分。适合刚接触栅格数据的学生、要出业务图的分析师以及被要求把人口叠到地块上的开发。这篇笔记按实操顺序走看懂元数据裁剪重投影可视化出图分区统计落表踩过的坑都写在对应章节里。2. 读懂这份 1km 人口栅格格式、坐标系与像元含义拿到 tif 格式的人口文件第一件事不是打开看图而是弄清三件事这个 GeoTIFF 是怎么组织的、像元值是什么口径、坐标系是哪一套。三件事里任何一件弄错后面统计和可视化都会白做这一章把三个问题一次讲透。2.1 GeoTIFF 的组织方式TIF 不只是一张图普通 TIF 只有像素值和压缩方式GeoTIFF 在文件头里多写了一套地理标签记录坐标原点、像元尺寸、投影定义和 NoData 值。所以同样一张 tif 格式的图用看图软件打开可能正常用 GIS 软件打开却要求你指定坐标系——这就是地理标签缺失或读不出来的表现。人口栅格一般用 Float32 单波段存储因为人口是一个连续数值一个格网算出 3.7 人这种值很正常用整数存就会丢总量。我拿到文件后的第一个动作永远是跑一遍 gdalinfo不猜gdalinfo CHN_ppp_v2b_2020.tif输出大致长这样示例实际以你拿到的文件为准Driver: GTiff/GeoTIFF Size is 7440, 4320 Coordinate System is: GEOGCRS[WGS 84, DATUM[World Geodetic System 1984]] Origin (73.000000000000000,54.000000000000000) Pixel Size (0.008333333333333,-0.008333333333333) Data Type: Float32 NoData Value-99999逐行解释Size 是宽和高中国范围的 1km 网格切片大约在 7000×4000 量级Pixel Size 是 0.0083333° 的经纬度步长即 1/120 度这就是1km的来历NoData Value 是无效值标记常见的有 -99999、0、65535不同来源不一样必须先确认。如果 NoData 没设对后面所有统计都会把无效区当真实数值算进去。在 Python 里对应读元数据的写法是这样import rasterio with rasterio.open(CHN_ppp_v2b_2020.tif) as src: profile src.profile # 波段数、类型、压缩、NoData 等完整参数 print(CRS:, src.crs) print(分辨率(度):, src.res) print(尺寸:, src.width, x, src.height) print(NoData:, src.nodata) print(边界:, src.bounds)profile返回一个字典包含驱动、数据类型、NoData、压缩方式等全部写入参数src.crs是投影对象src.res是像元尺寸src.bounds是地理范围。这一步花不了一分钟但能省掉后面数小时的返工。2.2 像元值含义人口数还是人口密度一字之差差出几十倍这一步最容易翻车。人口格网数据有两类口径先看对比表口径像元值含义求和结果常见来源每像元人数这个格网内约住多少人区域加总即人口总数WorldPop ppp、LandScan人口密度每平方公里多少人需乘像元面积才是人口国内机构发布的 1km 密度格网判断方法有两种。第一看文档里写的是 persons per pixel 还是 persons per km²。第二实测把全省范围的像元加起来和七普常住人口对比。加总结果是千万到亿级那就是每像元人数结果只有几百万且像元均值明显偏大那很可能是密度。如果是密度还要看投影投影后的 1000m 格网每个像元面积接近 1km²密度和人数数值上几乎相等但经纬度坐标下的 0.0083333° 格网像元实际面积随纬度变化东北地区一个格网约 0.7km²直接拿密度当人数会明显偏低。还有一个容易被忽略的点很多 1km 人口数据不是整数。人口分布在网格上是平滑的一个格网里有 3.7 人很正常。读取时不要自作主张四舍五入用 Float32 原值直接计算四舍五入到整数全省加总可能损失 2%3%对比年鉴时会造成错觉。2.3 坐标系先确认WGS84、CGCS2000 还是 Albers人口栅格常见坐标系有三类WGS84 经纬度EPSG:4326、CGCS2000 经纬度EPSG:4490、投影坐标如 Krasovsky 1940 Albers 或 CGCS2000 Albers 等积投影。世界公开数据多用 WGS84国内科研机构发布的 1km 人口格网很多用 Albers 等积投影中央经线 105°E、双标准纬线 25°N 和 47°N。用等积投影是因为人口统计要算面积等积投影下每个像元面积恒定统计才准确。为什么必须在裁剪前确认坐标系如果你的省界矢量是 CGCS2000栅格是 WGS84两者经纬度数值相差很小直接叠加肉眼几乎看不出来但裁剪时切出来的范围会整体偏移人口统计也随之偏差。更麻烦的是自定义 Albers 投影没有标准 EPSG 编码src.crs可能显示一长串 PROJCS 字符串这时最稳妥的做法是把矢量边界转换到栅格坐标系而不是反过来重投影整幅大栅格import geopandas as gpd prov gpd.read_file(province_boundary.geojson) print(矢量 CRS:, prov.crs) with rasterio.open(CHN_ppp_v2b_2020.tif) as src: print(栅格 CRS:, src.crs)两行代码就能对比。不一致时用prov.to_crs(src.crs)把边界转到栅格坐标系再做裁剪避免了对整个中国范围的大文件重投影速度和安全性都更好。3. 把人口 TIF 变成可用成果读取、裁剪、重投影一气呵成先读出数组、再裁剪、再重投影这三步是任何人口格网分析的地基。每一步都会遇到看着对、细算错的问题所以每一步的校验方法我一起写出来。3.1 用 rasterio 读取把第一波段读成二维数组读取是整个流程里最没技术含量、但最容易出错的一步错在 NoData 处理和数据类型上import rasterio import numpy as np with rasterio.open(CHN_ppp_v2b_2020.tif) as src: band src.read(1) # 读第 1 波段得到 (height, width) 的二维 ndarray nodata src.nodata # 人口数据不可能为负把负值统一处理掉 band_clean np.where(band 0, 0, band) # 粗校验全国加总应接近当次普查人口 total band_clean.sum() print(全国格网加总:, int(total))src.read(1)返回 Float32 数组内存可以估算宽 7440 × 高 4320 × 4 字节约 128MB普通机器完全放得下如果是更大范围的分块文件要分窗口读。把负值改成 0 而不是删掉是为了后面按多边形统计时数组形状不变量无人区像元值为 0不参与求和即可。粗校验加总如果和普查人口差一个数量级停在这一步回去查 2.2 节的口径问题不要继续往下做。对大文件还有一个省内存的写法用src.read(1, out...)复用缓冲区或者用src.windows()分窗口读取。对 1km 全国数据通常用不上但你要叠加多份变量一起算的时候这个方法能救急。3.2 按行政区裁剪先对齐坐标系再下刀裁剪最常见的翻车点是坐标系没对齐。我一般把裁剪写成固定流程每步都校验import geopandas as gpd import rasterio from rasterio.mask import mask city gpd.read_file(city_boundary.geojson) with rasterio.open(CHN_ppp_v2b_2020.tif) as src: # 第一步坐标系不一致就转换矢量 if city.crs ! src.crs: city city.to_crs(src.crs) # 第二步裁剪多边形外部填 nodata0 out_img, out_transform mask( src, city.geometry, cropTrue, # 输出贴合多边形范围 nodata0, # 外部区域标记为 0 filledTrue # 把外部像素填成 nodata避免混入真实值 ) # 第三步更新元数据并写盘 meta src.meta.copy() meta.update(driverGTiff, heightout_img.shape[1], widthout_img.shape[2], transformout_transform, nodata0) with rasterio.open(city_pop_1km.tif, w, **meta) as dst: dst.write(out_img)cropTrue决定输出范围是否贴合多边形filledTrue保证多边形外部的像素全部填成 0不让周边真实人口值混进裁剪文件这是统计正确的关键。nodata0给输出文件打上外部无数据标记QGIS 里打开外部区域不会被当成真实 0 值上色。裁剪完之后做一次验证把裁剪结果加总理论上等于整幅数据落在这个市范围内的那部分。如果总和明显偏少但边界又没错多半是边缘像元归属问题回到第 5 章 5.4 的解法处理。3.3 重投影与重采样人口总量要守恒别用 nearest 硬转当你要把 WGS84 的人口栅格转成 CGCS2000或把 1km 转成 5km 时重采样方法选错人口总量会直接崩掉。nearest 适合分类数据bilinear 适合连续表面而人口这种总量型数据降到更粗分辨率时必须用 sum 重采样保证新网格上叠加的源像元值加起来等于原来的总量import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path CHN_ppp_v2b_2020.tif dst_path pop_wgs84_5km.tif dst_crs EPSG:4326 with rasterio.open(src_path) as src: # 按源范围计算目标网格的变换参数 transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) profile src.profile.copy() profile.update(crsdst_crs, transformtransform, widthwidth, heightheight) with rasterio.open(dst_path, w, **profile) as dst: reproject( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.sum # 人口总量守恒 )calculate_default_transform根据源范围和目标坐标系算出合适的宽高与仿射变换避免目标网格过大造成冗余Resampling.sum是 rasterio 1.3 才有的选项适合人口、降水这类总量数据。版本不支持 sum 时退而求其次用 average再把结果乘上源像元面积 / 目标像元面积做换算但这只是近似不如 sum 干净。重采样方法适用数据类型对人口总量的影响nearest分类、土地利用重复或丢失格网总量失真大bilinear连续表面抹平峰值总量略变average降尺度通用总量须按面积比例换算sum人口等总量型守恒推荐升尺度转更细分辨率时没有完美方案任何插值都是想象出的分布总量保持但空间格局被平滑输出时要在说明里标注重采样方法和日期否则后手会把你平滑过的东西当原始数据用。4. 人口数据可视化分级设色、百分位拉伸与边界叠加数据处理完到出图环节。1km 人口数据可视化的难点不在工具而在人口分布极端不均匀绝大多数格网接近 0少数超大城市一个格网就好几万。默认的最小值-最大值拉伸会把色带全压在城市上全国图几乎全白所以出图的核心是拉伸策略。4.1 QGIS 流程三分钟出一张全国人口分布图我常用 QGIS 做快速出图操作路径固定拖入 TIF右键图层属性Symbology 选 Singleband pseudocolor色带选 YlOrRd 反转深色表示高值Mode 选 Quantile 分成 8 类Min/Max 不勾自动改用手动范围并按 2%98% 百分位裁剪。参数推荐值理由渲染类型Singleband pseudocolor连续数值按色带上色比灰度可读色带YlOrRd 反转低值浅黄、高值深红符合直觉ModeQuantile8 类等频分级避免大量格网挤在一个色级拉伸范围2%98% 百分位排除极小/极大值干扰Clip out of range勾选超范围像元显示为透明不被截断人口数据高度右偏等间距分级会让 80% 的格网落在第一级里图没法看分位数分级保证每级格网数量大致相同读图时看到的是相对稠密区而不是绝对数值。想进一步压缩动态范围可以对数据做 log1p 变换后再分级但出图说明里要写清楚变换方式。4.2 Python 出图百分位拉伸 色带的完整脚本不依赖 QGIS 的自动化出图用 matplotlib 三件事就能说清读数组、定拉伸范围、画图。完整脚本如下import matplotlib.pyplot as plt import numpy as np import rasterio with rasterio.open(CHN_ppp_v2b_2020.tif) as src: band src.read(1) # imshow 的 extent 顺序是 left, right, bottom, top extent [src.bounds.left, src.bounds.right, src.bounds.bottom, src.bounds.top] crs src.crs # 无效值和非人口区域统一置为 nan让它们不上色 band np.where(band 0, np.nan, band) # 百分位拉伸避开少数大城市把色阶拉平 p2, p98 np.nanpercentile(band, [2, 98]) fig, ax plt.subplots(figsize(12, 9)) im ax.imshow( band, extentextent, cmapYlOrRd, vminp2, vmaxp98, interpolationnearest ) ax.set_title(中国人口空间分布1km 格网, fontsize14) cbar plt.colorbar(im, axax, fraction0.03, pad0.02) cbar.set_label(人 / 格网) # 若是密度数据请改成 人 / 平方公里 plt.savefig(china_pop_1km.png, dpi300, bbox_inchestight)extent必须按 left/right/bottom/top 顺序组装直接用src.bounds会把边界顺序弄混图像坐标就错了。interpolationnearest是关键matplotlib 默认插值会把人口栅格平滑成热力图误导人以为人口连续渐变。nanpercentile计算百分位时自动跳过 nanNoData 和负值不参与拉伸。vmin/vmax取 2% 和 98% 分位是为了出图对比度而不是为了数值准确——这跟统计是两回事。4.3 叠加边界与色标让图和业务对得上光有栅格没有行政界线业务没法用。接着上一段代码叠加省界import geopandas as gpd prov gpd.read_file(province_boundary.geojson) if prov.crs is not None and prov.crs ! crs: prov prov.to_crs(crs) prov.boundary.plot(axax, color#444444, linewidth0.4) plt.savefig(china_pop_1km_boundary.png, dpi300, bbox_inchestight)叠加前必须确认两侧坐标系一致否则在 1km 尺度上偏移几百米一般看不出但统计对不上。输出印刷级图建议直接在 QGIS 打印布局里做矢量边界和栅格一起导出图例、比例尺、指北针都齐备Python 脚本更适合批量出图或做多期对比。5. 常见问题与避坑五个人口栅格高频翻车点这几条是我和同事在实际项目里真实踩过、并且每条都赔过时间的坑。按现象 → 原因 → 解决的顺序写方便对号入座。5.1 现象TIF 拖进软件全黑或全白什么都看不见排第一的坑。1km 人口数据像元值从 0 到几万90% 的像元集中在低值区只有少数城市格网到几万。默认渲染用全幅最小最大值做线性拉伸色带全被极端高值占据剩下区域在色带里分不到一级看上去就是全黑或全白。原因不是数据坏是渲染策略不适合右偏分布。解决QGIS 里 Symbology 选 Singleband pseudocolorMode 选 QuantileMin/Max 手动改成 2%98% 百分位Python 出图用np.nanpercentile算vmin/vmax。想进一步拉开低值区差异可以做 log1p 变换后渲染但图例要标注变换方式。5.2 现象区域加总人口只有年鉴的一成甚至出现负数这个坑我印象最深。第一次用这份数据算一个地级市结果只有年鉴的 8%我一度以为是数据坏了。逐格排查才发现 NoData 是 -99999统计时没排除无效值几十万无效像元全被当成真实人口求和直接把结果拉崩。另一种常见原因是口径混用把密度数据当每像元人数直接加总或者把户籍口径数据拿去和常住人口年鉴对比差出 20% 是常态。解决分三步先 gdalinfo 确认 NoData 值和数据类型再全省加总和七普常住人口对比确定这份文件属于哪类口径最后统计时显式传nodata参数不要依赖函数默认值。5.3 现象裁剪结果和行政边界错位像元好像漂了半格矢量边界是 CGCS2000栅格是 WGS84两者经纬度数值只差零点几弧秒叠图在省级尺度几乎看不出放大到城市、区县尺度几百米偏移就非常明显。还有一类更隐蔽国内机构发布的人口栅格用自定义 Albers 投影没有标准 EPSG 编码读出来是一长串 PROJCS 字符串如果偷懒不处理偏移量会到公里级。解决每次处理前打印两侧 CRS 对比不一致就to_crs把矢量转到栅格坐标系不要凭文件名猜坐标系一切以元数据输出为准。5.4 现象山区或建成区边界出现一排空白像元裁剪用cropTrue时多边形边界是弧线像元是方块弧线切过像元时边缘像元归属哪边存在取舍如果边界精度高于栅格边界上只有小半落在多边形内的像元会被丢弃视觉上就是边界锯齿加一圈空白。解决统计环节用rasterstats时传all_touchedTrue把与多边形相交的所有像元纳入计算裁剪环节先给边界做缓冲0.01° 或一个像元尺寸裁完再用原始边界做可视化遮罩。对全省层面影响可忽略但对市区级别的小范围统计一排像元可能代表上万人。5.5 现象重投影后全国总人口变了 20% 以上WGS84 经纬度网格转成 CGCS2000 或 Albers 投影时目标网格和源网格像元位置不可能完全对齐重采样必须决定每个新像元从哪些源像元取值。nearest 取最近的一个会重复取某些像元、丢掉另一些人口总量瞬间失真bilinear 对四邻加权平均城市峰值被抹平形状对了但总量变。解决降尺度用Resampling.sum保证新网格加总等于原网格加总升尺度没有守恒方案必须标注插值来源涉及面积或密度计算先把栅格转入 Albers 等积投影再用真实像元面积算不要在经纬度坐标下用固定系数糊弄。注意重采样参数要在成果目录里留一份说明标明源分辨率、目标分辨率、重采样方法和日期防止后手把插值结果当实测数据。6. 进阶分区统计、总量校验与一张可用的表处理完栅格、出完图业务最终要的是这个省多少人、那个市多少人。分区统计和总量校验是收尾的关键两步用 rasterstats 一次跑完import geopandas as gpd from rasterstats import zonal_stats provinces gpd.read_file(province_boundary.geojson) stats zonal_stats( provinces, CHN_ppp_v2b_2020.tif, stats[sum, mean, max], nodata-99999, all_touchedTrue, geojson_outTrue ) df gpd.GeoDataFrame.from_features(stats) result df[[NAME, sum, mean, max]].sort_values(sum, ascendingFalse) result.to_csv(province_pop_1km.csv, indexFalse) print(result.head(10))stats[sum, mean, max]中 sum 才是人口总数mean 是平均每格网人数max 是最高密度格网nodata-99999必须和文件实际 NoData 一致否则负数参与求和all_touchedTrue把边界相交像元算进去避免区县级边界切漏人口。输出 GeoDataFrame 直接写 CSV这张表就能进 PPT 或对接业务系统。统计完必须做总量校验把各省 sum 加总和全国普查常住人口对比。世界公开数据有 UN adjusted 和 unadjusted 两个版本后者没按联合国分年龄性别数据校准低龄和老年人口偏差明显国内机构的格网还分常住、户籍、夜间灯光修正版。误差在 ±5% 以内可以放心用超过 20% 先回查 NoData 和坐标系再查口径。1km 格网在人口稀疏区本来就有空间平滑拿来比较城市间相对密度没问题但拿某个偏远县的绝对人口去做预算审批风险自负。从那以后我每次拿到人口 TIF 都强制走一遍固定三步gdalinfo 看元数据、全省加总对口径、边界 CRS 对齐截图确认然后才允许进入可视化和统计。这三分钟换来的是一次做对而不是反复返工。希望帮到你。本文还有配套的精品资源点击获取