2010年全国乡镇界线矢量数据处理实战:从读取、清洗到空间分析
简介这份资源是2010年全国乡镇界线矢量数据集面向GIS从业者、城市规划人员、社会科学研究者及地图制图爱好者可用于历史行政区划变迁研究、区域对比分析与专题制图等场景。压缩包共含7个文件以Shapefile矢量格式为核心涵盖shp几何数据、dbf属性表、prj投影信息、shx索引及xml元数据等配套文件整体约103.56MB构成一套可直接在ArcGIS、QGIS等软件中读取分析的完整数据集。目前已有170人学习下载。数据以点、线、面形式精确刻画乡镇行政边界支持无损缩放、距离计算与缓冲区分析等操作便于研究人口分布、交通网络与资源管理等问题。需注意数据时点为2010年部分乡镇划分可能已有调整建议结合最新官方区划资料校对更新以更准确反映当前格局。1. 从一份 2010 年全国乡镇界线数据说起矢量数据到底能拿来干什么手里拿到一份标注为「2010年全国乡镇界线」的压缩包第一反应往往是这玩意儿能直接丢进 GIS 软件里用吗答案取决于你要做什么。这份数据本质上是全国乡镇级行政边界的矢量数据通常以 Shapefile 或 GeoJSON 格式组织包含乡镇名称、行政代码、所属区县等属性字段几何类型以面Polygon为主。它能解决的核心问题是当你需要按乡镇粒度做空间统计、制图或者叠加分析时不用自己从零数字化边界。适合做区域规划、人口经济数据落图、选址分析、专题制图的从业者也适合需要乡镇级底图做可视化展示的开发同学。但要注意2010 年的行政区划和现在有差异乡镇撤并、更名、界线调整都会导致属性对不上这是后面要重点聊的坑。2. 矢量数据格式拆解Shapefile 的组成与坐标系判断2.1 Shapefile 不是单个文件别只拷贝 .shp很多人第一次接触矢量数据看到.shp就以为是一个完整文件结果拷到另一台机器上打开发现属性表空了、投影丢了。Shapefile 实际上是一组文件的集合缺一不可。常见做法是把整个文件夹一起打包而不是只拿.shp。文件后缀作用缺失后果.shp存储几何形状无法显示图形.shx几何索引部分软件无法定位要素.dbf属性表丢失乡镇名称、代码等字段.prj坐标系定义无法正确投影坐标可能错位.cpg字符编码中文属性可能乱码拿到压缩包后先解压到一个纯英文路径下检查文件是否齐全。如果缺.prj你需要自己判断坐标系。2010 年前后的全国性行政边界数据常见的是 WGS84 地理坐标系EPSG:4326或者西安80、北京54等旧坐标系。判断方法很简单用 QGIS 或 ArcGIS 加载后看要素的坐标数值范围。如果经纬度在 73135、353 之间基本是地理坐标系如果数值是几十万到几百万那就是投影坐标系通常是高斯-克吕格投影。2.2 用 Python 快速读取并检查数据完整性不依赖桌面 GIS 软件直接用geopandas就能把数据读进来做体检。下面这段代码是我每次拿到新矢量数据必跑一遍的流程。import geopandas as gpd import os # 替换为实际解压后的 shp 路径注意路径不要有中文 shp_path rD:\data\2010_township\2010_township.shp # 读取矢量数据 gdf gpd.read_file(shp_path) # 基本信息体检 print(要素数量:, len(gdf)) print(几何类型:, gdf.geom_type.unique()) print(坐标系:, gdf.crs) print(字段列表:, gdf.columns.tolist()) # 查看前 5 行属性确认中文是否乱码 print(gdf.head()) # 检查几何有效性 invalid gdf[~gdf.is_valid] print(无效几何数量:, len(invalid)) # 检查空几何 empty gdf[gdf.is_empty] print(空几何数量:, len(empty))逻辑说明gpd.read_file()会自动识别同目录下的.shx、.dbf、.prj等配套文件所以路径指向.shp即可。gdf.crs返回坐标系信息如果是None说明.prj缺失或未被识别需要手动指定。is_valid用来排查自相交、悬挂节点等几何错误这类问题在早期数字化数据里很常见。is_empty检查有没有只有属性没有图形的要素。参数方面如果读进来发现中文乱码可以在read_file里加encodinggbk或encodingutf-8试一下。2010 年的数据用 GBK 编码的概率不低。如果坐标系显示为None可以用gdf.set_crs(epsg4326, inplaceTrue)强制指定但前提是你确认它确实是 WGS84。提示不要在原文件上直接做修改先复制一份再做清洗和转换避免原始数据被覆盖后无法回溯。3. 从读取到可用乡镇矢量数据的清洗与投影转换3.1 属性表清洗处理空值、重复与编码问题原始数据里属性表往往不干净。常见情况是乡镇名称为空、行政代码位数不对、同一乡镇出现多条记录。下面这段代码做一轮基础清洗。import geopandas as gpd import pandas as pd gdf gpd.read_file(rD:\data\2010_township\2010_township.shp, encodinggbk) # 查看关键字段的空值情况假设字段名为 NAME 和 CODE print(gdf[[NAME, CODE]].isnull().sum()) # 删除名称为空的记录 gdf gdf[gdf[NAME].notnull()] # 去除重复记录按行政代码去重保留第一条 gdf gdf.drop_duplicates(subset[CODE], keepfirst) # 统一名称字段格式去除首尾空格 gdf[NAME] gdf[NAME].str.strip() # 检查行政代码长度分布 print(gdf[CODE].astype(str).str.len().value_counts()) # 保存清洗后的数据 gdf.to_file(rD:\data\2010_township\cleaned.shp, encodingutf-8)逻辑说明isnull().sum()快速定位哪些字段有缺失。drop_duplicates按行政代码去重因为一个乡镇理论上只应有一条边界记录。str.strip()处理名称里混入的空格这种问题在 Excel 编辑过的数据里特别常见。行政代码长度检查是为了发现异常值乡镇级代码通常是 9 位或 12 位如果出现 6 位或更短可能是数据录入错误。参数说明to_file的encoding建议用utf-8兼容性更好。如果后续要在 ArcGIS 里打开且出现乱码可以改回gbk。保存时如果提示字段名过长可以用gdf.rename(columns{旧字段名: 新字段名})缩短字段名Shapefile 对字段名长度有限制。3.2 投影转换从地理坐标到适合制图的投影坐标地理坐标系经纬度适合存储但不适合做面积量算和距离分析。要做乡镇级面积统计必须转到投影坐标系。全国范围常用 Albers 等面积投影下面是对应的转换代码。import geopandas as gpd gdf gpd.read_file(rD:\data\2010_township\cleaned.shp) # 如果 crs 为空先指定为 WGS84 if gdf.crs is None: gdf.set_crs(epsg4326, inplaceTrue) # 转换到 Albers 等面积投影适合全国范围面积计算 # 中央经线 105双标准纬线 25 和 47 gdf_albers gdf.to_crs( projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm no_defs ) # 计算每个乡镇的面积单位平方公里 gdf_albers[area_km2] gdf_albers.geometry.area / 1e6 print(gdf_albers[[NAME, area_km2]].head()) # 保存投影后的数据 gdf_albers.to_file(rD:\data\2010_township\albers.shp, encodingutf-8)逻辑说明to_crs完成坐标系转换参数是一个 PROJ 字符串。Albers 投影的参数中lat_1和lat_2是双标准纬线lon_0是中央经线。这套参数是国内全国性制图常用的配置。面积计算前必须确保数据在投影坐标系下否则geometry.area算出来的是平方度没有实际意义。参数调整如果只做某个省的分析可以把lon_0改成该省中央经线比如做广东的可以改成 113。unitsm表示输出单位是米。转换后建议抽查几个已知面积的乡镇做验证比如找一个面积数据已知的乡镇对比计算结果是否在合理范围内。注意投影转换不可逆地引入变形不同投影的面积计算结果会有差异。如果要做严格的面积统计建议用 Albers 等面积投影如果只是做展示制图Web MercatorEPSG:3857也可以但不适合量算。4. 避坑与排查乡镇矢量数据最常见的五个翻车现场4.1 中文属性乱码打开全是问号现象用 ArcGIS 或 QGIS 打开后乡镇名称字段显示为乱码或问号。原因Shapefile 的.dbf文件对编码支持有限2010 年前后的数据多用 GBK 编码而部分软件默认按 UTF-8 读取。解决在 QGIS 里可以通过「图层属性 → 源 → 数据源编码」切换为 GBK在 Python 里读取时指定encodinggbk。如果已经读进来乱码了用gdf[NAME] gdf[NAME].encode(latin1).decode(gbk)尝试修复。4.2 坐标系丢失图形跑到海里去了现象加载后图形位置明显不对比如全国数据跑到了赤道附近或者非洲。原因.prj文件缺失软件默认按无投影处理把经纬度当成了米。解决先判断数据本身的坐标数值范围如果是经纬度手动指定 EPSG:4326如果是投影坐标需要找到对应的投影参数。实在不确定可以拿几个已知地点的坐标做比对。4.3 乡镇名称对不上属性连接全失败现象想把自己的业务数据按乡镇名称关联到边界上结果匹配率极低。原因2010 年至今乡镇撤并频繁名称变更、代码升级都会导致对不上。解决优先用行政代码做关联代码比名称稳定。如果代码也变了需要找行政区划变更对照表做映射。常见做法是保留原始名称字段同时新增一列标准化名称手动处理差异较大的记录。4.4 几何无效导致空间分析报错现象做相交、合并、缓冲区分析时软件报错或结果异常。原因早期数字化数据存在自相交、重复节点、悬挂边等几何错误。解决用gdf.is_valid排查对无效几何用gdf.buffer(0)做修复或者用shapely的make_valid方法。修复后再次检查有效性确认无误再做后续分析。4.5 数据量太大读取和渲染卡顿现象全国乡镇级数据要素数量可能在三万到五万之间直接全量加载到内存或渲染会很慢。原因矢量数据节点密集尤其是边界复杂的乡镇。解决如果只关注某个区域先用空间筛选裁切出目标范围如果用于 Web 展示建议做简化simplify并转为 GeoJSON 或 TopoJSON减小文件体积。简化时注意容差不要设太大否则边界会明显变形。5. 进阶用法用乡镇矢量数据做空间叠加与专题制图5.1 把业务数据落到乡镇边界上假设你有一份乡镇级的人口或经济数据想做成专题图。核心思路是用行政代码或名称做属性连接然后按数值分级渲染。下面是用 Python 做连接和分级统计的示例。import geopandas as gpd import pandas as pd import matplotlib.pyplot as plt # 读取乡镇边界 gdf gpd.read_file(rD:\data\2010_township\albers.shp) # 模拟业务数据实际使用时替换为你的 CSV 或 Excel data pd.DataFrame({ CODE: [110101001, 110101002, 110101003], population: [35000, 42000, 28000] }) # 按行政代码左连接 gdf_merged gdf.merge(data, onCODE, howleft) # 检查匹配率 matched gdf_merged[population].notnull().sum() print(f匹配成功: {matched} / {len(gdf_merged)}) # 按人口分级制图 fig, ax plt.subplots(1, 1, figsize(16, 12)) gdf_merged.plot( columnpopulation, axax, legendTrue, cmapOrRd, missing_kwds{color: lightgrey, label: 无数据}, edgecolorgrey, linewidth0.3 ) ax.set_title(乡镇人口分布专题图) ax.axis(off) plt.savefig(rD:\data\township_population.png, dpi300, bbox_inchestight) plt.show()逻辑说明merge做属性连接howleft保证边界数据不丢业务数据匹配不上的显示为空。missing_kwds把无数据区域标成灰色避免误读。cmapOrRd是橙色到红色的渐变色带适合表示数值高低。dpi300保证出图精度适合打印或报告使用。参数调整如果匹配率低先检查两边的代码格式是否一致比如一边是字符串一边是数字或者有前导零丢失。可以在连接前统一转成字符串并补齐位数。分级方式可以用schemequantiles做分位数分级或者schemefisher_jenks做自然断点分级后者在专题制图里更常用。5.2 用空间叠加做乡镇与流域的交叉分析另一个高频场景是把乡镇边界和自然地理要素如流域、地形做叠加分析每个乡镇落在不同流域内的面积占比。这种分析在资源环境领域很常见。import geopandas as gpd # 读取乡镇边界和流域边界 towns gpd.read_file(rD:\data\2010_township\albers.shp) basins gpd.read_file(rD:\data\basins\basins.shp) # 确保两者坐标系一致 basins basins.to_crs(towns.crs) # 空间相交得到每个乡镇与每个流域的重叠部分 intersect gpd.overlay(towns, basins, howintersection) # 计算重叠面积 intersect[overlap_km2] intersect.geometry.area / 1e6 # 按乡镇和流域分组汇总 result intersect.groupby([CODE, BASIN_NAME])[overlap_km2].sum().reset_index() # 计算每个乡镇内各流域的面积占比 town_total result.groupby(CODE)[overlap_km2].sum().reset_index() town_total.columns [CODE, total_km2] result result.merge(town_total, onCODE) result[ratio] result[overlap_km2] / result[total_km2] print(result.head(10))逻辑说明gpd.overlay做空间相交howintersection保留两者重叠部分。分组汇总后得到每个乡镇在每个流域内的面积再除以乡镇总面积得到占比。这套流程可以套用到土地利用、生态分区等多种叠加分析场景。参数说明overlay对几何有效性要求较高如果报错先用buffer(0)修复。如果数据量大叠加运算可能较慢可以先用clip按研究区范围裁切再操作。5.3 一个我踩过的坑别用地理坐标系做面积统计刚接触矢量数据那会儿我直接拿 WGS84 经纬度数据算面积结果一个乡镇算出来零点几还以为是数据错了。后来才反应过来经纬度下的面积单位是平方度不是平方米。从那以后我每次做面积相关的分析都强制先转投影坐标系再算面积。这个习惯帮我省了很多返工的时间。如果你手里正好有这份 2010 年全国乡镇界线数据建议先按第 2 章的流程做一次完整体检再根据实际需求决定是否做投影转换和属性清洗。数据本身不复杂但细节决定它能不能真正用起来。希望帮到你。本文还有配套的精品资源点击获取