成渝城市群矢量数据实战：shp体检、POI空间连接与坐标处理

简介：资源聚焦成渝城市群空间数据应用场景，面向地理信息系统、城市规划与科研项目中需要学校、医院等POI及道路、建筑轮廓等基础矢量数据的开发者和研究人员。数据年限为2019年，POI信息经网络地图爬取并完成导出、裁剪等处理，统一采用WGS84椭球投影；建筑、铁路与道路数据则源自OpenStreetMap，同样经过爬取与裁剪，覆盖建筑轮廓、铁路和道路三类要素。压缩包内共155个文件，以shp、dbf、shx、prj等矢量数据格式为主，其中shp用于空间几何信息，dbf保存属性字段，prj定义坐标参考，整体包体约40.3MB，可直接导入ArcGIS或QGIS等工具使用。已有1191人学习下载，资源包含医院、幼儿园、中小学、高校、专科医院等细分POI层级以及建筑物和路网数据，适合作为城市群分析、设施可达性评估与制图底图等项目的原始数据支撑。

1. 成渝城市群矢量数据：底图、路网、POI 一次给全

做城市群分析时最耗时间的永远是拼数据底图。我最早做成渝城市群基础设施评估时，行政边界、道路、建筑轮廓、POI 分了好几个来源下载，结果坐标系对不上、路网点位跑到边界外，光修数据就修了三天。这套成渝城市群矢量数据把四类图层打包成了 shp：市/区县边界、道路、建筑轮廓、POI 数据齐全，下载后直接丢进 ArcGIS 或 QGIS 就能用。对做空间分析、选址评估、城市群可视化的人尤其省事——你不需要再到处找底图和路网，也不用自己拼接多源数据。做 WebGIS 开发的人也可以拿它当测试数据，验证加载性能、跑空间查询都够用。

2. 数据体检先行：坐标系、图层结构与属性表预检

任何 shp 拿到手，第一件事都不是急着做分析，而是先体检。坐标系错位、字段乱码、缺文件是 shp 数据最常见的三座大山。这一章我把拆包后的检查流程写清楚，跟着走一遍就能避掉八成低级问题。

2.1 拆包先看图层与配套文件

一个完整的 shapefile 不是单个.shp文件，而是一组配套文件：.shp存几何、.dbf存属性、.shx存索引，还有可选的.prj存坐标系、.cpg存编码。下载解压后如果发现.dbf或.prj缺失，这个数据基本不能直接信任——几何和属性对不上，坐标系也成了黑匣子。

我在命令行里常用 Python 快速列一下包里的内容：

from pathlib import Path data_dir = Path("./chengyu_shp") for shp in data_dir.rglob("*.shp"): # 统计配套文件是否齐全 shx = shp.with_suffix(".shx").exists() dbf = shp.with_suffix(".dbf").exists() prj = shp.with_suffix(".prj").exists() print(f"{shp.name} | shx={shx} dbf={dbf} prj={prj}")

shx是几何索引，少了它很多 GIS 软件打不开；dbf保存属性字段，少了它图层只剩几何；prj记录坐标系，少了它后面对接必翻车。我一般的习惯是拆包后先跑这段脚本，缺文件的先补，补不了的就换数据源。这一步不花时间，但能省下后面好几个小时的排查时间。

2.2 坐标系：先统一再叠加是硬前提

成渝城市群这套数据里，POI 点位大概率是经纬度坐标（WGS84 或 CGCS2000），而行政区边界如果做过投影，坐标值就会变成以米为单位的平面坐标。这两类数据直接叠加，点位会跑到完全无关的位置，而且很难通过肉眼察觉——因为范围看起来是“大致对得上”的，细看才发现全都偏了。

用 Python 检查坐标系是最直接的方式：

import geopandas as gpd poi = gpd.read_file("chengyu_poi.shp", encoding="utf-8") boundary = gpd.read_file("chengyu_boundary.shp", encoding="utf-8") print("POI 坐标系:", poi.crs) print("边界坐标系:", boundary.crs)

常见坐标系就三种：

如果边界是用 CGCS2000 投影坐标系做的，而 POI 是 WGS84 经纬度，叠加前一定要统一。我一般统一到 CGCS2000 的投影坐标系做空间计算，出图展示时再转 Web Mercator。注意boundary自带的.prj文件里写的是什么就用什么，不要凭文件名猜。

2.3 属性表预检：编码、字段与空值率

shp 的属性表存在.dbf文件里，编码问题是最常见的坑。国内很多数据用 GBK 或 CP936 编码保存，而geopandas默认按 UTF-8 读，结果就是 POI 名称显示成乱码。先读一下字段结构再决定后续处理：

gdf = gpd.read_file("chengyu_poi.shp", encoding="utf-8") print(gdf.columns.tolist()) print(gdf.head(3)) print(gdf.isnull().sum())

encoding="utf-8"是显式指定 dbf 的读取编码，如果你的数据显示乱码，换成encoding="gbk"再试一次。字段列表决定你能做什么分析——比如 POI 图层有没有“大类/小类”字段，边界图层有没有“市/区县”名称字段，这直接影响后面的空间挂接和统计。空值率则告诉你数据质量：如果一个字段空了大半，那统计类分析基本别指望它。

属性表这步特别值得认真看。字段名规范、编码正确、空值可控的图层，后面每一步都顺畅；反之，你在半路发现字段不对，回头改编码再读一遍，等于所有下游分析全部重跑。

3. 把 POI 挂到区县边界：空间连接与统计输出

POI 数据的核心价值是“落到空间里”。单独看一个 POI 点位没有意义，把它挂到行政区边界里，统计每个区县有多少设施、什么类型占主导，这才是城市群分析的基本功。这一章从清洗、空间连接到按区县统计，给出完整可复现的操作。

3.1 清洗 POI：去重、补空、归一类别

原始 POI 数据里最常见的问题是重复和空字段。同一个门店可能被采集了两条记录，或者“名称”字段有值但“类别”字段空着。不洗直接用，统计结果会虚高。

我的做法分三步：先按“名称+地址+经纬度”去重，再处理空值，最后把类别做归一。

import pandas as pd import geopandas as gpd from shapely.geometry import Point # 读原始 POI 表（可能是 Excel 导出） poi_df = pd.read_excel("poi_raw.xlsx") # 去重：名称+地址完全一致则保留第一条 poi_df = poi_df.drop_duplicates(subset=["名称", "地址", "经度", "纬度"]) # 空值处理：没有类别的补“未知” poi_df["类别"] = poi_df["类别"].fillna("未知") # 把 DataFrame 转成 GeoDataFrame poi_gdf = gpd.GeoDataFrame( poi_df, geometry=[Point(x, y) for x, y in zip(poi_df["经度"], poi_df["纬度"])], crs="EPSG:4326" )

Point(x, y)的经纬度顺序别搞反——x是经度、y是纬度。crs="EPSG:4326"是显式指定 WGS84 经纬度，但如果你的原始坐标用的是 GCJ-02 火星坐标，这里指定 EPSG:4326 会造成几百米偏移，属于坐标系的另一个故事了。去重逻辑里，经纬度完全一致才能判定为同一点，坐标差一位小数点都不建议合并，宁可保留两条，也不误删真实存在的不同门店。

3.2 空间连接：sjoin 让点位落在边界内

POI 坐标是经纬度点，区县边界是多边形。要把点位挂到所在的区县，核心操作是空间连接。做过一次就知道，这个操作的本质是“对每一个点，找到包含它的那个多边形”，性能和数据量直接挂钩——边界图层有几万个面时，全量相交会慢到让人怀疑电脑当机。

boundary = gpd.read_file("chengyu_boundary.shp", encoding="utf-8") # 统一坐标系再做连接 boundary = boundary.to_crs("EPSG:4326") # 空间连接 joined = gpd.sjoin( poi_gdf, boundary, how="inner", predicate="within" )

两个参数值得展开说。predicate="within"要求点必须完全落在面内，落在边界线上或者面外的点都会被丢弃；如果你希望保留所有点、没有匹配到的标记为 NaN，就用how="left"。我的习惯是先跑一次how="inner"看看命中率，如果丢掉的点超过 10%，多半是坐标偏移或边界图层缺少辖区范围，先排查再继续。

joined结果里会带上边界图层的所有字段——尤其是区县名称字段。这时候你可以直接验证一个点位是否挂对了地方：随便抽几条记录检查“名称”和“区县名”是否合理匹配。这一步是手动的，但非常必要，别把空间连接当作可以盲信的黑匣子。

3.3 按区县统计并导出制图

连接完成之后，按区县分组统计是顺手的事：

result = joined.groupby(["区县名"]).size().reset_index(name="poi_count") result = result.sort_values("poi_count", ascending=False) result.to_csv("poi_count_by_district.csv", index=False, encoding="utf-8-sig") print(result.head(10))

size()统计每个区县有多少个 POI 点。encoding="utf-8-sig"值得单独说明：带 BOM 的 UTF-8 编码，Excel 直接双击打开才不乱码。用默认的utf-8写出来，Excel 打开中文列名全是乱码——这是最容易“在最后一步翻车”的地方。

导出的 CSV 可以拖进 Tableau 或者直接在地图软件里做专题图。如果要按类别细分，把groupby的字段改成["区县名", "类别"]就能得到每个区县不同 POI 类别的计数。再往后，你可以用这些统计结果做设施密度分析、人均设施覆盖评估，或者叠加道路数据做可达性计算——POI 挂接这一步做完，下游分析的路就通了。

4. 道路与建筑轮廓二次加工：拓扑修复、裁剪与简化

行政区边界、POI 这些数据拿来就能用的情况比较多，道路和建筑轮廓则基本都需要二次加工。原因很简单：一套打包好的 shp 数据覆盖范围大，直接加载进工程文件会卡顿，拓扑错误也多。真实项目里几乎不会全量使用，而是按区域裁剪再局部处理。

4.1 道路拓扑：悬挂点与伪节点排查

道路网络的拓扑错误主要分两类：悬挂点（道路没连到其他道路、变成断头）和伪节点（一段完整道路被无意义地打断成多段）。这些错误在做网络分析时是致命的——最短路径计算可能因为断头路直接算不出结果。

在 QGIS 里跑 GRASS 的清理工具是最经济的方式：

v.clean input=chengyu_road output=chengyu_road_clean \ tool=rmdupl,break,rmdangle \ threshold=0.01

rmdupl删除重复线段，break在交叉处打断道路，rmdangle删除短于阈值的悬挂线。threshold=0.01的单位取决于数据坐标系——如果数据是经纬度，0.01 度约等于 1 公里，这个阈值会误删大量短线；如果你处理的是平面坐标系，这个值通常设置在 1~5 米比较合理。我的血泪经验是：先确认坐标系单位，再设阈值，否则清理完发现路网缺了一大块，后悔药都没得吃。

清理完成不等于万事大吉。我一直保留一个习惯：清理后随机抽几条道路手动检查连通性。工具只能处理规则错误，路网逻辑上的断裂它看不见。

4.2 按区县边界裁剪出你要的片区

成渝城市群范围很大，全量加载道路和建筑数据会拖慢分析效率。大部分项目只需要其中的几个区县，这时裁剪比筛选字段更可靠——裁剪直接按空间范围切片。

import geopandas as gpd road = gpd.read_file("chengyu_road_clean.shp", encoding="utf-8") boundary = gpd.read_file("chengyu_boundary.shp", encoding="utf-8") # 选出你关心的区县 target = boundary[boundary["区县名"].isin(["渝中区", "锦江区"])] # 按目标区县裁剪道路 road_clip = gpd.clip(road, target)

执行裁剪前必须确认两个图层的坐标系完全一致，否则clip的结果可能是一堆空数据或者残缺几何。target直接用了边界图层里的区县名做筛选，如果字段名不叫“区县名”，要改成实际字段名。裁剪完成后检查一下road_clip的要素数量，如果比预期少太多，去看是不是坐标系不一致导致的空间范围错位。

4.3 建筑轮廓简化：精度换顺畅

建筑轮廓通常细节极多，一个城市可能有几十万个面。直接整个加载到地图里，卡顿是必然的。取舍之道是先投影、后简化、再出图。

# 先投影到平面坐标系，得到米的单位 building = building.to_crs("EPSG:3857") # 简化几何：容差 5 米，保留拓扑结构 building_simple = building.simplify(tolerance=5, preserve_topology=True)

tolerance=5的含义是：简化后的几何与原始几何的最大偏差不超过 5 米。值设得越大、文件越小、细节丢失越多。在经纬度坐标系下直接做simplify是常见误用——单位是度，tolerance=5代表 5 度，基本等于把建筑轮廓画成一个点。先投影成米制坐标再简化，这个顺序不能反。

preserve_topology=True也很关键，它保证简化过程不会产生自相交、镂空等拓扑错误。设为False性能更快，但可能出现几何错误。对建筑轮廓这种密集面数据，我一般优先保拓扑，速度慢一点没关系。

5. 常见问题与排查：坐标错位、属性乱码与拓扑破洞

shp 数据用多了，翻车场景来来回回就那么几个。这一章把高频问题按“现象 → 原因 → 解决”梳理出来，你可以直接对照排查。

5.1 坐标错位：数据“飘”出边界外

现象：POI 点叠到区县边界图层上，点位整体偏向一侧几百米，有些点直接落在边界外。

原因：最常见的两种情况。一是 POI 用了 GCJ-02 火星坐标，而边界是 WGS84 或 CGCS2000，两者差几百米；二是两个图层的坐标系没有统一，直接把经纬度和投影坐标叠加了。

解决：先看两个图层的crs是否一致。GCJ-02 的问题没有完美的数学转换，常见做法是用公开的纠偏算法转成 WGS84 再做分析。坐标系不统一的话，用to_crs()统一到同一个 EPSG 代码。我拿到新数据的第一反应永远是先打印crs，这已经成了肌肉记忆。

5.2 属性表中文乱码：读出来全是“口口”

现象：POI 名称、区县名等中文字段在 QGIS 里显示正常，但用 Python 读取后全是“口口”或奇怪符号。

原因：dbf 文件的编码是 GBK，而你读取时用了 UTF-8。国内很多 shp 数据在制作时没有写.cpg文件，软件只能靠猜。

解决：

# 优先尝试 GBK gdf = gpd.read_file("chengyu_poi.shp", encoding="gbk")

如果还是乱码，就逐个试gb18030、utf-8。定下来之后，顺手把数据转成 GeoJSON 保存一份，GeoJSON 默认 UTF-8 编码，以后不会再出问题。

5.3 拓扑破洞：裁剪结果缺块

现象：用边界裁剪道路或建筑时，结果里有些区域明显缺失，或者生成的面要素中间多了个空洞。

原因：输入几何存在自相交、重合点等拓扑错误，裁剪工具遇到这种错误直接跳过或生成错误结果。

解决：裁剪前先做一次修复：

# buffer(0) 是最常用的拓扑修复方式 boundary_clean = boundary.buffer(0) road_clip = gpd.clip(road, boundary_clean)

buffer(0)会重建几何的边界关系，消除自相交和重叠。处理后的数据可以再跑一次is_valid检查：

print(boundary_clean.is_valid.all())

输出True说明拓扑没有问题，再去裁剪就稳了。

5.4 中文路径导致的“找不到文件”

现象：代码和文件都没问题，但read_file报错找不到文件，或者 ArcGIS 直接打不开图层。

原因：部分 GIS 组件和 Python 库对中文路径支持不友好，尤其是历史版本的 ArcGIS 工具箱。路径带“重庆”“成都”这类中文字符时，底层 C++ 组件读不了。

解决：把整套 shp 数据复制到纯英文路径下处理，比如D:\gisdata\chengyu\。这个操作很玄学，但确实是最稳的解法。项目交付时我会刻意要求所有数据路径只用英文和数字，从源头断掉这个坑。

6. 按需导出：shp 转 WKT、GeoJSON 与 3DTiles

shp 是 GIS 桌面端的通用格式，但换到 Web 端、数据库或者三维场景里就需要转格式了。常见需求是把 shp 导出成 WKT 文本给 PostGIS 用，或者转 GeoJSON 给前端加载。我一般按用途选转换方式——文本交换用 WKT，Web 可视化用 GeoJSON，三维场景用 3DTiles。

把 shp 转成 WKT 可以直接拿属性加几何拼文本：

gdf = gpd.read_file("chengyu_poi.shp", encoding="utf-8") with open("poi.wkt", "w", encoding="utf-8") as f: for _, row in gdf.iterrows(): f.write(f"{row.geometry.wkt}\t{row['名称']}\n")

每行一个点，几何和名称用制表符分隔。row.geometry.wkt是 Shapely 提供的 WKT 字符串输出，导入 PostgreSQL 的 PostGIS 时可以直接用ST_GeomFromText读取。转 GeoJSON 更省事：

gdf.to_file("chengyu_poi.geojson", driver="GeoJSON")

GeoJSON 自带坐标系信息，Web 端用 MapLibre 或 Leaflet 加载非常顺畅。三维场景里把 shp 转 3DTiles 是另一类需求——我通常先把数据按渔网分割成小块再单独转换，避免一次转换整个城市导致生成的文件过大加载崩溃。工具上可以选支持批处理的转换器，转换前确认源数据是 WGS84，因为 3DTiles 的标准坐标系就是 WGS84。

从那以后我每次导出 shp 之前，都强制自己先确认三件事：坐标系是什么、编码是什么、目标格式要什么。这三件事确认完，基本不会再出现导出后数据对不上的情况。希望帮到你。

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