GPW v4人口密度栅格数据分析：Python实现乡镇级裁剪与人口统计

简介北京市乡镇级2000-2020年人口密度栅格数据基于约4万个行政单位生成空间分辨率约1km坐标系为WGS84人口计数已调整并匹配联合国国家总数修订版可直接用于GIS制图、城乡规划与人口分布研究。数据覆盖2000、2005、2010、2015、2020五个年份适合从宏观区域到乡镇尺度的连续人口密度分析尤其匹配国内以乡镇为最小统计单元的普查公布方式。资源包共16个文件包含5个TIFF栅格数据、5个TFW世界文件、5个XML元数据及1个PNG概览图压缩包仅360KB轻量易用。已有470人学习/下载配套TFW与XML文件便于地理配准与数据溯源可直接叠加矢量边界提取任意区域人口密度适配ArcGIS、QGIS等常用平台。适合地理信息相关从业者、科研人员及高校师生作为研究底图或教学案例使用。1. GPW v4人口栅格与北京市乡镇级分析的价值拿到这份“北京市乡镇级2000-2020年人口密度.rar”时我先确认了一下文件清单里面是五个年份的栅格数据文件名带着gpw_v4_density前缀从2000到2020每五年一条。这不是一份矢量面数据而是世界人口网格化项目 GPW v4 的北京裁剪结果单位是“每平方公里人数”空间分辨率约 1 公里。对做区域研究、城市规划或数据分析的人来说它的价值在于乡镇级行政单元往往找不到长时间序列的公开人口密度而这套栅格可以从 2000 年一直铺到 2020 年统一了口径还能回溯历史。缺点也很明确1 公里网格对面积较小的乡镇来说不够细边缘误差需要自己校正读到坑的时候你会有感觉。2. 数据文件结构与坐标系统先看清 tfw 和 aux.xml 再动手2.1 文件清单与命名规律解压后你会看到类似下面这样的文件列表gpw_v4_density_2000年_北京市_北京市.tif gpw_v4_density_2005年_北京市_北京市.tif gpw_v4_density_2010年_北京市_北京市.tif gpw_v4_density_2015年_北京市_北京市.tif gpw_v4_density_2020年_北京市_北京市.tif每个 tif 旁边还有同名的.tfw和.tif.aux.xml。tfw 是 ESRI World File记录栅格像元的地图坐标变换参数aux.xml 是 GDAL 在读取或生成 GeoTIFF 时自动写的辅助元数据。很多 GIS 新手只盯着 tif 看忽略这两个文件结果在老软件里打开后位置完全偏移其实问题就出在这里。2.1.1 用 gdalinfo 快速验证坐标系打开终端对 2020 年那个 tif 执行gdalinfo gpw_v4_density_2020年_北京市_北京市.tif重点看这几行输出Driver: GTiff/GeoTIFF Size is 54, 44 Coordinate System is: GEOGCRS[WGS 84, ...] Origin (115.5,41.5) Pixel Size (0.008333333333333, -0.008333333333333)说明四件事。第一坐标系是 WGS84不是投影坐标系。第二像元大小是 0.0083333 度即 1/120 度等于 30 弧秒。第三图像只有 54×44 个像元正好覆盖北京范围。第四Origin 左上角坐标在 115.5°E、41.5°N往下扫。这个横向尺寸意味着在北京纬度约 40°N附近每个像元的真正地表宽度大约是0.0083333 * 111.32 * cos(40°)算下来约 0.71 公里而南北方向约 0.93 公里所以它不是标准正方形网格计算面积时别直接用 1 公里×1 公里。2.2 tfw 文件的六个参数如果你用文本编辑器打开.tfw会看到六行数字。以 2020 年文件为例常见内容如下0.008333333333333 0.000000000000000 0.000000000000000 -0.008333333333333 115.500000000000000 41.500000000000000这六个数依次是X 方向像元大小、旋转项、旋转项、Y 方向像元大小负值表示从上往下排列、左上角 X 坐标、左上角 Y 坐标。从这里可以直接算出任意像元中心坐标第 i 行第 j 列的中心 X 第 5 个数字 j × 第 1 个数字 i × 第 3 个数字中心 Y 第 6 个数字 j × 第 2 个数字 i × 第 4 个数字。旋转项都是 0说明没有旋转是标准的北向上栅格。如果某个软件不识别 GeoTIFF 内部的坐标就会自动读 tfw 来完成地理配准理解这一点可以帮你排查坐标偏移问题。2.3 相邻年份文件命名差异与版本一致性资源包里文件名带有中文“年”并且年份顺序不是递增排列而是 2000、2005、2010、2015、2020 五个年份各一份。这类批量处理时最容易出问题的是文件通配符。如果你在 Linux 下用glob.glob(*.tif)返回顺序可能按拼音或字节排序导致年份错乱。我的习惯是先用正则把年份提取出来再按年份显示排列import re tifs [gpw_v4_density_2000年_北京市_北京市.tif, gpw_v4_density_2005年_北京市_北京市.tif, ...] tifs_sorted sorted(tifs, keylambda x: int(re.search(r(\d{4})年, x).group(1)))这样可以避免后续年份匹配错位。GPW v4 的数据版本本身是统一的所有年份基于同一套行政边界和人口估计方法所以做时间序列是安全的。3. 用 Python 读取人口密度栅格并裁剪到乡镇边界3.1 环境准备与依赖处理这套数据建议用 Python 3.9 以上安装 rasterio、geopandas、shapely、pyproj。如果你已经有 conda可以直接建一个新环境conda create -n gpw-analysis python3.10 conda activate gpw-analysis pip install rasterio geopandas pyprojrasterio 负责栅格读取和裁剪geopandas 负责矢量边界操作pyproj 负责坐标转换。注意不要用 PIL 读 tif因为 PIL 不会读取地理坐标信息直接把空间参考丢了。3.2 用 rasterio.mask 裁剪单个乡镇假设你有一份北京市乡镇边界线数据beijing_towns.shp里面包含乡镇名称字段name。下面这段代码把某个乡镇范围裁剪出来得到该乡镇的人口密度栅格import geopandas as gpd import rasterio from rasterio.mask import mask # 读取乡镇边界 towns gpd.read_file(beijing_towns.shp) town towns[towns[name] 怀柔镇].geometry.iloc[0] # 打开 2020 年人口密度栅格 src_path gpw_v4_density_2020年_北京市_北京市.tif with rasterio.open(src_path) as src: out_image, out_transform mask(src, [town], cropTrue, all_touchedTrue) out_meta src.meta.copy() out_meta.update({ driver: GTiff, height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(town_density_2020.tif, w, **out_meta) as dst: dst.write(out_image)这段代码核心是mask(src, [town], cropTrue, all_touchedTrue)。第一个参数[town]必须是几何对象列表cropTrue表示把范围裁剪到该几何的外接矩形all_touchedTrue表示凡是被边界碰到的像元都保留。对低分辨率栅格做乡镇级裁剪时all_touchedTrue非常重要因为它能减少边缘像元被排除的概率代价是边缘会多出几个像元。裁剪后的 out_image 是三维数组第一维是波段数这里是 1第二三维是行列。out_transform是根据裁剪窗口重新计算的仿射变换写入 GeoTIFF 时必须用它覆盖原元数据否则新文件的位置还是错的。3.3 批量裁剪全部乡镇单乡镇裁剪可以跑通但实际分析要处理几十上百个乡镇用循环批量处理import os import numpy as np output_dir towns_2000_2020 os.makedirs(output_dir, exist_okTrue) years [2000, 2005, 2010, 2015, 2020] towns gpd.read_file(beijing_towns.shp) for year in years: tif fgpw_v4_density_{year}年_北京市_北京市.tif with rasterio.open(tif) as src: for _, t in towns.iterrows(): out_path f{output_dir}/{t[name]}_{year}.tif out_image, out_transform mask(src, [t.geometry], cropTrue, all_touchedTrue) out_meta src.meta.copy() out_meta.update({height: out_image.shape[1], width: out_image.shape[2], transform: out_transform}) with rasterio.open(out_path, w, **out_meta) as dst: dst.write(out_image)这里注意嵌套循环每打开一个年份的 tif遍历所有乡镇裁剪。如果乡镇数量上百、年份有五个一会儿就能跑完因为原始栅格本身很小。如果换到全国数据这种做法就太慢了需要改用分块或先栅格化矢量再分区统计。3.4 裁剪失败时的检查点裁剪后最好立刻检查结果别急着算人口。常见问题是输出栅格全是 0 或缺失。原因有三类乡镇边界坐标是投影坐标系比如 CGCS2000 / 高斯-克吕格而栅格是 WGS84二者空间参考不一致。先执行towns towns.to_crs(src.crs)把乡镇边界转到 WGS84。乡镇边界离栅格范围太远没有相交此时 mask 返回的 out_image 全是 nodata 或面积为 0。用town.bounds和src.bounds对比一下范围。all_touchedFalse导致狭窄乡镇被漏掉栅格像元中心全不在多边形内。这种情况把all_touched改为 True 即可。提示北京市面积约 1.6 万平方公里GPW 栅格只有 54×44 像元裁剪后每个乡镇可能只有几个到几十个像元完全不构成“高清”。所以裁剪只是第一步真正的功夫在统计指标的计算和误差修正上。4. 乡镇级人口统计与时间序列对比从栅格到表格4.1 密度栅格不等于人口数第一次处理这种数据的人最容易犯的错误直接把像元值和面积相乘。GPW v4 的 density 栅格每个像元的值是“每平方公里的人数”不是该像元内的总人口。要得到某个乡镇的总人口需要把每个有效像元的密度值乘以该像元在实地对应的面积单位平方公里再求和。因为 WGS84 的像元并不是规则正方形你需要先做投影转换或者按纬度逐像元计算面积。最稳妥的做法把裁剪好的乡镇密度栅格重投影到等积投影比如针对北京地区使用 WGS84 / Albers 投影EPSG: 9844 或者自定义 Albers 中央经线 105°E、标准纬线 25°N 和 47°N。等积投影下像元面积是常数总人口 密度值 × 像元面积 × 像元数量的加权和。4.2 用 rasterstats 快速统计乡镇人口另一种更直接的方式是不手工裁剪每个乡镇而是用 zonalstats 对全部乡镇做区域统计。rasterio 搭配rasterstats库能省很多事pip install rasterstats统计五年所有乡镇的平均密度和总人口代码如下import geopandas as gpd import rasterstats as rs import pandas as pd towns gpd.read_file(beijing_towns.shp) years [2000, 2005, 2010, 2015, 2020] rows [] for year in years: tif fgpw_v4_density_{year}年_北京市_北京市.tif # 注意rasterstats 默认读取栅格要求矢量与栅格 CRS 一致 town_stats rs.zonnal_stats( towns, tif, stats[mean, sum, count], all_touchedTrue, nodata-9999 ) # sum 在密度栅格里没有意义必须加权计算 # 但 rasterstats 不能直接按面积加权所以这里临时用 mean 再乘乡镇面积这里我写了一半就发现问题stats[sum]是对像元密度值求和不是人口。正确做法是先把密度栅格转换成人口计数栅格人口数 密度 × 像元面积平方千米然后再用 sum。或者直接对每个乡镇计算平均密度 × 乡镇实际面积。考虑到 1km 网格的问题我更推荐先把整个北京市的密度栅格重投影到等积投影再用 zonnal_stats 的 sum 来统计人口。完整流程如下import rasterio from rasterio.warp import calculate_default_transform, reproject, Resampling src_path gpw_v4_density_2020年_北京市_北京市.tif dst_crs EPSG:9844 # WGS 84 / Albers North America 示例实际北京建议自定义 Albers with rasterio.open(src_path) as src: dst_transform, dst_width, dst_height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds, resolution800) dst_meta src.meta.copy() dst_meta.update({ crs: dst_crs, transform: dst_transform, width: dst_width, height: dst_height }) with rasterio.open(beijing_density_albers.tif, w, **dst_meta) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformdst_transform, dst_crsdst_crs, resamplingResampling.bilinear)这段代码把密度栅格重投影到 Albers 等积投影分辨率设为 800 米也就是 0.8 公里。在等积投影下每个像元面积就是800 × 800平方米即 0.64 平方公里。重投影之后每个像元的值依然是“每平方公里人数”所以单个像元人口数 值 × 0.64。此时再对乡镇做区域统计统计sum就有意义了。4.3 生成乡镇级人口时间序列表等投影就绪后用 rasterstats 统计每个乡镇的总人口和平均密度并生成宽表towns_albers towns.to_crs(EPSG:9844) pixel_area_km2 0.8 * 0.8 # 0.64 km² pop_records [] for year in years: tif_albers fbeijing_density_albers_{year}.tif # 假设每年都重投影过 stats rs.zonnal_stats(towns_albers, tif_albers, stats[sum, mean, count], all_touchedTrue) for town, s in zip(towns_albers[name], stats): pop_records.append({town: town, year: year, total_pop: s[sum] * PIXEL_AREA_KM2, mean_density: s[mean]}) df pd.DataFrame(pop_records) pivot df.pivot(indextown, columnsyear, valuestotal_pop) pivot.to_csv(beijing_town_population_2000_2020.csv)这里s[sum]是所有有效像元的密度值之和乘以每个像元对应的 0.64 平方公里后就是该乡镇总人口。高斯投影下的 Albers 投影面积误差很小对于乡镇级统计足够。mean_density也可以直接用但要注意它是按像元数量平均而不是按乡镇面积加权平均所以如果栅格边缘像元覆盖不完全会和年鉴值有偏差。4.4 对比同一乡镇五年变化把 pivot 表读出来后可以画简单的时间序列。用 matplotlib 画折线import matplotlib.pyplot as plt pivot.loc[怀柔镇].plot(markero) plt.ylabel(Total population) plt.title(Population change of Huairou Town (2000-2020)) plt.grid(alpha0.3) plt.show()五个年份的点即可看出乡镇人口增长趋势。但注意这五个年份不是等间隔同步普查数据2000、2005、2010、2015、2020 是 GPW v4 的估计年份中间年份没有。做趋势分析时把它当作“每五年估计值”看待不要解读成精确普查数据。5. 分辨率限制与边缘校正乡镇级应用的两个关键坑1 公里网格对北京市乡镇意味着什么北京市城六区乡镇面积平均约 2050 平方公里映射到 GPW 栅格上可能只有 2080 个像元。像元内人口被视为均匀分布这完全掩盖了人口沿街道、河流、地铁线集聚的事实。更麻烦的是乡镇边界切割栅格时边缘像元可能只覆盖了乡镇的一小部分如果用all_touchedTrue会把整个像元面积算进乡镇导致人口高估如果用all_touchedFalse又会漏掉真正的边缘人口。我的做法用all_touchedTrue保证人口总数不漏然后对边缘像元做面积加权修正——即对每个被边界切割的像元用实际与乡镇多边形相交的面积比例乘上该像元的人口数。计算面积加权需要做一步栅格化。用 rasterio 的rasterize把乡镇边界生成一个 0/1 掩膜掩膜像元大小与重投影后的密度栅格完全一致再与密度栅格相乘就可以得到每个像元实际被乡镇覆盖的面积比例。具体实现from rasterio.features import rasterize import numpy as np town_geom towns_albers.iloc[0].geometry # 创建一个和密度栅格尺寸相同的掩膜 mask rasterize([(town_geom, 1)], out_shape(dst_height, dst_width), transformdst_transform, fill0, all_touchedTrue) # 读取密度数组 with rasterio.open(tif_albers) as src: density src.read(1) # 有效像元掩膜nodata 处理 valid density 0 # 边缘像元比例估算这里用掩膜值作为面积比例0或1如果要更精确可使用 antialiasingTrue # 真实的比例计算更精细需要把像元多边形化。简单场景下先用掩膜过滤 adjusted_density density * mask total_pop adjusted_density.sum() * pixel_area_km2这里mask只是 0/1不是精确比例。更精确的做法是把每个像元转为 Polygon然后与乡镇求交面积比但计算成本高。对于大多数分析我建议用 0/1 掩膜再结合乡镇官方人口做比例校正调一个系数使得 2020 年统计值与第七次人口普查的多镇人口对上然后用同一系数修正其他年份。这样至少能消除系统性的面积计算偏差。另一个坑是行政区域边界变化。2000 年到 2020 年北京市乡镇撤并、改名频繁。如果你用的乡镇边界是基于 2020 年的那 2000 年的人口统计对应的是 2020 年的边界两地差异会混入“人口变化”里。解决办法是把边界统一到同一套乡镇编码或者在分析时合并小乡镇到更大尺度如区级比较。GPW v4 本身使用的是 2010 年左右的行政边界如果你用 2020 年乡镇边界裁剪它要注意边界不匹配会导致边缘误差。最后给你一个验证技巧把统计出的乡镇总人口除以北京统计年鉴公布的常住人口看比值是否在 0.91.1 之间。如果偏差超过 20%优先检查 nodata 处理——GPW v4 的密度值有时会出现 0 值表示无人区但偶尔水域是负值或 255 这种特殊标记。用numpy.unique(density)检查数值范围过滤掉所有不合理的像元再统计。用这套思路你不仅能用这份资源复现乡镇级 2000-2020 年人口密度时间序列还能把误差控制到可解释范围内。真正要动工时先从 2020 年一个区开始做通再扩展到全市不要一上来就全量跑。本文还有配套的精品资源点击获取