简介这份中国自然保护区空间分布数据以shp矢量格式整理面向从事GIS制图、生态研究与国土空间规划的人员可用于自然保护区分布格局分析、专题地图制作及空间叠加统计等场景。资源包共16个文件约7.95MB涵盖国家级与地方级两类保护区图层每个图层均配套shp主文件及dbf属性表、prj投影信息、shx索引、cpg编码说明、sbn与sbx空间索引、xml元数据等完整组件可直接在ArcGIS、QGIS等平台加载无需额外转换即可参与制图与空间分析。目前已有948人学习下载说明该数据在生态与地理信息领域具有较高的实用参考价值。借助这套数据读者能够快速获取全国自然保护区的位置与范围信息用于生态保护红线比对、区域生态格局研究或教学演示省去自行收集与数字化整理的时间成本适合作为GIS课程作业、科研论文配图及规划项目的基础底图数据。1. 中国自然保护区空间分布数据shp从一份边界文件到可复现的空间分析手里拿到一份“中国自然保护区空间分布数据shp”第一反应往往不是画图而是先确认它到底能不能用。这份数据通常以面要素记录各级自然保护区的边界范围属性表里带名称、级别、类型、所属省份等字段是做生态评估、国土空间规划、保护空缺分析的基础底图。它解决的核心问题是把“保护区在哪、有多大、什么级别”变成可计算、可叠加、可统计的矢量对象。适合生态遥感、GIS 分析、自然资源调查方向的从业者也适合需要把保护区边界和自身业务图层做叠加的人。但 shp 是上世纪九十年代的老格式字段名限 10 个字符、单文件 2GB 上限、编码容易乱拿到手先别急着入库后面几章会把选型、清洗、叠加和踩坑一次讲透。2. 先搞懂 shp 的脾气为什么保护区边界总在属性表上翻车2.1 shp 不是单个文件而是一组同名文件很多人从压缩包里解压出一个.shp就以为拿到了全部数据结果在软件里打开只有图形没有属性或者坐标直接飞到南极。shp 是 Esri 定义的矢量数据交换格式一个完整数据集至少包含三个必需文件.shp存几何、.shx存索引、.dbf存属性表。缺了.shx多数软件还能靠重建索引打开缺了.dbf属性就全丢了。实际分发时还常带.prj投影定义、.cpg字符编码声明、.sbn/.sbx空间索引。保护区数据因为要带名称、级别、面积等字段.dbf尤其关键。我一般拿到数据先做一次完整性检查用命令行比在图形界面里点更可靠# 列出所有同名文件确认 shp/shx/dbf/prj/cpg 是否齐全 ls -l 自然保护区.shp 自然保护区.shx 自然保护区.dbf 自然保护区.prj 自然保护区.cpg 2/dev/null # 用 ogrinfo 看图层概览不加载图形也能读元数据 ogrinfo -so -al 自然保护区.shpogrinfo -so只输出摘要-al表示所有图层。重点看三处Feature Count是不是和预期一致Extent的经纬度范围是否落在中国境内大约东经 73° 到 135°、北纬 3° 到 54°Layer SRS是否给出了投影。如果Extent出现几十万、几百万这种数值说明它是投影坐标而不是经纬度别当成错误。提示.cpg文件里通常只写一行编码名比如UTF-8或GBK。没有它时中文属性极易变乱码这是保护区数据最高频的问题。2.2 字段名 10 字符限制决定了你拿到的属性有多“残”shp 的.dbf沿用 dBASE 规范字段名最多 10 个字符且不支持中文长名。所以原始数据里“保护区名称”可能被压成NAME“级别”变成LEVEL“批准年份”变成YEAR。更麻烦的是有些分发者为了塞进更多信息把多个含义拼在一个字段里比如TYPE里同时写“森林生态系统/国家级”。这不是数据错是格式逼的。处理办法是先把字段映射关系固定下来再决定要不要转成 GeoPackage 或 PostGIS。下面这段用 Python 读取字段并打印前几行确认每个字段的真实含义import geopandas as gpd # 读取 shpencoding 按 .cpg 或实际乱码情况调整 gdf gpd.read_file(自然保护区.shp, encodingUTF-8) # 打印字段名、类型、非空数量 print(gdf.dtypes) print(gdf.notna().sum()) # 看前 3 行的关键属性判断字段语义 print(gdf[[NAME, LEVEL, TYPE, PROVINCE]].head(3))gpd.read_file会自动读取.prj作为 CRS。如果中文乱码把encoding换成GBK或GB18030再试。dtypes里geometry是几何列其余是属性列。notna().sum()能快速暴露哪些字段大面积缺失——保护区数据里“面积”“批准年份”缺失很常见别默认它们完整。2.3 坐标系不统一是叠加分析前必须解决的问题中国境内的空间数据常见三种坐标WGS84 经纬度EPSG:4326、CGCS2000EPSG:4490、以及各类高斯-克吕格投影带如 EPSG:4547 等。保护区边界如果来自不同年份、不同部门很可能一部分是经纬度、一部分是投影坐标。直接叠加会错位面积统计也会失真。判断和统一的做法是先看.prj再用代码检查并转换import geopandas as gpd gdf gpd.read_file(自然保护区.shp) print(原始 CRS:, gdf.crs) # 统一到 CGCS2000 地理坐标便于全国范围叠加 if gdf.crs is None: gdf gdf.set_crs(EPSG:4490) else: gdf gdf.to_crs(EPSG:4490) # 面积计算前投影到等面积坐标系避免经纬度算面积失真 gdf_equal_area gdf.to_crs(EPSG:6933) # 全球等面积圆柱投影 gdf[area_km2] gdf_equal_area.geometry.area / 1e6 print(gdf[[NAME, area_km2]].head())to_crs负责坐标转换set_crs只声明不转换两者别混。算面积一定要先转到等面积投影直接在 EPSG:4490 下算出来的是“平方度”没有物理意义。EPSG:6933是全球等面积投影全国尺度够用如果只做省级精细统计可以换对应省份的 Albers 投影。3. 从原始 shp 到可分析数据清洗、裁剪与格式转换3.1 几何有效性检查自相交和空几何必须先修保护区边界多来自扫描矢量化或部门上报常见自相交、重复节点、空几何。这些在画图时看不出来一做叠加或面积统计就报错。GeoPandas 提供is_valid和is_empty两个检查import geopandas as gpd from shapely.validation import make_valid gdf gpd.read_file(自然保护区.shp) # 统计无效几何和空几何 invalid gdf[~gdf.geometry.is_valid] empty gdf[gdf.geometry.is_empty] print(无效几何数:, len(invalid), 空几何数:, len(empty)) # 用 make_valid 修复无效几何 gdf[geometry] gdf.geometry.apply( lambda geom: make_valid(geom) if geom is not None and not geom.is_valid else geom ) # 修复后再次确认 print(修复后无效几何数:, len(gdf[~gdf.geometry.is_valid]))make_valid会把自相交面拆成多部件或修正边界修复后几何类型可能从Polygon变成MultiPolygon这是正常的。空几何直接删除或单独导出核查不要留在主表里。修复前建议先备份原始文件几何修复不可逆。3.2 按研究区裁剪用掩膜提取目标省份或流域全国保护区数据动辄几千个面做省级分析时先裁剪能大幅提速。常见做法是用研究区边界做clip注意裁剪会切断跨省保护区统计时要么按裁剪后面积算要么保留完整面只做空间筛选。import geopandas as gpd protected gpd.read_file(自然保护区.shp).to_crs(EPSG:4490) study_area gpd.read_file(四川省边界.shp).to_crs(EPSG:4490) # 方式一空间筛选保留与四川相交的完整保护区 selected protected[protected.intersects(study_area.unary_union)] # 方式二裁剪只保留落在四川内的部分 clipped gpd.clip(protected, study_area) print(相交完整面:, len(selected), 裁剪后面:, len(clipped)) clipped.to_file(四川保护区_裁剪.shp, encodingUTF-8)intersects判断是否相交保留完整几何适合统计“涉及四川的保护区总数”。gpd.clip做几何裁剪适合算“四川境内保护区面积”。两者结果不同按业务选。导出时显式写encodingUTF-8并确认同时生成了.cpg。3.3 转成 GeoPackage绕开 shp 的字段和编码限制如果后续要做多表关联、长字段名、存中文备注shp 会一直拖后腿。GeoPackage.gpkg是 OGC 标准单文件、支持长字段名和 UTF-8、没有 2GB 限制是 shp 的现代替代。转换一行代码import geopandas as gpd gdf gpd.read_file(自然保护区.shp, encodingGBK) gdf gdf.to_crs(EPSG:4490) # 写入 GeoPackage图层名 protected_areas gdf.to_file(protected_areas.gpkg, layerprotected_areas, driverGPKG)layer参数决定表名一个 gpkg 可以放多个图层。转换后字段名不再被截断中文属性也稳定。如果团队还在用 ArcGISgpkg 从 10.2 起就支持用 QGIS 更没问题。转完之后建议用ogrinfo再验一次要素数和范围。3.4 导出为其他格式按下游工具选不盲目转热搜里常出现 shp 转 txt、shp 转 wkt、shp 转 3dtiles这些需求背后是不同下游。转 WKT 适合入库或做文本比对转 txt 适合给不支持矢量的程序读坐标转 3dtiles 是三维可视化。别为了转而在中间格式上反复折腾直接从 gpkg 或 shp 出发。import geopandas as gpd gdf gpd.read_file(protected_areas.gpkg, layerprotected_areas) # 导出 WKT 文本保留名称和几何 with open(protected_wkt.txt, w, encodingUTF-8) as f: for _, row in gdf.iterrows(): f.write(f{row[NAME]}\t{row.geometry.wkt}\n) # 导出质心经纬度 txt便于轻量程序读取 centroids gdf.geometry.centroid with open(protected_centroid.txt, w, encodingUTF-8) as f: for name, pt in zip(gdf[NAME], centroids): f.write(f{name}\t{pt.x}\t{pt.y}\n)geometry.wkt输出标准 WKT 字符串注意大面要素的 WKT 会很长。质心导出前要确认几何有效无效几何的质心可能落在面外。转 3dtiles 属于三维管线通常需要先拉伸高度或贴合地形不在本文展开但源数据几何有效、坐标系正确是共同前提。4. 避坑与排查保护区 shp 处理中最容易翻车的 5 个点4.1 中文属性全是问号或乱码现象打开属性表保护区名称显示为????或锟斤拷。原因.dbf没有.cpg声明或声明为UTF-8但实际是GBK。解决先用ogrinfo看编码再用 Python 指定encodingGBK或GB18030读取确认无误后导出为 gpkg编码问题一次性解决。4.2 面积统计结果大得离谱或小得离谱现象一个国家级保护区算出几万平方公里。原因在经纬度坐标系下直接算面积单位是平方度。解决先to_crs到等面积投影如 EPSG:6933 或省级 Albers再算geometry.area并除以 1e6 转平方公里。4.3 叠加分析报拓扑异常现象做相交或联合时提示TopologyException。原因面要素自相交或边界重叠。解决先跑is_valid检查用make_valid修复修复后仍异常的要素单独导出人工核查。不要用buffer(0)一把梭它可能改变几何形状。4.4 字段名被截断导致关联失败现象用NAME关联外部表时对不上或字段名变成NAME_1。原因shp 字段名 10 字符限制导出时自动改名。解决转成 gpkg 后再做关联或在导出前把字段名改成短且唯一的英文名并保留一份字段映射表。4.5 裁剪后要素数量对不上现象裁剪后要素数比预期少很多。原因clip会丢弃完全落在研究区外的要素且对跨边界要素只保留内部部分。解决先明确业务要“完整面”还是“裁剪面”前者用intersects筛选后者用clip并在日志里记录两种结果的数量差异。5. 进阶用渔网分割做保护区覆盖度统计与结果验证5.1 渔网分割 shp 的生成与叠加做保护区覆盖度、保护空缺分析时常需要把研究区切成规则网格再统计。这就是热搜里“渔网分割 shp”的典型场景。思路是先生成渔网再和保护区做叠加算出每个网格内保护区面积占比。import geopandas as gpd import numpy as np from shapely.geometry import box # 读取研究区和保护区 study gpd.read_file(四川省边界.shp).to_crs(EPSG:4490) protected gpd.read_file(protected_areas.gpkg, layerprotected_areas).to_crs(EPSG:4490) # 生成 0.1 度渔网 minx, miny, maxx, maxy study.total_bounds step 0.1 cells [] x minx while x maxx: y miny while y maxy: cells.append(box(x, y, x step, y step)) y step x step grid gpd.GeoDataFrame(geometrycells, crsEPSG:4490) grid gpd.clip(grid, study) # 只保留研究区内网格 # 叠加计算每个网格内保护区面积占比 grid[cell_area] grid.to_crs(EPSG:6933).geometry.area / 1e6 inter gpd.overlay(grid, protected, howintersection) inter[inter_area] inter.to_crs(EPSG:6933).geometry.area / 1e6 cover inter.groupby(inter.index)[inter_area].sum() grid[protected_ratio] (cover / grid[cell_area]).fillna(0) grid.to_file(四川保护区覆盖度_渔网.gpkg, layergrid_cover, driverGPKG)step是网格边长单位跟 CRS 一致这里是度。0.1 度约 11 公里全国尺度常用 0.05 到 0.5 度。overlay做相交groupby按网格汇总。protected_ratio为 0 表示该网格无保护区覆盖可用于识别保护空缺。注意grid.index在overlay后可能变化实际使用时建议先给网格加唯一 ID 再叠加。5.2 结果验证三个必须做的检查第一抽查几个已知保护区看覆盖度是否合理。比如某个国家级保护区核心区所在网格protected_ratio应接近 1。第二统计总面积把所有网格的保护区面积加总和直接算的保护区总面积对比差异应在裁剪和投影误差范围内。第三检查边界网格clip后边缘网格可能不完整占比会偏低做趋势分析时要标记或剔除。# 验证一总面积对比 total_direct protected.to_crs(EPSG:6933).geometry.area.sum() / 1e6 total_grid inter[inter_area].sum() print(f直接统计: {total_direct:.2f} km2, 网格汇总: {total_grid:.2f} km2) # 验证二覆盖率分布 print(grid[protected_ratio].describe()) # 验证三标记边缘不完整网格 grid[is_edge] grid.geometry.boundary.intersects(study.geometry.boundary.unary_union) print(边缘网格数:, grid[is_edge].sum())总面积差异一般控制在 1% 以内超过就要查投影或几何有效性。describe()看覆盖率的中位数和最大值如果最大值远大于 1说明面积计算或叠加有误。边缘网格单独标记后续出图或统计时说明。5.3 我踩过的坑和现在的习惯最早做保护区覆盖度时我直接在经纬度下建渔网、算面积结果高纬度网格被拉长覆盖率整体偏高图出来自己都觉得别扭。后来固定成“渔网用地理坐标生成、面积一律转等面积投影算”再没翻过车。另一个血泪经验是每次导出 shp 前先转 gpkg字段名和编码的后悔药就在这一步。现在拿到任何保护区 shp我的顺序是ogrinfo看元数据、Python 查几何有效性、转 gpkg、再进分析四步走完才动手做业务。这套流程不新鲜但能省下大量返工时间。希望帮到你。本文还有配套的精品资源点击获取