企业官网建设流程全解析

简介：珠江三角洲城市群区域开发密度数据集面向城市规划、经济地理与区域可持续发展研究者，汇集1998年、2006年、2012年三个时点的地理空间信息，可用于揭示近二十年间城市空间扩展、经济集聚与开发强度的时空演变规律。压缩包共44个文件、约10.83MB，以Shapefile矢量格式（.shp、.shx、.dbf、.prj等）和GeoTIFF栅格数据（.tif、.tfw、.ovr）为主，分层存储研究区域边界、经济密度、发展紧凑度、发展强度以及1公里网格市域GDP，在ArcGIS、QGIS等常用GIS平台中可直接打开和横向比较。已有74人浏览学习。借助其中包含的城市行政边界、河流海岸线、GDP空间分布、建成区紧凑度、建筑覆盖率与人口密度等指标，研究者能逐项对比不同年份的空间结构变化，识别经济热点与城乡差异，亦可延伸至城市可持续性评估、交通规划、环境保护及社会经济不平等分析，为区域发展政策制定提供量化支撑。

1. 这份珠三角密度数据集为什么值得拆开看

拿到一个名为“珠江三角洲城市群区域开发密度数据集（1998，2006，2012）.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 "珠江三角洲城市群区域开发密度数据集（1998，2006，2012）.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是元数据。

值得注意的是DevelopmentIntensity文件夹内，.cpg被写成了大写.CPG。在Linux下，文件系统区分大小写，如果从Windows复制到Linux，geopandas自动搜索.cpg时可能找不到，导致属性编码识别失败。我遇到这种情况时，一般会用一个for循环把所有文件改为小写扩展名，或者直接在读取时显式指定encoding。

2.3 先看投影：.prj与.tfw告诉你在哪种坐标下

投影是一切空间计算的基准。用文本编辑器打开StudyArea.prj，里面是一段WKT字符串。珠三角区域常见的坐标系有WGS84地理坐标、WGS84 UTM 50N投影，以及国家2000坐标系，具体要看文件内容。我判断投影时重点关注三个信息：基准面（datum）、投影名称（projection）、单位（units）。地理坐标系单位为度，投影坐标系单位为米。NAGDP_1km_1998.tfw是GeoTIFF的世界文件，用普通文本编辑器打开就能看到六行数字，前两行是像元尺寸和旋转项，后两行是栅格原点。更直观的方式是用gdalinfo：

gdalinfo NAGDP_1km_1998.tif | head -n 40

gdalinfo输出里会显示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编码声明，如果乱码就指定encoding='gbk' eco_1998 = gpd.read_file('./2c_EconomicDensity/EconomicDensity.shp', encoding='utf-8') print(eco_1998.columns) print(eco_1998.head()) print(eco_1998[['GDP_KM2', 'YEAR']].describe())

这段代码先用read_file加载 Shapefile，encoding参数匹配.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', encoding='utf-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, how='left', predicate='intersects') # 按城市名聚合 grouped = joined.groupby('CITY_NAME')['COMPACT'].mean().reset_index()

代码中sjoin的predicate='intersects'表示两个要素只要任意部分相交就算匹配，这能规避边界微小缝隙。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_table：

import 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(index='CITY_NAME', columns='year', values='GDP_KM2', aggfunc='mean') pivot['growth_06_98'] = pivot[2006] / pivot[1998] - 1 pivot['growth_12_06'] = pivot[2012] / pivot[2006] - 1

pivot_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), dtype=np.float32) reproject( source=data_1998, destination=dst, src_transform=src_transform, src_crs=src_crs, dst_transform=transform, dst_crs='EPSG:32650', resampling=Resampling.bilinear, src_nodata=nodata, 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_touched=False只保留完整落入边界的像元，并在读取时显式指定masked=True。

提示：如果某一年份的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_touched=False # 只统计完整覆盖的像元 ) df = pd.DataFrame(stats)

zonal_stats第一个参数接收矢量边界，第二个参数接收栅格路径，stats支持sum/mean/max/min/median等聚合统计。nodata参数能过滤栅格中的无值像元，避免把-9999当作真实数值。all_touched=False表示只统计多边形中心点位于内部的像元，这比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_touched=False ) return pd.DataFrame(st) result = pd.concat( [extract_year(1998), extract_year(2006), extract_year(2012)], axis=1 ) result.columns = ['GDP_sum_1998', 'GDP_mean_1998', 'GDP_sum_2006', 'GDP_mean_2006', 'GDP_sum_2012', 'GDP_mean_2012']

这个封装直接生成一个DataFrame，再配合study里的行政区名称列即可保存为CSV：

study = gpd.read_file('./1_studyarea/StudyArea.shp') result.insert(0, 'city_name', study['CITY_NAME']) result.to_csv('./prd_gdp_summary.csv', index=False)

一个小技巧：如果发现某年的sum异常低，先回到第5章检查nodata和重投影步骤，或者换个思路用面积加权平均——因为栅格像元面积可能随纬度变化而不同，尤其在投影坐标系里，不同纬度像元实际面积并不完全一致。遇到这种情况，可以先用area属性求出单位面积值再乘实际面积。这套流程跑通后，后面再加2018年数据只需要复制这个函数模板，替换年份即可。

本文还有配套的精品资源，点击获取