企业官网建设流程全解析

简介：这份资源是面向地理空间分析学习者的全国地理分区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", encoding="utf-8") basins = gpd.read_file("nine_basins.shp", encoding="utf-8") climates = gpd.read_file("climate_zones.shp", encoding="utf-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编码，直接读会乱码，要改成encoding="GBK"。

字段里通常还有面积属性，但要注意这个面积是数据生产方算的，单位可能是平方千米，也可能是不规则单位。正规数据集会在元数据文档里说明，没文档的话宁可自己用投影坐标重算一遍，也别直接用原始面积字段做统计。九大流域之间用名称字段区分，三大分区和气候区带靠各自的分类字段，这三个文件互相独立又可以用空间位置叠加。

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）或CGCS2000（EPSG:4490）。你要是只做查询和出图，用原始坐标系就行；做面积统计或裁剪，就必须转投影坐标系。切记不要拿经纬度去做面积计算，1度经度对应地面的距离随纬度变化，在黑龙江和海南算出来的结果差好几倍，这种错误在论文里时有发生。

2.3 坐标系统的选择：经纬度存储、投影计算的取舍

我一般习惯分两套坐标体系用：

• 存储和展示用原始地理坐标，也就是EPSG:4326或EPSG:4490
• 面积计算和栅格裁剪用Albers等积圆锥投影

Albers是处理中国全域面积统计最常用的投影，没有之一。国家基础地理数据经常用它，因为双标准纬线设置在25°N和47°N，在国界范围内面积变形控制在很小的范围内。代码里用GeoPandas重投影：

albers_crs = "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=CGCS2000 +units=m +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", encoding="utf-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的几何列表 crop=True, nodata=-9999, all_touched=False ) 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)

参数的坑先从这几个说起：

• crop=True把输出范围裁剪成和流域外接矩形一致，文件体积小很多
• nodata=-9999保持和原栅格一致，避免把无效值误当成0
• all_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_out=True ) 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, how="intersection") 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", driver="GeoJSON", encoding="utf-8")

转KML有两个路线。一个是GDAL内置的KML驱动，通过GeoPandas的to_file指定扩展名和驱动程序即可：

basins_simple = basins[["NAME", "geometry"]].copy() basins_simple.to_file("nine_basins.kml", driver="KML")

但GDAL的KML驱动有个麻烦：它只输出有限的字段，而且属性值里如果有非法XML字符，比如“&”“<”这些符号，写入会失败或生成的KML打不开。实测下来字段多了经常翻车。我一般用simplekml库生成KML，字段控制和样式定制更灵活：

import simplekml kml = simplekml.Kml() for _, row in basins.iterrows(): # 取第一个多边形环的坐标，舍弃空洞 coords = list(row.geometry.exterior.coords) kml.newpolygon( name=row["NAME"], outerboundaryis=coords ) kml.save("nine_basins.kml")

这里用exterior.coords取外环坐标，把内环空洞丢掉了。对流域这种通常没有内部空洞的面要素影响不大，但遇到带复杂的岛屿要素时信息会丢，建议确认原始几何类型再决定。带洞的面要素会把洞填成实心，看着倒无所谓，但极值统计时会把周围要素覆盖。热词里常见的“arcgis shp转kml”其实就是这种需求，在ArcGIS里另存为KML，在Python里用上面两段代码替换即可。

TXT导出更直接，适合生成坐标表：

with open("basins_coords.txt", "w", encoding="utf-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, geometry=geometry, crs="EPSG: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", encoding="utf-8")

注意两点：坐标范围过滤是我每次必做的，避免个别站点经纬度填反、变成坐标跑到国外或海洋里；dropna要先执行，因为导入时空值的坐标行会在zip阶段直接报错。

Excel里常见坑是“经纬度列名带空格”或“单位不是十进制度而是度分秒”。度分秒要先用公式换算：十进制度=度+分/60+秒/3600，没换算直接进GeoDataFrame，画出来的点会全部挤在坐标原点附近，是个非常经典的翻车现场。

Arcgis里手动做“arcgis excel点转shp”是菜单操作，但数据量大或需要反复更新时会想吐。Python方式可以从读取Excel到写SHP一个脚本跑完，换数据只要改路径。另外这套脚本稍加改动就能接CSV，pd.read_csv即可。

4.3 DWG/DXF与SHP互转：坐标系漂移的处理方法

DWG转SHP是另一个高频热搜词。直说结论：GDAL不原生支持DWG格式，别指望一行代码装就能转。常见做法分两步走：

1. 在CAD里用SAVEAS命令把DWG另存为DXF，选ASCII版本
2. 再用GDAL/OGR或GeoPandas读取DXF转SHP

from 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_type=ogr.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转SHP，GeoJSON直接gpd.read_file("xxx.geojson")就吃进来了，不需要“网站”转。浏览器端在线转小文件尚可，数据量大时浏览器内存直接爆掉

这些衍生工具链我都不在这篇里写全，核心思路不变：先确认源数据的坐标系和目标格式的限制，再选择转换驱动。

5. 避坑实录：矢量数据集使用中常见的五个翻车现场

5.1 中文属性乱码：读出来全是“锟斤拷”

现象：gpd.read_file读SHP后，NAME字段显示“锟斤拷”或“锟斤拷锟斤拷”，过滤条件怎么写都匹配不上。

原因：数据生产方用的编码是GBK或GB2312，GeoPandas默认按UTF-8解码，编码错位。

解决：读文件时显式传encoding="GBK"，如果不确定原始编码，用文本编辑器打开.dbf文件看头部信息或直接用Python循环尝试多种编码：

for enc in ["utf-8", "GBK", "GB18030"]: try: gdf = gpd.read_file("nine_basins.shp", encoding=enc) if "流域" in str(gdf["NAME"].iloc[0]): print("OK:", enc) break except UnicodeDecodeError: continue

顺便说一个我的习惯：自己导出SHP时一律用encoding="utf-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("+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +datum=WGS84 +units=m +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, how="intersection")

缓冲会扩大边界范围，但注意缓冲操作也在修改几何，量大的时候记得把缓冲后超出研究区的部分裁掉。还有一个做法是切换“省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), flags=re.UNICODE) basins_simple = basins[["NAME", "geometry"]].copy() basins_simple["NAME"] = basins_simple["NAME"].map(clean_xml_text)

清理字段名长度、过滤非法字符后，再调用to_file(..., driver="KML")。还有一个细节：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(ax=ax, color="lightblue", edgecolor="black") ax.set_title(name, fontsize=14) fig.savefig(f"basin_{name}.png", dpi=150, bbox_inches="tight") plt.close(fig)

如果需要按省界+流域双重分组，就先用空间连接把流域名称挂到省界SHP上，再按拼接字段分组。这样能一次性输出诸如“安徽-淮河流域”这类文件名完整的图。缺点是用Matplotlib画的符号不好看，但胜在自动化、速度快。追求美观就导出GeoJSON到QGIS里套模板批量打印。

6.2 每个新版本数据到手，先过这三个流程再动手分析

第一，核对坐标系和边界范围，打印crs、total_bounds，和已知参考SHP做一次叠加对比。第二，检查拓扑有效性和字段类型，过滤字段值里的空值、异常值；is_valid.sum()一眼看出问题。第三，重算面积并和官方公开数据进行对照，面积差率超过5%就要检查投影或边界是否正确。

那套“省1/省2”两个版本的省界线在字段组织和精度上有细微区别，不检查直接混用，后面做任何带省界的统计都会“差之毫厘失之千里”。这套数据同样如此，等到的版本越新，字段名、面积单位越可能有变化，我的应对是写一个自动校验函数，每次拿到压缩包先跑一遍格式检测，再进正式流程。从那以后，我处理任何SHP数据集都强制走一遍“读属性、查坐标系、验拓扑、算面积”这四个动作，希望帮到你。

本文还有配套的精品资源，点击获取