尧图网页设计 
简介安徽省乡镇级2000至2020年人口密度栅格数据面向地理信息系统、城乡规划、人口统计及环境科学等专业人员填补了高分辨率长时序乡镇尺度人口数据的空缺可为区域研究与决策提供基础数据支撑。数据涵盖2000、2005、2010、2015、2020年五期栅格基于约4万个行政单元的人口统计并调整匹配联合国人口总数修订版采用WGS84地理坐标系分辨率约为1公里有效支撑从省到乡镇多层级的人口分布制图与时空演变分析。压缩包整体约3.68MB内含栅格图像文件可直接放入ArcGIS、QGIS等主流平台使用便于地图制图与空间统计。目前已有175人学习下载结合人口普查和年鉴资料可快速生成人口密度图、识别人口集聚与疏散区域为区域规划、公共服务设施布局及资源配置提供量化依据也可用于应急管理与社会经济研究。1. 乡镇级人口密度栅格一份能直接拆出来用的时空数据拿到这个“安徽省乡镇级2000-2020年人口密度.rar”很多人第一反应是解压看里面有几种文件但真正值钱的是数据本身的设计它把2000、2005、2010、2015、2020五期人口密度做成约1km分辨率的WGS84栅格每个像元值代表该像元对应区域的人口数或人口密度。因为是基于约4万个行政单位的人口普查数据拟合到30弧秒网格所以既能在省、市、县层面做对比也能切到乡镇边界做统计。下面就以这个数据集为对象从RAR里的栅格文件一路拆到Python/GDAL批量提取乡镇人口密度再讨论如何做多期变化分析最后给出边界误差修正的验证方法。2. 解开RAR之后认识栅格文件、WGS84坐标与30弧秒分辨率2.1 解压RAR文件时先看这三件事拿到.rar压缩包别急着拖到GIS软件里。我一般做三件事第一用支持RAR5的7-Zip或WinRAR解压解压前先看压缩包内文件如果单个栅格超过200MB后续在Python里读取时就要考虑内存映射或分块读取第二确认是否带.prj、.tfw等辅助文件如果只有裸的.tif或.bil需要从元数据描述里补坐标系第三若压缩包带有密码这类公开数据很少见直接联系来源方确认不要用网上的“rar密码移除”类破解工具避免拿到被改过的文件也可能顺手装一堆垃圾软件。提示RAR包内如果包含同名.aux.xml别删里面通常是统计信息和无效值定义读取NoData时会用到。2.2 30弧秒网格到底是多少米原文说“30弧秒网格单元”精度约1km。这里的弧秒是角度单位1度3600弧秒30弧秒1/120度。在赤道附近1度约111km所以30弧秒约0.925km。随着纬度升高经度方向的实际距离会缩短但纬度方向仍约1.11km在“约1km”范围内。对安徽省而言地处北纬29°41′34°38′一个网格的东西跨度约0.8-1.0km落在乡镇级别刚好能区分集中居民点和外围农田。很多人在ArcGIS里看到WGS84坐标系的栅格就直接用“度”去做长度或面积计算这是常见坑。栅格像元大小显示为0.0083333333度实际上横向和纵向代表的地面距离不同。要做面积统计必须用等积投影或按纬度换算否则乡镇面积会偏小人口密度会被高估。2.3 用GDAL快速查看栅格元数据拿到TIF后第一个命令往往是gdalinfogdalinfo population_anhui_2020.tif关注这样几个字段Size是宽高Origin是左上角坐标Pixel Size是像元尺寸Coordinate System是坐标系描述NoData Value是无效值。例如Size is 2400, 3600 Origin (114.870833333333337,34.662500000000001) Pixel Size (0.008333333333333,-0.008333333333333) Coordinate System is: GEOGCRS[WGS 84, ...] NoData Value -9999这里的Pixel Size正负号表示行列方向x方向为正y方向为负意味着数据按北到南存储。NoData Value不是0统计时必须先过滤掉-9999。如果输出的Size和Origin与你预期不一致说明数据是经过裁剪或重采样的需要记录下这个源信息。另外gdalinfo -stats还能输出栅格的最小/最大/平均值。如果看到最小值为负数除了NoData外还要怀疑是否存在-1这类“占位值”。有些数据源会把水域或者边界外设为-1而不是标准的NoData。读取时要把这些值一并排除。2.4 人口计数与人口密度的区别原文提到“人口计数已调整并匹配联合国国家总数修订版”又说这是“人口密度栅格”。这里要留意有些版本存的是每个像元的人口数整数有些存的是人/平方公里浮点。判断方法是打开栅格直方图人口数的像元值通常在几十到几千密度的像元值可以到数万。如果目标是要汇总得到乡镇总人口用整数人口数直接累加更自然如果只是做密度等级图则直接用人口密度。两类值不能混在一起比较否则结果会差出几个数量级。属性含义影响分辨率30弧秒约1km能支撑乡镇级不能支撑村级坐标系WGS84地理坐标系计算面积、距离前需投影单位每像元人口数或人/km²决定后续统计口径NoData通常是-9999统计时必须过滤时间2000/2005/2010/2015/2020可做时间序列变化检测3. 用PythonGDAL把安徽乡镇级人口密度“切”出来3.1 环境准备不需要ArcGIS也能跑通如果只用ArcGIS的“分区统计”工具也能完成但对五期数据、一百多个乡镇的小批量场景PythonGDAL更可控也方便放进批量处理脚本。环境准备pip install gdal numpy pandas geopandas shapely rasterstatsWindows下gdal的pip安装经常找不到dll推荐用condaconda install -c conda-forge gdal。rasterstats是专门做栅格分区统计的库能少写很多循环。装完以后先确认GDAL版本from osgeo import gdal print(gdal.__version__)如果版本号是3.x注意gdal.Warp的某些参数行为和2.x略有不同但本文的用法在两者都适用。如果你的乡镇边矢量不是WGS84坐标系比如是EPSG:2383安徽地方坐标系需要先转成WGS84否则裁剪时会出现错位或统计结果为空。转换用geopandas很简单import geopandas as gpd towns gpd.read_file(anhui_towns.shp) towns towns.to_crs(EPSG:4326)3.2 用gdal.Warp裁剪到安徽边界原始栅格通常覆盖全国甚至全球直接裁剪到安徽省能减少后续计算的内存占用。裁剪建议用省界矢量而不是简单矩形范围因为省界是不规则多边形。代码from osgeo import gdal, gdalconst in_tif rD:\data\pop_chn_2020.tif out_tif rD:\data\pop_anhui_2020.tif shp rD:\data\anhui_boundary.shp # WGS84坐标系 gdal.Warp(out_tif, in_tif, formatGTiff, cutlineDSNameshp, cropToCutlineTrue, dstNodata-9999, resampleAlggdalconst.GRA_NearestNeighbour, options[TARGET_ALIGNED_YES]) ds gdal.Open(out_tif) print(ds.RasterXSize, ds.RasterYSize)cutlineDSName指定裁剪矢量cropToCutlineTrue让输出范围与矢量范围对齐。dstNodata-9999设置裁剪后边界外区域的无效值要和原始数据一致。resampleAlg用最近邻不改变像元值语义。如果后续做多期计算还需要统一网格就要在gdal.Warp中加xRes和yRes。裁剪后建议检查一下输出范围与省界范围是否一致可以再用gdalinfo看一次。3.3 按乡镇矢量分区统计平均人口密度乡镇级行政区划一般有几百个面要素。想得到每个乡镇的人口密度通常做法是“先用乡镇面切栅格再统计面内像元”。rasterstats的zonal_stats封装了这个过程import geopandas as gpd from rasterstats import zonal_stats towns gpd.read_file(rD:\data\anhui_towns.shp) # 乡镇边界WGS84 out_tif rD:\data\pop_anhui_2020.tif stats zonal_stats(towns, out_tif, stats[sum, mean, count], nodata-9999) towns[pop_sum_2020] [s[sum] for s in stats] towns[pop_mean_2020] [s[mean] for s in stats] towns[valid_count] [s[count] for s in stats] towns.to_file(rD:\data\anhui_towns_pop_2020.geojson, driverGeoJSON)stats[sum, mean, count]分别表示像元值求和、平均值、有效像元数量。如果栅格值代表每像元人口数sum就是乡镇总人口数如果栅格值已经是密度sum就不是总人口而是密度之和需要小心。nodata-9999告诉库跳过无效值。valid_count可以用来排查统计质量如果某个乡镇的有效像元数是0说明乡镇边界和栅格没有交集多半是坐标系不匹配或乡镇边界有问题。3.4 批量处理5期并统一网格2000到2020共五期写循环即可但要注意网格对齐years [2000, 2005, 2010, 2015, 2020] for year in years: in_tif rfD:\data\pop_chn_{year}.tif out_tif rfD:\data\pop_anhui_{year}.tif gdal.Warp(out_tif, in_tif, formatGTiff, cutlineDSNameshp, cropToCutlineTrue, dstNodata-9999, xRes0.008333333333333, yRes0.008333333333333, resampleAlggdalconst.GRA_Bilinear) stats zonal_stats(towns, out_tif, stats[sum, mean], nodata-9999) towns[fpop_sum_{year}] [s[sum] for s in stats] towns[fpop_mean_{year}] [s[mean] for s in stats]这里我使用了双线性重采样GRA_Bilinear而不是最近邻目的是把不同期数据的网格统一到同一个像元网格上。双线性重采样会平滑像元值但对于人口这种连续变量影响不大。如果你做的是精度要求高的变化检测可以用GRA_Cubic但计算会更慢。提示xRes和yRes需要和原始像素大小一致否则会改变分辨率。这里写的是30弧秒。实际使用前建议打印原始数据的Origin和Pixel Size确认没有差值。4. 构建2000-2020乡镇人口密度时间序列从栅格代数到变化分级4.1 为什么不能直接拿两期栅格相减拿到五期栅格后有人直接打开栅格计算器用pop2020 - pop2000得到的“变化值”看似直观但有三个问题一是两期数据可能经过不同裁剪或网格对齐必须先在Python里做完一致性处理二是人口密度栅格存在空间自相关逐像元相减会放大噪声比如某个像元2000年是空地2010年人口正好被分配到相邻像元结果出现巨大负数三是乡镇级分析应该先把像元值聚合成乡镇值再做变化分析而不是先做差值再聚合。正确的顺序是先完成上一章的分区统计生成每个乡镇一行、包含5个年份人口数的表然后基于行数据计算变化率。4.2 用Pandas计算乡镇级人口变化率上一章已经生成GeoJSON现在读取并计算import geopandas as gpd import pandas as pd df gpd.read_file(rD:\data\anhui_towns_pop_2020.geojson) pop_cols [pop_sum_2000, pop_sum_2005, pop_sum_2010, pop_sum_2015, pop_sum_2020] # 二十年总变化率 df[total_change] (df[pop_sum_2020] - df[pop_sum_2000]) / df[pop_sum_2000] # 分阶段变化率 df[growth_2000_2005] (df[pop_sum_2005] - df[pop_sum_2000]) / df[pop_sum_2000] df[growth_2010_2015] (df[pop_sum_2015] - df[pop_sum_2010]) / df[pop_sum_2010] df[growth_2015_2020] (df[pop_sum_2020] - df[pop_sum_2015]) / df[pop_sum_2015] print(df[total_change].describe())total_change是二十年总变化率后续三个字段用来观察阶段性差异。describe()输出的四分位数能帮你判断全省是普遍增长还是两极分化。如果最大最小差异太大先排查是否有个别乡镇边界缺失或者NoData值被当成了0计算。4.3 划分增长型、收缩型、稳定型乡镇只看连续数值不够直观我一般按业务含义分三类总增长率小于-10%为收缩型大于20%为增长型中间为稳定型。为什么不直接用0做阈值因为人口普查和网格化拟合本身有误差-5%到5%大概率是噪声。def classify(row): v row[total_change] if v -0.10: return 收缩型 elif v 0.20: return 增长型 else: return 稳定型 df[type] df.apply(classify, axis1) print(df[type].value_counts())写回GeoJSONdf.to_file(rD:\data\anhui_town_change.geojson, driverGeoJSON)在QGIS里用type字段做分类渲染会立刻看到空间分布皖北平原的收缩型乡镇往往连片合肥和沿江城市的增长型乡镇像星星一样散落。如果你看到的模式完全随机那就要回头检查数据匹配是否出错。分类阈值不是固定的。如果研究区域以合肥、芜湖等城市为主20%的增长阈值会漏掉许多“低增长但绝对量大”的郊区乡镇。此时可以把阈值提高到30%或者改用自然断点法。判断标准应该是分类结果能体现出空间集聚性而不是把全省切成棋盘格。4.4 结果表的字段设计与导出为了展示最终的表结构我列出常用字段设计字段名示例值说明name某镇乡镇名称pop_sum_2000452312000年人口总数pop_sum_2020603402020年人口总数total_change0.334二十年变化率type增长型分类结果这个表可以直接输出给业务方df.to_csv(rD:\data\anhui_town_pop_change.csv, indexFalse, encodingutf-8-sig)utf-8-sig能让Excel正确识别中文字段。后续如果需要做空间计量把这个表join回乡镇shp即可。注意如果df还是GeoDataFrameto_csv会保留geometry列但CSV里会变成WKT字符串导出前可以先把geometry列去掉。5. 边缘乡镇的边界误差面积加权与QGIS验证5.1 为什么统计值总是差几十人用zonal_stats默认的“像元中心在面内”规则时面积小的乡镇很容易漏掉边缘像元。比如某个乡镇边界刚好切过像元而像元中心落在外侧整个像元就被丢弃。对只有几平方公里的乡镇漏掉一个像元就是几百人。解决思路有两个一是用all_touchedTrue把与面有交集的像元全部纳入这样会多算但不会漏二是做面积加权按像元落在面内的比例分配人口。实际项目中我通常先跑一遍默认规则再和统计年鉴里已知的乡镇人口对比。如果差异超过5%再改用加权算法。5.2 用rasterstats做面积加权统计rasterstats提供了一个不太常见的开关stats_w zonal_stats(towns, out_tif, stats[mean], nodata-9999, all_touchedFalse, rasterizeTrue)rasterizeTrue会先把多边形栅格化再计算每个像元和多边形相交的面积权重得到的均值比默认方法更接近真实。注意这个参数会明显增加计算时间200个乡镇以上时建议只对差异较大的乡镇单独跑。另一种方式是干脆用all_touchedTrue做一遍粗算然后检查每个乡镇的有效像元数是否和乡镇面积匹配。比如一个10km²的乡镇在30弧秒分辨率下大约有10-15个有效像元如果只有2个说明边界很可能有缝隙。5.3 用QGIS独立验证一遍无论Python结果多好看我都会在QGIS里再验证一次。步骤如下把pop_anhui_2020.tif和乡镇shp拖进QGIS确认两者坐标系都是WGS84。打开“处理工具箱” → “栅格分析” → “栅格区域统计”。输入栅格选pop_anhui_2020.tif输入矢量选乡镇shp统计量选“求和”和“平均值”。运行后把结果和Python的pop_sum_2020对比差异超过1%就回查NoData设置和裁剪范围。QGIS的算法和rasterstats类似但实现细节不同正好用来互相验证。如果你发现QGIS能算出值而Python返回0基本可以断定是NoData或文件路径问题。5.4 把元数据写进数据目录最后建议在数据压缩包同目录放一个README.md记录所有处理步骤和数据状态至少包括原始数据的NoData值重采样方法和目标分辨率裁剪矢量文件路径五期数据是否网格对齐乡镇边界的坐标系与数据日期这类生成数据往往会在项目里流转很多年不写清楚后期根本没法复盘。本文还有配套的精品资源点击获取 乡镇级人口密度栅格处理与2000-2020变化分析实战
简介安徽省乡镇级2000至2020年人口密度栅格数据面向地理信息系统、城乡规划、人口统计及环境科学等专业人员填补了高分辨率长时序乡镇尺度人口数据的空缺可为区域研究与决策提供基础数据支撑。数据涵盖2000、2005、2010、2015、2020年五期栅格基于约4万个行政单元的人口统计并调整匹配联合国人口总数修订版采用WGS84地理坐标系分辨率约为1公里有效支撑从省到乡镇多层级的人口分布制图与时空演变分析。压缩包整体约3.68MB内含栅格图像文件可直接放入ArcGIS、QGIS等主流平台使用便于地图制图与空间统计。目前已有175人学习下载结合人口普查和年鉴资料可快速生成人口密度图、识别人口集聚与疏散区域为区域规划、公共服务设施布局及资源配置提供量化依据也可用于应急管理与社会经济研究。1. 乡镇级人口密度栅格一份能直接拆出来用的时空数据拿到这个“安徽省乡镇级2000-2020年人口密度.rar”很多人第一反应是解压看里面有几种文件但真正值钱的是数据本身的设计它把2000、2005、2010、2015、2020五期人口密度做成约1km分辨率的WGS84栅格每个像元值代表该像元对应区域的人口数或人口密度。因为是基于约4万个行政单位的人口普查数据拟合到30弧秒网格所以既能在省、市、县层面做对比也能切到乡镇边界做统计。下面就以这个数据集为对象从RAR里的栅格文件一路拆到Python/GDAL批量提取乡镇人口密度再讨论如何做多期变化分析最后给出边界误差修正的验证方法。2. 解开RAR之后认识栅格文件、WGS84坐标与30弧秒分辨率2.1 解压RAR文件时先看这三件事拿到.rar压缩包别急着拖到GIS软件里。我一般做三件事第一用支持RAR5的7-Zip或WinRAR解压解压前先看压缩包内文件如果单个栅格超过200MB后续在Python里读取时就要考虑内存映射或分块读取第二确认是否带.prj、.tfw等辅助文件如果只有裸的.tif或.bil需要从元数据描述里补坐标系第三若压缩包带有密码这类公开数据很少见直接联系来源方确认不要用网上的“rar密码移除”类破解工具避免拿到被改过的文件也可能顺手装一堆垃圾软件。提示RAR包内如果包含同名.aux.xml别删里面通常是统计信息和无效值定义读取NoData时会用到。2.2 30弧秒网格到底是多少米原文说“30弧秒网格单元”精度约1km。这里的弧秒是角度单位1度3600弧秒30弧秒1/120度。在赤道附近1度约111km所以30弧秒约0.925km。随着纬度升高经度方向的实际距离会缩短但纬度方向仍约1.11km在“约1km”范围内。对安徽省而言地处北纬29°41′34°38′一个网格的东西跨度约0.8-1.0km落在乡镇级别刚好能区分集中居民点和外围农田。很多人在ArcGIS里看到WGS84坐标系的栅格就直接用“度”去做长度或面积计算这是常见坑。栅格像元大小显示为0.0083333333度实际上横向和纵向代表的地面距离不同。要做面积统计必须用等积投影或按纬度换算否则乡镇面积会偏小人口密度会被高估。2.3 用GDAL快速查看栅格元数据拿到TIF后第一个命令往往是gdalinfogdalinfo population_anhui_2020.tif关注这样几个字段Size是宽高Origin是左上角坐标Pixel Size是像元尺寸Coordinate System是坐标系描述NoData Value是无效值。例如Size is 2400, 3600 Origin (114.870833333333337,34.662500000000001) Pixel Size (0.008333333333333,-0.008333333333333) Coordinate System is: GEOGCRS[WGS 84, ...] NoData Value -9999这里的Pixel Size正负号表示行列方向x方向为正y方向为负意味着数据按北到南存储。NoData Value不是0统计时必须先过滤掉-9999。如果输出的Size和Origin与你预期不一致说明数据是经过裁剪或重采样的需要记录下这个源信息。另外gdalinfo -stats还能输出栅格的最小/最大/平均值。如果看到最小值为负数除了NoData外还要怀疑是否存在-1这类“占位值”。有些数据源会把水域或者边界外设为-1而不是标准的NoData。读取时要把这些值一并排除。2.4 人口计数与人口密度的区别原文提到“人口计数已调整并匹配联合国国家总数修订版”又说这是“人口密度栅格”。这里要留意有些版本存的是每个像元的人口数整数有些存的是人/平方公里浮点。判断方法是打开栅格直方图人口数的像元值通常在几十到几千密度的像元值可以到数万。如果目标是要汇总得到乡镇总人口用整数人口数直接累加更自然如果只是做密度等级图则直接用人口密度。两类值不能混在一起比较否则结果会差出几个数量级。属性含义影响分辨率30弧秒约1km能支撑乡镇级不能支撑村级坐标系WGS84地理坐标系计算面积、距离前需投影单位每像元人口数或人/km²决定后续统计口径NoData通常是-9999统计时必须过滤时间2000/2005/2010/2015/2020可做时间序列变化检测3. 用PythonGDAL把安徽乡镇级人口密度“切”出来3.1 环境准备不需要ArcGIS也能跑通如果只用ArcGIS的“分区统计”工具也能完成但对五期数据、一百多个乡镇的小批量场景PythonGDAL更可控也方便放进批量处理脚本。环境准备pip install gdal numpy pandas geopandas shapely rasterstatsWindows下gdal的pip安装经常找不到dll推荐用condaconda install -c conda-forge gdal。rasterstats是专门做栅格分区统计的库能少写很多循环。装完以后先确认GDAL版本from osgeo import gdal print(gdal.__version__)如果版本号是3.x注意gdal.Warp的某些参数行为和2.x略有不同但本文的用法在两者都适用。如果你的乡镇边矢量不是WGS84坐标系比如是EPSG:2383安徽地方坐标系需要先转成WGS84否则裁剪时会出现错位或统计结果为空。转换用geopandas很简单import geopandas as gpd towns gpd.read_file(anhui_towns.shp) towns towns.to_crs(EPSG:4326)3.2 用gdal.Warp裁剪到安徽边界原始栅格通常覆盖全国甚至全球直接裁剪到安徽省能减少后续计算的内存占用。裁剪建议用省界矢量而不是简单矩形范围因为省界是不规则多边形。代码from osgeo import gdal, gdalconst in_tif rD:\data\pop_chn_2020.tif out_tif rD:\data\pop_anhui_2020.tif shp rD:\data\anhui_boundary.shp # WGS84坐标系 gdal.Warp(out_tif, in_tif, formatGTiff, cutlineDSNameshp, cropToCutlineTrue, dstNodata-9999, resampleAlggdalconst.GRA_NearestNeighbour, options[TARGET_ALIGNED_YES]) ds gdal.Open(out_tif) print(ds.RasterXSize, ds.RasterYSize)cutlineDSName指定裁剪矢量cropToCutlineTrue让输出范围与矢量范围对齐。dstNodata-9999设置裁剪后边界外区域的无效值要和原始数据一致。resampleAlg用最近邻不改变像元值语义。如果后续做多期计算还需要统一网格就要在gdal.Warp中加xRes和yRes。裁剪后建议检查一下输出范围与省界范围是否一致可以再用gdalinfo看一次。3.3 按乡镇矢量分区统计平均人口密度乡镇级行政区划一般有几百个面要素。想得到每个乡镇的人口密度通常做法是“先用乡镇面切栅格再统计面内像元”。rasterstats的zonal_stats封装了这个过程import geopandas as gpd from rasterstats import zonal_stats towns gpd.read_file(rD:\data\anhui_towns.shp) # 乡镇边界WGS84 out_tif rD:\data\pop_anhui_2020.tif stats zonal_stats(towns, out_tif, stats[sum, mean, count], nodata-9999) towns[pop_sum_2020] [s[sum] for s in stats] towns[pop_mean_2020] [s[mean] for s in stats] towns[valid_count] [s[count] for s in stats] towns.to_file(rD:\data\anhui_towns_pop_2020.geojson, driverGeoJSON)stats[sum, mean, count]分别表示像元值求和、平均值、有效像元数量。如果栅格值代表每像元人口数sum就是乡镇总人口数如果栅格值已经是密度sum就不是总人口而是密度之和需要小心。nodata-9999告诉库跳过无效值。valid_count可以用来排查统计质量如果某个乡镇的有效像元数是0说明乡镇边界和栅格没有交集多半是坐标系不匹配或乡镇边界有问题。3.4 批量处理5期并统一网格2000到2020共五期写循环即可但要注意网格对齐years [2000, 2005, 2010, 2015, 2020] for year in years: in_tif rfD:\data\pop_chn_{year}.tif out_tif rfD:\data\pop_anhui_{year}.tif gdal.Warp(out_tif, in_tif, formatGTiff, cutlineDSNameshp, cropToCutlineTrue, dstNodata-9999, xRes0.008333333333333, yRes0.008333333333333, resampleAlggdalconst.GRA_Bilinear) stats zonal_stats(towns, out_tif, stats[sum, mean], nodata-9999) towns[fpop_sum_{year}] [s[sum] for s in stats] towns[fpop_mean_{year}] [s[mean] for s in stats]这里我使用了双线性重采样GRA_Bilinear而不是最近邻目的是把不同期数据的网格统一到同一个像元网格上。双线性重采样会平滑像元值但对于人口这种连续变量影响不大。如果你做的是精度要求高的变化检测可以用GRA_Cubic但计算会更慢。提示xRes和yRes需要和原始像素大小一致否则会改变分辨率。这里写的是30弧秒。实际使用前建议打印原始数据的Origin和Pixel Size确认没有差值。4. 构建2000-2020乡镇人口密度时间序列从栅格代数到变化分级4.1 为什么不能直接拿两期栅格相减拿到五期栅格后有人直接打开栅格计算器用pop2020 - pop2000得到的“变化值”看似直观但有三个问题一是两期数据可能经过不同裁剪或网格对齐必须先在Python里做完一致性处理二是人口密度栅格存在空间自相关逐像元相减会放大噪声比如某个像元2000年是空地2010年人口正好被分配到相邻像元结果出现巨大负数三是乡镇级分析应该先把像元值聚合成乡镇值再做变化分析而不是先做差值再聚合。正确的顺序是先完成上一章的分区统计生成每个乡镇一行、包含5个年份人口数的表然后基于行数据计算变化率。4.2 用Pandas计算乡镇级人口变化率上一章已经生成GeoJSON现在读取并计算import geopandas as gpd import pandas as pd df gpd.read_file(rD:\data\anhui_towns_pop_2020.geojson) pop_cols [pop_sum_2000, pop_sum_2005, pop_sum_2010, pop_sum_2015, pop_sum_2020] # 二十年总变化率 df[total_change] (df[pop_sum_2020] - df[pop_sum_2000]) / df[pop_sum_2000] # 分阶段变化率 df[growth_2000_2005] (df[pop_sum_2005] - df[pop_sum_2000]) / df[pop_sum_2000] df[growth_2010_2015] (df[pop_sum_2015] - df[pop_sum_2010]) / df[pop_sum_2010] df[growth_2015_2020] (df[pop_sum_2020] - df[pop_sum_2015]) / df[pop_sum_2015] print(df[total_change].describe())total_change是二十年总变化率后续三个字段用来观察阶段性差异。describe()输出的四分位数能帮你判断全省是普遍增长还是两极分化。如果最大最小差异太大先排查是否有个别乡镇边界缺失或者NoData值被当成了0计算。4.3 划分增长型、收缩型、稳定型乡镇只看连续数值不够直观我一般按业务含义分三类总增长率小于-10%为收缩型大于20%为增长型中间为稳定型。为什么不直接用0做阈值因为人口普查和网格化拟合本身有误差-5%到5%大概率是噪声。def classify(row): v row[total_change] if v -0.10: return 收缩型 elif v 0.20: return 增长型 else: return 稳定型 df[type] df.apply(classify, axis1) print(df[type].value_counts())写回GeoJSONdf.to_file(rD:\data\anhui_town_change.geojson, driverGeoJSON)在QGIS里用type字段做分类渲染会立刻看到空间分布皖北平原的收缩型乡镇往往连片合肥和沿江城市的增长型乡镇像星星一样散落。如果你看到的模式完全随机那就要回头检查数据匹配是否出错。分类阈值不是固定的。如果研究区域以合肥、芜湖等城市为主20%的增长阈值会漏掉许多“低增长但绝对量大”的郊区乡镇。此时可以把阈值提高到30%或者改用自然断点法。判断标准应该是分类结果能体现出空间集聚性而不是把全省切成棋盘格。4.4 结果表的字段设计与导出为了展示最终的表结构我列出常用字段设计字段名示例值说明name某镇乡镇名称pop_sum_2000452312000年人口总数pop_sum_2020603402020年人口总数total_change0.334二十年变化率type增长型分类结果这个表可以直接输出给业务方df.to_csv(rD:\data\anhui_town_pop_change.csv, indexFalse, encodingutf-8-sig)utf-8-sig能让Excel正确识别中文字段。后续如果需要做空间计量把这个表join回乡镇shp即可。注意如果df还是GeoDataFrameto_csv会保留geometry列但CSV里会变成WKT字符串导出前可以先把geometry列去掉。5. 边缘乡镇的边界误差面积加权与QGIS验证5.1 为什么统计值总是差几十人用zonal_stats默认的“像元中心在面内”规则时面积小的乡镇很容易漏掉边缘像元。比如某个乡镇边界刚好切过像元而像元中心落在外侧整个像元就被丢弃。对只有几平方公里的乡镇漏掉一个像元就是几百人。解决思路有两个一是用all_touchedTrue把与面有交集的像元全部纳入这样会多算但不会漏二是做面积加权按像元落在面内的比例分配人口。实际项目中我通常先跑一遍默认规则再和统计年鉴里已知的乡镇人口对比。如果差异超过5%再改用加权算法。5.2 用rasterstats做面积加权统计rasterstats提供了一个不太常见的开关stats_w zonal_stats(towns, out_tif, stats[mean], nodata-9999, all_touchedFalse, rasterizeTrue)rasterizeTrue会先把多边形栅格化再计算每个像元和多边形相交的面积权重得到的均值比默认方法更接近真实。注意这个参数会明显增加计算时间200个乡镇以上时建议只对差异较大的乡镇单独跑。另一种方式是干脆用all_touchedTrue做一遍粗算然后检查每个乡镇的有效像元数是否和乡镇面积匹配。比如一个10km²的乡镇在30弧秒分辨率下大约有10-15个有效像元如果只有2个说明边界很可能有缝隙。5.3 用QGIS独立验证一遍无论Python结果多好看我都会在QGIS里再验证一次。步骤如下把pop_anhui_2020.tif和乡镇shp拖进QGIS确认两者坐标系都是WGS84。打开“处理工具箱” → “栅格分析” → “栅格区域统计”。输入栅格选pop_anhui_2020.tif输入矢量选乡镇shp统计量选“求和”和“平均值”。运行后把结果和Python的pop_sum_2020对比差异超过1%就回查NoData设置和裁剪范围。QGIS的算法和rasterstats类似但实现细节不同正好用来互相验证。如果你发现QGIS能算出值而Python返回0基本可以断定是NoData或文件路径问题。5.4 把元数据写进数据目录最后建议在数据压缩包同目录放一个README.md记录所有处理步骤和数据状态至少包括原始数据的NoData值重采样方法和目标分辨率裁剪矢量文件路径五期数据是否网格对齐乡镇边界的坐标系与数据日期这类生成数据往往会在项目里流转很多年不写清楚后期根本没法复盘。本文还有配套的精品资源点击获取