简介北京市乡镇级2000—2020年人口密度栅格数据包面向地理信息、城乡规划、人口与公共健康等领域的科研人员与在校学生。数据源自全球人口栅格模型人口计数已调整至与联合国人口估计一致并利用约4万个行政单位校准再分配到30弧秒网格形成精度约1km、WGS84坐标系的人口密度栅格可支撑乡镇尺度的人口时空演变分析、公共服务设施配置及灾害风险评估等场景。压缩包共16个文件包含2000、2005、2010、2015、2020五个年份的tif栅格主文件配套tfw坐标配准文件与xml元数据文件另有说明图片整体体积约360KB轻量易用。目前已有470人学习下载。数据取自GPWv4人口密度产品并裁剪出北京市域范围可直接在ArcGIS、QGIS中加载制图免去自行下载、裁剪、统一投影的步骤适合需要连续年份人口密度底图的研究任务。1. 北京市乡镇级2000-2020年人口密度.rar先看清这个包再动手拿到“北京市乡镇级2000-2020年人口密度.rar”这种文件第一反应别是找解压密码而是先想清楚里面大概率装了什么三个年份2000、2010、2020的乡镇街道行政区划矢量面每个面挂一条人口属性人口除以面积就是密度。乡镇级粒度比区县级细一个量级能暴露城市内部空间结构所以这个包常见于城市规划、公共卫生、商业选址和学术研究。但它不是开箱即用的干净数据我见过太多人卡在shapefile中文乱码、乡镇边界对不上、面积算出来离谱这三类事上。下面按我自己的处理顺序从解包一路讲到出图。2. 解包读取从RAR压缩包到GeoDataFrame的完整链路2.1 先列目录确认这份RAR里到底是什么格式拿到RAR后我会先列压缩包内容而不是直接解压。用unrar的list模式看内层文件unrar l 北京市乡镇级2000-2020年人口密度.rar如果系统没有unrar用归档管理器也可以。列出来的内容通常分为几类2000/、2010/、2020/三个年份目录每个目录下的.shp、.shx、.dbf、.prj—— 这是ESRI Shapefile的标准四件套偶尔有.xlsx或.csv说明附带统计表。看到Shapefile意味着接下来的读取工具是GDAL/geopandas。如果内层是GeoJSON或gpkg代码也差不多只是不需要处理编码。.dbf是属性表人口数字就存在这里.shx是几何索引.prj虽然小但别删里面写着坐标系。确认完结构再解压能避免解出来一堆不知道是什么格式的文件。2.2 用unrar或7z命令行解出数据确认格式后解压到一个独立的目录。两个主流程任选# 方案Aunrar直接解压 mkdir -p beijing_pop unrar x 北京市乡镇级2000-2020年人口密度.rar beijing_pop/ # 方案B7z同样能解RAR mkdir -p beijing_pop 7z x -obeijing_pop 北京市乡镇级2000-2020年人口密度.rar这段代码里x代表完整解压并保留内层目录结构-o指定输出目录。unrar x会自动创建目录所以先mkdir -p只是个保险动作。两个工具在主流Linux发行版上都能用apt install unrar或apt install p7zip-full装上。注意RAR的解压实现是有专利的部分发行版只提供unarlibarchive遇到解压失败时先换工具再怀疑数据损坏。Windows下没有这两个命令时用7-Zip图形界面选“解压到当前文件夹”效果一样只是不方便写进批处理。2.3 在Python脚本里用rarfile解压当项目要在批处理里处理多个同类压缩包时命令行解压不方便这时候用rarfile这个库import os import rarfile rar_path 北京市乡镇级2000-2020年人口密度.rar out_dir beijing_pop os.makedirs(out_dir, exist_okTrue) # 指定unrar可执行文件路径Windows下必须配这一步 rarfile.UNRAR_TOOL /usr/bin/unrar with rarfile.RarFile(rar_path) as rf: for info in rf.infolist(): if not info.is_dir(): rf.extract(info, out_dir) print(解压完成文件列表) for root, dirs, files in os.walk(out_dir): for f in files: print(os.path.join(root, f))这个循环里infolist()返回压缩包内全部文件元数据is_dir()用来过滤目录项。其实直接rf.extractall(out_dir)也能达成目的我习惯保留这个循环是为了在正式解压前打印文件清单方便核对数量。rarfile只是外壳底层必须调用unrar或unar可执行文件UNRAR_TOOL就是告诉它去哪里找。如果找不到报错通常是Couldnt find unrar tool这时候不是代码问题而是缺二进制。2.4 geopandas读取Shapefile编码参数一次设对解压完成后用geopandas读入。最需要注意的是encoding参数——Shapefile的.dbf属性表没有强制编码标准国内数据源常见两种GBK和UTF-8。先按UTF-8试字段名或乡镇名字段出现“鍖椾含”这类乱码时换成GBK再读import geopandas as gpd # 2020年的数据通常较新先按UTF-8读 gdf_2020 gpd.read_file(beijing_pop/2020/town_2020.shp, encodingutf-8) # 2000/2010年数据多半是GBK编码直接指定 gdf_2000 gpd.read_file(beijing_pop/2000/town_2000.shp, encodinggbk) # 看一眼字段结构确认人口字段名 print(gdf_2020.columns.tolist()) print(gdf_2020[[name, pop, geometry]].head())encoding只影响属性表不影响几何。字段名有时是英文缩写例如POP2020、DENSITY也可能是中文如“常住人口”。如果实在猜不准用pd.read_csv去读同目录下的.dbf文件把编码猜出来或者直接print列名看规律。读进来的数据是一个GeoDataFramegeometry列是Polygon/MultiPolygonprojection信息读自同名的.prj文件读完后用gdf_2020.crs确认坐标系是不是WGS84这直接影响后面的面积计算。3. 人口密度计算面积重算与口径核对3.1 为什么不能信任Shapefile自带的面积字段很多乡镇面数据自带一个area或SHAPE_Area字段。这个字段能不能直接用取决于它当时是用什么坐标系计算的。如果原始数据的坐标系是WGS84经纬度这个面积字段的单位就是平方度毫无物理意义如果用的是Web墨卡托EPSG:3857面积会比真实值偏大在北京这个纬度大约偏大35%以上。北京乡镇面积从几平方公里到两百平方公里不等密度计算里面积在分母上差三分之一结果就完全失实。3.2 用北京市范围的Albers等积投影重算乡镇面积我一般拿到面数据的第一件事就是投影重算面积。北京市范围在经度115.5E到117.5E、纬度39.5N到41.1N之间针对这个范围定制一套Albers等积投影参数from pyproj import CRS import geopandas as gpd # 北京专用Albers等积投影双标准纬线39.5°N和41°N中央经线117°E beijing_albers CRS.from_proj4( projaea lat_139.5 lat_241 lat_00 lon_0117 x_00 y_00 datumWGS84 unitsm no_defs ) gdf_2020_albers gdf_2020.to_crs(beijing_albers) # 面积单位是平方米除以1e6得到平方公里 gdf_2020[area_km2] gdf_2020_albers.geometry.area / 1e6 # 人口密度 常住人口 / 面积人/平方公里 gdf_2020[density] gdf_2020[pop] / gdf_2020[area_km2] print(gdf_2020[[name, pop, area_km2, density]].describe())参数说明aea是Albers等面积投影它的关键特征是面积不变形适合做密度和面积统计。lat_1和lat_2是两条标准纬线把北京的主要纬度包在中间这样整个市域内的面积误差都压到极小lon_0117是中央经线这个值不会影响面积精度但决定图形在平面上的位置取区域中心更直观。unitsm特别重要不写的话有些实现默认输出英尺面积算出来直接大一圈。3.3 三次普查的人口口径不同加总前先核对2000、2010、2020三个年份对应第五次、第六次、第七次全国人口普查。普查指标主要是常住人口但有些数据包里混入了“户籍人口”或“年末总人口”字段名后面不带口径说明时很容易错用。核对办法很直接按区县汇总和普查公报对比。常见的字段命名对应关系可以先查一遍字段可能命名常见含义单位/备注pop / population / rk常住人口人huji_pop / hkrk户籍人口人与常住差异大area / SHAPE_Area前手算好的面积单位取决于投影不能盲用density / rkmd可能已被前手计算过用前先核一遍总量# 按区县代码汇总人口 district_pop gdf_2020.groupby(district_code).agg( district_pop(pop, sum), district_area(area_km2, sum) ).reset_index() # 粗略核对2020年第七次普查北京常住人口约2189万 total_pop_2020 district_pop[district_pop].sum() print(f乡镇汇总人口: {total_pop_2020 / 1e4:.0f} 万)如果汇总值和1356.9万2000年、1961.2万2010年、2189.3万2020年这几个普查总量相差超过5%先怀疑字段口径选错了再怀疑有乡镇要素缺失。差值在1%以内是正常的因为个别街道的界线在普查登记时与实际民政边界存在微差而且有些“园区管委会”“开发区”不在乡镇统计台账里。千万不能用全市总量硬平分到乡镇那个误差会被空间分布放大。4. 2000-2020年跨期对比行政区划边界变动怎么处理4.1 行政区划变动是跨期分析的最大敌人三个年份的数据放在一起最头疼的不是格式而是乡镇边界根本对不上。从2000年到2020年的二十年里北京经历了大规模“撤乡并镇”和“乡改街道”具体到乡镇一级有过多次合并、拆分、更名和代管调整。同一个2020年的乡镇2000年可能对应两三块小区域反过来一个2000年的乡镇名字到今天可能已经不存在了。直接把两个年份的shp做join会丢掉一大半记录。跨期对比必须先把边界关系讲清楚。4.2 行政区划变动类型与处理策略对照表处理之前先判断每个乡镇属于哪一类变动再选对应策略不要一把梭用同一种join方式变动类型举例场景常见处理策略乡镇合并两个乡合为一个镇将旧年份的人口加总到新年份目标单元乡镇拆分一个镇拆成新的镇和街道按面积权重拆分人口或按居村人口比例拆分乡改街道乡镇建制改为街道办事处行政代码变化边界基本不动直接映射名称变更仅改名字、代码不变按行政区划代码关联边界微调相邻乡镇间界线小幅变动当作不动或按面积比例修正实际操作里拆分的情况权重最难定。没有居村粒度数据的时候面积加权是次优但可复现的做法旧乡镇的人口密度假设在该范围内均匀分布切出去的部分按面积比例分走人口。这个假设在山区和城区都不完全成立所以报告里要写明“按面积折算”。如果有居村委会级别的重心点或人口数据优先按那个粒度做比例分配误差会小很多。4.3 用质心空间连接对齐多期边界当两份数据的乡镇代码和名称都已乱掉最可靠的办法是空间关系兜底用旧年份乡镇面的质心打到最新年份的面上落在哪个2020乡镇里就归属谁。import geopandas as gpd import pandas as pd # 读取三个年份数据并统一到WGS84地理坐标 gdf_2000 gpd.read_file(beijing_pop/2000/town_2000.shp, encodinggbk).to_crs(EPSG:4326) gdf_2020 gpd.read_file(beijing_pop/2020/town_2020.shp, encodingutf-8).to_crs(EPSG:4326) # 用representative_point取面内代表点避免质心落在边界外 gdf_2000[centroid] gdf_2000.geometry.representative_point() points_2000 gdf_2000.set_geometry(centroid) # 空间连接找每个2000年乡镇代表点落在哪个2020年乡镇里 mapped gpd.sjoin( points_2000, gdf_2020[[town_code, town_name, geometry]], howleft, predicatewithin ) # 统计没有匹配上的旧乡镇 missed mapped[mapped[town_code_right].isna()] print(f未匹配乡镇数: {len(missed)}) print(missed[[town_name, pop]].head())这里绕开了两个容易翻车的细节。一是质心改用representative_point()而不是centroid——质心可能落在凹型乡镇的边界线外一个细长弯曲的山谷乡镇很容易发生这种情况落在外面就匹配不到。“代表点”保证点一定在面内。二是sjoin的predicate参数用within而不是默认的intersects。用intersects时边界线略微重叠会产生一对多匹配密度计算会重复计人。within要求点完整落在面内多期数据如果边界不完全重合个别代表点可能挂在裂缝里所以最后要打印missed检查。4.4 把三期数据合并成“乡镇-年份”长表匹配完成后把人口字段按年份和2020年乡镇代码透视得到可供分析的长表# 为2010年重复做一次同样的空间映射 gdf_2010 gpd.read_file(beijing_pop/2010/town_2010.shp, encodinggbk).to_crs(EPSG:4326) gdf_2010[centroid] gdf_2010.geometry.representative_point() points_2010 gdf_2010.set_geometry(centroid) mapped_2010 gpd.sjoin( points_2010, gdf_2020[[town_code, town_name, geometry]], howleft, predicatewithin ) # 给两期映射结果打上年份标签再纵向拼接 mapped[year] 2000 mapped_2010[year] 2010 combined gpd.GeoDataFrame( pd.concat([mapped, mapped_2010], ignore_indexTrue) ) # 以2020年乡镇为单位透视出三列人口 pivot_pop combined.pivot_table( indextown_code_right, columnsyear, valuespop, aggfuncsum ).reset_index() pivot_pop.columns [town_code] [fpop_{int(c)} for c in pivot_pop.columns[1:]] # 与2020年乡镇面合并附上面积与密度变化 result gdf_2020.merge(pivot_pop, ontown_code, howleft)pivot_table里aggfuncsum的语义是如果一个2020乡镇包含多个2000乡镇人口全部加总对应的是“合并”型变动。如果确认发生过拆分一个旧乡镇分到两个新乡镇pivot_table无法处理同一条旧记录被复制到多行的场景那种情况要先按面积权重拆行再透视。做拆分时给每一行加上weight列调pivot_table的values权重计算量不大但逻辑要单独写清楚。5. 验证数据质量并用分级设色图输出密度图5.1 三个快速检查空值、异常密度、总量对比出图之前先做三件小事。第一密度值有没有空值通常是有乡镇人口字段为空或面积为0导致的第二密度有没有离谱的极大极小值第三三期总量与公报对比。# 空值和异常值检查 print(空值数量, result[density].isna().sum()) print(result.sort_values(density, ascendingFalse)[[town_name, density]].head(5)) print(result.sort_values(density).head(5))北京人口密度的正常数量级城市核心区东城、西城部分街道每平方公里一两万人郊区浅山乡镇几百到一千人远郊深山区可能低于100人。数值太大通常是人口错成了全区总量数值为0或NaN要查面积是否被算成0投影失败导致的退化几何很常见。查法是把那几行的geometry打印出来看是不是面积非常小的退化多边形。5.2 用Quantile分级的matplotlib出图模板地图分级配色我直接给一个可复用的matplotlib模板。这里用mapclassify的Quantiles做五分级比等间距更适合密度这类长尾分布数据import matplotlib.pyplot as plt import mapclassify fig, ax plt.subplots(figsize(8, 10)) result.plot( axax, columndensity, cmapYlOrRd, schemequantiles, k5, legendTrue, edgecolorwhite, linewidth0.3, legend_kwds{ loc: lower right, bbox_to_anchor: (1.35, 0), fmt: {:.0f} } ) # 隐藏坐标轴突出图面 ax.set_axis_off() ax.set_title(Beijing township population density (2020), fontsize14) plt.tight_layout() plt.savefig(bj_density_2020.png, dpi300, bbox_inchestight)scheme参数由mapclassify提供quantiles把数据按分位数切成五段。对于大城市内部少数高密度街道与大量低密度山区并存的情况Quantiles的优点是每级乡镇数量基本相同图面不全是浅色高密度区能被看出层次。legend_kwds里的fmt控制图例数字格式密度往往几千上万带小数点反而分散注意力。如果三张图2000/2010/2020要放在一起对比三张图必须用同一个分级边界不能各自算分位数否则图与图之间色阶不对等视觉上会误判增长。固定边界的一个简单做法是取三期密度合并后的分位数作为统一的classification_kwds传入。5.3 一个重要数据导出技巧存成GeoPackage分析做完要学会把结果存成可再次加载的格式。shp文件有2GB体积上限和字段名10字符限制dBase III规范密度计算后字段多而且中文名很容易被截断。GeoPackage没有这两条限制字段名可长可短后缀一个文件全搞定result.to_file(beijing_pop_density_2000_2020.gpkg, layertown_density, driverGPKG)下次再读这个文件时gpd.read_file(beijing_pop_density_2000_2020.gpkg, layertown_density)直接拿到包含密度和面积字段的GeoDataFrame不再依赖三个年份各自的编码问题。如果团队里有人用ArcGIS或QGIS这个文件同样能直接打开。本文还有配套的精品资源点击获取