自然保护区shp空间分布数据处理：从坐标系检查到叠加分析实战

简介：这份资源为ESRI Shapefile格式的自然保护区空间分布矢量数据，面向从事GIS分析、生态保护研究、环境政策制定及地理信息教学的师生与研究人员，可用于保护区边界制图、空间统计与多源数据集成等场景。压缩包共8个文件，约533KB，包含shp几何数据、dbf属性表、shx索引、prj投影定义，以及cpg、sbn、sbx、xml等辅助文件，完整覆盖Shapefile读取与坐标解析所需组件。数据记录保护区多边形边界，属性表可能含名称、类型、面积、成立年份与管理机构等字段，便于匹配与统计。目前已有323人学习下载，适合作为GIS入门练手或环境研究的底图数据，结合Geopandas等工具可进一步计算保护区总面积、识别相邻区域并制作专题地图。

1. 拿到「自然保护区空间分布数据.zip」先别急着双击：这份 shp 到底能干什么

很多人第一次拿到「自然保护区空间分布数据.zip」这种压缩包，第一反应是双击解压，然后对着里面一堆.shp、.shx、.dbf、.prj发懵——明明只想要一张能看的图，怎么冒出来五六个同名文件。先说清楚：这份数据本质是一套矢量边界，用 shp 格式存储全国或某区域自然保护区的空间范围，每个保护区是一个多边形要素，属性表里通常带名称、级别（国家级/省级/市县级）、面积、所属省份这些字段。它能干的事很具体：做生态红线叠加、算保护区与工程项目的空间冲突、统计某流域内保护区数量、给规划图做底图。适合生态评估、国土空间规划、环境咨询、遥感地信方向的从业者，也适合要写论文需要边界数据的学生。zip 只是打包外壳，真正干活的是解压后的 shp 文件组，所以第一步不是打开 ArcGIS，而是先搞懂这套文件怎么组织、坐标系是什么、字段全不全。这一章先把「这是什么、能解决什么」讲透，后面再动手。

2. 解压之后先看什么：shp 文件组的构成与坐标系判断

2.1 shp 不是单个文件，是一组同名文件在协同工作

shp 格式从 1990 年代沿用至今，它的设计是「一个几何主体 + 若干辅助文件」，缺一个都可能打不开或者属性丢失。拿到「自然保护区空间分布数据.zip」解压后，你大概率会看到下面这些：

所以「shp 文件下载」下来如果只有一个.shp，那基本是残缺的。常见做法是：把这一组同名文件放在同一个文件夹里，不要改名、不要拆散，ArcGIS、QGIS、GeoPandas 都是按文件名前缀去自动找同组文件的。

判断坐标系是解压后的第一要务。用记事本打开.prj文件，如果看到GCS_WGS_1984或GEOGCS["GCS_WGS_1984"...]，说明是地理坐标系（经纬度，单位度）；如果看到CGCS2000_3_Degree_GK_Zone_xx或带Meter的投影坐标系，说明是投影坐标（单位米）。这一步决定了你后面算面积、做叠加时数字对不对。

2.2 用 QGIS 或 GeoPandas 快速体检数据

不装 ArcGIS 也能干活。QGIS 免费，拖进去就能看；如果要做批量处理，用 Python 的 GeoPandas 更顺手。先做一次「体检」，看要素数量、坐标系、字段、几何是否有效：

import geopandas as gpd # 读取解压后的 shp，注意路径指向 .shp 文件即可 gdf = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp") print("要素数量:", len(gdf)) print("坐标系:", gdf.crs) print("字段列表:", list(gdf.columns)) print("几何类型:", gdf.geom_type.unique()) # 检查几何有效性，无效几何会导致后续叠加分析报错 invalid = gdf[~gdf.is_valid] print("无效几何数量:", len(invalid)) # 看前几行属性，确认名称、级别字段是否存在 print(gdf.head())

逻辑说明：read_file会自动关联同组的 shx、dbf、prj；gdf.crs直接读出坐标系，比翻 prj 文件快；is_valid是排查「玄学报错」的第一道关，很多叠加失败都是因为个别多边形自相交。参数上，如果数据量大（几万个要素），read_file可以加rows=1000先抽样看结构，避免一次性吃满内存。

如果gdf.crs返回None，说明 prj 缺失或没被识别，这时候不要硬算面积，先补坐标系。常见做法是查数据说明或根据坐标数值范围判断：经纬度数值在 -180~180 之间，投影坐标数值通常是六位或七位数。

3. 把保护区边界用起来：从坐标转换到空间叠加的完整链路

3.1 坐标系不统一是叠加分析翻车的头号原因

假设你要把保护区边界和某个工程项目的红线做叠加，判断是否压占。如果保护区是 WGS84 经纬度，工程项目是 CGCS2000 投影坐标，直接叠加要么报错，要么结果错得离谱。正确做法是先统一到同一个投影坐标系，再算面积和相交。

import geopandas as gpd # 读取保护区和项目红线 pa = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp") project = gpd.read_file(r"D:\data\project\redline.shp") # 统一转到 CGCS2000 三度带投影，适合中国区域算面积 # EPSG:4547 是 CGCS2000 / 3-degree Gauss-Kruger CM 114E，按实际经度选带号 pa_proj = pa.to_crs(epsg=4547) project_proj = project.to_crs(epsg=4547) # 空间相交，找出与项目红线有重叠的保护区 overlap = gpd.overlay(pa_proj, project_proj, how="intersection") # 计算压占面积，单位平方米转公顷 overlap["压占面积_公顷"] = overlap.geometry.area / 10000 print(overlap[["名称", "级别", "压占面积_公顷"]])

逻辑说明：to_crs做坐标转换，epsg=4547只是示例，实际要按数据所在经度选对应带号，选错带号会导致东西方向偏移几百米。overlay的how="intersection"只保留相交部分，how="difference"则用来做「去掉一部分矢量」——这正是热搜里「如何在 shp 图中去掉一部分矢量」的典型操作。参数上，overlap后一定要检查结果要素数量，如果比预期多，可能是保护区之间有重叠或项目红线自相交。

提示：算面积前务必确认坐标系是投影坐标。经纬度坐标下geometry.area得到的是「平方度」，没有任何实际意义。

3.2 属性表清洗：级别字段不统一会让统计全错

自然保护区数据常见的坑是级别字段写法五花八门：「国家级」「国家」「National」「国家级保护区」混在一起。直接分组统计会得到一堆重复类别。清洗思路是先看唯一值，再映射归一。

import geopandas as gpd gdf = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp") # 查看级别字段的所有唯一值 print(gdf["级别"].value_counts(dropna=False)) # 建立映射规则，把各种写法归一到三类 level_map = { "国家级": "国家级", "国家": "国家级", "National": "国家级", "省级": "省级", "省": "省级", "市县级": "市县级", "市级": "市县级", "县级": "市县级" } gdf["级别_归一"] = gdf["级别"].map(level_map).fillna("其他") # 按归一后的级别统计数量和总面积 gdf_proj = gdf.to_crs(epsg=4547) gdf_proj["面积_公顷"] = gdf_proj.geometry.area / 10000 print(gdf_proj.groupby("级别_归一")["面积_公顷"].agg(["count", "sum"]))

逻辑说明：value_counts(dropna=False)把空值也列出来，避免漏掉缺失；map加fillna保证没匹配上的归到「其他」而不是变成 NaN；统计前转投影坐标，面积才可信。参数上，如果字段名不是「级别」，先用list(gdf.columns)确认实际字段名，中文列名在不同软件导出时可能带空格或 BOM 头。

3.3 导出成 kml 或 GeoJSON 给别人用

保护区数据经常要给非 GIS 同事看，或者嵌到网页地图里。shp 转 kml 用 QGIS 右键导出即可，命令行或脚本方式如下：

import geopandas as gpd gdf = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp") # 转成 WGS84，kml 和网页地图都要求经纬度 gdf_wgs = gdf.to_crs(epsg=4326) # 导出 GeoJSON，注意 ensure_ascii=False 保留中文 gdf_wgs.to_file(r"D:\data\output\protected_areas.geojson", driver="GeoJSON") # 导出 kml，需要 fiona 支持 KML driver gdf_wgs.to_file(r"D:\data\output\protected_areas.kml", driver="KML")

逻辑说明：kml 和 GeoJSON 都基于 WGS84 经纬度，所以必须先to_crs(epsg=4326)。driver="KML"依赖 GDAL 编译时带 KML 支持，如果报错，改用 QGIS 导出更稳。参数上，GeoJSON 文件会比 shp 大，因为文本格式冗余，几万个要素建议先简化几何（gdf.simplify(0.001)）再导出。

4. 避坑与排查：自然保护区 shp 处理中最容易翻车的 5 个点

4.1 现象：打开后一片空白，什么都看不到

原因：坐标系缺失或错误，数据被定位到地球另一端；或者几何范围极小，缩放到全图时看不见。解决：先看.prj是否存在，用gdf.crs确认；如果 crs 为 None，根据坐标数值判断是经纬度还是投影，手动set_crs后再to_crs。另外用gdf.total_bounds看范围，如果数值是几百几千，基本是投影坐标被当成了经纬度。

4.2 现象：中文属性全是乱码

原因：dbf 文件的字符编码没有声明，或者 cpg 文件缺失。解决：读取时指定编码，gpd.read_file(path, encoding="gbk")或encoding="utf-8"试；如果 cpg 缺失，手动创建一个同名.cpg文件，内容写UTF-8或GBK。导出时也要注意，GeoJSON 用ensure_ascii=False。

4.3 现象：叠加分析报拓扑错误，提示 invalid geometry

原因：个别多边形自相交、有重复点或环方向错误。解决：先gdf.is_valid定位，再用gdf.buffer(0)修复——这是最常用的「后悔药」，对绝大多数自相交有效。修复后重新检查有效性，再跑叠加。注意buffer(0)对极小的碎片可能产生空几何，修复后要过滤掉。

4.4 现象：面积算出来大得离谱或小得可怜

原因：在经纬度坐标系下直接算面积，得到的是平方度；或者投影带号选错，导致比例失真。解决：统一转到合适的投影坐标系，中国区域常用 CGCS2000 三度带或 Albers 等积投影。等积投影（如 EPSG:102025）更适合做面积统计，因为它保证面积不变形。

4.5 现象：解压后文件缺失，只有 shp 没有 dbf

原因：压缩包本身不完整，或者解压时被杀毒软件拦截了 dbf。解决：重新解压，关闭实时防护再试；如果确实只有 shp，属性无法恢复，只能重新获取完整数据。这也是为什么「shp 文件下载」后要第一时间检查文件组是否齐全，别等到分析到一半才发现。

5. 进阶：用渔网分割和批量裁剪把保护区数据接到自己的分析流程里

5.1 渔网分割 shp：把大边界切成规则网格做统计

做生态评估时经常要把保护区切成规则网格，统计每个网格内的保护区面积占比。这就是热搜里「渔网分割 shp」的典型场景。思路是先按研究区范围生成渔网，再和保护区做叠加。

import geopandas as gpd import numpy as np from shapely.geometry import box # 读取保护区并转投影 pa = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp").to_crs(epsg=4547) # 按保护区范围生成 10km x 10km 渔网 minx, miny, maxx, maxy = pa.total_bounds cell = 10000 # 10km，单位米 cols = np.arange(minx, maxx, cell) rows = np.arange(miny, maxy, cell) cells = [box(x, y, x + cell, y + cell) for x in cols for y in rows] grid = gpd.GeoDataFrame({"id": range(len(cells))}, geometry=cells, crs=pa.crs) # 叠加，计算每个网格内保护区面积 inter = gpd.overlay(grid, pa, how="intersection") inter["面积_公顷"] = inter.geometry.area / 10000 result = inter.groupby("id")["面积_公顷"].sum().reset_index() # 合并回渔网，没有保护区的网格面积为 0 grid = grid.merge(result, on="id", how="left").fillna({"面积_公顷": 0}) grid.to_file(r"D:\data\output\grid_stats.shp", encoding="utf-8")

逻辑说明：total_bounds拿到整体范围，box逐个生成方格，overlay求交后按网格 id 汇总。参数上，cell决定网格大小，10km 适合省级评估，1km 适合精细分析但要素数量会暴涨，几万个网格叠加会很慢，建议先按研究区裁剪保护区再生成渔网。fillna把没有交集的网格补 0，否则统计时会漏掉空白区域。

5.2 批量裁剪与批量导出：一次处理多个保护区

如果属性表里每个保护区要单独出一个文件，用循环按名称分组导出：

import geopandas as gpd import os gdf = gpd.read_file(r"D:\data\自然保护区空间分布\protected_areas.shp").to_crs(epsg=4326) out_dir = r"D:\data\output\by_name" os.makedirs(out_dir, exist_ok=True) for name, group in gdf.groupby("名称"): # 文件名去掉非法字符 safe = "".join(c for c in str(name) if c not in r'\/:*?"<>|') group.to_file(os.path.join(out_dir, f"{safe}.shp"), encoding="utf-8") print("已导出:", safe, "要素数:", len(group))

逻辑说明：groupby("名称")按保护区名分组，每组单独导出。参数上，文件名清洗很重要，保护区名称里常带「/」或「（）」，直接当文件名会报错。导出编码统一用 utf-8，避免中文乱码。如果名称字段有重复或空值，先清洗再分组，否则会覆盖或报错。

5.3 验证数据质量的一个习惯

我一般拿到任何 shp 数据，先跑三行：看要素数、看坐标系、看几何有效性。这三步能挡掉八成后续报错。自然保护区数据尤其要注意边界是否和行政区划对得上——有时候数据版本旧，保护区范围已经调整过，直接用于正式报告前最好和最新公告核对。这个习惯帮我省过很多返工，希望帮到你。

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