简介珠江三角洲城市群区域开发密度数据集面向城市规划、经济地理与区域可持续发展研究者汇集1998年、2006年、2012年三个时点的地理空间信息可用于揭示近二十年间城市空间扩展、经济集聚与开发强度的时空演变规律。压缩包共44个文件、约10.83MB以Shapefile矢量格式.shp、.shx、.dbf、.prj等和GeoTIFF栅格数据.tif、.tfw、.ovr为主分层存储研究区域边界、经济密度、发展紧凑度、发展强度以及1公里网格市域GDP在ArcGIS、QGIS等常用GIS平台中可直接打开和横向比较。已有74人浏览学习。借助其中包含的城市行政边界、河流海岸线、GDP空间分布、建成区紧凑度、建筑覆盖率与人口密度等指标研究者能逐项对比不同年份的空间结构变化识别经济热点与城乡差异亦可延伸至城市可持续性评估、交通规划、环境保护及社会经济不平等分析为区域发展政策制定提供量化支撑。1. 这份珠三角密度数据集为什么值得拆开看拿到一个名为“珠江三角洲城市群区域开发密度数据集199820062012.rar”的压缩包第一反应不是解压而是先问它能回答什么问题我给城市地理项目做前期数据梳理时常遇到这种混合了矢量边界和栅格GDP的“研究型数据集”。它不像标准测绘产品有统一规范往往带有一套自己的命名逻辑比如目录里的2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity以及后面的3_NAGDP_1km。这套数据覆盖广州、深圳、佛山、东莞、中山、珠海、江门、肇庆、惠州九市跨越三个时间断面用两种数据模型矢量面与1km栅格表达同一区域的发展演变。对规划师和计量研究者来说最直接的用途是重现“2006是拐点”之类的空间经济命题对GIS工程师来说它则是一份练习时空数据清洗、投影转换和栅格-矢量联动的标准样本。下面从解压这一步开始按实际分析链路把每个文件的作用、读取方式和常见坑逐层拆开尽量做到拿这份数据的人能直接照着操作。2. 先把.rar拆开文件结构、坐标系统与栅格/矢量格式核对拿到压缩包先不要急着在ArcGIS里拖进图层。数据集使用的.rar压缩格式在Windows上有WinRAR、7-Zip等常见工具在Linux服务器上则依赖unrar或7z命令行。解压前先确认目标目录不含中文路径否则后续用Python处理时容易在路径解析上踩编码坑。2.1 解压不是双击了事Linux命令与完整性校验在Linux下常见的做法是先安装unrar然后执行解压命令sudo apt install unrar # Debian/Ubuntu unrar x 珠江三角洲城市群区域开发密度数据集199820062012.rar ./prd_dataset/参数说明x表示保留压缩包内的目录结构./prd_dataset/是解压目标目录。解压后建议用unrar t 压缩包名.rar做一次完整性测试因为后续数据里的.sbn/.sbx如果损坏虽然不致命但会影响ArcGIS的空间索引读取。在Windows下我更推荐用7-Zip的右键“解压到当前文件夹”它对中文文件名和全角括号的处理比某些版本更稳定。解压完成后先看目录树确认是否存在1_studyarea、2a_DevelopmentIntensity等顶层文件夹以及每个文件夹内是否都有成对的.shp和.dbf文件。这一步虽然琐碎却能在正式分析前暴露文件缺失。2.2 Shapefile的十个零件哪些文件必须一起拷贝Shapefile不是单文件而是多个文件集合。本数据集中每个矢量图层都有至少8个同名文件例如EconomicDensity.shp、.shx、.dbf、.prj、.cpg、.sbn、.sbx、.shp.xml。其中.shp存几何.shx存索引.dbf存属性三者缺一不可.prj存投影信息.cpg存属性表编码.sbn/.sbx是ArcGIS维护的空间索引.shp.xml是元数据。扩展名作用缺失后的影响.shp几何要素无法读取.shx几何索引部分库会报错.dbf属性表没有任何属性字段.prj坐标系定义GIS软件可能不识别坐标.cpg字符集声明中文字段名或值乱码.sbn/.sbx空间索引仅影响查询性能.shp.xml元数据可忽略值得注意的是DevelopmentIntensity文件夹内.cpg被写成了大写.CPG。在Linux下文件系统区分大小写如果从Windows复制到Linuxgeopandas自动搜索.cpg时可能找不到导致属性编码识别失败。我遇到这种情况时一般会用一个for循环把所有文件改为小写扩展名或者直接在读取时显式指定encoding。2.3 先看投影.prj与.tfw告诉你在哪种坐标下投影是一切空间计算的基准。用文本编辑器打开StudyArea.prj里面是一段WKT字符串。珠三角区域常见的坐标系有WGS84地理坐标、WGS84 UTM 50N投影以及国家2000坐标系具体要看文件内容。我判断投影时重点关注三个信息基准面datum、投影名称projection、单位units。地理坐标系单位为度投影坐标系单位为米。NAGDP_1km_1998.tfw是GeoTIFF的世界文件用普通文本编辑器打开就能看到六行数字前两行是像元尺寸和旋转项后两行是栅格原点。更直观的方式是用gdalinfogdalinfo NAGDP_1km_1998.tif | head -n 40gdalinfo输出里会显示Coordinate System is、Origin和Pixel Size。如果三个年份片的坐标系不一致不要贸然做差值运算必须先统一投影。对于这份数据1_studyarea提供的边界可以看作所有后续操作的空间基准。3. 三种密度指标的定义与属性表读取用geopandas把字段读出来目录里三个以数字开头的子文件夹分别叫2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity。这个“2c→2a”的字母顺序有讲究经济密度是结果紧凑度是形态开发强度是过程入口。摘要里解释得很清楚经济密度度量单位面积经济产值紧凑度度量建成区集中程度开发强度则反映建设密度和人口密度。三者放在一起才能解释城市扩张的机理。3.1 三个shapefile的字段差异与命名启示我用ogrinfo或geopandas快速列出字段总能发现一些规律。通常EconomicDensity.dbf持有类似GDP_KM2、VAL1998这样的字段DevelopmentCompactness.dbf里有类似SHAPE_INDEX的值DevelopmentIntensity.dbf里则可能是BUILD_PCT或POP_KM2。虽然名字不同但数值都是区域级统计量。命名前缀“2c、2b、2a”更像是一个生产流水线先算出经济产出密度再评估形态紧凑度最后落到开发强度。这样设计的好处是属性表里可以不冗余存储几何信息靠FID关联。3.2 geopandas读取示例import geopandas as gpd # 注意cpg编码声明如果乱码就指定encodinggbk eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(eco_1998.columns) print(eco_1998.head()) print(eco_1998[[GDP_KM2, YEAR]].describe())这段代码先用read_file加载 Shapefileencoding参数匹配.cpg文件声明的字符集。然后打印所有列名和统计摘要。describe()会输出count、mean、std等字段能快速判断GDP_KM2是否存在异常负值或大范围空值。拿到属性表后最好顺带检查YEAR字段因为一个shp文件里可能存储多期数据靠YEAR区分。3.3 紧凑度与强度的计算口径紧凑度计算公式在文献里常用C 2√(πA)/P其中A是建成区面积P是周长结果落在0到1之间。本数据集里的DevelopmentCompactness如果直接是数值大概率就是这类形状指数或紧凑比。至于DevelopmentIntensity我见过更合理的字段是“单位面积建设用地面积”或“人口密度”的网格化统计结果。比较年份时要注意如果行政边界本身发生了调整比如2006年某镇并入街道则必须先使用1_studyarea里统一后的边界做裁剪或交集不然紧凑度会被区域面积变化干扰。4. 跨年份空间分析如何提取“开发密度”的变化信号三个年份的数据放一起最直接的分析是看同一行政单元在不同指标上的走势。但前提是边界一致。1_studyarea提供了研究区边界而各年份shp的几何可能来自不同来源必须先处理掉缝隙、重叠与坐标系差异。4.1 统一坐标系与属性对齐先检查三个文件的坐标系是否一致import geopandas as gpd study gpd.read_file(./1_studyarea/StudyArea.shp) eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(study.crs) print(eco_1998.crs) if eco_1998.crs ! study.crs: eco_1998 eco_1998.to_crs(study.crs)如果.prj文件缺失或定义不完整to_crs会抛错。应急方案是从一个已知正确的GeoDataFrame复制crs对象例如eco_1998.crs study.crs但这样做的前提是你确认两套数据的椭球和基准面一致否则会出现几米的偏移。属性对齐更麻烦。如果2006年的字段结构变了比如新增了“空间GDP”字段而1998年没有直接用pd.concat会生成大量NaN。我常用的做法是先只保留几个核心字段统一命名为region、year、value再按区域名进行merge或groupby。# 读取紧凑度数据并做空间连接 compact gpd.read_file(./2b_DevelopmentCompactness/DevelopmentCompactness.shp) joined gpd.sjoin(study, compact, howleft, predicateintersects) # 按城市名聚合 grouped joined.groupby(CITY_NAME)[COMPACT].mean().reset_index()代码中sjoin的predicateintersects表示两个要素只要任意部分相交就算匹配这能规避边界微小缝隙。groupby取均值是因为一个城市可能由多个面要素组成。如果希望面积加权则需要先计算相交面积joined[inter_area] joined.geometry.intersection(study.unary_union).area weighted joined.groupby(CITY_NAME).apply( lambda df: (df[COMPACT] * df[inter_area]).sum() / df[inter_area].sum() )注意apply里的加权公式当区域面积差异大时简单均值会偏向面积小但数值高的单元加权均值能反映整体状态。4.2 生成三年对比表把三个年份的经济密度整合到一张宽表可以直接用pivot_tableimport pandas as pd # 假设三个gdf都有YEAR和VALUE字段 eco_1998[year] 1998 eco_2006[year] 2006 eco_2012[year] 2012 frames [eco_1998[[CITY_NAME, year, GDP_KM2]], eco_2006[[CITY_NAME, year, GDP_KM2]], eco_2012[[CITY_NAME, year, GDP_KM2]]] all_df pd.concat(frames) pivot all_df.pivot_table(indexCITY_NAME, columnsyear, valuesGDP_KM2, aggfuncmean) pivot[growth_06_98] pivot[2006] / pivot[1998] - 1 pivot[growth_12_06] pivot[2012] / pivot[2006] - 1pivot_table的aggfunc默认mean它会处理同一城市下多个面要素的重复记录。用增长率列能直接识别哪几年扩张最快。4.3 变化可视化可视化可以先用matplotlib画折线趋势再用geopandas给每个城市填充增长率颜色。需要注意用颜色分级图时最好将增长率分位数分成5类避免极端值把配色拉平。5. 处理1km GDP栅格用rasterio读取、重投影与统计栅格部分是这份数据里最硬的骨头。NAGDP_1km_1998.tif这类文件是GeoTIFF除了.tif本体还有.tfw世界文件、.aux.xml辅助信息、.ovr金字塔。这套文件对新手不太友好但处理思路很固定先看元数据再处理空值最后做统计或重投影。5.1 栅格文件的组成.tif、.tfw、.ovr各干什么.tfw是文本世界文件记录像元尺寸和左上角坐标.ovr是金字塔文件如果缺失大范围缩放变慢但不影响数值精度.aux.xml包含统计信息和nodata标记。读取时我习惯用rasterio而不是gdal命令行因为能直接返回numpy数组import rasterio with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data src.read(1) nodata src.nodata profile src.profile print(nodata , nodata) print(数据类型 , data.dtype) print(shape , data.shape)src.read(1)读出第一波段nodata是无值标记。GDP数据通常用-9999或0表示无数据拿到数组后第一步要把这些值屏蔽掉否则计入均值会让结果严重偏低import numpy as np masked np.ma.masked_where(data 0, data) print(masked.mean())5.2 做年度差值并重投影三个年份的tif如果像元尺寸和原点不完全一致不要直接做数组减法。正确做法是先统一到同一个格网。可以用reproject完成from rasterio.warp import reproject, Resampling, calculate_default_transform with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data_1998 src.read(1) src_crs src.crs src_transform src.transform nodata src.nodata # 目标坐标系如果源是WGS84地理坐标这里用EPSG:4326 transform, width, height calculate_default_transform( src_crs, EPSG:32650, src.width, src.height, *src.bounds ) dst np.zeros((height, width), dtypenp.float32) reproject( sourcedata_1998, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformtransform, dst_crsEPSG:32650, resamplingResampling.bilinear, src_nodatanodata, dst_nodata-9999 )上述代码中calculate_default_transform根据源范围和目标投影计算新尺寸bilinear重采样适用于连续型GDP值如果处理分类栅格则改用nearest。src_nodata与dst_nodata保持一致保证重投影后空洞仍是同一标记。如果不想写这么长的代码优先用rioxarray它将坐标和投影封装得更好一行就能做差值import rioxarray ds_1998 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_1998.tif) ds_2006 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_2006.tif) diff ds_2006 - ds_1998但要先确认两个对象的空间参考一致若不一致则调用ds_1998.rio.reproject_match(ds_2006)对齐。5.3 高分辨率GDP数据的局限与坑1km格网GDP并不是真实观测值它一般基于夜间灯光、土地利用和统计年鉴做空间化插值因此“某个像元上的GDP”并不代表那个位置真有这么多产出。解读时更适合用在区域总量或相对比上。我在做行政单元统计时会先用行政区边界裁剪然后逐像元累加。这里最典型的坑有两个一是边缘像元被边界切成窄条如果连通面积小于一个像元直接累加会高估边界地带数值二是.ovr金字塔文件版本不一致可能让某些软件读到错误nodata。解决办法是用rasterstats的all_touchedFalse只保留完整落入边界的像元并在读取时显式指定maskedTrue。提示如果某一年份的GDP均值明显异常优先检查nodata值是否被默认为0再检查重投影时是否把边缘的nodata像元插值成小数值。6. 收尾技巧用rasterstats一键把1km GDP汇总到行政区前面的操作分别处理了矢量和栅格实际交付成果时常需要“一张Excel表”每一行是一个区县列是三个年份的GDP密度与总量。这时用rasterstats包可以省掉大量手写循环。它专门用来把栅格数据按矢量多边形进行分区统计pip install rasterstats然后写一个极简函数import pandas as pd from rasterstats import zonal_stats stats zonal_stats( ./1_studyarea/StudyArea.shp, # 矢量边界 ./3_NAGDP_1km/NAGDP_1km_1998.tif, # 栅格文件 stats[sum, mean, min, max], nodata-9999, # 与tif元数据一致 all_touchedFalse # 只统计完整覆盖的像元 ) df pd.DataFrame(stats)zonal_stats第一个参数接收矢量边界第二个参数接收栅格路径stats支持sum/mean/max/min/median等聚合统计。nodata参数能过滤栅格中的无值像元避免把-9999当作真实数值。all_touchedFalse表示只统计多边形中心点位于内部的像元这比True更严格也更能避免边缘噪声。三个年份只需包装成一个循环def extract_year(year): raster f./3_NAGDP_1km/NAGDP_1km_{year}.tif st zonal_stats( ./1_studyarea/StudyArea.shp, raster, stats[sum, mean], nodata-9999, all_touchedFalse ) return pd.DataFrame(st) result pd.concat( [extract_year(1998), extract_year(2006), extract_year(2012)], axis1 ) result.columns [GDP_sum_1998, GDP_mean_1998, GDP_sum_2006, GDP_mean_2006, GDP_sum_2012, GDP_mean_2012]这个封装直接生成一个DataFrame再配合study里的行政区名称列即可保存为CSVstudy gpd.read_file(./1_studyarea/StudyArea.shp) result.insert(0, city_name, study[CITY_NAME]) result.to_csv(./prd_gdp_summary.csv, indexFalse)一个小技巧如果发现某年的sum异常低先回到第5章检查nodata和重投影步骤或者换个思路用面积加权平均——因为栅格像元面积可能随纬度变化而不同尤其在投影坐标系里不同纬度像元实际面积并不完全一致。遇到这种情况可以先用area属性求出单位面积值再乘实际面积。这套流程跑通后后面再加2018年数据只需要复制这个函数模板替换年份即可。本文还有配套的精品资源点击获取 珠三角开发密度数据集全解析：从Shapefile到栅格处理实战
简介珠江三角洲城市群区域开发密度数据集面向城市规划、经济地理与区域可持续发展研究者汇集1998年、2006年、2012年三个时点的地理空间信息可用于揭示近二十年间城市空间扩展、经济集聚与开发强度的时空演变规律。压缩包共44个文件、约10.83MB以Shapefile矢量格式.shp、.shx、.dbf、.prj等和GeoTIFF栅格数据.tif、.tfw、.ovr为主分层存储研究区域边界、经济密度、发展紧凑度、发展强度以及1公里网格市域GDP在ArcGIS、QGIS等常用GIS平台中可直接打开和横向比较。已有74人浏览学习。借助其中包含的城市行政边界、河流海岸线、GDP空间分布、建成区紧凑度、建筑覆盖率与人口密度等指标研究者能逐项对比不同年份的空间结构变化识别经济热点与城乡差异亦可延伸至城市可持续性评估、交通规划、环境保护及社会经济不平等分析为区域发展政策制定提供量化支撑。1. 这份珠三角密度数据集为什么值得拆开看拿到一个名为“珠江三角洲城市群区域开发密度数据集199820062012.rar”的压缩包第一反应不是解压而是先问它能回答什么问题我给城市地理项目做前期数据梳理时常遇到这种混合了矢量边界和栅格GDP的“研究型数据集”。它不像标准测绘产品有统一规范往往带有一套自己的命名逻辑比如目录里的2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity以及后面的3_NAGDP_1km。这套数据覆盖广州、深圳、佛山、东莞、中山、珠海、江门、肇庆、惠州九市跨越三个时间断面用两种数据模型矢量面与1km栅格表达同一区域的发展演变。对规划师和计量研究者来说最直接的用途是重现“2006是拐点”之类的空间经济命题对GIS工程师来说它则是一份练习时空数据清洗、投影转换和栅格-矢量联动的标准样本。下面从解压这一步开始按实际分析链路把每个文件的作用、读取方式和常见坑逐层拆开尽量做到拿这份数据的人能直接照着操作。2. 先把.rar拆开文件结构、坐标系统与栅格/矢量格式核对拿到压缩包先不要急着在ArcGIS里拖进图层。数据集使用的.rar压缩格式在Windows上有WinRAR、7-Zip等常见工具在Linux服务器上则依赖unrar或7z命令行。解压前先确认目标目录不含中文路径否则后续用Python处理时容易在路径解析上踩编码坑。2.1 解压不是双击了事Linux命令与完整性校验在Linux下常见的做法是先安装unrar然后执行解压命令sudo apt install unrar # Debian/Ubuntu unrar x 珠江三角洲城市群区域开发密度数据集199820062012.rar ./prd_dataset/参数说明x表示保留压缩包内的目录结构./prd_dataset/是解压目标目录。解压后建议用unrar t 压缩包名.rar做一次完整性测试因为后续数据里的.sbn/.sbx如果损坏虽然不致命但会影响ArcGIS的空间索引读取。在Windows下我更推荐用7-Zip的右键“解压到当前文件夹”它对中文文件名和全角括号的处理比某些版本更稳定。解压完成后先看目录树确认是否存在1_studyarea、2a_DevelopmentIntensity等顶层文件夹以及每个文件夹内是否都有成对的.shp和.dbf文件。这一步虽然琐碎却能在正式分析前暴露文件缺失。2.2 Shapefile的十个零件哪些文件必须一起拷贝Shapefile不是单文件而是多个文件集合。本数据集中每个矢量图层都有至少8个同名文件例如EconomicDensity.shp、.shx、.dbf、.prj、.cpg、.sbn、.sbx、.shp.xml。其中.shp存几何.shx存索引.dbf存属性三者缺一不可.prj存投影信息.cpg存属性表编码.sbn/.sbx是ArcGIS维护的空间索引.shp.xml是元数据。扩展名作用缺失后的影响.shp几何要素无法读取.shx几何索引部分库会报错.dbf属性表没有任何属性字段.prj坐标系定义GIS软件可能不识别坐标.cpg字符集声明中文字段名或值乱码.sbn/.sbx空间索引仅影响查询性能.shp.xml元数据可忽略值得注意的是DevelopmentIntensity文件夹内.cpg被写成了大写.CPG。在Linux下文件系统区分大小写如果从Windows复制到Linuxgeopandas自动搜索.cpg时可能找不到导致属性编码识别失败。我遇到这种情况时一般会用一个for循环把所有文件改为小写扩展名或者直接在读取时显式指定encoding。2.3 先看投影.prj与.tfw告诉你在哪种坐标下投影是一切空间计算的基准。用文本编辑器打开StudyArea.prj里面是一段WKT字符串。珠三角区域常见的坐标系有WGS84地理坐标、WGS84 UTM 50N投影以及国家2000坐标系具体要看文件内容。我判断投影时重点关注三个信息基准面datum、投影名称projection、单位units。地理坐标系单位为度投影坐标系单位为米。NAGDP_1km_1998.tfw是GeoTIFF的世界文件用普通文本编辑器打开就能看到六行数字前两行是像元尺寸和旋转项后两行是栅格原点。更直观的方式是用gdalinfogdalinfo NAGDP_1km_1998.tif | head -n 40gdalinfo输出里会显示Coordinate System is、Origin和Pixel Size。如果三个年份片的坐标系不一致不要贸然做差值运算必须先统一投影。对于这份数据1_studyarea提供的边界可以看作所有后续操作的空间基准。3. 三种密度指标的定义与属性表读取用geopandas把字段读出来目录里三个以数字开头的子文件夹分别叫2c_EconomicDensity、2b_DevelopmentCompactness、2a_DevelopmentIntensity。这个“2c→2a”的字母顺序有讲究经济密度是结果紧凑度是形态开发强度是过程入口。摘要里解释得很清楚经济密度度量单位面积经济产值紧凑度度量建成区集中程度开发强度则反映建设密度和人口密度。三者放在一起才能解释城市扩张的机理。3.1 三个shapefile的字段差异与命名启示我用ogrinfo或geopandas快速列出字段总能发现一些规律。通常EconomicDensity.dbf持有类似GDP_KM2、VAL1998这样的字段DevelopmentCompactness.dbf里有类似SHAPE_INDEX的值DevelopmentIntensity.dbf里则可能是BUILD_PCT或POP_KM2。虽然名字不同但数值都是区域级统计量。命名前缀“2c、2b、2a”更像是一个生产流水线先算出经济产出密度再评估形态紧凑度最后落到开发强度。这样设计的好处是属性表里可以不冗余存储几何信息靠FID关联。3.2 geopandas读取示例import geopandas as gpd # 注意cpg编码声明如果乱码就指定encodinggbk eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(eco_1998.columns) print(eco_1998.head()) print(eco_1998[[GDP_KM2, YEAR]].describe())这段代码先用read_file加载 Shapefileencoding参数匹配.cpg文件声明的字符集。然后打印所有列名和统计摘要。describe()会输出count、mean、std等字段能快速判断GDP_KM2是否存在异常负值或大范围空值。拿到属性表后最好顺带检查YEAR字段因为一个shp文件里可能存储多期数据靠YEAR区分。3.3 紧凑度与强度的计算口径紧凑度计算公式在文献里常用C 2√(πA)/P其中A是建成区面积P是周长结果落在0到1之间。本数据集里的DevelopmentCompactness如果直接是数值大概率就是这类形状指数或紧凑比。至于DevelopmentIntensity我见过更合理的字段是“单位面积建设用地面积”或“人口密度”的网格化统计结果。比较年份时要注意如果行政边界本身发生了调整比如2006年某镇并入街道则必须先使用1_studyarea里统一后的边界做裁剪或交集不然紧凑度会被区域面积变化干扰。4. 跨年份空间分析如何提取“开发密度”的变化信号三个年份的数据放一起最直接的分析是看同一行政单元在不同指标上的走势。但前提是边界一致。1_studyarea提供了研究区边界而各年份shp的几何可能来自不同来源必须先处理掉缝隙、重叠与坐标系差异。4.1 统一坐标系与属性对齐先检查三个文件的坐标系是否一致import geopandas as gpd study gpd.read_file(./1_studyarea/StudyArea.shp) eco_1998 gpd.read_file(./2c_EconomicDensity/EconomicDensity.shp, encodingutf-8) print(study.crs) print(eco_1998.crs) if eco_1998.crs ! study.crs: eco_1998 eco_1998.to_crs(study.crs)如果.prj文件缺失或定义不完整to_crs会抛错。应急方案是从一个已知正确的GeoDataFrame复制crs对象例如eco_1998.crs study.crs但这样做的前提是你确认两套数据的椭球和基准面一致否则会出现几米的偏移。属性对齐更麻烦。如果2006年的字段结构变了比如新增了“空间GDP”字段而1998年没有直接用pd.concat会生成大量NaN。我常用的做法是先只保留几个核心字段统一命名为region、year、value再按区域名进行merge或groupby。# 读取紧凑度数据并做空间连接 compact gpd.read_file(./2b_DevelopmentCompactness/DevelopmentCompactness.shp) joined gpd.sjoin(study, compact, howleft, predicateintersects) # 按城市名聚合 grouped joined.groupby(CITY_NAME)[COMPACT].mean().reset_index()代码中sjoin的predicateintersects表示两个要素只要任意部分相交就算匹配这能规避边界微小缝隙。groupby取均值是因为一个城市可能由多个面要素组成。如果希望面积加权则需要先计算相交面积joined[inter_area] joined.geometry.intersection(study.unary_union).area weighted joined.groupby(CITY_NAME).apply( lambda df: (df[COMPACT] * df[inter_area]).sum() / df[inter_area].sum() )注意apply里的加权公式当区域面积差异大时简单均值会偏向面积小但数值高的单元加权均值能反映整体状态。4.2 生成三年对比表把三个年份的经济密度整合到一张宽表可以直接用pivot_tableimport pandas as pd # 假设三个gdf都有YEAR和VALUE字段 eco_1998[year] 1998 eco_2006[year] 2006 eco_2012[year] 2012 frames [eco_1998[[CITY_NAME, year, GDP_KM2]], eco_2006[[CITY_NAME, year, GDP_KM2]], eco_2012[[CITY_NAME, year, GDP_KM2]]] all_df pd.concat(frames) pivot all_df.pivot_table(indexCITY_NAME, columnsyear, valuesGDP_KM2, aggfuncmean) pivot[growth_06_98] pivot[2006] / pivot[1998] - 1 pivot[growth_12_06] pivot[2012] / pivot[2006] - 1pivot_table的aggfunc默认mean它会处理同一城市下多个面要素的重复记录。用增长率列能直接识别哪几年扩张最快。4.3 变化可视化可视化可以先用matplotlib画折线趋势再用geopandas给每个城市填充增长率颜色。需要注意用颜色分级图时最好将增长率分位数分成5类避免极端值把配色拉平。5. 处理1km GDP栅格用rasterio读取、重投影与统计栅格部分是这份数据里最硬的骨头。NAGDP_1km_1998.tif这类文件是GeoTIFF除了.tif本体还有.tfw世界文件、.aux.xml辅助信息、.ovr金字塔。这套文件对新手不太友好但处理思路很固定先看元数据再处理空值最后做统计或重投影。5.1 栅格文件的组成.tif、.tfw、.ovr各干什么.tfw是文本世界文件记录像元尺寸和左上角坐标.ovr是金字塔文件如果缺失大范围缩放变慢但不影响数值精度.aux.xml包含统计信息和nodata标记。读取时我习惯用rasterio而不是gdal命令行因为能直接返回numpy数组import rasterio with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data src.read(1) nodata src.nodata profile src.profile print(nodata , nodata) print(数据类型 , data.dtype) print(shape , data.shape)src.read(1)读出第一波段nodata是无值标记。GDP数据通常用-9999或0表示无数据拿到数组后第一步要把这些值屏蔽掉否则计入均值会让结果严重偏低import numpy as np masked np.ma.masked_where(data 0, data) print(masked.mean())5.2 做年度差值并重投影三个年份的tif如果像元尺寸和原点不完全一致不要直接做数组减法。正确做法是先统一到同一个格网。可以用reproject完成from rasterio.warp import reproject, Resampling, calculate_default_transform with rasterio.open(./3_NAGDP_1km/NAGDP_1km_1998.tif) as src: data_1998 src.read(1) src_crs src.crs src_transform src.transform nodata src.nodata # 目标坐标系如果源是WGS84地理坐标这里用EPSG:4326 transform, width, height calculate_default_transform( src_crs, EPSG:32650, src.width, src.height, *src.bounds ) dst np.zeros((height, width), dtypenp.float32) reproject( sourcedata_1998, destinationdst, src_transformsrc_transform, src_crssrc_crs, dst_transformtransform, dst_crsEPSG:32650, resamplingResampling.bilinear, src_nodatanodata, dst_nodata-9999 )上述代码中calculate_default_transform根据源范围和目标投影计算新尺寸bilinear重采样适用于连续型GDP值如果处理分类栅格则改用nearest。src_nodata与dst_nodata保持一致保证重投影后空洞仍是同一标记。如果不想写这么长的代码优先用rioxarray它将坐标和投影封装得更好一行就能做差值import rioxarray ds_1998 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_1998.tif) ds_2006 rioxarray.open_rasterio(./3_NAGDP_1km/NAGDP_1km_2006.tif) diff ds_2006 - ds_1998但要先确认两个对象的空间参考一致若不一致则调用ds_1998.rio.reproject_match(ds_2006)对齐。5.3 高分辨率GDP数据的局限与坑1km格网GDP并不是真实观测值它一般基于夜间灯光、土地利用和统计年鉴做空间化插值因此“某个像元上的GDP”并不代表那个位置真有这么多产出。解读时更适合用在区域总量或相对比上。我在做行政单元统计时会先用行政区边界裁剪然后逐像元累加。这里最典型的坑有两个一是边缘像元被边界切成窄条如果连通面积小于一个像元直接累加会高估边界地带数值二是.ovr金字塔文件版本不一致可能让某些软件读到错误nodata。解决办法是用rasterstats的all_touchedFalse只保留完整落入边界的像元并在读取时显式指定maskedTrue。提示如果某一年份的GDP均值明显异常优先检查nodata值是否被默认为0再检查重投影时是否把边缘的nodata像元插值成小数值。6. 收尾技巧用rasterstats一键把1km GDP汇总到行政区前面的操作分别处理了矢量和栅格实际交付成果时常需要“一张Excel表”每一行是一个区县列是三个年份的GDP密度与总量。这时用rasterstats包可以省掉大量手写循环。它专门用来把栅格数据按矢量多边形进行分区统计pip install rasterstats然后写一个极简函数import pandas as pd from rasterstats import zonal_stats stats zonal_stats( ./1_studyarea/StudyArea.shp, # 矢量边界 ./3_NAGDP_1km/NAGDP_1km_1998.tif, # 栅格文件 stats[sum, mean, min, max], nodata-9999, # 与tif元数据一致 all_touchedFalse # 只统计完整覆盖的像元 ) df pd.DataFrame(stats)zonal_stats第一个参数接收矢量边界第二个参数接收栅格路径stats支持sum/mean/max/min/median等聚合统计。nodata参数能过滤栅格中的无值像元避免把-9999当作真实数值。all_touchedFalse表示只统计多边形中心点位于内部的像元这比True更严格也更能避免边缘噪声。三个年份只需包装成一个循环def extract_year(year): raster f./3_NAGDP_1km/NAGDP_1km_{year}.tif st zonal_stats( ./1_studyarea/StudyArea.shp, raster, stats[sum, mean], nodata-9999, all_touchedFalse ) return pd.DataFrame(st) result pd.concat( [extract_year(1998), extract_year(2006), extract_year(2012)], axis1 ) result.columns [GDP_sum_1998, GDP_mean_1998, GDP_sum_2006, GDP_mean_2006, GDP_sum_2012, GDP_mean_2012]这个封装直接生成一个DataFrame再配合study里的行政区名称列即可保存为CSVstudy gpd.read_file(./1_studyarea/StudyArea.shp) result.insert(0, city_name, study[CITY_NAME]) result.to_csv(./prd_gdp_summary.csv, indexFalse)一个小技巧如果发现某年的sum异常低先回到第5章检查nodata和重投影步骤或者换个思路用面积加权平均——因为栅格像元面积可能随纬度变化而不同尤其在投影坐标系里不同纬度像元实际面积并不完全一致。遇到这种情况可以先用area属性求出单位面积值再乘实际面积。这套流程跑通后后面再加2018年数据只需要复制这个函数模板替换年份即可。本文还有配套的精品资源点击获取