简介：面向地理信息系统分析、城市与区域规划及社会科学研究人员，这份资源提供了中国人口密度公里格网栅格数据。数据按一公里乘一公里的网格组织，包含覆盖全国的人口密度栅格图层及配套矢量边界与投影定义文件，可用于人口分布可视化、公共服务设施布局、生态承载力评估等场景。压缩包共二百四十个文件，以矢量文件承载边界与属性，以栅格文件存储密度值，并附有投影与元数据说明，整体约六十兆字节，解压后可直接导入主流地理信息软件。目前已有两千七百三十一人学习下载。数据经预处理与标准化，用户可快速进行叠加分析、密度分级渲染或与其他地理数据整合，为研究人口流动、城乡差异及区域发展提供扎实数据底座。
1. 人口密度公里格网栅格数据：把“哪个区住多少人”摊平到一公里格网
做商业选址、应急疏散评估或者基站规划时，最大的尴尬是手里只有区级统计年鉴——全市五百万人，可你压根不知道这五百万人散在哪几条街道。这份人口密度公里格网栅格数据，把人口估算放到一公里见方的格网上，每个像元直接对应一个区域的常住人口规模或密度，配合 Python 和 GIS 就能裁出任意多边形范围内的人口总量。它特别适合城市规划、零售网点评估、交通流量预测这类需要“落到地面上看人口”的场景。跟传统普查数据相比，它的优势是空间分辨率高到能看出人口梯度变化，代价是得先熟悉栅格数据的坐标系、NoData 值和重采样这些基础概念。
2. 读懂一平方公里格网：分辨率、坐标系与像元值的门道
2.1 公里格网数据到底长什么样
人口密度公里格网数据最常见载体是 GeoTIFF，本质是一个二维数组加一套地理配准信息。数组的每个元素是一个像元，代表着一平方公里的地面范围；配准信息则告诉你这个数组左上角在哪、每个像元跨多少经纬度或米。理解这个结构，后面所有操作都有方向了。
像元值这块要格外小心。不同来源的数据，像元值可能是“这个格网里的常住人口总数”，也可能是“每平方公里的人口密度”。两者只相差一个面积因子。如果数据坐标是地理坐标系（WGS84 之类），一平方公里的像元跨 0.008983 度左右；如果是投影坐标系（比如 CGCS2000 / UTM），像元就直接是 1000 米乘 1000 米。取值是总数还是密度，直接决定后续汇总算法，拿到数据第一件事就是看数据说明或者用 gdalinfo 确认。
| 像元值语义 | 常见表示 | 求和方式 | 典型来源 |
|---|---|---|---|
| 总数 | 每格网总人口 | 直接求和 | 部分成品数据 |
| 密度 | 每平方公里人数 | 乘像元面积后求和 | 国际主流 1km 格网产品 |
2.2 动手前先 gdalinfo：给 GeoTIFF 拍个 X 光
我拿到一份人口密度公里格网栅格，从来不急着写代码，先跑一条最简单的命令：
gdalinfo pop_1km.tif输出里的关键字段有几个：Driver 是格式类型，通常 GTiff；Size 是列数乘行数，对应影像范围；Coordinate System at the origin 这一段是投影定义；Origin 是左上角坐标；Pixel Size 是两个数——第一个是东西方向分辨率，第二个是南北方向分辨率，注意它通常是负数，表示 y 轴从上往下递增。还有 Band 1 里的 NoData Value——无效区域怎么编码，是 0 还是 -9999，这个直接决定后面统计时要不要掩膜。
这段输出的信息量很大。举例来说，如果 Pixel Size 显示 (0.0083333, -0.0083333) 而 Origin 是 (73.5, 53.5)，那这份数据是 WGS84 经纬度下的近似公里格网——0.0083333 度在地球中纬度大约就是 700 到 900 米，跟 1 公里有差距，但产品通常仍按公里格网命名。如果 Pixel Size 显示的是 (1000, -1000)，那就是标准投影坐标系的公里格网，栅格本身就已经摊平到平面上了。两种坐标系下裁剪、面积计算、重采样的处理方式不一样，先认清这一点，能少踩很多坑。
2.3 用 GDAL 把栅格读进 NumPy 数组
确认了头信息后，读取这一步用 Python 的 GDAL 绑定就很顺。传统做法是直接用 osgeo.gdal，新版环境我更习惯用 rasterio——API 更简洁，而且和矢量库 geopandas 配合得很自然。下面这段是读取并同时取出坐标信息的代码：
from osgeo import gdal import numpy as np ds = gdal.Open('pop_1km.tif', gdal.GA_ReadOnly) band = ds.GetRasterBand(1) nodata = band.GetNoDataValue() data = band.ReadAsArray() # 二维NumPy数组，形状是(rows, cols) gt = ds.GetGeoTransform() # gt = (x左上角, x方向分辨率, x旋转项, y左上角, y旋转项, y方向分辨率) west = gt[0] north = gt[3] pixel_width = abs(gt[1]) pixel_height = abs(gt[5]) east = west + pixel_width * ds.RasterXSize south = north - pixel_height * ds.RasterYSize读取 RasterBand 这一步，GDAL 返回的 data 形状是 rows 行、cols 列，数值类型由原数据决定——常见的是 Int16 或 Float32，人口密度常带小数，所以一般是 Float32。GetGeoTransform 返回的六个参数里，索引 0 和 3 是左上角坐标，索引 1 和 5 是分辨率，索引 2 和 4 通常是 0，代表没有旋转。用 west、east、south、north 把范围算出来，后续和矢量边界做对比时能快速判断投影是否吻合。
读进来之后建议马上做两个校验：一是打印 data.dtype，确认精度没有丢失；二是用 np.unique(data) 看看有没有奇怪的极值和 NoData 编码。我一般会再顺手检查一下 data.min() 和 data.max()，如果最小值是 -9999 而 NoData Value 又恰好没设置，后面统计就会把这个值当真实人口算进去，总量直接失真。这一步相当于数据体检，前后不到半分钟，能避免后续所有分析全错的局面。
3. 从读取到出结果：裁剪、提取与统计的实战操作
3.1 按行政区裁剪：用 rasterio.mask 一步到位
最常见的需求就是“我要某座城市、某个区的总人口”，这时候第一步是矢量边界裁剪栅格。Geopandas 读进 shp，rasterio 做 mask，一气呵成：
import rasterio from rasterio.mask import mask import geopandas as gpd input_tif = 'pop_1km.tif' city_shp = 'city_boundary.shp' output_tif = 'city_pop.tif' with rasterio.open(input_tif) as src: # 关键：矢量必须和栅格坐标系一致，不一致时先转换 shp = gpd.read_file(city_shp) if shp.crs != src.crs: shp = shp.to_crs(src.crs) out_image, out_transform = mask( src, shp.geometry, crop=True, filled=True, nodata=src.nodata if src.nodata is not None else 0 ) out_meta = src.meta.copy() out_meta.update({ "height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform }) with rasterio.open(output_tif, 'w', **out_meta) as dst: dst.write(out_image)mask 函数接受矢量几何对象作为第二个参数，src 是已打开的栅格数据集，crop=True 表示把输出范围裁剪到矢量的外接矩形，filled=True 会把边界外区域填充成 nodata，这样后续统计时有效像元和无效像元分得很清楚。这里的踩点在于 shp.crs 和 src.crs 不一致时程序不会报错，只会输出一张范围对不上的图，所以一定要先做 to_crs。src.nodata 如果为空，要显式指定一个值，不然 mask 输出里没有无效区标记。
裁剪完，建议用 gdalinfo 或 rasterio 重新打开输出文件，检查尺寸和范围的合理性。比如原栅格是几千列乘几千行，裁到单个城市后行列数应该大幅缩小；边界范围应该落在行政区内。如果发现宽高和原图几乎一样，八成是矢量几何没被正确匹配，回查坐标系和几何有效性。
3.2 用密度阈值提取高密度建成区
很多场景不需要看全部格网，只需要把“人多的地方”拎出来。城市蔓延研究里常用 5000 人/km² 作为城市核心区的经验阈值。从密度栅格里提取高密度区，本质上就是一次条件重分类：
import numpy as np # 原数据里无效值编码为 -9999，先确认再替换 density = np.where(data == -9999, 0, data) threshold = 5000 # 人/km² dense_zone = np.where(density >= threshold, 1, 0) from osgeo import gdal driver = gdal.GetDriverByName('GTiff') rows, cols = dense_zone.shape out_ds = driver.Create('dense_zone.tif', cols, rows, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band = out_ds.GetRasterBand(1) out_band.WriteArray(dense_zone) out_band.SetNoDataValue(-1) out_ds.FlushCache()阈值不是拍脑袋定的。一线做城市规划时，5000 人/km² 常用来识别城市核心建成区，小于这个数就属于一般建成区或城乡结合部。但如果对象是乡镇街道，这个阈值可能太高——农村乡镇中心常年在 2000 到 3000 人/km² 之间。我一般先跑一个分位数统计：np.percentile(density, [50, 75, 90, 99])，看看数据集里的数值分布，再用分位数反推阈值，比直接套经验值稳得多。
这段代码里驱动创建输出栅格时，我显式指定了 GDT_Byte（8 位整型），因为分类结果只有 0 和 1，不需要 32 位浮点的精度。GeoTransform 和 Projection 直接复用输入栅格的，保证空间参考不丢。NoData 设为 -1 而不是 0，是为了区分“格网外区域”（-1）和“非高密度区”（0），这个区分在后续做面积统计时非常关键。
3.3 算区域总人口：先把单位与 NoData 对齐
统计任意区域内总人口，逻辑很简单：把该区域内所有格网的人口值加起来。麻烦在于两个地方：一是像元值到底是密度还是总数，二是 NoData 的存在。先看密度版本的统计代码：
# 像元值为密度（人/km²），像元面积为 1 km² # 直接对有效像元求和即可 def summarize_population(data, nodata, pixel_area_km2=1.0): if nodata is None: valid = np.ones_like(data, dtype=bool) else: valid = data != nodata if data[valid].size == 0: return 0, 0, 0 total_pop = (data[valid] * pixel_area_km2).sum() total_area = valid.sum() * pixel_area_km2 mean_density = total_pop / total_area return total_pop, total_area, mean_density # 使用：先打开裁剪后的tif读取数据 total, area_km2, avg_density = summarize_population(out_image[0], src.nodata) print(f"总人口约 {total:.0f} 人，面积约 {area_km2:.0f} km²，平均密度 {avg_density:.0f} 人/km²")如果像元值是密度（人/km²），每个像元面积又是 1 km²，那单个像元的人口贡献就是密度乘以 1，所以求和前必须乘上 pixel_area_km2。如果数据是标准投影坐标系且像元确实是 1000 米乘 1000 米，pixel_area_km2 就固定是 1.0；如果数据是经纬度的 0.0083333 度格网，面积不是恒定值，在中高纬度要按纬度做余弦校正，否则到北方城市误差会扩大。这个细节是很多人漏掉的。
另一个常见情况是像元值已经是“该格网内总人数”。这时不能乘面积，直接求和就是总人口。判断方式很简单：把全图总值和当地官方发布的统计数据对一下，如果量级差 10 倍左右，大概率就是面积因子搞反了。我用 rasterio 处理完数据后，通常会把结果和统计口径做一次对照落位，对不上就回头检查单位与 NoData 这一步。
4. 栅格数据四大经典翻车现场：避坑排查手册
4.1 现象一：裁剪出来的结果全是空值
症状：用 shp 边界做 mask 裁剪，输出 tif 打开后要么全黑，要么只有零星几个像元有值。
原因：栅格和矢量坐标系不一致。最常见的是栅格是 WGS84 经纬度坐标，矢量是 CGCS2000 高斯投影坐标或者 UTM 坐标，rasterio.mask 按相同数值范围的坐标匹配，二者 X、Y 量纲完全不同，裁剪结果自然对不上。另一个隐蔽原因是 shp 自身坐标系定义错乱，比如某个图层的 prj 文件丢了或写错，geopandas 读进来显示的是一个错误 CRS。
解决：在 mask 前强制做一次坐标对齐。代码里加上 if shp.crs != src.crs: shp = shp.to_crs(src.crs)，并打印转换后的边界范围与栅格 extent 做人工核对。如果 shp.crs 是 None，优先看 prj 文件修复，不要赌它和栅格一定一致。这其实不是玄学，就是坐标系对齐流程没走严格。
4.2 现象二：算出来的总人口比官方数据少了三分之一
症状：用 valid = data != nodata 做掩膜后求和，结果显著低于官方统计数据，有些格网明明在建成区却是空值。
原因：NoData 编码不一致。同一份数据里可能有多种无效方式：比较规范的数据用 -9999 表示陆地外区域，但海域和湖泊有时被编码为 0，个别数据源还会用 -1。只用一种 NoData 值去掩膜，其他无效值会被当成真实人口加进统计结果里。另外，如果数据说明写着“无人口区域为 0”，那 0 确实应该保留，而不是剔除——剔除 0 会把大量无人居住的真实格网也丢掉。
解决：拿到数据后先看波段统计信息。用 np.unique(data, return_counts=True) 把唯一值和频次打印出来，确认哪些值是无效标记、哪些值是合法 0，再按语义设置掩膜。我一般会把“水域/境外区域”和“真实无人区”分开处理：前者设 NoData，后者保留为 0。区分方法就是交叉对比土地利用数据，或者直接看数值分布直方图里是否有孤立尖峰。
4.3 现象三：重采样到 250m 后总人口缩水
症状：把 1km 栅格用双线性插值重采样到 250m，叠加统计后总人口少了百分之三十多。
原因：重采样方法选错。人口密度是空间统计量，双线性或三次卷积会把高密度格网和周边低密度格网做平滑，边缘信息丢失，结果就是高值区被压低、低值区被抬升，总人口必然缩水。这跟高程 DEM 重采样的思路不一样——DEM 允许平滑，人口数据平滑后就是失真。
解决：如果源数据是密度图，重采样前先乘以像元面积换算成总量图，重采样用 sum 聚合，再除以新面积恢复密度。顺序不能反。用 rasterio 实现如下：
import rasterio from rasterio.enums import Resampling from rasterio.warp import reproject, calculate_default_transform # 目标分辨率：250米，坐标系沿用源数据的投影 with rasterio.open('pop_1km.tif') as src: transform, width, height = calculate_default_transform( src.crs, src.crs, src.width, src.height, *src.bounds, resolution=250 ) dst_kwargs = src.meta.copy() dst_kwargs.update({ 'transform': transform, 'width': width, 'height': height }) with rasterio.open('pop_250m.tif', 'w', **dst_kwargs) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=src.crs, resampling=Resampling.sum )这段代码的核心是 Resampling.sum，它表示输出像元值是落在该像元范围内的所有输入像元值之和。对“总量图”来说，sum 重采样才是保面积的正确做法。用 nearest 则只是把旧值直接塞过去，面积变了但值没变，密度图可以这么干，总量图不能。反过来，如果直接对密度图用 sum，每个 250m 格网的密度值会被缩小大约 16 倍，因为贡献面积只有原来的四分之一。
4.4 现象四：高密度区提取结果锯齿严重
症状：阈值提取出来的 dense_zone.tif，边界跟测试区域边缘参差不齐，放大看全是锯齿。
原因：公里格网数据的天然颗粒度。1km 格网本身就很粗，从栅格出来的等高线式边界必然带锯齿。更隐蔽的坑是，提取前没有做 NoData 掩膜，导致大块无人区里孤立的异常高值被当成高密度区。
解决：两种策略。其一是提取后做形态学开运算，去掉小噪点和断开区，但注意开运算会同时削减真实高密度区的边缘。其二是从矢量出发，先栅格化矢量边界再做提取，保证行政区边界外不会有像元残存。形态学处理用 scipy.ndimage 的 binary_opening，结构元用 3x3 即可：
from scipy import ndimage dense_clean = ndimage.binary_opening( dense_zone, structure=np.ones((3, 3)) ).astype(np.uint8)binary_opening 先腐蚀再膨胀，能把孤立像素点和小碎片清掉。结构元尺寸越大，清除的对象就越大，但在 1km 格网上用 5x5 会把真实高密度区的边缘吃掉一圈——这一圈少则几个平方公里，多则几十个平方公里，对后续面积统计影响相当大。我一般只用 3x3 做一遍，然后人工对比原密度图和提取结果，确认边缘位置基本吻合。
5. 出图与进阶用法：从密度栅格到能汇报的成果
5.1 用对数色标快速可视化密度分布
密度分布天然偏态：绝大多数格网在 0 到 200 人/km² 之间，少数核心区上万。用线性色标出图，整张图只有市中心几点亮。我一般用 LogNorm 做色标，同时把 0 值映射为透明，这样底图能透出来：
import matplotlib.pyplot as plt from matplotlib.colors import LogNorm import numpy as np fig, ax = plt.subplots(figsize=(12, 8)) masked = np.ma.masked_where(data <= 0, data) valid_max = data[data > 0].max() norm = LogNorm(vmin=1, vmax=valid_max) im = ax.imshow(masked, cmap='YlOrRd', norm=norm, interpolation='nearest') plt.colorbar(im, ax=ax, label='人/km²', extend='max') ax.set_title('人口密度分布（对数色标）') ax.set_axis_off() plt.tight_layout() plt.savefig('pop_density_map.png', dpi=300, bbox_inches='tight')LogNorm 把数值压缩到对数空间，色带区分度就出来了。imshow 里 interpolation 设成 nearest 是为了避免 matplotlib 在低分辨率格网上做插值平滑，否则会看到模糊糊的一片。加一个 extend='max' 让色条顶部的超范围值显示成箭头，比 clamp 后颜色失真更直观。如果希望输出有地理坐标背景，可以叠加 geopandas 读入的县界边界线，但注意要先 to_crs 统一坐标系再用 transform 对应。
5.2 密度栅格与其他数据的关联分析
进阶用法是把人口密度图和其他空间变量做相关性分析。比如评估“餐饮 POI 密度和常住人口密度是否同步”：
import geopandas as gpd import numpy as np import rasterio from rasterio.features import rasterize from rasterio.enums import MergeAlg import scipy.stats as stats with rasterio.open('pop_1km.tif') as src: data = src.read(1) transform = src.transform out_shape = (src.height, src.width) crs = src.crs poi = gpd.read_file('poi_points.shp') if poi.crs != crs: poi = poi.to_crs(crs) counts = rasterize( [(geom, 1) for geom in poi.geometry], out_shape=out_shape, transform=transform, fill=0, merge_alg=MergeAlg.add ) valid = (data > 0) & (counts > 0) log_density = np.log10(data[valid]) log_poi = np.log10(counts[valid]) r, p = stats.pearsonr(log_density, log_poi) print(f"Pearson r = {r:.3f}, p = {p:.3g}")rasterize 的 out_shape 和 transform 必须与密度栅格完全一致，这样才能保证每个格网对齐。merge_alg=MergeAlg.add 表示对落在同一格网内的多个 POI 点做累加，而不是直接覆盖；fill=0 保证没有 POI 的格网值为 0。取对数后再算相关，是因为人口和 POI 的分布都极度右偏，直接算 Pearson 相关系数容易被异常高值带偏，对数变换能压制偏态对线性相关的干扰。
做双变量制图时，把人口密度分五档，把夜间灯光或设施密度也分五档，交叠色块生成网格图，是给规划部门汇报时很受欢迎的一种形式。这个技巧不复杂，但需要确保两套栅格的坐标系和 NoData 处理完全一致，否则横竖对不上，出来的图就没法解释了。
其实无论做哪种进阶分析，我的第一步永远是重新检查投影与单位的文档记录。这份人口密度公里格网栅格数据的空间精度是 1 公里，适合宏观尺度判断，不适合精确到街区内部的微观问题。把数据能力用到边界附近时，宁可花五分钟做一次 gdalinfo，也不要等到结果对不上官方数据再返工。想起前几年有一次我在某个项目里把像元值直接当总数做了统计，结果和官方数据差了一个数量级，后来排查才知道是密度不是总数。从那以后，我拿到任何栅格的第一件事就是把 NoData 值、像元值语义、坐标系三个字段写进项目说明文档，无论多急都先走一遍这三步。希望帮到你。
本文还有配套的精品资源，点击获取
