南京矢量数据包实战：SHP图层解析、坐标系处理与空间分析全流程

简介：这份南京市区县乡镇街道及交通设施矢量数据包，面向城市规划、地理信息研究与GIS教学实训场景，解决本地基础空间数据获取零散、图层不全的问题。资源覆盖南京市全域，包含市、区、县、镇、村多级行政边界，国道、省道、县道、高速、地铁轻轨、铁路及乡村道路等线状交通要素，并整合水系、绿地、岛屿等面状自然要素，以及公交站、地铁出入口、停车场、收费站、学校、医院、银行、超市、餐饮、药店、宾馆酒店、政府机构、公园、大厦等二十余类点状POI。压缩包共184个文件，以shp、shx、dbf、prj、xml等标准SHP配套文件为主，另有mxd工程、py脚本与说明文本，整体约11.17MB，坐标系为WGS84或CGCS2000，投影信息完整，可直接在ArcGIS、QGIS中叠加分析与制图。目前已有84人学习下载，适合需要快速搭建南京空间分析底图、开展城市研究与课程实训的读者。

1. 南京市区县乡镇街道及交通设施矢量数据包：从一份SHP到一套可用的空间底图

如果你做过南京相关的空间分析、地图可视化或者选址评估，大概率遇到过这种局面：底图能拿到，但行政边界只到区一级，乡镇街道这一层要么缺失，要么年份对不上，道路和水系更是各说各话。一份包含区县、乡镇街道、道路、水系、POI 的 SHP 矢量数据包，解决的正是这个断层——它把多层级行政单元和交通设施要素统一到同一套坐标体系下，让你不用再花两三天做数据对齐。这类数据包适合做城市研究、物流路径预判、商业网点覆盖分析，也适合需要快速搭建本地空间数据库的开发者。但拿到压缩包只是起点，真正决定成败的是坐标系、字段结构和拓扑质量这三件事。

2. 先搞清楚数据包里到底有什么：图层清单与字段结构

2.1 五类核心图层的用途与取舍

一份典型的南京矢量数据包，通常按要素类型拆成多个 SHP 文件。SHP 格式本身是 Esri 定义的一种矢量数据交换格式，一个完整图层至少包含 .shp、.shx、.dbf 三个文件，缺一不可。下面这张表是我拿到任何一份行政+交通数据包时，第一遍会核对的图层清单：

这里有个容易翻车的点：乡镇街道边界在不同数据源里，有的按民政区划，有的按实际管理范围，两者在城乡结合部差异明显。我一般会先看字段里有没有“区划代码”和“统计用代码”两套编码，如果有，说明数据制作者考虑过口径问题；如果只有名称，那就要警惕同名街道跨区的情况。

2.2 用 QGIS 快速检查图层完整性

拿到数据后不要急着导入数据库，先用 QGIS 做一次目视检查。QGIS 是开源桌面 GIS，对 SHP 的支持最直接。打开后把 .shp 文件拖进去，重点看三件事：坐标系是否统一、几何有没有明显偏移、属性表字段是否完整。

# 用 geopandas 批量读取并检查坐标系与字段 import geopandas as gpd import os data_dir = "./nanjing_shp" layers = ["districts", "townships", "roads", "water", "poi"] for name in layers: path = os.path.join(data_dir, f"{name}.shp") if not os.path.exists(path): print(f"[缺失] {name}.shp 不存在") continue gdf = gpd.read_file(path) print(f"[{name}] 要素数={len(gdf)}, CRS={gdf.crs}, 字段={list(gdf.columns)}") # 检查几何有效性 invalid = (~gdf.geometry.is_valid).sum() if invalid > 0: print(f" 警告：{invalid} 个无效几何")

这段代码的逻辑很直白：逐个读取图层，打印要素数量、坐标参考系和字段名，同时统计无效几何。参数上，data_dir换成你解压后的实际路径；layers列表按压缩包内文件名调整。如果某个图层 CRS 显示为 None，说明缺少 .prj 文件，这时候不要硬猜坐标系，先找数据说明。

提示：SHP 的 .prj 文件只记录坐标系定义，不记录投影转换参数。如果 .prj 丢失，即使你知道是 CGCS2000，也需要确认是地理坐标还是投影坐标。

2.3 坐标系判断：为什么这一步不能省

南京地区常见坐标系有三种：WGS84 地理坐标（EPSG:4326）、CGCS2000 地理坐标（EPSG:4490）、以及各类高斯克吕格投影坐标。很多数据包为了兼容 Web 地图，会直接给 WGS84；但做面积统计和距离计算时，地理坐标的单位是度，直接算面积会得到荒谬结果。

判断方法很简单：看坐标数值范围。如果经纬度数值在 118~119 和 31~33 之间，那是地理坐标；如果数值是六位数甚至七位数，那是投影坐标。我一般会统一转成 EPSG:4326 做展示，转成 EPSG:3857 或当地投影做量算。转换用 geopandas 一行搞定：

# 统一转为 WGS84 地理坐标，便于后续叠加 gdf_roads = gpd.read_file("./nanjing_shp/roads.shp") if gdf_roads.crs is None: gdf_roads = gdf_roads.set_crs("EPSG:4490") # 按数据说明设定 gdf_roads_wgs = gdf_roads.to_crs("EPSG:4326") gdf_roads_wgs.to_file("./output/roads_wgs84.shp", encoding="utf-8")

set_crs是“声明”坐标系，to_crs是“转换”坐标系，两者不能混用。如果数据本身有 .prj，直接to_crs即可；如果没有，先用set_crs声明再转换。这个顺序搞反，数据会飞到几千公里外。

3. 把 SHP 用起来：入库、裁剪与路网分析的最小流程

3.1 导入 PostGIS 做空间查询

SHP 适合交换，不适合频繁查询。要把数据用起来，常见做法是导入 PostGIS。PostGIS 是 PostgreSQL 的空间扩展，支持空间索引和空间函数。导入用 shp2pgsql 命令行工具最稳：

# 将 roads.shp 导入 PostGIS，SRID 按实际坐标系填写 shp2pgsql -s 4326 -I -W UTF-8 -D ./nanjing_shp/roads.shp public.nj_roads | \ psql -h localhost -U postgres -d gis_db

参数说明：-s 4326指定 SRID；-I创建空间索引；-W UTF-8解决中文乱码；-D使用 dump 模式，比默认的 insert 模式快很多。导入后可以用一条 SQL 查某个区县内的道路长度：

-- 统计某区县内各等级道路总长度 SELECT r.road_class, SUM(ST_Length(r.geom::geography)) / 1000 AS km FROM nj_roads r JOIN nj_districts d ON ST_Intersects(r.geom, d.geom) WHERE d.name = '某区' GROUP BY r.road_class ORDER BY km DESC;

这里::geography是把几何转为地理类型，让长度单位变成米。如果直接用ST_Length(geom)，在 4326 下得到的是度，没有意义。ST_Intersects配合空间索引，比用属性字段关联快一个量级。

3.2 按行政边界裁剪出研究区

做具体项目时，往往只需要一个区或几个街道。用 geopandas 做裁剪比在 GIS 软件里点选更可复现：

import geopandas as gpd districts = gpd.read_file("./nanjing_shp/districts.shp") roads = gpd.read_file("./nanjing_shp/roads.shp") poi = gpd.read_file("./nanjing_shp/poi.shp") # 统一坐标系后再裁剪 target_crs = districts.crs roads = roads.to_crs(target_crs) poi = poi.to_crs(target_crs) study_area = districts[districts["name"] == "某区"] roads_clip = gpd.clip(roads, study_area) poi_clip = gpd.clip(poi, study_area) roads_clip.to_file("./output/roads_study.shp", encoding="utf-8") poi_clip.to_file("./output/poi_study.shp", encoding="utf-8")

gpd.clip会保留落在研究区内的几何，边界处自动切断。注意裁剪前必须统一 CRS，否则 clip 结果为空或报错。如果研究区是多个不相邻的街道，先把它们合并成一个 GeoDataFrame 再裁剪。

3.3 路网连通性检查与最短路径

道路数据用于分析时，最大的坑是“断头路”——几何上相交但属性上没有节点。做最短路径前，需要先构建拓扑。常见做法是用 networkx 把道路中心线转成图结构：

import networkx as nx import geopandas as gpd from shapely.geometry import LineString roads = gpd.read_file("./output/roads_study.shp") G = nx.Graph() for idx, row in roads.iterrows(): line = row.geometry if not isinstance(line, LineString): continue coords = list(line.coords) # 用首尾点作为节点，中间点忽略（简化处理） start, end = coords[0], coords[-1] length = line.length G.add_edge(start, end, weight=length, road_id=idx) # 检查连通分量 components = list(nx.connected_components(G)) print(f"连通分量数：{len(components)}")

这段代码把每条道路的首尾点作为图节点，道路长度作为权重。connected_components返回的连通分量数如果远大于 1，说明路网存在大量断头。实际项目中，我一般会先用ST_Node或shapely.ops.unary_union做节点打断，再构建图。参数上，weight用几何长度还是通行时间，取决于分析目标；如果做应急响应，应该用限速字段换算时间。

注意：SHP 字段名有 10 字符限制，导入 PostGIS 或转 GeoJSON 时可能被截断。如果原始字段名较长，建议在导入前用 ogr2ogr 重命名。

4. 避坑与排查：SHP 数据包最常见的五个翻车现场

4.1 中文乱码：属性表打开全是问号

现象：用 QGIS 打开正常，用 Python 读取后中文变成乱码或问号。原因：SHP 的 .dbf 文件默认编码是系统区域设置，中文 Windows 下常为 GBK，而 geopandas 默认按 UTF-8 读取。解决：读取时显式指定编码，gpd.read_file(path, encoding="gbk")；如果写入后乱码，在to_file时加encoding="utf-8"。更彻底的办法是用 ogr2ogr 转成 GeoPackage，GeoPackage 原生支持 UTF-8。

4.2 坐标系“看起来对，算起来错”

现象：地图上位置正确，但面积统计结果离谱。原因：数据是地理坐标（度），直接参与了几何计算。解决：所有量算前先投影到合适的投影坐标系。南京地区可用 EPSG:4547（CGCS2000 3 度带）或 EPSG:3857 做近似。判断方法：gdf.crs.is_geographic返回 True 就是地理坐标，需要to_crs后再算。

4.3 乡镇街道边界与区县边界不套合

现象：乡镇街道拼起来后，边缘和区县边界有缝隙或重叠。原因：两套数据来源不同，或制图综合程度不一致。解决：以区县边界为准，用gpd.overlay做相交修正；如果缝隙很小（小于容差），可以用buffer(0)或snap处理。但要注意，行政边界调整涉及口径问题，修正后要在报告里说明。

4.4 POI 重复与坐标偏移

现象：同一家店出现多次，或 POI 整体偏移几十米。原因：POI 数据多源采集，去重不彻底；坐标偏移常见于不同坐标系混用。解决：按名称+坐标距离去重，距离阈值设 50 米左右；偏移问题用已知控制点做配准，或直接与道路数据叠加目视检查。如果偏移是系统性的，考虑数据是否经过加密处理。

4.5 道路等级字段缺失或口径不一

现象：想做路网密度分析，发现道路等级字段为空或取值混乱。原因：数据包只保留了几何，属性简化。解决：用道路宽度或车道数反推等级，或与开源路网数据做空间连接补全。如果都没有，只能按几何长度做加权，在结论里注明局限性。

5. 进阶技巧：用空间连接把 POI 挂到街道上

前面把数据读进来、裁出来、跑通路径，最后一步往往是把 POI 和行政单元关联起来，做分街道的设施统计。这个操作叫空间连接，核心是判断点落在哪个面内。geopandas 的sjoin一行就能完成，但有几个参数决定结果对不对。

import geopandas as gpd townships = gpd.read_file("./nanjing_shp/townships.shp") poi = gpd.read_file("./nanjing_shp/poi.shp") # 统一坐标系 poi = poi.to_crs(townships.crs) # 空间连接：点落在面内 poi_with_town = gpd.sjoin(poi, townships, how="left", predicate="within") # 统计每个街道的 POI 数量 count_by_town = poi_with_town.groupby("name_right").size().reset_index(name="poi_count") print(count_by_town.sort_values("poi_count", ascending=False).head(10))

how="left"保证所有 POI 都保留，落在边界外的会标记为 NaN；predicate="within"表示点完全在面内。如果有些 POI 正好落在边界线上，within会漏掉，可以改用intersects，但那样边界点可能被两个街道同时统计。我一般先用within，再单独检查未匹配的点，看是数据问题还是边界问题。

另一个实用技巧是计算每个街道的 POI 密度。直接用数量除以面积，面积要从投影坐标算：

# 投影后计算面积，再算密度 townships_proj = townships.to_crs("EPSG:4547") townships_proj["area_km2"] = townships_proj.geometry.area / 1e6 # 合并统计结果 result = townships_proj.merge(count_by_town, left_on="name", right_on="name_right", how="left") result["poi_density"] = result["poi_count"] / result["area_km2"] print(result[["name", "poi_count", "area_km2", "poi_density"]].sort_values("poi_density", ascending=False).head())

这里EPSG:4547是 CGCS2000 的 3 度带投影，单位是米，面积除以 1e6 得到平方公里。密度排序能快速看出哪些街道商业设施集中，但要注意面积口径——如果街道边界包含大片水域或山地，密度会被低估，这时候需要先扣除不可达区域。

最后说一个我踩过的坑：空间连接前一定要检查两个图层的坐标系是否一致。有一次我直接 sjoin，结果所有 POI 都没匹配上，排查半天才发现一个图层是 4326，另一个是 3857。现在我的习惯是，任何空间操作前先打印gdf.crs，确认一致再往下走。希望帮到你。

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