简介本资源为2022年7月整理的四川省地理空间数据集面向GIS专业人员、城市规划与交通研究人员及高校师生用于行政区域分析、交通网络规划与空间可视化等场景。压缩包共33个文件约66.64MB以shp、dbf、prj、shx等Shapefile核心格式为主辅以sbn、sbx空间索引及xml元数据文件覆盖四川省、市、县三级行政边界以及道路网与铁路网数据属性表与投影信息完整可在ArcGIS、QGIS中直接加载。借助这些数据读者可开展人口分布与经济发展差异比较、道路铁路布局评估、交通瓶颈识别、选址与服务覆盖分析并结合其他数据构建预测模型。目前已有2630人学习下载适合需要四川省基础地理底图与交通网络数据的中高级GIS用户参考使用。1. 从一份 2022.07 的四川省 SHP 数据说起省、市、县、路网到底怎么用如果你手头拿到一份标注为「2022.07」的四川省行政区划 SHP 文件里面同时包含省、市、县三级边界外加道路网和公路网两层线数据第一反应大概率不是「太好了」而是「这几层怎么叠、坐标系对不对、字段能不能直接拿来分组」。我见过太多人把这类数据直接拖进 GIS 软件里出图结果县级边界和路网错位几百米或者属性表里全是拼音缩写根本没法做统计。这份数据的价值不在于「有」而在于你能不能把它拆成可复用的分析底座行政区划用来做空间聚合和裁剪道路网用来做可达性、缓冲区、连通性分析公路网用来区分等级和通行能力。适合谁做国土空间、物流选址、应急调度、区域经济分析的人以及需要一份干净底图做可视化的人。下面我按「先验数据、再跑流程、最后避坑」的顺序把这份 SHP 从打开到出结果的全过程拆开讲。2. 拿到 SHP 先别急着画图三层数据的字段与坐标系核对2.1 省、市、县三级边界的属性表该看哪几列行政区划 SHP 通常是一个面图层但省、市、县三级可能分三个文件也可能合在一个图层里用层级字段区分。我一般先看字段名里有没有NAME、ADCODE、LEVEL、PARENT这类列。如果没有就得靠面积或包含关系反推层级。常见做法是省级面面积最大市级次之县级最小但四川有甘孜、阿坝、凉山这些面积巨大的县单靠面积会翻车。更稳的办法是用ADCODE的位数判断省级通常是 6 位前两位非零后四位零市级是前四位非零后两位零县级是 6 位全非零。如果数据里没有ADCODE那就用空间包含关系一个面被另一个面完全包含且面积更小大概率是下级。import geopandas as gpd # 读取行政区划面数据注意编码中文属性常出现乱码 gdf gpd.read_file(sc_admin_202207.shp, encodingutf-8) # 查看字段名和前几行确认层级字段 print(gdf.columns.tolist()) print(gdf[[NAME, ADCODE, LEVEL]].head(10) if LEVEL in gdf.columns else gdf.head(10)) # 如果没有 LEVEL 字段用 ADCODE 位数生成层级 if LEVEL not in gdf.columns and ADCODE in gdf.columns: gdf[LEVEL] gdf[ADCODE].astype(str).apply( lambda x: province if x[2:] 0000 else (city if x[4:] 00 else county) )这段代码先解决「有没有层级字段」的问题。encodingutf-8不是万能的如果报UnicodeDecodeError换成gbk或gb18030再试。ADCODE转字符串后判断后四位是否为零是区分省市县的硬规则比面积靠谱。参数上read_file默认读第一个图层如果 SHP 是多图层数据集得用layer参数指定。2.2 道路网和公路网的线数据怎么区分等级道路网和公路网经常被混在一起但它们的用途不同。道路网通常包含城市道路、乡村道路、匝道等公路网则偏向国道、省道、县道带有ROAD_TYPE、GRADE、WIDTH这类字段。我拿到线数据后先看GRADE或TYPE字段的取值分布再决定怎么分类。如果字段里是「高速」「一级」「二级」这种中文直接分组如果是数字编码得找对应的编码表。没有等级字段的话用WIDTH或LANES反推但这条路容易踩坑后面避坑章节细说。# 读取公路网线数据 roads gpd.read_file(sc_road_202207.shp, encodingutf-8) # 查看等级字段的取值分布 if GRADE in roads.columns: print(roads[GRADE].value_counts()) elif TYPE in roads.columns: print(roads[TYPE].value_counts()) else: print(无等级字段检查 WIDTH/LANES 列) print(roads.columns.tolist()) # 按等级筛选高速和国道用于后续缓冲区分析 highway roads[roads[GRADE].isin([高速, 国道])] if GRADE in roads.columns else roadsvalue_counts()能快速暴露字段里的脏数据比如「高速」和「高速公路」混用。isin筛选时要把所有可能的写法列全否则会漏。如果数据里等级字段是数字比如 1 代表高速、2 代表国道那就得先建映射字典再筛选。2.3 坐标系不统一是错位的头号原因SHP 文件自带.prj文件里面写着坐标系。行政区划常用CGCS2000或WGS84道路网可能用GCJ02或者投影坐标系。如果两个图层的.prj不一致直接叠图就会错位。我一般先读.crs属性看是不是地理坐标系单位度还是投影坐标系单位米。做面积统计和缓冲区必须用投影坐标系否则算出来的面积是平方度没有意义。# 检查各图层坐标系 print(行政区划 CRS:, gdf.crs) print(公路网 CRS:, roads.crs) # 如果行政区划是地理坐标系转为投影坐标系再做面积统计 if gdf.crs and gdf.crs.is_geographic: # 四川常用 CGCS2000 3度带中央经线根据范围选这里用 EPSG:4544 示例 gdf_proj gdf.to_crs(epsg4544) gdf_proj[area_km2] gdf_proj.geometry.area / 1e6 print(gdf_proj[[NAME, area_km2]].head())to_crs是转换坐标系的核心方法epsg4544是 CGCS2000 3 度带的一个示例实际用哪个带号要看数据覆盖的经度范围。四川跨多个 3 度带如果做全省分析建议统一用EPSG:4490地理坐标系做展示用EPSG:4544或EPSG:4545做面积量算。area / 1e6是把平方米转成平方公里这个除数别写错。3. 用 GeoPandas 跑通裁剪、叠加与路网缓冲区的最小流程3.1 用县级边界裁剪路网得到每个县的公路密度公路密度是区域分析里最常用的指标之一算法很简单县内公路总长度除以县面积。但前提是先把路网裁剪到县级边界内。裁剪用clip或overlayclip更快适合线被面裁。裁剪后要重新计算长度因为一条路可能跨多个县裁剪后每个县只保留自己那一段。# 确保两个图层坐标系一致 roads_proj roads.to_crs(gdf_proj.crs) # 用县级边界裁剪路网 county gdf_proj[gdf_proj[LEVEL] county] roads_clip gpd.clip(roads_proj, county) # 计算每个县内的公路总长度 roads_clip[length_km] roads_clip.geometry.length / 1000 density roads_clip.groupby(ADCODE)[length_km].sum().reset_index() density density.merge(county[[ADCODE, NAME, area_km2]], onADCODE) density[road_density] density[length_km] / density[area_km2] print(density[[NAME, length_km, area_km2, road_density]].sort_values(road_density, ascendingFalse).head())gpd.clip的第二个参数是裁剪面返回的是线被面切割后的结果。groupby(ADCODE)按行政区代码汇总长度再和县级面积表合并。road_density单位是公里每平方公里数值越大路网越密。注意clip会保留线的属性如果路网数据里有重复线段长度会偏大后面避坑章节会讲怎么去重。3.2 道路网缓冲区分析以高速出入口为例缓冲区是路网分析里最直观的操作。比如要分析高速出入口 5 公里范围内覆盖了多少人口或多少县就得先对高速点做缓冲区再和行政区划做叠加。如果只有线数据没有出入口点可以用线数据的端点或节点近似但精度会差一些。我一般用buffer生成面再用sjoin做空间连接。# 筛选高速线转成投影坐标系后做 5 公里缓冲区 highway_proj highway.to_crs(gdf_proj.crs) highway_buffer highway_proj.copy() highway_buffer[geometry] highway_proj.geometry.buffer(5000) # 单位米 # 用缓冲区去匹配县级面看哪些县被覆盖 joined gpd.sjoin(county, highway_buffer, howinner, predicateintersects) covered_counties joined[NAME].unique() print(高速 5 公里缓冲覆盖的县数量:, len(covered_counties)) print(covered_counties[:10])buffer(5000)里的 5000 是米因为已经转成投影坐标系。如果没转投影5000 就是度结果完全不对。sjoin的predicateintersects表示面与缓冲区有交集就算命中。howinner只保留匹配上的记录。这个流程可以扩展到服务区、收费站只要把点数据换成对应的图层。3.3 省、市、县三级聚合的两种写法做统计时经常需要从县级汇总到市级再从市级汇总到省级。如果数据里已经有PARENT字段直接dissolve就行如果没有得用空间关系或代码前缀匹配。代码前缀匹配最快市级ADCODE前四位等于县级ADCODE前四位省级前两位等于市级前两位。# 方法一用 ADCODE 前缀做聚合 county[city_code] county[ADCODE].astype(str).str[:4] city_stats county.groupby(city_code).agg( total_road(length_km, sum), total_area(area_km2, sum) ).reset_index() # 方法二用 dissolve 按市级面合并 city gdf_proj[gdf_proj[LEVEL] city] city_dissolved county.dissolve(bycity_code, aggfuncsum)方法一适合纯属性统计速度快方法二适合需要生成新几何体的场景比如把多个县合并成一个市的面。dissolve的aggfuncsum会对数值列求和但几何体合并后面积可能因为边界重叠而偏大需要事后校验。4. 避坑与排查SHP 处理中最容易翻车的 5 个点4.1 中文属性乱码现象是字段值显示为问号或方块现象打开属性表NAME列全是????或锟斤拷。原因SHP 的.dbf文件默认编码是GBK或GB18030但读取时用了UTF-8。解决读文件时显式指定encodinggb18030如果还不行用chardet检测编码后再读。写入时也要指定encodingutf-8否则导出给别人的文件又会乱码。4.2 路网重复线段导致长度翻倍现象某个县的公路密度算出来是其他县的十几倍明显不合理。原因路网数据里同一条路被拆成多段或者上下行分开存储裁剪后每段都算了一次长度。解决先按ROAD_ID或NAME去重或者用unary_union合并后再算长度。如果上下行是两条线做密度分析时应该只取一条或者除以 2。4.3 坐标系误判以为都是 WGS84 结果错位现象行政区划和路网叠图后路网整体偏移几百米。原因一个图层是WGS84另一个是GCJ02两者之间的偏移是非线性的。解决先看.prj文件如果缺失用已知地物比对。GCJ02转WGS84需要专门的转换库不能简单用to_crs解决。做分析前统一转成CGCS2000或WGS84。4.4 裁剪后面积变小但长度没变现象用县级边界裁剪路网后length_km和裁剪前一样。原因clip只裁剪几何体但属性表里的长度字段是裁剪前算好的不会自动更新。解决裁剪后重新用geometry.length计算长度覆盖旧字段。这个坑很隐蔽因为属性表看起来「有值」但值是错的。4.5 空几何和无效几何导致分析中断现象buffer或sjoin报错TopologyException或返回空结果。原因SHP 里存在自相交、重复点、空几何。解决先跑gdf.is_valid检查对无效几何用buffer(0)修复空几何直接dropna或过滤。修复后再做空间分析能省掉大量调试时间。5. 进阶技巧把这份 SHP 变成可复用的分析底座如果你不止做一次分析建议把这份 2022.07 的数据整理成一个标准化的 GeoPackage而不是每次读 SHP。GeoPackage 支持多图层、字段类型更规范、中文编码问题少而且读写速度比 SHP 快。我一般会建三个图层admin省、市、县三级带层级字段、road道路网带等级和长度、highway公路网带等级和宽度。建好之后后续做任何区域分析都直接从这个库取数不用再重复处理编码和坐标系。# 导出为 GeoPackage多图层存储 gdf_proj.to_file(sc_202207.gpkg, layeradmin, driverGPKG) roads_proj.to_file(sc_202207.gpkg, layerroad, driverGPKG) highway_proj.to_file(sc_202207.gpkg, layerhighway, driverGPKG) # 后续读取指定图层 admin gpd.read_file(sc_202207.gpkg, layeradmin) road gpd.read_file(sc_202207.gpkg, layerroad)另一个技巧是给行政区划加一个bbox字段存每个面的最小外接矩形做空间查询时先用bbox过滤再用几何精确判断速度能快几倍。对于路网可以预先算好每段的start_point和end_point做连通性分析时直接用点匹配不用每次解析几何。# 给行政区划加 bbox 字段加速空间查询 admin[bbox] admin.geometry.apply(lambda geom: geom.bounds) # 给路网加起终点字段 road[start_pt] road.geometry.apply(lambda geom: geom.interpolate(0)) road[end_pt] road.geometry.apply(lambda geom: geom.interpolate(geom.length))bounds返回(minx, miny, maxx, maxy)存成元组后可以用str或list类型写入 GeoPackage。interpolate(0)取线的起点interpolate(geom.length)取终点。这两个字段在做路径匹配和网络构建时非常有用省去每次重新计算。最后说一个我自己的习惯拿到任何 SHP先跑一遍gdf.info()和gdf.geometry.geom_type.value_counts()确认几何类型和空值情况再决定后续操作。这个动作花不到十秒但能避免后面几小时的调试。希望帮到你。本文还有配套的精品资源点击获取