简介这份南京市区县乡镇街道及交通设施矢量数据包面向城市规划、地理信息教学与空间分析从业者解决本地基础地理数据获取零散、图层不全的问题。包内共184个文件以SHP、SHX、DBF、PRJ、XML等标准GIS格式为主另有MXD工程文件与Python脚本压缩包约11.17MB适配ArcGIS、QGIS等主流软件。数据覆盖市、区、县、镇、村多级行政边界包含国道、省道、县道、高速、地铁轻轨、铁路及乡村道路等线状要素并整合水系、绿地、岛屿等面状自然要素以及公交站、地铁出入口、停车场、学校、医院、银行、超市、餐饮、药店、酒店、政府机构、公园等二十余类点状POI。坐标系为WGS84或CGCS2000投影信息完整属性字段命名规范可直接用于空间分析、专题制图与城市研究。目前已有84人学习下载适合需要快速搭建南京本地GIS实验环境、开展叠加分析与教学实训的读者。1. 南京市区县乡镇街道及交通设施矢量数据包从一份 SHP 到能跑的空间分析底图手上拿到一份标着“南京市区县乡镇街道及交通设施矢量数据包”的 SHP 集合第一反应往往不是兴奋而是犯嘀咕坐标系是哪套、字段编码对不对、道路和水系能不能叠上、POI 去重了没有。这份数据包的价值不在于“有数据”而在于它把行政边界、道路、水系、POI 这几层常用底图一次性凑齐省掉了从多个来源分别抓取再对齐的麻烦。它适合做城市尺度空间分析的人——比如做选址评估、做通勤可达性、做行政区划统计、做地图可视化底图。但 SHP 是上世纪留下来的格式字段名限 10 个字符、编码容易乱、多文件成套直接拖进 GIS 软件能打开不代表能用。这篇笔记按“先认清数据、再跑通加载、再处理坑、最后落到分析”的顺序讲重点放在坐标系对齐、字段清洗和拓扑检查这三件最容易翻车的事上。2. 先认清数据包里到底有什么图层、字段与坐标系2.1 行政边界三级的层级关系与字段设计区县、乡镇街道这两级行政边界是整份数据包的骨架。常见做法是区县一层、乡镇街道一层各自独立成 SHP靠行政代码字段关联。区县层的字段一般包含行政区名称、行政区代码、上级市代码乡镇街道层多一个上级区县代码。这里的关键是代码字段的位数和编码规则——行政区划代码是 6 位区县 3 位乡镇共 9 位或 12 位含村级不同来源截取长度不一样。拿到手先看字段别急着画图。import geopandas as gpd # 读取区县与乡镇街道两层先只看字段和坐标系 county gpd.read_file(nanjing_county.shp, encodingutf-8) town gpd.read_file(nanjing_town.shp, encodingutf-8) print(区县层字段:, list(county.columns)) print(乡镇层字段:, list(town.columns)) print(区县坐标系:, county.crs) print(乡镇坐标系:, town.crs) print(区县要素数:, len(county), 乡镇要素数:, len(town))这段代码的作用是先做“体检”而不是“动手术”。encodingutf-8是因为 SHP 的.dbf属性表中文编码经常是 GBK 或 UTF-8 混用读出来乱码就换gbk再试。crs打印出来如果是EPSG:4326说明是经纬度如果是EPSG:4526这类带带号的说明是投影坐标。两者混用会导致面积、距离全错。要素数用来判断数据完整度——南京下辖区县数量是固定的乡镇街道数量级在百位如果读出来只有个位数多半是文件没读全或路径指错了。2.2 道路、水系、POI 三类设施的几何类型差异交通设施层通常拆成道路线、水系线/面、POI 点三类。道路是 LineString水系可能是 LineString河流中心线也可能是 Polygon水面POI 是 Point。几何类型不统一意味着后续做空间连接时要用不同的谓词——点到线用最近距离面到面用相交。POI 层还要看有没有分类字段比如餐饮、教育、医疗没有分类字段的 POI 基本只能当密度点用。# 分别检查几何类型避免后续空间操作报错 for name, path in [(道路, nanjing_road.shp), (水系, nanjing_water.shp), (POI, nanjing_poi.shp)]: g gpd.read_file(path, encodingutf-8) geom_types g.geometry.geom_type.unique() print(f{name}: 几何类型{geom_types}, 要素数{len(g)}) # 检查是否有空几何空几何会让空间连接直接崩 print(f{name}: 空几何数量{g.geometry.is_empty.sum()})geom_type.unique()能一眼看出这层是不是纯线或纯点。如果道路层里混进了 Polygon说明数据源里把立交桥面也塞进来了做网络分析前得先过滤。is_empty统计空几何空几何在sjoin或buffer时会抛异常属于典型的“打开没事、一算就崩”的坑。这一步花两分钟能省后面半小时排错。2.3 坐标系判断为什么不能直接叠加坐标系是这份数据包最容易被忽略、后果最严重的一环。经纬度坐标EPSG:4326下算距离一度约 111 公里但经度方向随纬度收缩南京在北纬 32 度附近经度一度约 94 公里。如果直接拿 4326 算缓冲区圆形会变成椭圆面积统计偏差能到百分之十几。正确做法是统一投影到适合南京的平面坐标系常见的是 CGCS2000 的 3 度带投影南京大致落在带号 39 或 40 的范围。# 统一投影先确认各层原始 CRS再统一转到投影坐标 target_crs EPSG:4526 # CGCS2000 / 3度带中央经线117度附近 layers {} for name, path in [(county, nanjing_county.shp), (town, nanjing_town.shp), (road, nanjing_road.shp), (water, nanjing_water.shp), (poi, nanjing_poi.shp)]: g gpd.read_file(path, encodingutf-8) if g.crs is None: print(f{name} 缺少 CRS 定义需手动指定) else: layers[name] g.to_crs(target_crs) print(f{name} 已转到 {target_crs})to_crs是重投影不是简单改标签。如果某层crs是None说明文件里没写坐标系信息这时候不能瞎转得先根据数据范围判断原始坐标系再手动set_crs。判断方法看经纬度范围如果 x 在 118 到 119 之间、y 在 31 到 32 之间那就是 4326如果 x 是六位数、y 是七位数那是投影坐标。这一步判断错了后面所有距离和面积都是错的而且错得很隐蔽。3. 把 SHP 跑通加载、裁剪与空间连接的完整流程3.1 用 GeoPandas 批量读取并统一投影实际项目里不会一层层手写读取通常写个批量函数把路径、编码、目标坐标系作为参数传进去。这样换一批数据只改配置不改逻辑。批量读取时要注意 SHP 的伴生文件——.shp、.shx、.dbf、.prj必须同名同目录缺一个都读不进来或读进来没坐标系。import geopandas as gpd from pathlib import Path def load_layers(base_dir, layer_map, target_crsEPSG:4526): 批量读取 SHP 并统一投影 base_dir: 数据目录 layer_map: {逻辑名: 文件名} 映射 result {} for name, fname in layer_map.items(): path Path(base_dir) / fname # 先试 utf-8失败回退 gbk中文属性表常见 try: g gpd.read_file(path, encodingutf-8) except UnicodeDecodeError: g gpd.read_file(path, encodinggbk) if g.crs is not None: g g.to_crs(target_crs) result[name] g print(f{name}: {len(g)} 要素, CRS{g.crs}) return result layer_map { county: nanjing_county.shp, town: nanjing_town.shp, road: nanjing_road.shp, water: nanjing_water.shp, poi: nanjing_poi.shp, } data load_layers(./nanjing_shp, layer_map)layer_map把逻辑名和文件名解耦后面代码里用data[road]就行不用记文件名。编码回退逻辑是血泪经验——很多从旧系统导出的 SHP 属性表是 GBK直接 utf-8 读会抛UnicodeDecodeError捕获后换 gbk 基本能解决。to_crs只在crs不为空时执行避免对无坐标系数据做错误转换。3.2 按行政区裁剪出研究区做分析很少用全市数据通常先框定一个区或几个街道。裁剪用clip前提是裁剪边界和被裁图层坐标系一致。这里用区县层里某个区去裁道路和水系得到研究区内的设施。# 假设区县层有名称字段取某个区作为研究区 study_area data[county][data[county][名称] 某区] # 用研究区边界裁剪道路和水系 road_clip gpd.clip(data[road], study_area) water_clip gpd.clip(data[water], study_area) print(裁剪后道路数:, len(road_clip)) print(裁剪后水系数:, len(water_clip)) # 裁剪后重新计算长度注意单位是投影坐标的米 road_clip[length_m] road_clip.geometry.length print(道路总长(km):, road_clip[length_m].sum() / 1000)clip保留的是几何上落在边界内的部分跨边界的线会被切断这是符合预期的。geometry.length在投影坐标系下单位是米直接除以 1000 得公里。如果忘了投影还在 4326这个长度就是“度”数值小得离谱且没有物理意义。裁剪后要素数应该比裁剪前少如果反而多了多半是边界本身有问题或坐标系没对齐。3.3 POI 与行政区的空间连接统计POI 最常见的用法是统计每个区或每个街道里有多少个某类设施。用sjoin做空间连接把 POI 落到行政区上再分组计数。这里要注意sjoin的predicate参数点在面内用within面包含点用contains两者结果一样但方向相反。# 把 POI 落到区县上统计每个区的 POI 数量 poi_join gpd.sjoin(data[poi], data[county], howleft, predicatewithin) # 按区县名称分组计数 poi_count poi_join.groupby(名称).size().reset_index(namepoi_count) print(poi_count.sort_values(poi_count, ascendingFalse).head(10)) # 如果 POI 有分类字段可以按类别透视 if 类别 in poi_join.columns: pivot poi_join.pivot_table(index名称, columns类别, valuesgeometry, aggfunccount) print(pivot.fillna(0).astype(int))sjoin的howleft保证落在边界外的 POI 也保留但行政区字段为空方便排查边界问题。groupby后计数是最基础的密度指标。如果 POI 有类别字段pivot_table能直接出交叉表aggfunccount对几何列计数等价于数点。注意sjoin在大数据量下会慢POI 上万条时可以先做边界框粗筛再精确连接。4. 避坑与排查SHP 数据包最容易翻车的五个地方4.1 中文属性表乱码现象、原因与解决现象是打开属性表看到“å京”这类字符或者直接报编码错误。原因是 SHP 的.dbf文件头里有个语言驱动标识但很多导出工具不写或写错读取时默认按系统编码解析。解决办法是读取时显式指定编码utf-8 不行就 gbk再不行用chardet探测。如果已经读进来是乱码可以尝试用g.to_file重新导出时指定encodingutf-8覆盖原文件。更稳妥的做法是读进来后立刻把字段值转成 Python 字符串再处理避免后续操作反复触发编码问题。4.2 坐标系缺失或标错最隐蔽的错误现象是数据能打开、能画图但一算距离面积就不对或者和其他图层叠不上。原因是.prj文件缺失或者内容写的是 WGS84 但实际是别的坐标系。排查方法是看坐标数值范围经纬度量级在几十到一百多投影坐标在几十万到几百万。如果范围对不上标注的坐标系就以数值范围为准手动set_crs。另一个办法是拿一个已知坐标的点比如某个明显地标去比对偏差在几十米内说明坐标系基本对偏差几公里说明错了。4.3 几何无效自相交与重复点现象是buffer、intersection等操作报拓扑异常或者面积算出来是负的。原因是多边形自相交、线有重复点、环没闭合。解决办法是用make_valid修复GeoPandas 新版本有g.geometry.make_valid()老版本用buffer(0)这个经典偏方。修复后要重新检查面积自相交的多边形修复后面积会变。对于线数据重复点用simplify或remove_repeated_points处理。# 修复无效几何 invalid ~data[county].geometry.is_valid print(无效几何数:, invalid.sum()) if invalid.sum() 0: data[county][geometry] data[county].geometry.make_valid() # 修复后重新验证 print(修复后无效数:, (~data[county].geometry.is_valid).sum())is_valid返回布尔序列取反统计无效数。make_valid是 Shapely 2.0 之后的标准做法比buffer(0)更可控不会把多边形炸成多块。修复后一定要重新验证因为有些几何修复后仍无效需要进一步处理。4.4 字段名截断SHP 的 10 字符限制现象是字段名变成行政区名_1、POP_2020_这种或者两个字段名撞车。原因是 SHP 格式规定字段名最多 10 个字符超长会被截断中文按字节算更容易超。解决办法是读进来后立刻重命名字段用简短英文或拼音在代码里维护一份映射表。如果必须保留长字段名考虑转成 GeoPackage 或 GeoJSON 格式这两种没有字段名长度限制。4.5 道路拓扑断裂网络分析前的必修课现象是做路径分析时找不到连通路径或者最短路径绕远。原因是道路线在交叉口没有打断两条路交叉但几何上不相交或者相交但没共享节点。解决办法是用unary_union合并所有线再polygonize或line_merge然后按节点打断。这一步比较重通常用专门的网络分析库处理。如果只是做可视化或长度统计可以跳过一旦涉及可达性、服务区分析拓扑必须处理干净。5. 从底图到分析三个能直接复用的进阶技巧5.1 用缓冲区算设施覆盖率覆盖率是选址和评估里最常用的指标。做法是给设施点做缓冲区再和行政区相交算覆盖面积占比。缓冲区半径按设施类型定比如学校按 1 公里医院按 3 公里。注意缓冲区要在投影坐标系下做否则半径单位是度完全不对。# 给 POI 做 1000 米缓冲区算每个区的覆盖率 poi_buf data[poi].copy() poi_buf[geometry] poi_buf.geometry.buffer(1000) # 缓冲区与区县相交算覆盖面积 covered gpd.overlay(poi_buf, data[county], howintersection) covered_area covered.dissolve(by名称, aggfuncsum).geometry.area # 各区总面积 total_area data[county].set_index(名称).geometry.area coverage (covered_area / total_area).dropna() print(覆盖率:) print(coverage.sort_values(ascendingFalse))buffer(1000)在投影坐标下就是 1000 米。overlay做相交得到覆盖部分dissolve按区合并避免同一区多个缓冲区重叠重复计算。covered_area / total_area得到覆盖率。这里有个细节如果缓冲区超出区边界overlay会自动裁掉所以分子不会超过分母。覆盖率超过 1 说明坐标系或面积计算有问题。5.2 道路密度分级与网格统计道路密度是衡量区域开发强度的常用指标。做法是把研究区切成规则网格算每个网格内的道路长度再除以网格面积。网格大小按分析尺度定城市尺度用 500 米或 1 公里。用geopandas生成网格再和道路相交。import numpy as np # 生成 1 公里网格 bounds data[county].total_bounds cell 1000 # 米 x np.arange(bounds[0], bounds[2], cell) y np.arange(bounds[1], bounds[3], cell) from shapely.geometry import box grid gpd.GeoDataFrame( geometry[box(xi, yi, xicell, yicell) for xi in x for yi in y], crsdata[county].crs ) # 网格与道路相交算每个网格的道路长度 road_grid gpd.overlay(data[road], grid, howintersection) road_grid[len] road_grid.geometry.length density road_grid.groupby(road_grid.index_right)[len].sum() / (cell * cell) print(道路密度分布:) print(density.describe())total_bounds拿到研究区外接矩形np.arange按 1 公里步长生成网格左下角坐标box生成方格。overlay相交后每条道路被网格切断groupby按网格索引汇总长度。除以cell*cell得密度单位是米每平方米乘 10000 换成米每公顷更直观。describe看分布密度特别高的网格通常是主干道交汇处。5.3 数据导出与格式选择分析完要导出SHP 不是唯一选择。如果字段名长、要存中文、要单文件用 GeoPackage.gpkg更合适如果要给前端用GeoJSON 更通用但体积大如果只是内部交换Parquet 加几何列性能最好。导出时注意坐标系——给别人用的数据最好同时保留投影和经纬度两个版本投影用于计算经纬度用于展示。# 导出为 GeoPackage保留长字段名和中文 data[county].to_file(output_county.gpkg, layercounty, driverGPKG) data[road].to_file(output_road.gpkg, layerroad, driverGPKG) # 导出经纬度版本用于展示 data[county].to_crs(EPSG:4326).to_file( output_county_wgs84.geojson, driverGeoJSON)to_file的driver参数指定格式GPKG 支持多图层共存一个文件放多层。to_crs(EPSG:4326)转回经纬度再导出 GeoJSON给前端地图库用。导出前建议先drop掉中间计算字段只留原始字段和关键结果避免文件臃肿。我自己的习惯是拿到任何一份 SHP 数据包先花十分钟做三件事打印每层的 CRS 和几何类型、统计空几何和无效几何、看字段名有没有被截断。这三步做完后面百分之八十的报错都能提前避开。数据清洗没有捷径但把检查做成固定流程就能少熬几个夜。希望帮到你。本文还有配套的精品资源点击获取
