简介这份国家基础地理信息系统数据集面向GIS学习者、城市规划与交通研究人员以及需要中国基础地理要素做空间分析的用户提供一套可直接加载的SHP矢量数据解决行政边界、路网、水系等底图数据获取零散的问题。压缩包共93个文件约10.63MB以shp、shx、dbf、prj、sbn、sbx等Shapefile配套格式为主另有少量xml元数据其中shp存储几何要素dbf保存属性字段prj定义坐标投影shx与sbn/sbx支撑索引与空间查询便于在ArcGIS、QGIS等平台直接打开。内容覆盖主要公路、铁路、河流、湖泊、国界与省界、县界、地州界、省会及县城驻地、经纬网和县级统计数据等要素可支撑交通规划、物流分析、行政区划统计与生态研究。已有1943人学习下载适合作为课程设计、论文制图与空间统计的底图素材。1. 拿到一份 foreste75 的 shp 主要公路数据先搞清楚它到底能干什么你从某个渠道拿到一个压缩包名字叫「国家基础地理信息系统数据.zip」解压后里面有个 foreste75 目录再往下翻是一堆 .shp、.shx、.dbf、.prj 文件其中一组文件名里带着「主要公路」。这不是什么新鲜事但很多人第一次面对它时会卡在同一个地方这些文件怎么打开、坐标系是什么、主要公路图层跟普通道路图层有什么区别、能不能直接拿来做路网分析或者出图。这份数据的核心价值在于「主要公路」这个图层本身。它通常记录的是国道、省道及以上等级的干线公路几何类型以 LineString 为主属性表里一般带公路编号、名称、等级等字段。foreste75 这个前缀暗示它可能来自某个分幅编号体系但不管来源如何你拿到手之后要做的第一件事不是急着加载而是先确认三件事坐标系、属性字段、几何完整性。这三件事决定了你后面能不能把它跟其他数据叠在一起、能不能做空间查询、能不能导出成其他格式。适合读这篇的人做 GIS 实习的学生、需要快速出路网底图的规划从业者、想把 shp 转成其他格式做可视化或分析的人。如果你只是想把主要公路画到地图上看看那很简单但如果你要拿它做缓冲区分析、路径规划、或者跟行政区划做叠加统计那下面这些步骤和坑你就绕不开。2. 主要公路 shp 的坐标系、属性表和几何检查动手前必须过的三关2.1 先看 .prj 文件别等叠加时才发现坐标对不上.shp 本身不存坐标系信息坐标系写在同名的 .prj 文件里。很多人拿到数据后直接拖进 ArcGIS 或 QGIS软件不报错就以为没问题结果跟另一份数据叠在一起时发现一个在赤道附近、一个在北极圈。常见做法是先用文本编辑器打开 .prj看里面的 PROJCS 或 GEOGCS 字段。如果 .prj 缺失或者内容明显不对你就需要根据数据范围来判断。foreste75 这类分幅数据国内常见的是 CGCS2000 或西安80 坐标系投影方式多为高斯-克吕格 3 度带或 6 度带。判断方法看 shp 的坐标值。如果 X 坐标在 6 到 7 位数、Y 坐标在 7 到 8 位数基本是投影坐标如果 X 在 70 到 140、Y 在 10 到 60那是经纬度。# 用 ogrinfo 查看 shp 的坐标系和基本信息 ogrinfo -so -al foreste75_highway.shp # 输出里重点关注这几行 # Geometry: Line String # Feature Count: 1247 # Extent: (xxx, xxx) - (xxx, xxx) # Layer SRS WKT: PROJCS[...-so表示只输出摘要-al表示所有图层。输出里的 Extent 能帮你快速判断坐标范围是否合理Layer SRS WKT 就是坐标系定义。如果 SRS 显示为 unknown那就得手动指定。提示不要依赖软件自动识别坐标系。QGIS 和 ArcGIS 的自动识别有时会把 CGCS2000 认成 WGS84虽然差异不大但做精确叠加时会有米级偏移。2.2 属性表里哪些字段真正有用用 ogrinfo 看完坐标系后接着看属性字段。主要公路图层的属性表通常不会太复杂但字段名可能是英文缩写或拼音需要你对照着判断。# 查看属性表的字段定义和前几条记录 ogrinfo -al -geomNO foreste75_highway.shp | head -60常见字段及含义字段名示例含义是否关键NAME / ROADNAME公路名称出图标注用CODE / ROADCODE公路编号如 G318分类和查询用CLASS / GRADE公路等级符号化用LENGTH路段长度统计用但可能是错的TYPE道路类型辅助判断重点看 CLASS 或 GRADE 字段的取值分布。如果取值是「高速」「国道」「省道」这种中文说明数据整理过如果是数字代码你需要找对应的编码表。没有编码表怎么办按公路编号前缀判断G 开头是国道S 开头是省道X 开头是县道。如果连编号都没有那就只能按几何长度和连通性来推断重要性。LENGTH 字段要特别小心。很多 shp 的 LENGTH 字段是数字化时自动计算的单位可能是度而不是米直接拿来统计会闹笑话。正确做法是用投影后的几何重新计算长度。2.3 几何检查断线、重复线、自相交怎么发现主要公路数据最常见的几何问题是断线。一条国道在分幅边界处被切断或者数字化时漏了一段导致路网不连通。如果你要做路径分析断线是致命的。import geopandas as gpd from shapely.validation import explain_validity gdf gpd.read_file(foreste75_highway.shp) # 检查几何有效性 invalid gdf[~gdf.is_valid] print(f无效几何数量: {len(invalid)}) for idx, row in invalid.iterrows(): print(f FID {idx}: {explain_validity(row.geometry)}) # 检查重复几何 dup gdf[gdf.duplicated(subsetgeometry, keepFalse)] print(f重复几何数量: {len(dup)}) # 检查零长度线 zero_len gdf[gdf.geometry.length 0] print(f零长度线数量: {len(zero_len)})这段代码做了三件事is_valid检查几何是否自相交或环方向错误duplicated找出完全重复的线length 0找出退化的零长度线。如果无效几何数量不为零用gdf.geometry gdf.geometry.buffer(0)可以修复大部分自相交问题但 buffer(0) 会改变几何形状修复后要重新检查。断线检测没有一行代码的解决方案。实用做法是把线图层转成端点然后看有多少端点是悬空的只连接一条线。悬空端点就是断线处。# 提取所有线的端点统计每个端点的连接数 from shapely.geometry import Point import pandas as pd endpoints [] for geom in gdf.geometry: coords list(geom.coords) endpoints.append(Point(coords[0])) endpoints.append(Point(coords[-1])) ep_gdf gpd.GeoDataFrame(geometryendpoints) ep_gdf[x] ep_gdf.geometry.x.round(6) ep_gdf[y] ep_gdf.geometry.y.round(6) counts ep_gdf.groupby([x, y]).size() dangling counts[counts 1] print(f悬空端点数量: {len(dangling)})悬空端点数量除以 2 大致就是断线处数量。如果这个数字很大说明数据不适合直接做路网分析需要先做拓扑修复。修复工具可以用 ArcGIS 的「拓扑检查」或 QGIS 的「v.clean」但更快的办法是找一份更完整的路网数据做参考手动补上关键断点。3. 把主要公路 shp 用起来从格式转换到空间分析的完整链路3.1 shp 转 GeoJSON、KML、PostGIS 的命令与参数拿到 shp 之后不同场景需要不同格式。做网页可视化用 GeoJSON做 Google Earth 展示用 KML做后端空间查询用 PostGIS。转换工具首选 GDAL 的 ogr2ogr一条命令搞定。# shp 转 GeoJSON指定坐标系为 WGS84 ogr2ogr -f GeoJSON -t_srs EPSG:4326 highway.geojson foreste75_highway.shp # shp 转 KML注意 KML 只支持 WGS84 ogr2ogr -f KML -t_srs EPSG:4326 highway.kml foreste75_highway.shp # shp 导入 PostGIS ogr2ogr -f PostgreSQL PG:hostlocalhost dbnamegis userpostgres passwordxxx \ foreste75_highway.shp -nln highway -lco GEOMETRY_NAMEgeom -lco FIDgid-t_srs EPSG:4326是重投影参数把源数据的坐标系转到 WGS84。如果源数据已经是 WGS84这个参数可以省略但加上更保险。-nln指定导入后的表名-lco是图层创建选项GEOMETRY_NAME 指定几何字段名FID 指定主键字段名。转 KML 时有个坑如果属性表里有中文KML 的编码可能出问题。加-lco ENCODINGUTF-8可以解决大部分情况。另外 KML 对线要素的节点数有限制节点太多的线会被截断必要时先用-simplify抽稀。# 抽稀后再转 KML容差 0.0001 度 ogr2ogr -f KML -t_srs EPSG:4326 -simplify 0.0001 highway_simple.kml foreste75_highway.shp-simplify的容差单位跟源坐标系一致。如果是经纬度数据0.0001 度大约是 10 米这个精度对公路展示足够了。3.2 用 GeoPandas 做缓冲区分析和叠加统计GeoPandas 是 Python 里做矢量分析最顺手的库。假设你要统计每个行政区划内主要公路的里程或者做公路两侧 500 米的缓冲区下面这套流程可以直接抄。import geopandas as gpd # 读取主要公路和行政区划 highway gpd.read_file(foreste75_highway.shp) districts gpd.read_file(districts.shp) # 统一坐标系到投影坐标系方便算长度和面积 # 根据数据所在位置选择合适的投影这里以 CGCS2000 3度带为例 highway_proj highway.to_crs(epsg4547) districts_proj districts.to_crs(epsg4547) # 计算每条公路的长度米 highway_proj[length_m] highway_proj.geometry.length # 空间叠加把公路按行政区划分组 joined gpd.sjoin(highway_proj, districts_proj, howleft, predicateintersects) # 统计每个区划内的公路总里程 stats joined.groupby(district_name)[length_m].sum().reset_index() stats.columns [区划名称, 公路总里程(米)] stats[公路总里程(公里)] stats[公路总里程(米)] / 1000 print(stats.sort_values(公路总里程(公里), ascendingFalse))关键参数说明to_crs(epsg4547)是 CGCS2000 3 度带第 39 带适用于东经 114 到 120 度区域。如果你不知道用哪个 EPSG可以用highway.estimate_utm_crs()让 GeoPandas 自动推荐。sjoin的predicateintersects表示只要公路与区划有交集就算如果只想统计完全在区划内的公路改成predicatewithin。缓冲区分析同理# 公路两侧 500 米缓冲区 buffer highway_proj.copy() buffer[geometry] highway_proj.geometry.buffer(500) # 合并所有缓冲区避免重叠 buffer_union buffer.unary_union # 统计缓冲区覆盖了哪些区划 covered districts_proj[districts_proj.intersects(buffer_union)] print(f500米缓冲区影响到的区划数量: {len(covered)})buffer(500)的单位是米因为已经投影了。unary_union把所有缓冲区合并成一个几何体避免后续叠加时重复计算。3.3 路网连通性修复与最短路径的实操步骤如果你拿主要公路做路径分析断线问题必须先解决。下面是一个基于 NetworkX 的简化流程适合中小规模路网。import networkx as nx from shapely.geometry import LineString import geopandas as gpd highway gpd.read_file(foreste75_highway.shp).to_crs(epsg4547) # 构建图每条线的端点是节点线是边 G nx.Graph() for idx, row in highway.iterrows(): coords list(row.geometry.coords) start (round(coords[0][0], 3), round(coords[0][1], 3)) end (round(coords[-1][0], 3), round(coords[-1][1], 3)) length row.geometry.length G.add_edge(start, end, weightlength, fididx) # 找到最大的连通分量 components list(nx.connected_components(G)) largest max(components, keylen) print(f最大连通分量包含 {len(largest)} 个节点总节点 {G.number_of_nodes()} 个) # 在最大连通分量内做最短路径 nodes list(largest) source, target nodes[0], nodes[-1] try: path nx.shortest_path(G, source, target, weightweight) print(f最短路径经过 {len(path)} 个节点) except nx.NetworkXNoPath: print(两点之间不连通需要修复断线)节点坐标用round(..., 3)是为了让几乎重合的端点合并成同一个节点。3 位小数在投影坐标系下大约是 1 毫米精度足够合并数字化误差。如果连通分量很多说明断线严重需要回到 2.3 节的悬空端点检测找到断点位置手动补线。注意NetworkX 的图是内存中的节点数超过百万时会很慢。主要公路数据通常几千到几万条线用 NetworkX 没问题。如果数据量更大换 igraph 或 PostGIS 的 pgRouting。4. 主要公路 shp 处理中的避坑与排查5 个血泪教训4.1 坑一坐标系缺失导致面积和长度全错现象用 GeoPandas 算公路长度结果每条路只有零点几明显不对。原因数据是经纬度坐标EPSG:4326直接算geometry.length得到的是度数不是米。1 度大约 111 公里但随纬度变化。解决先to_crs()转到投影坐标系再算长度。不知道用哪个投影就用estimate_utm_crs()。如果数据本身没有 .prj先用set_crs(epsg4326)指定再to_crs()转换。4.2 坑二中文属性字段乱码现象用 ogrinfo 或 GeoPandas 读 shp属性表里的中文显示成乱码。原因shp 的 .dbf 文件默认编码是 GBK 或 Latin-1而 GDAL 默认按 UTF-8 读。解决读的时候指定编码。GeoPandas 用gpd.read_file(file.shp, encodinggbk)。ogr2ogr 转换时加-lco ENCODINGUTF-8。如果还是乱码用encodinglatin1试试或者用 QGIS 打开后重新导出。4.3 坑三shp 转 KML 后线要素消失现象转出来的 KML 在 Google Earth 里打开部分公路线看不到。原因KML 对单个要素的节点数有限制超过约 1000 个节点的线会被截断或丢弃。另外 KML 只支持 WGS84如果没加-t_srs EPSG:4326坐标会跑到奇怪的地方。解决转换前先抽稀-simplify 0.0001。如果抽稀后还是丢把长线拆成多段。用-t_srs EPSG:4326确保坐标系正确。4.4 坑四sjoin 叠加后统计结果翻倍现象用gpd.sjoin把公路和行政区划叠加后统计每个区划的公路里程发现总数比实际多。原因一条公路跨越两个区划时sjoin 会生成两条记录每条都带完整的公路长度。直接 groupby 求和会把跨界的公路算两次。解决先做几何切割把公路按区划边界打断再统计。或者用gpd.overlay(highway, districts, howintersection)得到切割后的路段再算长度。overlay 会自动处理边界切割。4.5 坑五PostGIS 导入后查询慢现象shp 导入 PostGIS 后用 ST_Intersects 查询要等好几秒。原因没有建空间索引。PostGIS 默认不会自动为导入的几何字段建索引。解决导入后手动建索引。CREATE INDEX idx_highway_geom ON highway USING GIST (geom);建完索引后空间查询速度通常能提升几十倍。另外确认 SRID 设置正确SELECT UpdateGeometrySRID(highway, geom, 4326);如果导入时没指定 SRIDPostGIS 会默认设为 0导致空间函数无法正确使用索引。5. 从主要公路 shp 到可复用的路网底图一个我常用的验证习惯处理完坐标系、属性、几何这三关之后我一般会做一件事把主要公路叠加到一份已知正确的底图上目视检查。底图可以用在线瓦片也可以用另一份可信的行政区划数据。这一步花不了几分钟但能发现很多数据层面的问题——比如公路整体偏移、部分路段缺失、属性字段错位。具体做法在 QGIS 里加载主要公路 shp再加一个 XYZ Tiles 底图比如 OpenStreetMap把公路图层设为红色半透明线。如果公路跟底图上的道路基本重合说明坐标系没问题如果整体偏移几百米说明坐标系选错了如果部分路段对不上可能是数据本身的问题。# 用 contextily 快速给 GeoPandas 加底图做目视检查 import geopandas as gpd import contextily as ctx import matplotlib.pyplot as plt highway gpd.read_file(foreste75_highway.shp).to_crs(epsg3857) fig, ax plt.subplots(figsize(12, 12)) highway.plot(axax, colorred, linewidth1.5, alpha0.7) ctx.add_basemap(ax, sourcectx.providers.OpenStreetMap.Mapnik) ax.set_axis_off() plt.savefig(highway_check.png, dpi150, bbox_inchestight)to_crs(epsg3857)是 Web Mercator 投影contextily 的底图默认用这个坐标系。alpha0.7让公路半透明方便看底图上的道路是否对齐。保存成图片后放大看细节比在 GIS 软件里来回缩放快得多。另一个习惯是每次处理完 shp都把关键参数记在一个文本文件里——源坐标系、目标坐标系、使用的 EPSG 代码、处理日期、数据行数。下次再用这份数据时不用重新猜。这个习惯帮我省了很多重复排查的时间。如果你拿到的 foreste75 主要公路数据只是整个国家基础地理信息系统数据的一部分那它可能还有配套的铁路、水系、居民地图层。处理思路是一样的先查坐标系再看属性最后做几何检查。把这一套流程跑顺了以后拿到任何 shp 数据都能快速上手。希望帮到你。本文还有配套的精品资源点击获取
