ARTICLE DETAIL

资讯详情

深耕编程入门与网站建设的一线实战洞察。

中国人口密度公里格网栅格数据:从读取、裁剪到分区统计的完整流程

中国人口密度公里格网栅格数据:从读取、裁剪到分区统计的完整流程 简介这份资源是中国人口密度公里格网栅格数据包面向地理信息系统分析、城市规划、人口研究与社会科学领域的从业者和学生用于识别人口密集区与稀疏区、支撑公共服务布局与生态承载力评估等空间分析场景。压缩包共240个文件约60.58MB以shp、dbf、prj、shx、cpg等Shapefile配套文件为主另有adf栅格文件与xml元数据覆盖矢量边界、属性表、投影信息与栅格像元值可直接导入ArcGIS、QGIS等软件使用。目前已有2732人学习下载。数据以1平方公里格网组织每个像元对应人口密度值便于进行空间查询、叠加分析与可视化渲染读者可据此绘制人口分布专题图、开展城乡差异与人口迁移趋势研究或作为政策制定与科研建模的基础数据省去自行清洗与格式转换的环节。1. 中国人口密度公里格网栅格数据从压缩包到可分析图层的落地路径拿到「中国人口密度公里格网栅格数据.zip」这个压缩包时多数人的第一反应是解压、拖进 GIS 软件、出图。但真正做过空间分析的人都知道从压缩包到能参与建模的图层中间隔着坐标系、无效值、量纲和分辨率匹配四道坎。这份数据本质上是把全国人口按 1km×1km 网格做面积加权分配后形成的栅格表面每个像元值代表该平方公里内的人口数单位通常是人/km²。它能解决的核心问题是当你手头只有行政区划级别的统计人口却需要做站点选址、资源可达性、环境暴露评估这类需要连续空间粒度的工作时公里格网就是最直接的底图。适合谁用做城乡规划、公共卫生、物流选址、遥感反演验证的从业者以及需要把人口作为协变量接入机器学习模型的数据工程师。这篇笔记按「先搞懂数据长什么样、再跑通读取与裁剪、最后处理坑」的顺序展开每一步都给出可复现的命令和参数。2. 公里格网人口数据的结构、量纲与坐标系先看清再动手2.1 栅格值到底代表什么密度还是总量很多人第一次打开这类数据看到像元值在几百到几万之间会误以为是该网格的总人口。实际上公里格网人口数据的标准做法是密度值即人/km²。因为 1km×1km 的网格面积恰好是 1km²所以密度值在数值上等于该网格的总人口这也是公里格网被广泛采用的原因——密度和总量在数值上统一了。但要注意如果数据经过了投影变换或重采样网格面积不再是精确的 1km²这个等价关系就会失效。判断方法很简单找一个已知人口约 10 万的乡镇看覆盖它的网格像元值之和是否接近 10 万。如果差了一个数量级说明你拿到的可能是总量栅格而非密度栅格。常见做法是先用 GDAL 读取元数据确认像元大小和坐标参考。下面这段 Python 用 rasterio 打开数据打印关键信息import rasterio import numpy as np # 打开压缩包解压后的 tif 文件 with rasterio.open(china_pop_1km.tif) as src: print(CRS:, src.crs) # 坐标参考系 print(分辨率:, src.res) # 像元大小单位与 CRS 一致 print(行列数:, src.height, src.width) print(波段数:, src.count) print(无效值:, src.nodata) # 读取第一波段 arr src.read(1) # 统计有效像元的分布 valid arr[arr ! src.nodata] if src.nodata is not None else arr.ravel() print(有效像元数:, valid.size) print(像元值范围:, valid.min(), valid.max()) print(像元值均值:, valid.mean())逻辑说明src.res返回的像元大小单位取决于 CRS。如果 CRS 是地理坐标系如 EPSG:4326单位是度1km 约等于 0.0083 度如果是投影坐标系如 Albers 等积投影单位是米像元大小应接近 1000。参数上src.nodata可能是 -9999、-3.4e38 或 None这直接决定后续统计是否会把无效值算进去。我一般会先跑一遍这段把 CRS、分辨率、无效值三个信息记下来再决定要不要重投影。2.2 坐标系选择为什么等积投影不能省地理坐标系下的栅格每个像元的实际面积随纬度变化。在中国范围内纬度从 18°N 到 53°N1km×1km 的网格在高纬度地区的实际面积会明显缩小。如果你直接用 EPSG:4326 的数据做面积统计或密度计算北方地区的人口会被系统性低估。正确做法是转到等积投影国内常用 Albers 等积投影参数为中央经线 105°E、双标准纬线 25°N 和 47°N。用 gdalwarp 重投影的命令如下gdalwarp -t_srs projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm \ -tr 1000 1000 \ -r bilinear \ -srcnodata -9999 -dstnodata -9999 \ china_pop_1km.tif china_pop_1km_albers.tif参数说明-t_srs指定目标投影这里用 PROJ 字符串写 Albers-tr 1000 1000强制输出像元为 1000m×1000m保证网格面积精确为 1km²-r bilinear是重采样方法人口密度属于连续表面双线性比最近邻更平滑但如果你的分析要求总量守恒应改用-r sum并配合面积权重这一点后面避坑章节会展开。-srcnodata和-dstnodata保持一致避免无效值在重采样时被当成真实值参与计算。2.3 数据分块与压缩包内的文件组织这类压缩包通常包含一个主 tif 文件外加 .tfw 世界文件、.prj 投影文件、.aux.xml 辅助文件有时还有一份 readme 或元数据 xml。解压后先确认文件是否齐全尤其是 .prj 缺失时GDAL 会无法识别 CRS读出来是 None。如果压缩包内按省份或年份分了多个 tif建议先用 Python 批量检查每个文件的 CRS 和分辨率是否一致不一致的先统一再合并。合并用 gdal_merge 或 rasterio 的 merge 都可以但要注意合并前所有文件的 nodata 值必须相同否则无效值会污染拼接结果。3. 用 Python 跑通读取、裁剪与分区统计的最小流程3.1 按行政区裁剪mask 与 clip 的差别拿到全国数据后绝大多数分析只需要一个省或一个流域。裁剪有两种思路一是用矢量边界做 mask保留边界内像元边界外设为 nodata二是用矩形范围做 clip速度快但不精确。做人口统计必须用 mask否则边界外的像元会被错误计入。下面用 rasterio 的 mask 功能按矢量边界裁剪import geopandas as gpd import rasterio from rasterio.mask import mask # 读取行政区矢量确保与栅格 CRS 一致 gdf gpd.read_file(boundary.shp) with rasterio.open(china_pop_1km_albers.tif) as src: gdf gdf.to_crs(src.crs) # 矢量转到栅格坐标系 geoms [geom for geom in gdf.geometry] # 裁剪cropTrue 会缩小输出范围 out_image, out_transform mask(src, geoms, cropTrue, nodata-9999) out_meta src.meta.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform, nodata: -9999 }) # 写出裁剪结果 with rasterio.open(province_pop.tif, w, **out_meta) as dest: dest.write(out_image)逻辑说明mask函数返回裁剪后的数组和新的仿射变换。cropTrue会把输出范围收缩到矢量边界的外接矩形减少文件体积。参数nodata-9999要和你原始数据的 nodata 一致否则边界外像元会保留原值。注意gdf.to_crs(src.crs)这一步不能省矢量和栅格 CRS 不一致时 mask 会报错或产生错位。3.2 分区统计用 zonal stats 算每个县的人口裁剪出省份后下一步往往是按县统计人口总量。因为像元值是密度每个像元面积是 1km²所以直接对县边界内的像元值求和就是该县总人口。用 rasterstats 可以一行搞定from rasterstats import zonal_stats # 按县矢量统计人口总和 stats zonal_stats( county.shp, province_pop.tif, stats[sum, mean, count], nodata-9999, geojson_outFalse ) # stats 是一个列表每个元素对应一个县 for i, s in enumerate(stats[:5]): print(f县 {i}: 总人口{s[sum]:.0f}, 均值{s[mean]:.1f}, 有效像元{s[count]})参数说明stats[sum,mean,count]分别返回总和、均值和有效像元数。nodata-9999告诉 rasterstats 忽略无效值。count是参与统计的像元数如果某个县的 count 为 0说明该县边界与栅格没有重叠通常是 CRS 不一致或边界文件有问题。这里 sum 的单位是人因为密度×面积1km²总量。如果之前重采样改变了像元面积sum 就不再等于总人口需要乘以实际像元面积。3.3 把栅格转成表格接入机器学习流程做选址模型或回归分析时往往需要把栅格值提取到点或面上。用 rasterio 的 sample 可以批量提取点位人口密度import rasterio import pandas as pd # 读取点位数据 points pd.read_csv(sites.csv) # 含 lon, lat 列 coords list(zip(points[lon], points[lat])) with rasterio.open(china_pop_1km_albers.tif) as src: # 如果点位是经纬度需要先转到栅格 CRS from pyproj import Transformer transformer Transformer.from_crs(EPSG:4326, src.crs, always_xyTrue) coords_proj [transformer.transform(lon, lat) for lon, lat in coords] # 采样 values [v[0] for v in src.sample(coords_proj)] points[pop_density] values points.to_csv(sites_with_pop.csv, indexFalse)逻辑说明src.sample接收投影后的坐标列表返回每个点的像元值。如果点位本身已经是投影坐标跳过 Transformer 那一步。参数上always_xyTrue保证经纬度顺序是 lon, lat避免 x/y 颠倒。采样得到的值是密度如果模型需要总量还得结合点位所在网格的面积但公里格网下面积就是 1km²直接用密度即可。4. 避坑与排查公里格网人口数据最常见的五个翻车点4.1 无效值被当成真实人口参与统计现象某个县统计出来的人口比实际多了几十万或者全国求和后远超官方总人口。原因数据中的 nodata 可能是 -9999但你在统计时没有排除负值被计入求和或者 nodata 被设成了 0而 0 在人口密度里是合法值无人区导致无法区分。解决先用src.nodata确认无效值标记统计时显式传入nodata参数。如果 nodata 是 0 且你确定无人区也用 0 表示那就要接受这个模糊性或者用掩膜矢量排除水域和无人区。4.2 重采样后总量不守恒现象重投影到 Albers 后全省人口总和比原始数据少了或多了百分之几。原因双线性重采样是加权平均不保证总量守恒最近邻虽然保总量但会引入块状伪影。解决如果分析目标是总量用-r sum重采样并确保输出像元面积与输入一致如果目标是密度表面双线性可以接受但要在报告中说明总量可能有微小偏差。我一般会在重投影后跑一遍分区统计和官方统计对一下偏差超过 2% 就回去检查参数。4.3 矢量与栅格 CRS 不一致导致裁剪错位现象裁剪出来的省界和栅格上的实际边界对不上整体偏移几百米到几公里。原因矢量是 EPSG:4326栅格是 Albers直接 mask 时 GDAL 不会自动转换而是按数值强行套用。解决裁剪前统一用to_crs把矢量转到栅格 CRS。如果矢量没有 CRS 定义先用set_crs指定再转换。这个坑在跨省分析时尤其常见因为不同来源的矢量 CRS 可能不同。4.4 压缩包内文件缺失 .prj 导致 CRS 丢失现象rasterio 打开文件后src.crs返回 None后续所有投影操作报错。原因压缩包在打包时漏掉了 .prj 文件或者解压工具没有正确释放隐藏文件。解决如果知道数据来源手动指定 CRS用rasterio.open的crs参数或 gdal_edit 补上。不确定来源时用数据范围反推中国范围经度 73°E-135°E纬度 18°N-53°N如果像元值范围对应这个经纬度大概率是 EPSG:4326。4.5 像元值量纲误判导致模型输入偏差现象把密度值直接当成总人口输入回归模型系数解释完全错误。原因没有确认像元面积是否为 1km²或者数据本身是总量栅格而非密度栅格。解决用src.res确认像元大小如果是 1000m 且 CRS 是投影坐标系面积就是 1km²密度等于总量如果是 500m面积是 0.25km²总量等于密度×0.25。最稳妥的办法是找一个已知总人口的行政区用分区统计求和和官方数据对比反推量纲。5. 进阶技巧用多分辨率金字塔和分区并行加速全国分析全国公里格网数据量不小一个波段约 1.5 亿像元直接全量读取做分区统计会吃满内存。我一般会先建金字塔再按省并行。建金字塔用 gdaladdogdaladdo -r average --config COMPRESS_OVERVIEW DEFLATE china_pop_1km_albers.tif 2 4 8 16 32参数说明-r average对密度数据用平均值降采样保持密度语义2 4 8 16 32是金字塔层级对应 2×2、4×4 等聚合。建完金字塔后缩放浏览和低精度统计会快很多。并行分区统计可以用 Python 的 multiprocessing按省拆分矢量每个进程独立读取栅格窗口。注意 rasterio 的多进程要每个进程单独 open不能共享句柄。我习惯把全国数据按省预裁剪成小文件再并行跑 zonal_stats这样单省内存占用可控整体速度比单进程快 4 到 6 倍。最后提醒一句所有统计结果落地前务必和官方行政区人口对一遍偏差超过 5% 就回去查 nodata 和重采样参数。希望帮到你。本文还有配套的精品资源点击获取
返回列表