粤港澳大湾区shp数据全攻略：从坐标系到可视化

简介：这份粤港澳大湾区shp数据面向从事GIS制图、空间分析与区域规划的研究者、学生及从业者，提供可直接加载的矢量底图，解决大湾区跨城市边界数据获取零散、格式不统一的问题。资源包共20个文件，以shp、shx、dbf、prj等Shapefile核心组件为主，辅以sbx、sbn空间索引及xml元数据，压缩包约272KB，体量轻便，便于快速导入ArcGIS、QGIS等平台进行可视化与叠加分析。目前已有291人学习下载，适合需要快速搭建大湾区制图底图或开展空间统计的初中级用户。数据覆盖粤港澳大湾区整体范围及香港、澳门等局部线要素，属性表与投影信息完整，可直接用于专题图制作、行政区划展示与空间关系验证，省去自行搜集与坐标校正的环节，为区域研究提供规范、即用的基础地理数据支撑。

1. 拿到粤港澳大湾区 shp 数据后，先搞清楚它能干什么

做 GIS 可视化的人多半有过这种经历：底图配好了、配色调顺了，结果卡在边界数据上——要么行政区划对不上，要么海岸线缺一块，要么投影一叠加就飘到几百公里外。粤港澳大湾区 shp 数据解决的正是这个前置问题：它把广州、深圳、珠海、佛山、惠州、东莞、中山、江门、肇庆九个城市以及香港、澳门两个特别行政区的行政边界、海岸线、岛屿轮廓整理成标准矢量图层，直接拖进 QGIS、ArcGIS 或 GeoPandas 就能用。适合做区域经济热力图、交通可达性分析、城市扩张监测、专题地图出图的人。它不是一张图，而是一套可编辑、可裁剪、可叠加的几何底料，你后续所有空间统计和可视化都建立在它之上。

2. 拆开这份 shp 数据：图层结构、坐标系与字段设计

2.1 常见图层构成与属性字段

一套完整的粤港澳大湾区 shp 数据，通常不是一个孤零零的.shp文件，而是按行政层级拆成多个图层。常见做法是分成市级边界、区县级边界，有的还会附带海岸线和主要岛屿。每个图层对应一组同名文件：.shp存几何、.dbf存属性、.shx存索引、.prj存投影信息。少一个.prj，坐标系就得靠猜，这是新手最容易翻车的地方。

属性表里一般会有这些字段：

字段命名不同版本会有差异，但NAME和ADCODE基本是标配。ADCODE是后续做数据关联的关键——你要把统计年鉴的 GDP、人口数据挂到地图上，靠的就是它。

2.2 坐标系：WGS84 还是 CGCS2000

这是整个流程里最容易被忽略、后果最严重的一环。粤港澳大湾区 shp 数据常见的坐标系有两种：WGS84（EPSG:4326）和 CGCS2000（EPSG:4490）。两者在地图上看差别不大，但做面积计算或叠加其他数据时，差之毫厘谬以千里。

如果你只是做展示出图，WGS84 够用。但如果要算真实面积、做缓冲区分析，必须投影到平面坐标系。大湾区常用的投影是 UTM 49N（EPSG:32649）或 CGCS2000 3度带高斯克吕格投影。判断当前数据是什么坐标系，最直接的办法是看.prj文件内容，或者在 QGIS 里右键图层属性查看。

# 用 ogrinfo 快速查看 shp 的坐标系和字段信息 ogrinfo -al -so gba_city.shp # 输出里重点关注这几行： # Geometry: Polygon # Feature Count: 11 # Extent: (111.xxx, 21.xxx) - (115.xxx, 24.xxx) # Layer SRS WKT: GEOGCS["GCS_WGS_1984", ...]

ogrinfo是 GDAL 自带的命令行工具，-al表示列出所有要素，-so表示只输出摘要不输出几何。如果你看到Layer SRS WKT里写着WGS_1984，那就是地理坐标系，单位是度，不能直接用来算面积。

2.3 用 Python 读取并检查数据完整性

实际动手前，我习惯先用 GeoPandas 把数据读一遍，确认几何有效、字段完整、没有空值。

import geopandas as gpd # 读取市级边界 gdf = gpd.read_file("gba_city.shp") # 基本信息 print(f"要素数量: {len(gdf)}") print(f"坐标系: {gdf.crs}") print(f"字段列表: {list(gdf.columns)}") # 检查几何有效性 invalid = gdf[~gdf.geometry.is_valid] print(f"无效几何数量: {len(invalid)}") # 检查空几何 empty = gdf[gdf.geometry.is_empty] print(f"空几何数量: {len(empty)}") # 查看每个城市的面积（度为单位，仅作对比） gdf["area_deg"] = gdf.geometry.area print(gdf[["NAME", "area_deg"]].sort_values("area_deg", ascending=False))

这段代码做了四件事：确认要素数量是否符合预期（大湾区应该是11个城市）、检查坐标系是否已定义、排查无效几何、对比各城市面积量级。如果len(gdf)只有9，说明香港和澳门可能被单独放在另一个图层里。如果gdf.crs返回None，说明.prj文件缺失，需要手动指定坐标系。

几何无效是 shp 数据的常见病，表现为自相交、悬挂节点、重复顶点。GeoPandas 提供了buffer(0)这个经典偏方来修复大部分自相交问题：

# 修复无效几何 gdf["geometry"] = gdf.geometry.buffer(0) print(f"修复后无效几何: {len(gdf[~gdf.geometry.is_valid])}")

buffer(0)的原理是把多边形先向外缓冲0距离再收回，这个过程会重建拓扑关系，把自相交的环理顺。注意它不保证100%修复，极端情况还是得回 QGIS 用拓扑检查工具手动处理。

3. 从 shp 到可视化：投影转换、裁剪与出图全流程

3.1 投影转换：把度变成米

拿到地理坐标系的数据后，第一步是投影转换。不转投影就做缓冲区分析，等于用尺子量空气。

import geopandas as gpd gdf = gpd.read_file("gba_city.shp") # 如果 crs 为空，先手动指定为 WGS84 if gdf.crs is None: gdf = gdf.set_crs("EPSG:4326") # 投影到 UTM 49N（适合大湾区） gdf_proj = gdf.to_crs("EPSG:32649") # 验证：面积单位变成平方米 gdf_proj["area_km2"] = gdf_proj.geometry.area / 1e6 print(gdf_proj[["NAME", "area_km2"]].sort_values("area_km2", ascending=False)) # 保存投影后的数据 gdf_proj.to_file("gba_city_utm.shp", encoding="utf-8")

to_crs是 GeoPandas 的投影转换方法，参数可以是 EPSG 代码或 proj 字符串。UTM 49N 覆盖东经111°到114°，大湾区大部分在这个范围内，边缘区域变形可接受。如果你要做跨区域分析，可以考虑 Albers 等积投影。保存时加encoding="utf-8"是为了防止中文属性乱码，这个坑后面还会细说。

3.2 用 matplotlib 快速出图

投影转好后，出图就是几行代码的事。但要让图能看，还得调几个参数。

import matplotlib.pyplot as plt import geopandas as gpd gdf = gpd.read_file("gba_city_utm.shp") fig, ax = plt.subplots(1, 1, figsize=(12, 10)) # 绘制城市边界，按面积填色 gdf.plot( column="area_km2", cmap="YlOrRd", edgecolor="black", linewidth=0.8, ax=ax, legend=True, legend_kwds={"label": "面积（平方公里）", "orientation": "horizontal"} ) # 标注城市名 for idx, row in gdf.iterrows(): centroid = row.geometry.centroid ax.annotate( row["NAME"], xy=(centroid.x, centroid.y), ha="center", fontsize=9, fontproperties="SimHei" # 中文字体 ) ax.set_title("粤港澳大湾区各城市面积分布", fontsize=14) ax.set_axis_off() plt.tight_layout() plt.savefig("gba_area_map.png", dpi=300, bbox_inches="tight") plt.show()

column参数指定填色依据的字段，cmap是色带，legend=True自动生成图例。中文标注需要指定fontproperties，否则会显示成方框。set_axis_off()去掉坐标轴，出图更干净。dpi=300保证打印精度。

3.3 裁剪：只保留研究区范围

做具体项目时，往往只需要大湾区的一部分。比如研究深圳都市圈，就得把其他城市裁掉。

import geopandas as gpd from shapely.geometry import box gdf = gpd.read_file("gba_city_utm.shp") # 定义研究区范围（UTM坐标，单位米） # 这里以深圳为中心的大致范围 study_area = box(100000, 2480000, 250000, 2600000) # 用矩形裁剪 gdf_clip = gpd.clip(gdf, study_area) # 或者按属性筛选 gdf_sz = gdf[gdf["NAME"] == "深圳市"] # 保存 gdf_clip.to_file("gba_clip.shp", encoding="utf-8")

gpd.clip是几何裁剪，保留与裁剪框相交的部分。box的四个参数是(minx, miny, maxx, maxy)，单位跟数据坐标系一致。如果你只知道经纬度范围，得先转成投影坐标再裁剪。按属性筛选更简单，但要求字段值精确匹配，"深圳市"和"深圳"是两个结果。

3.4 叠加统计数据的完整链路

光有边界没意义，真正有用的是把统计数据挂上去。假设你有一份 CSV，包含各城市2023年GDP。

import pandas as pd import geopandas as gpd gdf = gpd.read_file("gba_city.shp") df = pd.read_csv("gdp_2023.csv") # 查看两边用于关联的字段 print(gdf[["NAME", "ADCODE"]].head()) print(df[["city_name", "gdp"]].head()) # 按城市名关联 merged = gdf.merge(df, left_on="NAME", right_on="city_name", how="left") # 检查是否有未匹配的记录 unmatched = merged[merged["gdp"].isna()] if len(unmatched) > 0: print(f"未匹配城市: {unmatched['NAME'].tolist()}") # 出图 merged.plot(column="gdp", cmap="Blues", legend=True, edgecolor="gray")

merge是 pandas 的操作，how="left"保证保留所有几何要素，匹配不上的 GDP 为空。关联字段两边名称可以不同，用left_on和right_on指定。关联后一定要检查未匹配记录，常见原因是城市名有空格、简称全称不一致、或者 CSV 里压根缺了某个城市。

4. 避坑指南：shp 数据处理的五个血泪教训

4.1 中文乱码：dbf 文件的编码陷阱

现象：用 GeoPandas 读取 shp 后，NAME字段显示为广å·å¸这样的乱码。

原因：shp 的.dbf文件默认用系统编码，Windows 下中文常是 GBK，而 GeoPandas 默认按 UTF-8 读取。

解决：读取时指定编码，或者用 QGIS 另存为 UTF-8。

# 方法一：读取时指定编码 gdf = gpd.read_file("gba_city.shp", encoding="gbk") # 方法二：用 ogr2ogr 转换编码 # ogr2ogr -f "ESRI Shapefile" output.shp input.shp -lco ENCODING=UTF-8

如果乱码已经发生且原始数据没备份，基本没有后悔药。我现在的习惯是拿到任何 shp 先ogrinfo看一眼编码，再决定怎么读。

4.2 坐标系丢失：.prj 文件缺失

现象：数据能读进来，但gdf.crs返回None，做任何空间操作都报错或结果离谱。

原因：拷贝数据时只复制了.shp、.shx、.dbf，漏掉了.prj。或者数据来源本身就没带投影文件。

解决：根据数据范围判断坐标系，手动指定。

# 如果确定是 WGS84 经纬度 gdf = gdf.set_crs("EPSG:4326") # 如果确定是 CGCS2000 gdf = gdf.set_crs("EPSG:4490") # 指定后再转投影 gdf_proj = gdf.to_crs("EPSG:32649")

判断依据：看坐标值范围。经纬度的话 x 在 111-115 之间，y 在 21-24 之间。如果是投影坐标，数值会是几十万到几百万。指定错了坐标系，后续所有分析都是错的。

4.3 几何无效：自相交与悬挂节点

现象：做intersection或union时抛出TopologyException，或者结果多边形出现奇怪的空洞。

原因：数据生产过程中数字化误差导致多边形自相交，或者相邻多边形共享边不完全重合。

解决：先用buffer(0)修复，再检查。

# 修复 gdf["geometry"] = gdf.geometry.buffer(0) # 如果还有问题，用 shapely 的 make_valid from shapely.validation import make_valid gdf["geometry"] = gdf.geometry.apply(make_valid)

make_valid是 shapely 1.8+ 提供的更严谨的修复方法，比buffer(0)更可靠，但速度慢一些。数据量大的话建议先buffer(0)批量处理，再对残余问题用make_valid。

4.4 面积计算单位错误

现象：算出广州面积是 0.12，明显不对。

原因：数据还是地理坐标系，geometry.area返回的是平方度，不是平方米。

解决：先投影再算面积。

# 错误做法 gdf = gpd.read_file("gba_city.shp") # WGS84 area_wrong = gdf.geometry.area # 单位是平方度 # 正确做法 gdf_proj = gdf.to_crs("EPSG:32649") area_correct = gdf_proj.geometry.area / 1e6 # 平方公里

平方度到平方公里的换算不是固定系数，因为不同纬度上1度对应的距离不同。所以必须投影，没有捷径。

4.5 文件体积过大导致读取缓慢

现象：区县级边界数据加载要等十几秒，出图卡顿。

原因：高精度边界包含大量冗余顶点，一个区县可能有几万个坐标点。

解决：做几何简化，在精度和体积之间取平衡。

# 简化几何，tolerance 单位与坐标系一致 gdf_simplified = gdf_proj.copy() gdf_simplified["geometry"] = gdf_proj.geometry.simplify(tolerance=100, preserve_topology=True) # 对比顶点数量 print(f"原始顶点数: {sum(len(g.exterior.coords) for g in gdf_proj.geometry)}") print(f"简化后顶点数: {sum(len(g.exterior.coords) for g in gdf_simplified.geometry)}")

tolerance=100表示100米，对于市级展示够用。preserve_topology=True防止简化后多边形自相交。如果是做精细分析，tolerance 要调小到10米甚至更小。

5. 进阶技巧：用 shp 数据做空间连接与热力图

5.1 空间连接：把点数据聚合到面

假设你有一批企业注册地址的经纬度点，想统计每个城市的企业数量。这就是典型的点面空间连接。

import geopandas as gpd import pandas as pd from shapely.geometry import Point # 读取城市边界 cities = gpd.read_file("gba_city.shp").to_crs("EPSG:4326") # 模拟企业点数据 points_df = pd.DataFrame({ "name": ["企业A", "企业B", "企业C"], "lng": [113.26, 114.06, 113.55], "lat": [23.13, 22.54, 22.19] }) points = gpd.GeoDataFrame( points_df, geometry=[Point(xy) for xy in zip(points_df["lng"], points_df["lat"])], crs="EPSG:4326" ) # 空间连接：每个点落在哪个城市 joined = gpd.sjoin(points, cities, how="left", predicate="within") # 统计各城市企业数 count = joined.groupby("NAME").size().reset_index(name="企业数量") print(count)

sjoin的predicate="within"表示点在城市多边形内部。how="left"保留所有点，落在边界外的点城市名为空。如果点正好落在两个城市交界线上，可能匹配到两个城市，需要去重。

5.2 热力图：用点密度反映空间聚集

有了点数据，热力图比单纯的分级填色更能反映聚集特征。

import matplotlib.pyplot as plt import geopandas as gpd import numpy as np from scipy.stats import gaussian_kde cities = gpd.read_file("gba_city_utm.shp") points = gpd.read_file("enterprises.shp").to_crs("EPSG:32649") # 提取坐标 x = points.geometry.x.values y = points.geometry.y.values # 核密度估计 kde = gaussian_kde(np.vstack([x, y])) xi, yi = np.mgrid[x.min():x.max():500j, y.min():y.max():500j] zi = kde(np.vstack([xi.flatten(), yi.flatten()])).reshape(xi.shape) # 出图 fig, ax = plt.subplots(figsize=(12, 10)) cities.plot(ax=ax, facecolor="none", edgecolor="gray", linewidth=0.5) ax.contourf(xi, yi, zi, levels=15, cmap="hot_r", alpha=0.7) ax.set_axis_off() plt.savefig("heatmap.png", dpi=300, bbox_inches="tight")

gaussian_kde是 scipy 的核密度估计，500j表示生成500个等间距点。contourf画填充等值线，alpha=0.7让底图边界透出来。带宽参数bw_method可以调整平滑程度，默认是 Scott 规则，点密集时可能过平滑，可以手动设小一点。

5.3 验证数据质量的三个习惯

从那以后我每次拿到新的 shp 数据，都强制走一遍这三步：先ogrinfo看坐标系和字段，再is_valid查几何，最后用total_bounds确认范围是否合理。这三步花不了两分钟，但能省下后面几小时的排查时间。尤其是total_bounds，如果返回的坐标范围跟预期差了一个数量级，那肯定是坐标系出了问题，后面所有操作都白搭。

希望帮到你。

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