简介这份资源是2017年中国县级行政区划的矢量边界数据集基于1:100万基础地理数据库制作面向从事GIS分析、城市规划、人口统计与灾害管理等工作的技术人员及地理信息专业学习者。压缩包共7个文件约46.33MB以SHP格式为核心配套DBF属性表、PRJ坐标系定义、SHX/SBN/SBX索引及XML元数据完整覆盖几何、属性与空间参考信息可直接在ArcGIS、QGIS等软件中加载使用。坐标系采用CGS_WGS_1984兼顾国内定位精度与全球兼容性便于与其他数据源叠加分析。目前已有615人学习下载。借助这份数据读者可快速获取精细到县级的行政边界用于专题制图、空间统计、区域规划及多源数据融合省去自行数字化与坐标校正的繁琐流程是开展大范围地理研究时较为实用的基础底图资料。1. 1100万基础地理数据库里的县级行政区一份2017年shp文件到底能干什么如果你手头正好拿到一份标注为“1100万基础地理数据库”的2017年县级行政区shp文件第一反应大概率不是兴奋而是犯嘀咕这玩意儿精度够不够坐标系是什么属性表里有没有我需要的行政区代码能不能直接叠到我的业务底图上我最初接触这类数据时也踩过同样的坑——兴冲冲加载进GIS软件结果发现和项目里用的在线底图偏了几百米排查半天才发现是地理坐标系和投影坐标系混用。1100万比例尺意味着全国范围一张图县级边界在这种尺度下必然经过综合取舍它不是测绘级精度但胜在覆盖完整、属性规整、年份明确特别适合做全国尺度的行政区划统计、专题图底图、空间连接时的行政单元挂接。这份2017年的县级行政区shp核心价值在于“全国县级单元一次性拿齐且带标准行政区划代码”对于做人口、经济、资源环境类数据空间化的人来说它是一个可靠的骨架。适合谁用做全国性数据可视化、需要按县聚合统计、或者要给点位数据反查行政归属的从业者。不适合谁做地块级权属确认、需要亚米级边界精度的工程测量。先把这个定位搞清楚后面所有操作才不会跑偏。2. 拿到shp先别急着画图坐标系、字段与几何自检2.1 为什么1100万县级行政区shp的坐标系必须第一个确认很多人拿到shp文件双击打开看到图形就以为万事大吉结果一做面积统计数值离谱到怀疑人生。1100万基础地理数据库的县级行政区数据常见做法是采用地理坐标系如CGCS2000或WGS84存储单位是度不是米。你直接用它算面积得到的是“平方度”没有任何实际意义。更隐蔽的问题是有些分发版本为了显示方便会预先投影成某种等积投影或兰勃特投影但属性表里不写清楚。我一般会先做三件事看.prj文件内容、在GIS软件里查图层属性、用Python读一遍CRS。如果.prj缺失或内容为空那就要警惕了这份数据可能被“洗”过坐标系信息丢了你得靠经验判断或者找元数据。import geopandas as gpd # 读取shp文件路径按实际修改 gdf gpd.read_file(rD:\data\county_2017\county_2017.shp) # 查看坐标系定义 print(CRS:, gdf.crs) # 查看几何类型和要素数量 print(几何类型:, gdf.geom_type.unique()) print(要素数量:, len(gdf)) # 查看属性表字段和前几行 print(字段列表:, gdf.columns.tolist()) print(gdf.head(3))这段代码的逻辑很直接先确认CRS是否存在再确认几何类型是否统一县级行政区通常是Polygon或MultiPolygon最后看字段名。参数方面gpd.read_file默认用fiona引擎如果遇到中文路径报错可以改用open方式或把数据放到纯英文路径下。如果gdf.crs返回None说明坐标系丢失这时候不要盲目设置先根据数据来源判断1100万全国数据2017年版本大概率是CGCS2000地理坐标系EPSG代码4490。但注意这只是常见情况不是绝对你得结合数据分发说明来定。2.2 属性表里哪些字段真正有用行政区代码、名称与层级县级行政区shp的属性表通常包含行政区代码、行政区名称、上级行政区代码、上级行政区名称、备注等字段。但不同分发版本字段名差异很大有的用拼音缩写有的用中文有的甚至只有代码没有名称。我见过最坑的一种情况字段名是“ADCODE99”这种带年份后缀的你以为是2017年代码结果里面混了更早的编码。所以第一步是字段清洗和核对。行政区代码是连接其他统计数据的关键6位县级代码如110101是标准但有些版本会给12位或带省代码前缀的变体。你需要确认代码位数和编码规则否则和统计年鉴挂接时会对不上。# 假设字段名为ADCODE和NAME实际按需替换 # 检查行政区代码位数分布 gdf[code_len] gdf[ADCODE].astype(str).str.len() print(gdf[code_len].value_counts()) # 检查名称字段是否有空值 print(名称空值数量:, gdf[NAME].isna().sum()) # 检查几何有效性 invalid_geom gdf[~gdf.geometry.is_valid] print(无效几何数量:, len(invalid_geom))这里的关键参数是code_len正常县级代码应该是6位。如果出现4位或12位说明编码体系不同需要做映射。几何有效性检查也很重要1100万数据经过综合偶尔会有自相交或悬挂节点虽然不影响显示但做空间分析时可能报错。如果发现无效几何可以用buffer(0)做一次修复但要注意buffer(0)可能改变微小边界对精度要求高的场景慎用。2.3 几何自检面积异常与飞地问题县级行政区里有一类特殊单元飞地。比如某些县在邻市有一块不相连的辖区。1100万数据里飞地通常被处理成MultiPolygon但有些版本会拆成多条记录导致同一个县出现两次。你按县聚合统计时如果不做dissolve就会重复计算。另一个常见问题是面积异常某个县的面积明显偏大或偏小可能是边界综合时把相邻县合并了或者漏掉了某个岛屿。我一般会先算一遍面积投影到等积坐标系后然后看分布把面积排名前后几位的县挑出来对照行政区划简册核对。# 投影到等积坐标系计算面积这里以Albers为例 gdf_proj gdf.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm) gdf_proj[area_km2] gdf_proj.geometry.area / 1e6 # 查看面积最大和最小的10个县 print(gdf_proj.nlargest(10, area_km2)[[NAME, area_km2]]) print(gdf_proj.nsmallest(10, area_km2)[[NAME, area_km2]])参数说明Albers投影的参数lat_1和lat_2是中国常用的双标准纬线lon_0取105度。如果你做的是分省研究可以改用该省中央经线。面积单位换算成平方公里后正常县的面积范围大概在几百到十几万平方公里之间如果出现个位数或超过20万基本可以判定有问题。这一步不是吹毛求疵而是提前把数据质量问题暴露出来避免后面统计结果被质疑。3. 从shp到可用底图裁剪、投影与空间连接实操3.1 按省裁剪县级行政区clip与intersection的选择全国1100万县级行政区shp文件如果你只研究某个省直接全量加载会拖慢速度也没必要。常见做法是按省裁剪。这里有两个函数容易混淆clip和intersection。clip是裁剪保留被裁剪要素在裁剪框内的部分属性表不变intersection是求交会生成新的几何属性表可能合并。对于县级行政区如果你只是想要某省的县用clip把省界作为裁剪框得到的是完整的县边界因为县界不会跨省除非飞地。但如果你要处理跨省飞地就得用intersection并且注意属性字段的保留。# 假设已有省界shp读取后取某个省 province gpd.read_file(rD:\data\province_2017\province_2017.shp) target_prov province[province[NAME] 某省] # 用clip裁剪 county_clip gpd.clip(gdf, target_prov) # 保存结果 county_clip.to_file(rD:\data\county_2017\county_target.shp, encodingutf-8)逻辑说明gpd.clip的第一个参数是被裁剪图层第二个是裁剪框。参数encodingutf-8是为了避免中文属性保存后乱码尤其在Windows环境下。如果裁剪后要素数量比预期少检查省界和县界的坐标系是否一致不一致先统一。另外clip对MultiPolygon的处理是保留整体不会拆分所以如果某个县部分在省界外clip后该县会保留省界内的部分但属性还是原来的县这可能导致面积偏小。这时候你要决定是接受这个近似还是用intersection精确切割。3.2 投影转换从地理坐标到Web墨卡托的坑做Web地图展示时你需要把地理坐标转成Web墨卡托EPSG:3857。这一步看似简单但有个经典坑Web墨卡托在高纬度地区面积变形极大如果你在Web地图上直接算面积会得到荒谬的结果。所以正确流程是展示用3857统计用等积投影。转换时用to_crs参数写EPSG代码或proj字符串。注意有些老版本geopandas对EPSG:3857支持有问题建议用proj字符串。# 转Web墨卡托用于展示 gdf_web gdf.to_crs(epsg3857) # 转等积投影用于统计 gdf_equal gdf.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm)参数说明epsg3857是Web墨卡托的标准代码projaea是Albers等积投影。转换后建议抽查几个已知点的坐标比如某县政府的经纬度转换后看是否落在合理范围。如果转换后图形跑到地球另一边多半是坐标系定义搞反了比如把投影坐标当成了地理坐标。3.3 空间连接把点位数据挂到县级行政区这是县级行政区shp最高频的用法你有一批带经纬度的点位比如企业、监测站想知道每个点属于哪个县。用spatial join但要注意点位数据的坐标系必须和县界一致。如果点位是WGS84经纬度县界也是地理坐标可以直接做如果县界已经投影点位也要先投影。连接方式用within即点在县内。# 假设点位数据为points已读取为GeoDataFrame # 确保坐标系一致 points points.to_crs(gdf.crs) # 空间连接 joined gpd.sjoin(points, gdf[[NAME, ADCODE, geometry]], howleft, predicatewithin) # 查看结果 print(joined[[point_name, NAME, ADCODE]].head())逻辑说明gpd.sjoin的predicatewithin表示点完全在县界内。如果点落在边界上可能匹配不到这时候可以用intersects但会带来重复匹配一个点落在两个县边界上会匹配两次。参数howleft保留所有点即使没匹配到县也保留方便排查。常见问题是点位在县界外一点点比如落在水域或邻县导致匹配为空。这时候可以加一个缓冲区但缓冲区大小要谨慎一般不超过100米。4. 避坑与排查县级行政区shp的5个血泪教训4.1 现象面积统计结果比官方数据大了一倍原因几何重复或未dissolve解决按行政区代码去重并dissolve我最早用这份数据做全国县域面积汇总时发现总数比官方公布的国土面积大了将近一倍。排查后发现数据里每个县有多条记录原因是飞地被拆成了独立要素但行政区代码相同。正确做法是先按行政区代码dissolve把同一个县的所有多边形合并成一个MultiPolygon再算面积。dissolve时注意保留名称和代码字段。# 按行政区代码合并 gdf_dissolved gdf.dissolve(byADCODE, aggfuncfirst).reset_index() # 重新计算面积 gdf_dissolved_proj gdf_dissolved.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm) gdf_dissolved_proj[area_km2] gdf_dissolved_proj.geometry.area / 1e6参数说明aggfuncfirst表示其他属性取第一条记录的值因为同一个县的名称应该一致。dissolve后几何类型可能变成MultiPolygon这是正常的。4.2 现象空间连接后部分点位没有匹配到县原因点位坐标系与县界不一致或点位在边界外解决统一坐标系并检查边界有一次帮同事处理一批监测点连接后发现有几十个点没有县名。检查发现点位数据用的是GCJ02火星坐标系而县界是WGS84两者偏差几百米导致点落在县界外。解决办法是把点位转成WGS84或者把县界转成GCJ02但后者涉及加密算法不建议。统一到WGS84后重新连接问题解决。另外如果点位确实在县界外比如在省界附近的飞地需要人工核对。4.3 现象裁剪后县的数量不对原因省界与县界的坐标系或拓扑不一致解决先统一坐标系再用空间查询核对裁剪时如果省界和县界坐标系不同clip会报错或结果为空。即使坐标系相同如果省界有拓扑错误比如自相交裁剪结果也会异常。我一般会先用gdf.crs province.crs判断不一致先转。然后检查省界的几何有效性无效的用buffer(0)修复。裁剪后用len(county_clip)和该省实际县数对比差一两个可能是省界飞地差很多就是数据问题。4.4 现象属性表中文乱码原因编码不一致解决读取和保存时指定encodingshp文件的属性表编码是个老问题。用geopandas读取时如果遇到乱码可以尝试encodinggbk或encodingutf-8。保存时也要指定否则在ArcGIS里打开可能乱码。我一般统一用utf-8但如果对方用ArcGIS可能需要gbk。最稳妥的办法是读进来后检查字段值如果乱码就换编码重读。# 尝试不同编码读取 try: gdf gpd.read_file(path, encodingutf-8) except: gdf gpd.read_file(path, encodinggbk)4.5 现象数据加载慢或内存溢出原因全量加载全国1100万数据解决按需读取或使用空间索引全国县级行政区有近3000个要素1100万比例尺下几何节点不算多但如果你同时加载多个图层内存还是吃紧。我一般会先只读需要的字段用columns参数或者用bbox参数只读某个范围。另外做空间连接前给县界建空间索引geopandas默认会建但如果你手动操作shapely对象记得用STRtree。# 只读必要字段 gdf gpd.read_file(path, columns[ADCODE, NAME, geometry]) # 按范围读取 gdf_sub gpd.read_file(path, bbox(100, 20, 120, 40))参数说明bbox是(minx, miny, maxx, maxy)单位与数据坐标系一致。如果数据是地理坐标就是经纬度范围。5. 进阶用法用2017县级行政区做时间序列对比与自动化出图5.1 跨年份对比2017版与其他年份的匹配与差异检测如果你手头还有别的年份的县级行政区shp可以做行政区划变更分析。常见做法是拿2017版和最新版做空间连接找出边界变化的县。但要注意行政区划调整包括撤县设区、合并、拆分等单纯靠空间叠加不一定能完全识别。我一般会先用行政区代码做左连接代码相同的直接匹配代码不同的再用空间相交面积比来判断。面积比超过90%的可能是同一区域改了代码面积比在50%到90%之间的可能是边界调整低于50%的可能是拆分或合并。# 假设gdf_2017和gdf_2020都有ADCODE和geometry # 先按代码合并 merged gdf_2017.merge(gdf_2020, onADCODE, howouter, suffixes(_2017, _2020)) # 对于代码不匹配的用空间连接找重叠 # 这里省略具体代码思路是计算相交面积参数说明空间相交面积比的计算需要先投影到等积坐标系。这一步比较耗时建议先按省筛选再逐省处理。5.2 自动化出图批量生成县级专题图做全国县级专题图手动一张张出图不现实。我一般用geopandas加matplotlib批量出图按省循环每个省一张图标注县名。关键参数是figsize和dpi全国图建议dpi不低于300省图可以200。中文标注需要设置字体否则显示方块。import matplotlib.pyplot as plt from matplotlib import font_manager # 设置中文字体 plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False # 按省循环出图 for prov_name in province[NAME].unique(): subset gdf[gdf[省名称] prov_name] fig, ax plt.subplots(figsize(10, 8), dpi200) subset.plot(axax, column指标值, legendTrue, cmapOrRd) ax.set_title(f{prov_name}县级专题图) plt.savefig(foutput/{prov_name}.png, bbox_inchestight) plt.close()逻辑说明column指定要映射的数值字段cmap是色带。bbox_inchestight去掉多余白边。如果县名标注重叠可以用adjustText库或者只标注重点县。这一步的坑在于有些省的县太多图例会挡住地图建议把图例放到图外。5.3 一个具体技巧用空间索引加速点位挂接当点位数据量达到几十万时直接sjoin会很慢。我一般会先给县界建空间索引然后用geopandas.tools.sjoin的rtree方式或者直接用shapely.STRtree。实测下来STRtree比默认的sjoin快3到5倍。具体做法是把县界的几何取出来建树然后遍历点位查询。from shapely.strtree import STRtree # 建索引 tree STRtree(gdf.geometry.values) # 查询每个点 for point in points.geometry: idx tree.query(point) if len(idx) 0: # 取第一个匹配的县 county gdf.iloc[idx[0]]参数说明tree.query返回的是几何对象的索引数组。注意STRtree查询是粗略的可能返回多个候选需要再用within精确判断。这个技巧在处理百万级点位时特别有用能省下大量时间。最后说个我自己的习惯每次拿到新的行政区划shp先跑一遍坐标系、字段、几何有效性、面积分布这四项检查把问题记在一个checklist里。这份2017年县级行政区数据我用了好几年最大的教训就是别假设它“应该没问题”而是假设它“一定有问题”然后逐个排除。希望帮到你。本文还有配套的精品资源点击获取