尧图网站建设 
简介这份资源是面向地理空间分析学习者的全国地理分区SHP矢量数据集以标准矢量格式覆盖中国九大流域、三大地理分区、六大地理区域、气候区划、农业区划及林业工程空间分布适合水资源管理、生态区划、农业规划等教学与科研场景。压缩包共65个文件总大小40.48MB主要由shp、dbf、prj、sbn、sbx、shx等矢量配套文件构成并配有jpg预览图、说明文档及将SHP转为JSON的Python脚本这些格式分别承载几何图形、属性表、投影坐标与空间索引可直接在ArcGIS、QGIS等软件中打开分析。目前已有207人在线学习与浏览。除矢量数据外附带的Python转换脚本能批量输出JSON便于Web地图展示与二次开发随包说明简明标注图层含义附赠辅助资料压缩包则提供额外参考价值。整套资料适合GIS、遥感、环境科学方向的学生、教师和研究人员既能用于地图制图与区域对比也可支撑流域治理、气候区划等专题建模。1. 全国地理分区数据集SHP矢量数据为什么水文气候分析都先拿它开刀一套全国地理分区数据集SHP矢量数据把三大分区、九大流域和气候区带拆成了可直接读取的边界文件再配上一段能落地的Python代码——这几年我做水文气象方面的空间统计最常用的底图就是这套东西。过去处理塔里木河流域的降水插值、给南盘江流域的气象站做分组最头疼的不是算法而是手上没有标准边界经费有限时只能拿省界硬凑。这套数据能解决的就是这个问题三大自然区、九大流域、气候区带三套边界叠加使用配合Python脚本做筛选、裁剪、统计和格式转换。适合水文、气象、环保、规划方向的从业者和学生尤其是那些有GIS基础、但不想在ArcGIS里反复点鼠标的人。2. 底子先摸清三大自然区、九大流域与气候区带的边界结构与坐标系2.1 这套数据集里到底有什么从文件名到字段组织拿到压缩包后先别急着跑代码把文件结构理一遍。常见目录一般是这样组织的three_regions.shp三大自然分区东部季风区、西北干旱半干旱区、青藏高寒区nine_basins.shp九大流域片松辽、海河、淮河、黄河、长江、珠江、东南诸河、西南诸河、西北诸河climate_zones.shp气候区带按温度带和干湿区组合如中温带半干旱区、亚热带湿润区每个SHP文件都带着配套的.dbf、.shx、.prj。.prj里写的是坐标系这个文件很容易被忽略但它是后面所有空间计算的地基。这套数据里字段命名不一定统一——有的叫NAME有的叫名称有的叫REGION——我每次拿到数据的第一件事就是打印字段列表看结构。import geopandas as gpd regions gpd.read_file(three_regions.shp, encodingutf-8) basins gpd.read_file(nine_basins.shp, encodingutf-8) climates gpd.read_file(climate_zones.shp, encodingutf-8) for name, gdf in [(regions, regions), (basins, basins), (climates, climates)]: print(name, gdf.columns.tolist()) print(gdf.head(3))这段代码用GeoPandas把三个SHP读进来先看列名和头几行。为什么先做这一步属性字段名决定了后面怎么写筛选条件。比如九大流域里“长江”字段写的是“长江流域”还是“长江”直接关系到SQL查询写法和通配符的使用。编码问题也在这里暴露有些SHP的.dbf属性表是GBK编码直接读会乱码要改成encodingGBK。字段里通常还有面积属性但要注意这个面积是数据生产方算的单位可能是平方千米也可能是不规则单位。正规数据集会在元数据文档里说明没文档的话宁可自己用投影坐标重算一遍也别直接用原始面积字段做统计。九大流域之间用名称字段区分三大分区和气候区带靠各自的分类字段这三个文件互相独立又可以用空间位置叠加。2.2 用GeoPandas快速预览SHP属性表、投影与边界检查读进来之后要先确认几何类型和边界范围。有的SHP看着是面实际里面混了线要素或者GeometryCollection里藏着多部件面这些都会在后续空间操作里埋雷。快速检查方法如下:print(regions.crs) # 查看原始坐标系 print(regions.total_bounds) # 查看整体边界范围 [minx, miny, maxx, maxy] print(regions.geometry.type.value_counts()) # 几何类型统计 print(basins.geometry.is_valid.sum(), valid /, len(basins), total)total_bounds是个非常实用的检查如果范围不是中国国界附近说明投影或裁剪有问题。常见的中国边界范围差不多在73°E-135°E、18°N-54°N附近如果你看到的是几百上千的坐标值说明这是投影坐标而不是经纬度。几何类型统计能筛出多部件面is_valid检查的是拓扑合法自相交的面在后续裁剪时结果会悬空这一行就能筛出来。拿到手的数据如果有几个面是invalid的后面第5章会讲怎么修。这套数据默认坐标系大概率是WGS84经纬度EPSG:4326或CGCS2000EPSG:4490。你要是只做查询和出图用原始坐标系就行做面积统计或裁剪就必须转投影坐标系。切记不要拿经纬度去做面积计算1度经度对应地面的距离随纬度变化在黑龙江和海南算出来的结果差好几倍这种错误在论文里时有发生。2.3 坐标系统的选择经纬度存储、投影计算的取舍我一般习惯分两套坐标体系用存储和展示用原始地理坐标也就是EPSG:4326或EPSG:4490面积计算和栅格裁剪用Albers等积圆锥投影Albers是处理中国全域面积统计最常用的投影没有之一。国家基础地理数据经常用它因为双标准纬线设置在25°N和47°N在国界范围内面积变形控制在很小的范围内。代码里用GeoPandas重投影albers_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumCGCS2000 unitsm no_defs basins_albers basins.to_crs(albers_crs) basins_albers[area_km2] basins_albers.geometry.area / 1e6 print(basins_albers[[NAME, area_km2]])to_crs重投影后geometry.area返回的单位是平方米除以1e6就是平方千米。这比直接使用属性表里的面积字段可靠因为投影基准、舍入方式全部由你自己控制。如果你处理的只是北京市这种较小范围也可以用UTM分区投影但代际差异不大。要注意的是重投影会改变几何精度如果到后面还要和其他地理坐标数据做空间连接记得把结果再转回经纬度。这三个文件之间的空间关系是嵌套的——同一地理位置既是某个气候区带又是某个流域的一部分。所以做叠加分析时用哪个图层做底都行但统计口径必须一致。最稳妥的做法是把流域作为基础图层气候区带作为属性挂接上去这样最后每个流域多边形都携带了气候区带信息不需要反复做空间连接。3. 用Python做流域级空间分析从读SHP到掩膜裁剪的完整流程3.1 筛选指定流域属性查询与要素抽取典型场景你想提取塔里木河流域边界用它裁剪降水栅格再统计整个流域的年降水量。第一步是属性筛选。basins_albers gpd.read_file(nine_basins.shp, encodingutf-8).to_crs(albers_crs) # 模糊匹配兼容字段值带流域后缀的情况 target 塔里木河 basin basins_albers[basins_albers[NAME].str.contains(target)].copy() print(basin[NAME].tolist())用str.contains而不是等于判断是为了兼容字段值可能是“塔里木河流域”或“塔里木河”的情况实测下来这种写法容错率最高。筛选结果可能不止一个要素——如果数据集把塔里木河流域拆成了多个子流域tolist()会把所有匹配名称显示出来这时就需要人工确认是取全流域还是只取其中一个子流域。全流域就用gpd.unary_union把多个要素合并成一个多边形。这一步最容易被忽略的是筛选前一定要copy()否则后续操作可能触发链式赋值的警告。还有一点如果你之后要把边界坐标点导出来写论文插图直接操作的是投影后的坐标单位是米不是经纬度输出前记得转回去。basin_union basin.geometry.unary_union print(Target basin area:, basin_union.area / 1e6, km2)unary_union把同属一个流域的所有多边形融合成单个面。融合后算一次面积该数字可作为后续裁剪结果的参照——如果裁剪出来的有效像元总面积和这个数字差太远说明裁剪过程中出了问题。3.2 用流域边界裁剪栅格Rasterio掩膜裁剪和参数设置拿到流域边界后下一步裁降水或气温栅格。常见的数据是NetCDF转的GeoTIFF或直接从气象共享平台下载的降水GeoTIFF。import rasterio from rasterio.mask import mask as rio_mask with rasterio.open(precipitation.tif) as src: out_image, out_transform rio_mask.mask( src, [basin_union], # 必须是GeoJSON-like的几何列表 cropTrue, nodata-9999, all_touchedFalse ) out_meta src.profile.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(tarin_basin_precip.tif, w, **out_meta) as dst: dst.write(out_image)参数的坑先从这几个说起cropTrue把输出范围裁剪成和流域外接矩形一致文件体积小很多nodata-9999保持和原栅格一致避免把无效值误当成0all_touched默认False意味着只有像元中心落在流域内才算有效。栅格分辨率低时细长的流域边界的边缘像元会被切掉这时可考虑设True但边缘会变得更粗糙[basin_union]必须放在列表里因为mask函数是按图层要素列表逐要素处理的直接传Polygon它会报错裁剪后的栅格范围是矩形流域外区域仍然是nodata。后边做统计的时候要先用有效值掩膜把流域外的nodata滤掉不然最小值算出来一定是-9999平均降水也会被拉低。3.3 分区统计每个气候带/流域的均值与总量计算裁剪只是中间步骤最终要落到分区统计。最常用的方案是rasterstats这个库它内置了zonal_stats函数对每个多边形计算栅格像元的统计量from rasterstats import zonal_stats stats zonal_stats( basins_albers, # 多边形图层 tarin_basin_precip.tif, stats[mean, sum, count], nodata-9999, geojson_outTrue ) print({f[properties][NAME]: f[properties][mean] for f in stats})zonal_stats内部会先做栅格化和重投影。多边形建议用与栅格相同的投影否则它内部会做一次矢量重投影结果可能不准。stats参数里mean是像元均值sum是全部有效像元的累加常用于计算区域降水总量count统计有多少个像元参与了计算可辅助判断结果是否合理——如果某个流域的count异常少大概率是边界与栅格没对齐。rasterstats的一个缺点是对非常大、要素非常多的矢量图层处理偏慢几千个面可能会跑几分钟。如果做全国尺度的气候带统计建议先把栅格裁成和矢量边界同一范围再跑别直接喂整张图。血泪经验我一开始图省事对全国1km分辨率降水栅格做31个省级行政区统计跑了快半个钟头后来改成先裁剪再统计几分钟搞定。如果要按气候区带分组统计各流域的平均气温场景就是空间叠加加属性分组。做法是先把气候区带和流域两个图层做空间连接再用groupby聚合overlay_result gpd.overlay(basins_albers, climates_albers, howintersection) grouped overlay_result.groupby(CLIMATE)[area_km2].sum()gpd.overlay会割开相交的面输出的是两块区域叠加后的细分多边形每个多边形仍保留两个父图层的属性字段。此时每个“流域-气候带”组合的面积就可以直接用groupby求总和。这套流程一步到位比在ArcGIS里做两轮相交、联接快得多。4. 格式转换与格式互转SHP转KML、JSON、TXT与点表生成4.1 SHP转KML/JSON几种驱动方式的参数差异用途不同选不同的转换路线丢到Google Earth或手机地图里看转KML做Web前端或传给后端接口转GeoJSON做简单属性导出或交给不装GIS软件的人转TXT表格GeoJSON最简单一行代码解决:basins.to_file(nine_basins.geojson, driverGeoJSON, encodingutf-8)转KML有两个路线。一个是GDAL内置的KML驱动通过GeoPandas的to_file指定扩展名和驱动程序即可basins_simple basins[[NAME, geometry]].copy() basins_simple.to_file(nine_basins.kml, driverKML)但GDAL的KML驱动有个麻烦它只输出有限的字段而且属性值里如果有非法XML字符比如“”“”这些符号写入会失败或生成的KML打不开。实测下来字段多了经常翻车。我一般用simplekml库生成KML字段控制和样式定制更灵活import simplekml kml simplekml.Kml() for _, row in basins.iterrows(): # 取第一个多边形环的坐标舍弃空洞 coords list(row.geometry.exterior.coords) kml.newpolygon( namerow[NAME], outerboundaryiscoords ) kml.save(nine_basins.kml)这里用exterior.coords取外环坐标把内环空洞丢掉了。对流域这种通常没有内部空洞的面要素影响不大但遇到带复杂的岛屿要素时信息会丢建议确认原始几何类型再决定。带洞的面要素会把洞填成实心看着倒无所谓但极值统计时会把周围要素覆盖。热词里常见的“arcgis shp转kml”其实就是这种需求在ArcGIS里另存为KML在Python里用上面两段代码替换即可。TXT导出更直接适合生成坐标表with open(basins_coords.txt, w, encodingutf-8) as f: for _, row in basins.iterrows(): xy ,.join([f{x:.4f},{y:.4f} for x, y in row.geometry.exterior.coords]) f.write(f{row[NAME]}|{xy}\n)注意TXT里保留的是经纬度如果你想在CAD里直接打开做成线稿得先转投影再输出坐标。转TXTSHP转TXT的好处是能在Excel里打开后排序、筛选、美化表格。坐标精度保留4位小数约合10米做概览图足够做精确放样没必要。4.2 Excel/CSV点表转SHP从经纬度列到矢量点的完整代码这是另一个高频需求手里有一份Excel表格两列是经纬度一列是站点名/测值要转成SHP做空间插值或叠加分析。GeoPandas提供了非常直接的路径import pandas as pd import geopandas as gpd from shapely.geometry import Point df pd.read_excel(stations.xlsx) # 必须包含 lon, lat 两列 geometry [Point(x, y) for x, y in zip(df[lon], df[lat])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) gdf_valid gdf.dropna(subset[lon, lat]) gdf_valid gdf_valid[(gdf_valid[lon].between(73, 135)) (gdf_valid[lat].between(18, 54))] gdf_valid.to_file(stations.shp, encodingutf-8)注意两点坐标范围过滤是我每次必做的避免个别站点经纬度填反、变成坐标跑到国外或海洋里dropna要先执行因为导入时空值的坐标行会在zip阶段直接报错。Excel里常见坑是“经纬度列名带空格”或“单位不是十进制度而是度分秒”。度分秒要先用公式换算十进制度度分/60秒/3600没换算直接进GeoDataFrame画出来的点会全部挤在坐标原点附近是个非常经典的翻车现场。Arcgis里手动做“arcgis excel点转shp”是菜单操作但数据量大或需要反复更新时会想吐。Python方式可以从读取Excel到写SHP一个脚本跑完换数据只要改路径。另外这套脚本稍加改动就能接CSVpd.read_csv即可。4.3 DWG/DXF与SHP互转坐标系漂移的处理方法DWG转SHP是另一个高频热搜词。直说结论GDAL不原生支持DWG格式别指望一行代码装就能转。常见做法分两步走在CAD里用SAVEAS命令把DWG另存为DXF选ASCII版本再用GDAL/OGR或GeoPandas读取DXF转SHPfrom osgeo import ogr ogr.UseExceptions() src_ds ogr.Open(boundary.dxf) if src_ds is None: raise RuntimeError(DXF打开失败检查文件是否被CAD占用) out_ds ogr.GetDriverByName(ESRI Shapefile).CreateDataSource(boundary.shp) src_layer src_ds.GetLayer(0) out_layer out_ds.CreateLayer(boundary, geom_typeogr.wkbPolygon) out_layer.CreateField(ogr.FieldDefn(layer, ogr.OFTString)) for feat in src_layer: geom feat.GetGeometryRef() lyr_name feat.GetField(Layer) out_feat ogr.Feature(out_layer.GetLayerDefn()) out_feat.SetGeometry(geom) out_feat.SetField(layer, lyr_name) out_layer.CreateFeature(out_feat)这段代码的最大意义是保留CAD图纸里的图层名。DXF文件里不同的图层通常代表不同的要素类型——有的是边界、有的是建筑、有的是标注——转SHP时如果不把Layer字段带出来所有要素混在一起没法用。转换后还有一步必须做坐标系核对。CAD图纸普遍不带投影信息转出来的要素如果和全国边界叠不上大概率是坐标系没对齐。可以用几个已知地理坐标点做参照反向校准或者把转换结果和这套地理分区数据用gpd.overlay做一次叠加检查跨度超过几百米就相当可疑。同理“shp文件转mapgis线文件”这类需求本质也是格式转换但是要注意SHP的线要素在MAPGIS里对应线文件面要素对应区文件属性字段也可能丢失不要指望无损往返。SHP转MAPGIS线文件的方法通常是先导出DXF再导入MAPGIS转换后检查要素个数和图层。这里不做展开但记得一个核心原则转换格式永不回头导出前先把属性该合并的合并该删的删。4.4 其他衍生场景SHP转3D Tiles、渔网分割与JSON互转的边界提醒热搜词里还有“shp转3dtiles”“渔网分割shp”“json转shp网站”这些实际都是同一类问题的变体SHP转3D Tiles常用在Cesium展示三维场景。常见做法是先把SHP要素转为带高度的GeoJSON或glTF再用CesiumLab等工具瓦片化。那个流程吃显卡和CPU对要素数量有上限要求几十万以上级别的面要素直接卡死先化简顶点再转渔网分割SHP是把研究区切成等尺寸网格再算网格内属性。QGIS的创建网格工具或者ArcGIS的渔网工具都能生成网格面然后用交集分割矢量。这里有一个常见误区渔网做成后网格面要素和原始SHP相交需要先用空间索引过滤不然几千个网格对几十万条边界的交集运算能跑上一个下午JSON转SHPGeoJSON直接gpd.read_file(xxx.geojson)就吃进来了不需要“网站”转。浏览器端在线转小文件尚可数据量大时浏览器内存直接爆掉这些衍生工具链我都不在这篇里写全核心思路不变先确认源数据的坐标系和目标格式的限制再选择转换驱动。5. 避坑实录矢量数据集使用中常见的五个翻车现场5.1 中文属性乱码读出来全是“锟斤拷”现象gpd.read_file读SHP后NAME字段显示“锟斤拷”或“锟斤拷锟斤拷”过滤条件怎么写都匹配不上。原因数据生产方用的编码是GBK或GB2312GeoPandas默认按UTF-8解码编码错位。解决读文件时显式传encodingGBK如果不确定原始编码用文本编辑器打开.dbf文件看头部信息或直接用Python循环尝试多种编码for enc in [utf-8, GBK, GB18030]: try: gdf gpd.read_file(nine_basins.shp, encodingenc) if 流域 in str(gdf[NAME].iloc[0]): print(OK:, enc) break except UnicodeDecodeError: continue顺便说一个我的习惯自己导出SHP时一律用encodingutf-8这是最通用的方案但给别人发数据前我会提醒一句对方如果有乱码就试试GBK。编码问题解决不了的时候就直接读dbf属性表单独处理用dbfread把dbf转成pandas DataFrame再和几何表合并。5.2 面积计算离谱没投影就统计的锅现象用原始SHP直接算几何面积结果显示塔里木河流域面积数千亿平方公里明显是地级市级别数据却算出比地球还大的结果。原因几何坐标是经纬度geometry.area在未投影的GeoDataFrame里返回的是度数单位的“面积”也算得出来但没有任何物理意义。解决无论数据源是什么坐标系统计面积前一律转Albers等积投影。我的检验标准是流域级面积算完和官方公报对一下三大流域中长江流域面积约180万平方公里黄河流域约75万平方公里大概超过合理范围的10%就说明投影有问题。def safe_area_km2(gdf): from copy import deepcopy tmp deepcopy(gdf) tmp tmp.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs) return tmp.geometry.area / 1e6用deepcopy是为了不修改原始GeoDataFrame的几何对象。这个函数在项目里反复用沉淀成了我的“后悔药函数”。5.3 流域与省界套不上叠加分析结果为空现象用流域边界剪省界数据结果输出图层是空的或者相交面积比预想小了两个数量级。原因两套SHP的坐标系不同或者边界精度版本差异很大——县级边界和流域边界在边缘有几十米的偏差小面要素相交时可能完全落空。南盘江流域边界比较曲折如果用的省界是最新版高精度数据两边误差叠加部分沿岸的窄条相交结果为0。解决先统一坐标系再叠加不要直接拿两套原始SHP做overlay。对边界接缝做一次拓扑容差缓冲basins_buffered basins.copy() basins_buffered.geometry basins_buffered.buffer(0.01) # 0.01度约1km overlay_result gpd.overlay(basins_buffered, provinces, howintersection)缓冲会扩大边界范围但注意缓冲操作也在修改几何量大的时候记得把缓冲后超出研究区的部分裁掉。还有一个做法是切换“省1”和“省2”两种版本的数据集所谓省1、省2通常是不同精度版本节点密度和属性字段组织有差异就看哪个版本和你的流域边界在坐标基准上更一致实际以你对交测试为准不要凭文件名猜。5.4 SHP自动拓扑修复shapechk与Geometry Collection重建现象读取SHP报错“Corrupt geometry”或者某些要素无法输出shapefile还打不开。原因SHP文件本质由多个子文件组成.shp、.dbf、.shx任何一个损坏或非法都会引发读取失败。拓扑错误则常在多边形自相交、重复节点、闭合环起点不一致时发生。解决先用shapechk检测SHP文件头的完整性修复后再进入GeoPandas处理。对一个存有拓扑错误的面用buffer(0)重置几何是经典处理bad_rows gdf[~gdf.geometry.is_valid].index print(Invalid rows:, len(bad_rows), bad_rows[:10].tolist()) gdf.loc[bad_rows, geometry] gdf.loc[bad_rows].geometry.buffer(0)buffer(0)的原理是对自相交的面做一次零宽度的缓冲操作让GIS引擎重建几何节点消除自相交和重复点。它能修复大部分拓扑问题但如果面本身有大量细碎的狭长裂缝最好先做一次simplify再buffer。修复后重新检查is_valid还不行就放弃而只保留大面主体的算法见很多GIS教材我遇到极端情况直接改用R语言的sf包做相同处理有时GDAL修不了的它能修好。5.5 KML写入报错或文件打不开字段限制带来的玄学问题现象SHP转KML后Google Earth打不开或者打开后要素名显示blank、属性表全空。原因GDAL的KML驱动对属性字段数量有限制超出后会被踢掉字段值里如果有非UTF-8字符或XML非法字符、、引号写入直接报错。解决转换前精简字段并做字符过滤import re def clean_xml_text(value): if value is None: return return re.sub(r[^\w\s\u4e00-\u9fff], , str(value), flagsre.UNICODE) basins_simple basins[[NAME, geometry]].copy() basins_simple[NAME] basins_simple[NAME].map(clean_xml_text)清理字段名长度、过滤非法字符后再调用to_file(..., driverKML)。还有一个细节KML驱动写入时会对字段名自动截断所以输出前把字段名改成短英文是值得的比如NM、AREA。从那以后我做任何转换第一步都是“精简字段过滤字符”这成了肌肉记忆。6. 进阶技巧按省界流域双重分组批量出图与数据验收习惯6.1 一个脚本生成多个分幅图避免手动重复劳动当天底下数据拿到手最占时间的是挨个区域出图。比如给领导汇报淮河、黄河、长江各流域的降水分布图逐个手动选范围出图很费劲。我通常的做法是按流域分组遍历生成地图import matplotlib.pyplot as plt for name, grp in basins.groupby(NAME): fig, ax plt.subplots(figsize(6, 6)) grp.plot(axax, colorlightblue, edgecolorblack) ax.set_title(name, fontsize14) fig.savefig(fbasin_{name}.png, dpi150, bbox_inchestight) plt.close(fig)如果需要按省界流域双重分组就先用空间连接把流域名称挂到省界SHP上再按拼接字段分组。这样能一次性输出诸如“安徽-淮河流域”这类文件名完整的图。缺点是用Matplotlib画的符号不好看但胜在自动化、速度快。追求美观就导出GeoJSON到QGIS里套模板批量打印。6.2 每个新版本数据到手先过这三个流程再动手分析第一核对坐标系和边界范围打印crs、total_bounds和已知参考SHP做一次叠加对比。第二检查拓扑有效性和字段类型过滤字段值里的空值、异常值is_valid.sum()一眼看出问题。第三重算面积并和官方公开数据进行对照面积差率超过5%就要检查投影或边界是否正确。那套“省1/省2”两个版本的省界线在字段组织和精度上有细微区别不检查直接混用后面做任何带省界的统计都会“差之毫厘失之千里”。这套数据同样如此等到的版本越新字段名、面积单位越可能有变化我的应对是写一个自动校验函数每次拿到压缩包先跑一遍格式检测再进正式流程。从那以后我处理任何SHP数据集都强制走一遍“读属性、查坐标系、验拓扑、算面积”这四个动作希望帮到你。本文还有配套的精品资源点击获取 全国地理分区SHP矢量数据与Python空间分析：从边界处理到栅格裁剪的完整实践
简介这份资源是面向地理空间分析学习者的全国地理分区SHP矢量数据集以标准矢量格式覆盖中国九大流域、三大地理分区、六大地理区域、气候区划、农业区划及林业工程空间分布适合水资源管理、生态区划、农业规划等教学与科研场景。压缩包共65个文件总大小40.48MB主要由shp、dbf、prj、sbn、sbx、shx等矢量配套文件构成并配有jpg预览图、说明文档及将SHP转为JSON的Python脚本这些格式分别承载几何图形、属性表、投影坐标与空间索引可直接在ArcGIS、QGIS等软件中打开分析。目前已有207人在线学习与浏览。除矢量数据外附带的Python转换脚本能批量输出JSON便于Web地图展示与二次开发随包说明简明标注图层含义附赠辅助资料压缩包则提供额外参考价值。整套资料适合GIS、遥感、环境科学方向的学生、教师和研究人员既能用于地图制图与区域对比也可支撑流域治理、气候区划等专题建模。1. 全国地理分区数据集SHP矢量数据为什么水文气候分析都先拿它开刀一套全国地理分区数据集SHP矢量数据把三大分区、九大流域和气候区带拆成了可直接读取的边界文件再配上一段能落地的Python代码——这几年我做水文气象方面的空间统计最常用的底图就是这套东西。过去处理塔里木河流域的降水插值、给南盘江流域的气象站做分组最头疼的不是算法而是手上没有标准边界经费有限时只能拿省界硬凑。这套数据能解决的就是这个问题三大自然区、九大流域、气候区带三套边界叠加使用配合Python脚本做筛选、裁剪、统计和格式转换。适合水文、气象、环保、规划方向的从业者和学生尤其是那些有GIS基础、但不想在ArcGIS里反复点鼠标的人。2. 底子先摸清三大自然区、九大流域与气候区带的边界结构与坐标系2.1 这套数据集里到底有什么从文件名到字段组织拿到压缩包后先别急着跑代码把文件结构理一遍。常见目录一般是这样组织的three_regions.shp三大自然分区东部季风区、西北干旱半干旱区、青藏高寒区nine_basins.shp九大流域片松辽、海河、淮河、黄河、长江、珠江、东南诸河、西南诸河、西北诸河climate_zones.shp气候区带按温度带和干湿区组合如中温带半干旱区、亚热带湿润区每个SHP文件都带着配套的.dbf、.shx、.prj。.prj里写的是坐标系这个文件很容易被忽略但它是后面所有空间计算的地基。这套数据里字段命名不一定统一——有的叫NAME有的叫名称有的叫REGION——我每次拿到数据的第一件事就是打印字段列表看结构。import geopandas as gpd regions gpd.read_file(three_regions.shp, encodingutf-8) basins gpd.read_file(nine_basins.shp, encodingutf-8) climates gpd.read_file(climate_zones.shp, encodingutf-8) for name, gdf in [(regions, regions), (basins, basins), (climates, climates)]: print(name, gdf.columns.tolist()) print(gdf.head(3))这段代码用GeoPandas把三个SHP读进来先看列名和头几行。为什么先做这一步属性字段名决定了后面怎么写筛选条件。比如九大流域里“长江”字段写的是“长江流域”还是“长江”直接关系到SQL查询写法和通配符的使用。编码问题也在这里暴露有些SHP的.dbf属性表是GBK编码直接读会乱码要改成encodingGBK。字段里通常还有面积属性但要注意这个面积是数据生产方算的单位可能是平方千米也可能是不规则单位。正规数据集会在元数据文档里说明没文档的话宁可自己用投影坐标重算一遍也别直接用原始面积字段做统计。九大流域之间用名称字段区分三大分区和气候区带靠各自的分类字段这三个文件互相独立又可以用空间位置叠加。2.2 用GeoPandas快速预览SHP属性表、投影与边界检查读进来之后要先确认几何类型和边界范围。有的SHP看着是面实际里面混了线要素或者GeometryCollection里藏着多部件面这些都会在后续空间操作里埋雷。快速检查方法如下:print(regions.crs) # 查看原始坐标系 print(regions.total_bounds) # 查看整体边界范围 [minx, miny, maxx, maxy] print(regions.geometry.type.value_counts()) # 几何类型统计 print(basins.geometry.is_valid.sum(), valid /, len(basins), total)total_bounds是个非常实用的检查如果范围不是中国国界附近说明投影或裁剪有问题。常见的中国边界范围差不多在73°E-135°E、18°N-54°N附近如果你看到的是几百上千的坐标值说明这是投影坐标而不是经纬度。几何类型统计能筛出多部件面is_valid检查的是拓扑合法自相交的面在后续裁剪时结果会悬空这一行就能筛出来。拿到手的数据如果有几个面是invalid的后面第5章会讲怎么修。这套数据默认坐标系大概率是WGS84经纬度EPSG:4326或CGCS2000EPSG:4490。你要是只做查询和出图用原始坐标系就行做面积统计或裁剪就必须转投影坐标系。切记不要拿经纬度去做面积计算1度经度对应地面的距离随纬度变化在黑龙江和海南算出来的结果差好几倍这种错误在论文里时有发生。2.3 坐标系统的选择经纬度存储、投影计算的取舍我一般习惯分两套坐标体系用存储和展示用原始地理坐标也就是EPSG:4326或EPSG:4490面积计算和栅格裁剪用Albers等积圆锥投影Albers是处理中国全域面积统计最常用的投影没有之一。国家基础地理数据经常用它因为双标准纬线设置在25°N和47°N在国界范围内面积变形控制在很小的范围内。代码里用GeoPandas重投影albers_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumCGCS2000 unitsm no_defs basins_albers basins.to_crs(albers_crs) basins_albers[area_km2] basins_albers.geometry.area / 1e6 print(basins_albers[[NAME, area_km2]])to_crs重投影后geometry.area返回的单位是平方米除以1e6就是平方千米。这比直接使用属性表里的面积字段可靠因为投影基准、舍入方式全部由你自己控制。如果你处理的只是北京市这种较小范围也可以用UTM分区投影但代际差异不大。要注意的是重投影会改变几何精度如果到后面还要和其他地理坐标数据做空间连接记得把结果再转回经纬度。这三个文件之间的空间关系是嵌套的——同一地理位置既是某个气候区带又是某个流域的一部分。所以做叠加分析时用哪个图层做底都行但统计口径必须一致。最稳妥的做法是把流域作为基础图层气候区带作为属性挂接上去这样最后每个流域多边形都携带了气候区带信息不需要反复做空间连接。3. 用Python做流域级空间分析从读SHP到掩膜裁剪的完整流程3.1 筛选指定流域属性查询与要素抽取典型场景你想提取塔里木河流域边界用它裁剪降水栅格再统计整个流域的年降水量。第一步是属性筛选。basins_albers gpd.read_file(nine_basins.shp, encodingutf-8).to_crs(albers_crs) # 模糊匹配兼容字段值带流域后缀的情况 target 塔里木河 basin basins_albers[basins_albers[NAME].str.contains(target)].copy() print(basin[NAME].tolist())用str.contains而不是等于判断是为了兼容字段值可能是“塔里木河流域”或“塔里木河”的情况实测下来这种写法容错率最高。筛选结果可能不止一个要素——如果数据集把塔里木河流域拆成了多个子流域tolist()会把所有匹配名称显示出来这时就需要人工确认是取全流域还是只取其中一个子流域。全流域就用gpd.unary_union把多个要素合并成一个多边形。这一步最容易被忽略的是筛选前一定要copy()否则后续操作可能触发链式赋值的警告。还有一点如果你之后要把边界坐标点导出来写论文插图直接操作的是投影后的坐标单位是米不是经纬度输出前记得转回去。basin_union basin.geometry.unary_union print(Target basin area:, basin_union.area / 1e6, km2)unary_union把同属一个流域的所有多边形融合成单个面。融合后算一次面积该数字可作为后续裁剪结果的参照——如果裁剪出来的有效像元总面积和这个数字差太远说明裁剪过程中出了问题。3.2 用流域边界裁剪栅格Rasterio掩膜裁剪和参数设置拿到流域边界后下一步裁降水或气温栅格。常见的数据是NetCDF转的GeoTIFF或直接从气象共享平台下载的降水GeoTIFF。import rasterio from rasterio.mask import mask as rio_mask with rasterio.open(precipitation.tif) as src: out_image, out_transform rio_mask.mask( src, [basin_union], # 必须是GeoJSON-like的几何列表 cropTrue, nodata-9999, all_touchedFalse ) out_meta src.profile.copy() out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(tarin_basin_precip.tif, w, **out_meta) as dst: dst.write(out_image)参数的坑先从这几个说起cropTrue把输出范围裁剪成和流域外接矩形一致文件体积小很多nodata-9999保持和原栅格一致避免把无效值误当成0all_touched默认False意味着只有像元中心落在流域内才算有效。栅格分辨率低时细长的流域边界的边缘像元会被切掉这时可考虑设True但边缘会变得更粗糙[basin_union]必须放在列表里因为mask函数是按图层要素列表逐要素处理的直接传Polygon它会报错裁剪后的栅格范围是矩形流域外区域仍然是nodata。后边做统计的时候要先用有效值掩膜把流域外的nodata滤掉不然最小值算出来一定是-9999平均降水也会被拉低。3.3 分区统计每个气候带/流域的均值与总量计算裁剪只是中间步骤最终要落到分区统计。最常用的方案是rasterstats这个库它内置了zonal_stats函数对每个多边形计算栅格像元的统计量from rasterstats import zonal_stats stats zonal_stats( basins_albers, # 多边形图层 tarin_basin_precip.tif, stats[mean, sum, count], nodata-9999, geojson_outTrue ) print({f[properties][NAME]: f[properties][mean] for f in stats})zonal_stats内部会先做栅格化和重投影。多边形建议用与栅格相同的投影否则它内部会做一次矢量重投影结果可能不准。stats参数里mean是像元均值sum是全部有效像元的累加常用于计算区域降水总量count统计有多少个像元参与了计算可辅助判断结果是否合理——如果某个流域的count异常少大概率是边界与栅格没对齐。rasterstats的一个缺点是对非常大、要素非常多的矢量图层处理偏慢几千个面可能会跑几分钟。如果做全国尺度的气候带统计建议先把栅格裁成和矢量边界同一范围再跑别直接喂整张图。血泪经验我一开始图省事对全国1km分辨率降水栅格做31个省级行政区统计跑了快半个钟头后来改成先裁剪再统计几分钟搞定。如果要按气候区带分组统计各流域的平均气温场景就是空间叠加加属性分组。做法是先把气候区带和流域两个图层做空间连接再用groupby聚合overlay_result gpd.overlay(basins_albers, climates_albers, howintersection) grouped overlay_result.groupby(CLIMATE)[area_km2].sum()gpd.overlay会割开相交的面输出的是两块区域叠加后的细分多边形每个多边形仍保留两个父图层的属性字段。此时每个“流域-气候带”组合的面积就可以直接用groupby求总和。这套流程一步到位比在ArcGIS里做两轮相交、联接快得多。4. 格式转换与格式互转SHP转KML、JSON、TXT与点表生成4.1 SHP转KML/JSON几种驱动方式的参数差异用途不同选不同的转换路线丢到Google Earth或手机地图里看转KML做Web前端或传给后端接口转GeoJSON做简单属性导出或交给不装GIS软件的人转TXT表格GeoJSON最简单一行代码解决:basins.to_file(nine_basins.geojson, driverGeoJSON, encodingutf-8)转KML有两个路线。一个是GDAL内置的KML驱动通过GeoPandas的to_file指定扩展名和驱动程序即可basins_simple basins[[NAME, geometry]].copy() basins_simple.to_file(nine_basins.kml, driverKML)但GDAL的KML驱动有个麻烦它只输出有限的字段而且属性值里如果有非法XML字符比如“”“”这些符号写入会失败或生成的KML打不开。实测下来字段多了经常翻车。我一般用simplekml库生成KML字段控制和样式定制更灵活import simplekml kml simplekml.Kml() for _, row in basins.iterrows(): # 取第一个多边形环的坐标舍弃空洞 coords list(row.geometry.exterior.coords) kml.newpolygon( namerow[NAME], outerboundaryiscoords ) kml.save(nine_basins.kml)这里用exterior.coords取外环坐标把内环空洞丢掉了。对流域这种通常没有内部空洞的面要素影响不大但遇到带复杂的岛屿要素时信息会丢建议确认原始几何类型再决定。带洞的面要素会把洞填成实心看着倒无所谓但极值统计时会把周围要素覆盖。热词里常见的“arcgis shp转kml”其实就是这种需求在ArcGIS里另存为KML在Python里用上面两段代码替换即可。TXT导出更直接适合生成坐标表with open(basins_coords.txt, w, encodingutf-8) as f: for _, row in basins.iterrows(): xy ,.join([f{x:.4f},{y:.4f} for x, y in row.geometry.exterior.coords]) f.write(f{row[NAME]}|{xy}\n)注意TXT里保留的是经纬度如果你想在CAD里直接打开做成线稿得先转投影再输出坐标。转TXTSHP转TXT的好处是能在Excel里打开后排序、筛选、美化表格。坐标精度保留4位小数约合10米做概览图足够做精确放样没必要。4.2 Excel/CSV点表转SHP从经纬度列到矢量点的完整代码这是另一个高频需求手里有一份Excel表格两列是经纬度一列是站点名/测值要转成SHP做空间插值或叠加分析。GeoPandas提供了非常直接的路径import pandas as pd import geopandas as gpd from shapely.geometry import Point df pd.read_excel(stations.xlsx) # 必须包含 lon, lat 两列 geometry [Point(x, y) for x, y in zip(df[lon], df[lat])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) gdf_valid gdf.dropna(subset[lon, lat]) gdf_valid gdf_valid[(gdf_valid[lon].between(73, 135)) (gdf_valid[lat].between(18, 54))] gdf_valid.to_file(stations.shp, encodingutf-8)注意两点坐标范围过滤是我每次必做的避免个别站点经纬度填反、变成坐标跑到国外或海洋里dropna要先执行因为导入时空值的坐标行会在zip阶段直接报错。Excel里常见坑是“经纬度列名带空格”或“单位不是十进制度而是度分秒”。度分秒要先用公式换算十进制度度分/60秒/3600没换算直接进GeoDataFrame画出来的点会全部挤在坐标原点附近是个非常经典的翻车现场。Arcgis里手动做“arcgis excel点转shp”是菜单操作但数据量大或需要反复更新时会想吐。Python方式可以从读取Excel到写SHP一个脚本跑完换数据只要改路径。另外这套脚本稍加改动就能接CSVpd.read_csv即可。4.3 DWG/DXF与SHP互转坐标系漂移的处理方法DWG转SHP是另一个高频热搜词。直说结论GDAL不原生支持DWG格式别指望一行代码装就能转。常见做法分两步走在CAD里用SAVEAS命令把DWG另存为DXF选ASCII版本再用GDAL/OGR或GeoPandas读取DXF转SHPfrom osgeo import ogr ogr.UseExceptions() src_ds ogr.Open(boundary.dxf) if src_ds is None: raise RuntimeError(DXF打开失败检查文件是否被CAD占用) out_ds ogr.GetDriverByName(ESRI Shapefile).CreateDataSource(boundary.shp) src_layer src_ds.GetLayer(0) out_layer out_ds.CreateLayer(boundary, geom_typeogr.wkbPolygon) out_layer.CreateField(ogr.FieldDefn(layer, ogr.OFTString)) for feat in src_layer: geom feat.GetGeometryRef() lyr_name feat.GetField(Layer) out_feat ogr.Feature(out_layer.GetLayerDefn()) out_feat.SetGeometry(geom) out_feat.SetField(layer, lyr_name) out_layer.CreateFeature(out_feat)这段代码的最大意义是保留CAD图纸里的图层名。DXF文件里不同的图层通常代表不同的要素类型——有的是边界、有的是建筑、有的是标注——转SHP时如果不把Layer字段带出来所有要素混在一起没法用。转换后还有一步必须做坐标系核对。CAD图纸普遍不带投影信息转出来的要素如果和全国边界叠不上大概率是坐标系没对齐。可以用几个已知地理坐标点做参照反向校准或者把转换结果和这套地理分区数据用gpd.overlay做一次叠加检查跨度超过几百米就相当可疑。同理“shp文件转mapgis线文件”这类需求本质也是格式转换但是要注意SHP的线要素在MAPGIS里对应线文件面要素对应区文件属性字段也可能丢失不要指望无损往返。SHP转MAPGIS线文件的方法通常是先导出DXF再导入MAPGIS转换后检查要素个数和图层。这里不做展开但记得一个核心原则转换格式永不回头导出前先把属性该合并的合并该删的删。4.4 其他衍生场景SHP转3D Tiles、渔网分割与JSON互转的边界提醒热搜词里还有“shp转3dtiles”“渔网分割shp”“json转shp网站”这些实际都是同一类问题的变体SHP转3D Tiles常用在Cesium展示三维场景。常见做法是先把SHP要素转为带高度的GeoJSON或glTF再用CesiumLab等工具瓦片化。那个流程吃显卡和CPU对要素数量有上限要求几十万以上级别的面要素直接卡死先化简顶点再转渔网分割SHP是把研究区切成等尺寸网格再算网格内属性。QGIS的创建网格工具或者ArcGIS的渔网工具都能生成网格面然后用交集分割矢量。这里有一个常见误区渔网做成后网格面要素和原始SHP相交需要先用空间索引过滤不然几千个网格对几十万条边界的交集运算能跑上一个下午JSON转SHPGeoJSON直接gpd.read_file(xxx.geojson)就吃进来了不需要“网站”转。浏览器端在线转小文件尚可数据量大时浏览器内存直接爆掉这些衍生工具链我都不在这篇里写全核心思路不变先确认源数据的坐标系和目标格式的限制再选择转换驱动。5. 避坑实录矢量数据集使用中常见的五个翻车现场5.1 中文属性乱码读出来全是“锟斤拷”现象gpd.read_file读SHP后NAME字段显示“锟斤拷”或“锟斤拷锟斤拷”过滤条件怎么写都匹配不上。原因数据生产方用的编码是GBK或GB2312GeoPandas默认按UTF-8解码编码错位。解决读文件时显式传encodingGBK如果不确定原始编码用文本编辑器打开.dbf文件看头部信息或直接用Python循环尝试多种编码for enc in [utf-8, GBK, GB18030]: try: gdf gpd.read_file(nine_basins.shp, encodingenc) if 流域 in str(gdf[NAME].iloc[0]): print(OK:, enc) break except UnicodeDecodeError: continue顺便说一个我的习惯自己导出SHP时一律用encodingutf-8这是最通用的方案但给别人发数据前我会提醒一句对方如果有乱码就试试GBK。编码问题解决不了的时候就直接读dbf属性表单独处理用dbfread把dbf转成pandas DataFrame再和几何表合并。5.2 面积计算离谱没投影就统计的锅现象用原始SHP直接算几何面积结果显示塔里木河流域面积数千亿平方公里明显是地级市级别数据却算出比地球还大的结果。原因几何坐标是经纬度geometry.area在未投影的GeoDataFrame里返回的是度数单位的“面积”也算得出来但没有任何物理意义。解决无论数据源是什么坐标系统计面积前一律转Albers等积投影。我的检验标准是流域级面积算完和官方公报对一下三大流域中长江流域面积约180万平方公里黄河流域约75万平方公里大概超过合理范围的10%就说明投影有问题。def safe_area_km2(gdf): from copy import deepcopy tmp deepcopy(gdf) tmp tmp.to_crs(projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs) return tmp.geometry.area / 1e6用deepcopy是为了不修改原始GeoDataFrame的几何对象。这个函数在项目里反复用沉淀成了我的“后悔药函数”。5.3 流域与省界套不上叠加分析结果为空现象用流域边界剪省界数据结果输出图层是空的或者相交面积比预想小了两个数量级。原因两套SHP的坐标系不同或者边界精度版本差异很大——县级边界和流域边界在边缘有几十米的偏差小面要素相交时可能完全落空。南盘江流域边界比较曲折如果用的省界是最新版高精度数据两边误差叠加部分沿岸的窄条相交结果为0。解决先统一坐标系再叠加不要直接拿两套原始SHP做overlay。对边界接缝做一次拓扑容差缓冲basins_buffered basins.copy() basins_buffered.geometry basins_buffered.buffer(0.01) # 0.01度约1km overlay_result gpd.overlay(basins_buffered, provinces, howintersection)缓冲会扩大边界范围但注意缓冲操作也在修改几何量大的时候记得把缓冲后超出研究区的部分裁掉。还有一个做法是切换“省1”和“省2”两种版本的数据集所谓省1、省2通常是不同精度版本节点密度和属性字段组织有差异就看哪个版本和你的流域边界在坐标基准上更一致实际以你对交测试为准不要凭文件名猜。5.4 SHP自动拓扑修复shapechk与Geometry Collection重建现象读取SHP报错“Corrupt geometry”或者某些要素无法输出shapefile还打不开。原因SHP文件本质由多个子文件组成.shp、.dbf、.shx任何一个损坏或非法都会引发读取失败。拓扑错误则常在多边形自相交、重复节点、闭合环起点不一致时发生。解决先用shapechk检测SHP文件头的完整性修复后再进入GeoPandas处理。对一个存有拓扑错误的面用buffer(0)重置几何是经典处理bad_rows gdf[~gdf.geometry.is_valid].index print(Invalid rows:, len(bad_rows), bad_rows[:10].tolist()) gdf.loc[bad_rows, geometry] gdf.loc[bad_rows].geometry.buffer(0)buffer(0)的原理是对自相交的面做一次零宽度的缓冲操作让GIS引擎重建几何节点消除自相交和重复点。它能修复大部分拓扑问题但如果面本身有大量细碎的狭长裂缝最好先做一次simplify再buffer。修复后重新检查is_valid还不行就放弃而只保留大面主体的算法见很多GIS教材我遇到极端情况直接改用R语言的sf包做相同处理有时GDAL修不了的它能修好。5.5 KML写入报错或文件打不开字段限制带来的玄学问题现象SHP转KML后Google Earth打不开或者打开后要素名显示blank、属性表全空。原因GDAL的KML驱动对属性字段数量有限制超出后会被踢掉字段值里如果有非UTF-8字符或XML非法字符、、引号写入直接报错。解决转换前精简字段并做字符过滤import re def clean_xml_text(value): if value is None: return return re.sub(r[^\w\s\u4e00-\u9fff], , str(value), flagsre.UNICODE) basins_simple basins[[NAME, geometry]].copy() basins_simple[NAME] basins_simple[NAME].map(clean_xml_text)清理字段名长度、过滤非法字符后再调用to_file(..., driverKML)。还有一个细节KML驱动写入时会对字段名自动截断所以输出前把字段名改成短英文是值得的比如NM、AREA。从那以后我做任何转换第一步都是“精简字段过滤字符”这成了肌肉记忆。6. 进阶技巧按省界流域双重分组批量出图与数据验收习惯6.1 一个脚本生成多个分幅图避免手动重复劳动当天底下数据拿到手最占时间的是挨个区域出图。比如给领导汇报淮河、黄河、长江各流域的降水分布图逐个手动选范围出图很费劲。我通常的做法是按流域分组遍历生成地图import matplotlib.pyplot as plt for name, grp in basins.groupby(NAME): fig, ax plt.subplots(figsize(6, 6)) grp.plot(axax, colorlightblue, edgecolorblack) ax.set_title(name, fontsize14) fig.savefig(fbasin_{name}.png, dpi150, bbox_inchestight) plt.close(fig)如果需要按省界流域双重分组就先用空间连接把流域名称挂到省界SHP上再按拼接字段分组。这样能一次性输出诸如“安徽-淮河流域”这类文件名完整的图。缺点是用Matplotlib画的符号不好看但胜在自动化、速度快。追求美观就导出GeoJSON到QGIS里套模板批量打印。6.2 每个新版本数据到手先过这三个流程再动手分析第一核对坐标系和边界范围打印crs、total_bounds和已知参考SHP做一次叠加对比。第二检查拓扑有效性和字段类型过滤字段值里的空值、异常值is_valid.sum()一眼看出问题。第三重算面积并和官方公开数据进行对照面积差率超过5%就要检查投影或边界是否正确。那套“省1/省2”两个版本的省界线在字段组织和精度上有细微区别不检查直接混用后面做任何带省界的统计都会“差之毫厘失之千里”。这套数据同样如此等到的版本越新字段名、面积单位越可能有变化我的应对是写一个自动校验函数每次拿到压缩包先跑一遍格式检测再进正式流程。从那以后我处理任何SHP数据集都强制走一遍“读属性、查坐标系、验拓扑、算面积”这四个动作希望帮到你。本文还有配套的精品资源点击获取