
简介中国第七次人口普查网格化人口数据集空间分辨率100米、时间口径为2020年是面向GIS、人口地理、城市规划及灾害风险评估研究者的全国人口栅格产品。基于第七次人口普查数据生成不但适合区域人口分布特征分析和人口密度制图也能与WorldPop、LandScan等国际人口产品在2020年中国范围内进行空间对比与精度验证。压缩包采用rar格式总大小126.04MB共包含8个文件核心为GeoTIFF格式的全国100m人口栅格同时提供tfw坐标配准参数、ovr金字塔用于快速缩放显示、xml元数据、jpg预览图以及dbf和cpg栅格属性表可直接导入ArcGIS和QGIS分析与制图。该资源已有171人学习适合用于人口密度分析、区域人口估算、城市扩张监测、人居环境与灾害风险评估等教学科研任务。借助这套数据用户能快速获取全国范围100m分辨率的七普人口栅格可对任意区域进行分区统计、多尺度聚合并可提取栅格值表开展专题制图为与其他人口栅格产品进行一致性比较和人口空间自相关研究提供可靠数据基础。1. 100m 栅格人口数据集从普查公报到网格的落地做区域人口分析时最缺的往往不是矢量边界而是能直接和道路、水系叠加的全国人口栅格数据。中国第七次人口普查网格化人口数据集把 2020 年普查结果重新分配到 100m 分辨率网格上每个像元值代表该格网内的估计人口适用于省、市、县到街区的多尺度分析。它和 WorldPop、LandScan 一样都是把统计人口“散开”到空间里但底层约束和建模思路有明显差异。对 GIS 工程师而言先搞清楚文件包中每个文件的作用再按统一方法做区域统计和跨数据比较才能避免数值误用。下面直接从前述文件包的实际内容说起逐步给出可复现的 GDAL 处理路径。2. 文件包拆解GeoTIFF、TFW、VAT 与金字塔之间的关系2.1 文件清单中每个后缀的角色压缩包解开后会看到一批同名不同后缀的文件很多人只认得 .tif拷数据时漏掉 .tfw 或 .ovr到了 Linux 服务器上就出现位置漂移或缩放缓慢。主文件PopSE_China2020_100m.tif是真正的栅格数据内部通常已经带有地理标签但有时软件也依赖外部的 .tfw 世界文件来定位。.tfw文本里保存六个仿射参数比如第 3、6 行是地图坐标原点的横纵坐标第 2、4 行是像元旋转和分辨率一旦缺失ArcGIS 和 QGIS 会尝试按“无坐标”方式打开。.tif.ovr是金字塔概览用于加速全图显示.vat.dbf与.vat.cpg是值属性表记录每一个像元值出现的次数和统计信息方便你在 ArcCatalog 里连接。.aux.xml是 GDAL 辅助信息可能记录 NoData、色彩映射和统计直方图。.jpg是预览图不是数据本身。这些文件共同描述一个完整的栅格数据集缺一不可。这个文件列表最能反映一个点国内的数据包往往带 VAT 与 TFW 组合而 WorldPop 通常只有裸 TIF。处理前要用gdalinfo检查 TIF 内嵌的元数据是否完整如果不完整建议保留辅助文件。因为 100m 全国范围跨度大TIF 文件经常超过 10GB.ovr对缩略图与在线发布极其关键。删除它虽然不影响数据精度但会让 web map 服务在缩放时反复读取原始像元拖慢响应。文件后缀内容说明工程上的作用.tif主栅格数据地理分析、统计、建模.tfw世界文件存储仿射变换参数定位栅格角点.tif.ovr金字塔概览保证大图缩放与预览性能.tif.aux.xmlGDAL 辅助元数据记录 NoData、坐标系统等.vat.dbf值属性表说明每个像元值与数量的关系便于分类显示.jpg预览图人工快速查看大致效果2.2 用 gdalinfo 验证坐标参考系与 NoData直接打开 TIF 不是不行但我更习惯先在终端执行gdalinfo PopSE_China2020_100m.tif这样能快速获得 Size、Coordinate System、GEOGCS 信息以及 Band 1 的数据类型。如果输出末尾出现NoData Value-99999说明无数据的像元已经定义如果没出现就要用gdalinfo -stats PopSE_China2020_100m.tif让 GDAL 扫描一遍数据输出最小值和最大值。注意最小值若是 0不能马上断定无数据是 0因为大量山区确实没人住像元值本来就应该是 0。区分方法很简单查看 .vat.dbf 或原数据说明。就这个数据集而言明确基于普查人口空间化无数据通常为 -99999而非 0。Coordinate System通常显示为 EPSG:4326 或 EPSG:4490两者都属于经纬度坐标但 4490 是 CGCS2000 地理坐标系4490 与 4326 的差值虽在米级但做高精度影像叠加时建议统一到 CGCS2000 或 WGS84 其中一个。2.3 理解人口网格的统计口径避免数值误读第七次人口普查公布的是“常住人口”WorldPop 和 LandScan 虽然也用普查数据做约束但引入其他地理因子重新分配最终结果并不是普查数本身。这个数据集的网格值应该理解为一个空间插值结果每一个像元值可能是通过权重模型从乡镇街道层面离散到 100m 格网的。因此你无法把一个像元和实际家庭一一对应只能做区域聚合分析。比如在乡镇边界聚合后人口总和与公报一致但在单个像元抽出某个具体数值放到现场核对可能偏差巨大。这提醒我们处理栅格时不要用 Point Sampling 直接去挑某个坐标的人口值而应先做区域求和。VAT.dbf正是服务于分类图与直方图而不是精确人口表的。后面所有统计步骤都应该以像元聚合为基础。3. 用 Python 与 GDAL 批量读取人口栅格并做区域统计3.1 使用 Rasterio 分块读取避免内存爆炸全国 100m 分辨率栅格的尺寸通常有数万列、数万行直接src.read(1)会把成亿级像素一次性读入内存导致机器卡死或 OOM。常见做法是使用rasterio.windows.Window分块读取同时获取窗口内的仿射变换这样既能控制内存也能把读取结果与地理坐标对应起来。import rasterio from rasterio.windows import Window with rasterio.open(PopSE_China2020_100m.tif) as src: profile src.profile print(分辨率:, profile[width], x, profile[height]) nodata src.nodata # 从左上角读取 4096x4096 像元 w Window(0, 0, 4096, 4096) data src.read(1, windoww) transform src.window_transform(w)Window(0, 0, 4096, 4096)的参数分别是列偏移、行偏移、宽、高顺序不能反。读取到的data是一个二维 numpy 数组transform是Affine对象用于把行列坐标转成地图坐标。如果想遍历全部数据可以用src.block_shapes获得每个分块的大小然后逐块处理。nodata变量用于判断无效像元通常是一个负数比如 -99999。注意人口栅格的像元值可能超过 2^31读取时最好用src.read(..., out_dtypefloat64)或data.astype(float64)避免 32 位整数求和溢出。3.2 按行政区边界裁剪并汇总人口拿到区县边界 GeoJSON通常可以采用rasterio.mask函数批量裁剪。这个函数要求输入面要素的坐标参考系统与栅格一致。如果不一致先重投影。以下代码遍历每个区县统计每个区域的人口总数import geopandas as gpd import numpy as np from rasterio.mask import mask import rasterio adm gpd.read_file(districts.geojson) result [] for _, row in adm.iterrows(): geom row.geometry with rasterio.open(PopSE_China2020_100m.tif) as src: out_image, out_transform mask( src, [geom], cropTrue, nodata-99999 ) out_image out_image[0].astype(float64) out_image[out_image -99999] np.nan if np.isnan(out_image).all(): continue population np.nansum(out_image) result.append({name: row[name], population: population}) print(pd.DataFrame(result).head())mask函数的第一个参数是栅格对象第二个参数是几何列表cropTrue会在裁剪之前计算最小外接范围降低 IO 读取量。nodata-99999告诉函数哪些像元是无数据避免把它们当作 0 参与统计。代码中将无数据替换为np.nan是因为np.nansum可以安全跳过无效值。如果adm包含多个县但人口字段类型是浮点也可以直接对裁剪后的每个数组做加权平均。这里population的单位取决于 TIF 文件像元值如果像元值是“每像元人数”直接求和如果像元值是“每平方公里人口密度”就要乘上像元面积100×100 米 0.01 km²。这一步骤必须从数据提供方确认不能靠猜。若在公文里发现某县人口少了一半多半就是这里出了问题。频繁打开栅格文件会拖慢速度。可以先把整个 TIF 打开一次然后在循环内调用 mask但 mask 函数内部会重新打开src。如果数据量很大建议使用rasterio.vrt.WarpedVRT将文件整体虚拟重投影到统一 CRS再用 global 的窗口算子但这需要更多代码。此处给出的是最直观的做法。3.3 坐标系不一致时的重投影与注意点行政边界如果来自高德或百度地图坐标系可能是 GCJ-02 或 BD-09绝不能直接用于掩膜。全域都会偏移几百米。建议用geopandas转换到与栅格一致但要注意转换参数是七参数还是简化三参数。对于全国尺度的统计误差影响较小但对于街道级精细裁剪需要确认原 CRS。# 先统一到 WGS84再转到栅格自身 CRS adm_wgs84 adm.to_crs(EPSG:4326) with rasterio.open(PopSE_China2020_100m.tif) as src: adm_proj adm_wgs84.to_crs(src.crs)上述代码里adm的原始坐标系若不是 EPSG:4326会以原始坐标为准先做一次转换。由于国内常见的 GeoJSON 可能是 GCJ02单纯调用to_crs并不会进行火星坐标纠正必须先用coord_transform函数施加密钥。这个点经常被忽略GDAL 不支持 GCJ-02 转换所以如果你发现边界与人口网格偏移几百米直接放弃to_crs改用高德官方接口或内部转换算法。对比这些全球数据集时最好统一使用 ESPG:4490 或 4326这样比较才有意义。4. 与 WorldPop、LandScan 对比分辨率、口径与验证方法4.1 三套数据的建模差异这个数据集与 WorldPop、LandScan 最大的不同在于其底层约束是第七次人口普查而 WorldPop 用的是多个年份的统计组合LandScan 则分布到每天 24 小时的平均人口。因此三类数据很难在栅格像元级别直接对比。WorldPop 通常使用随机森林从道路、建筑、夜间灯光、植被指数等数据中学习人口分布LandScan 使用“环境模式”将人口权重分配给合适的土地利用类型。100m 分辨率不代表它可以反映单个建筑的真实人数最多是空间化后的平滑人口。理解这一点之后你就不会再问“为什么某个小区像元值是 0但那里明显有房子”了。这些模型输出的是“估计概率分布”不是人口普查登记名册。指标公式说明绝对误差 AE预测值 - 真值对区域总量直接相减相对误差 RE预测 - 真值 / 真值评价偏差比例MAE平均绝对误差聚合到单元的平均误差RMSE均方根误差放大离群单元的影响R²1 - SSE/SST评估两套数据线性相关程度4.2 构建一个可重复的两两对比流程首选把多套数据聚合到同一网格单元。比如把全国按 1km×1km 划分然后分别求每个网格内的总人口。使用rasterstats的zonal_stats可以同时完成区域统计对比流程建议如下from rasterstats import zonal_stats import pandas as pd grid_stats zonal_stats( china_grid_1km.shp, PopSE_China2020_100m.tif, statssum, nodata-99999, geojson_outTrue, ) worldpop_stats zonal_stats( china_grid_1km.shp, worldpop_2020.tif, statssum, nodata-99999, geojson_outTrue, ) df pd.DataFrame(grid_stats) df[worldpop_sum] [f[properties][sum] for f in worldpop_stats] df[popse_sum] [f[properties][sum] for f in grid_stats] df[rel_diff] (df[popse_sum] - df[worldpop_sum]) / (df[worldpop_sum] 0.01) print(df[df[worldpop_sum] 0].describe())这里的zonal_stats参数简洁但要注意china_grid_1km.shp必须与 TIF 坐标系一致否则会使用慢速的经纬度转换。使用statssum时如果 TIF 内有 NoData函数会自动跳过。输出的字段名默认放在 properties 里所以使用 geojson_out 才能合并数据。后续可以从 DataFrame 里筛选相对误差大于 20% 的网格看看它们分布在哪个区域。4.3 注意分辨率差异带来的人工噪声100m 与 1km 数据重叠时即使做了网格聚合尺度效应仍然存在。100m 栅格在山区和农田地区变化剧烈相邻像元可能从 0 跳到 20 人而 1km 数据在聚合时已经平滑了这些波动。让两套数据在同一尺度比较时必须先将原始人口栅格用gdalwarp -r sum聚合到目标分辨率再开始统计。如果直接在原始分辨率下做差值会出现大量斑块噪声。正确的验证方法是把 100m 聚合到 1km然后计算“相对偏差”并绘制空间分布图。至于聚合后的结果一般要求有 95% 以上的网格人口误差在正负 30% 以内。如果达不到多半是数据结构有问题或坐标系偏移。这一步骤不需要改原始文件只用于建模评估。5. 重投影、掩膜提取与制图出图的进阶技巧5.1 重投影时保留人口总量用gdalwarp对人口栅格做重投影时千万不要用默认的最近邻重采样这可能改变每个像元的统计总量。一般我会用gdalwarp -t_srs EPSG:3857 -tr 100 100 -r sum -overwrite PopSE_China2020_100m.tif PopSE_web.tif-r sum确保在重采样过程中将覆盖的多个像元相加而不是取均值或最近邻。输出的像元分辨率设定为 100m。这里有一个限制如果目标坐标系是 3857100m 在纬度 60 度附近代表实际距离约 200m人口会因为像元面积变大而失真。因此做人口统计时应优先选择 Albers 等积投影EPSG:102025而不是 Web 墨卡托。如果只用于 Web 地图显示那可以保留 3857但后续任何数值统计都要重新投影回等积投影。5.2 提取圆形缓冲区人口用于选址分析商业选址时经常需要知道一个 POI 周边 2km 范围内的总人口。这时可以在 Python 中生成圆再执行 mask 提取。from shapely.geometry import Point import geopandas as gpd point Point(116.391, 39.905) gdf gpd.GeoDataFrame(geometry[point], crsEPSG:4326) # 转到投影坐标系后再做 2km 缓冲区 gdf_proj gdf.to_crs(EPSG:32650) gdf_proj.geometry gdf_proj.geometry.buffer(2000) # 再转回原始经纬度用于后续 rasterio.mask gdf_wgs84 gdf_proj.to_crs(EPSG:4326)更严谨的写法是先在投影坐标系中建立缓冲再转回地理坐标。比如用gdf.to_crs(EPSG:32650)然后gdf.geometry[0] gdf.geometry[0].buffer(2000)这样缓冲区半径 2000 米精确。之后再对 gdf 做mask提取。这个过程能看到人口热力分布与核心路网的关系很适合网格化人口数据落到商业场景。5.3 制图时用对数拉伸表现低密度区域出图时直接用线性拉伸会把绝大多数面积渲染成纯黑或纯白因为中国大量地区每平方公里只有几十人少数城市核心却超过上万人。建议在 QGIS 符号系统中选择“单波段伪彩色”渲染类型选择“对数”或自定义公式log10(x1)。同时将色带设定为高透明度的 viridis防止遮盖底图。如果使用gdal_translate导出 PNG可以使用-scale加上自定义-outsize但不建议修改原始 TIF。显示层的调整始终与统计层分离这是一个工程准则。这样发布出的底图既能表现全国人口密度形态又不会给决策者造成“到处都是高密度”的错觉。本文还有配套的精品资源点击获取