简介：这份2021年中国自然保护区矢量面数据包，面向GIS从业者、生态科研人员、政府规划部门及环保组织，用于绘制保护区分布地图、开展生态评估与政策规划。压缩包共5个文件，约2.9MB，以shp、dbf、prj、shx、xml等标准GIS格式为主：shp存储保护区边界几何形状，dbf承载名称、级别、类型、面积、设立时间等属性，prj定义坐标系统，shx提供空间索引，xml记录元数据说明。数据可被ArcGIS、QGIS等软件直接读取，支持保护区管理、保护成效与生态价值评估、人类活动影响分析等研究场景。目前已有6432人学习下载，适合需要中国自然保护区空间边界数据的中高级GIS用户，帮助快速搭建分析底图、开展空间统计与制图输出。
1. 从一份 2021 年自然保护区边界数据说起
如果你做过生态评估、国土空间规划、生物多样性分析或者环境类课题，大概率遇到过这样的场景：手头有物种分布点、遥感影像、土地利用栅格，但缺一层权威的保护区边界，没法做叠加统计，也没法算保护区内的地表覆盖变化。这份「2021 年中国自然保护区矢量面数据」解决的正是这个问题——它把全国各级自然保护区的范围整理成矢量面，直接能拖进 GIS 里用。适合生态、地理、环境、规划方向的从业者和学生，也适合需要做空间叠加分析的开发者。它不是什么高深算法，但属于那种「没有就卡住、有了就顺畅」的基础底图数据，能不能用、准不准、怎么接进现有流程，才是真正要花时间的地方。
2. 这份矢量面数据到底是什么：字段、坐标系与选型判断
2.1 矢量面数据的结构长什么样
自然保护区矢量面，本质是用多边形（Polygon）或多多边形（MultiPolygon）把每个保护区的边界圈出来，每个面附带属性信息。常见的字段包括保护区名称、级别（国家级/省级/市县级）、类型（森林生态、野生动物、湿地、荒漠等）、所在省份，有的版本还会带批复面积、建立年份。这份 2021 年的数据，时间节点比较关键——2021 年前后正好是自然保护地体系整合优化的阶段，很多保护区经历了范围调整、功能区划从「核心区/缓冲区/实验区」向「核心保护区/一般控制区」过渡，所以拿到手第一件事不是急着叠加，而是先确认它的区划口径属于哪一套。
从格式上看，矢量面数据通常以 Shapefile（.shp/.shx/.dbf/.prj）或 GeoPackage（.gpkg）交付。Shapefile 兼容性最好，几乎所有 GIS 软件和 Python 库都认，但字段名有 10 字符限制，中文容易乱码；GeoPackage 是单文件、支持长字段名和 UTF-8，QGIS、ArcGIS Pro、GeoPandas 都能读。如果你后续要做自动化处理，我更推荐先把 Shapefile 转成 GeoPackage，省掉一堆编码麻烦。
2.2 坐标系与投影：别在第一步就埋雷
拿到任何矢量数据，先看.prj文件或元数据里的坐标系。国内这类数据常见两种：地理坐标系 WGS84（EPSG:4326）或 CGCS2000（EPSG:4490），单位是度；投影坐标系常见 Albers 等积投影或 Web Mercator（EPSG:3857），单位是米。这里有个血泪经验：做面积统计必须用等积投影，不能用 WGS84 直接算。在经纬度坐标下算面积，越往高纬误差越大，东北的保护区能给你算出离谱的结果。
判断方法很简单，用 GeoPandas 读进来打印 CRS：
import geopandas as gpd # 读取矢量面数据，注意中文路径和编码 gdf = gpd.read_file("nature_reserves_2021.gpkg") print("当前坐标系:", gdf.crs) print("要素数量:", len(gdf)) print("字段列表:", list(gdf.columns)) print(gdf.head(3))这段代码做了三件事：确认坐标系、看有多少个保护区面、检查字段名有没有乱码。如果gdf.crs是EPSG:4326，说明是地理坐标，后面算面积前必须投影转换。len(gdf)返回的是要素数，注意一个保护区可能是 MultiPolygon（比如被道路切开的飞地），要素数不等于多边形个数。
2.3 为什么选矢量面而不是栅格
有人会问，直接用栅格化的保护区掩膜不行吗？行，但矢量面有三个不可替代的优势：一是边界精度高，栅格化必然有锯齿和面积损失；二是属性可查询，能按级别、类型、省份筛选；三是可编辑，如果发现某个保护区边界和最新批复不符，能直接改节点。代价是叠加分析时计算量大，尤其和全国高分辨率栅格做分区统计时，矢量转栅格那一步很吃内存。常见做法是：小范围分析直接用矢量，全国尺度先按需裁剪再转栅格。
提示：如果数据来源没有明确标注坐标系，不要默认它是 WGS84。用 QGIS 叠加在线底图看一眼，边界和影像对不上就说明坐标系有问题。
3. 把数据接进分析流程：裁剪、投影与分区统计
3.1 按研究区裁剪，别一上来就全国跑
全国保护区面数据动辄几百 MB，直接和你的研究区做叠加，大部分算力浪费在无关区域。第一步永远是裁剪。假设你研究的是某个省或流域，用一个边界面去裁：
import geopandas as gpd # 读取全国保护区数据和研究区边界 reserves = gpd.read_file("nature_reserves_2021.gpkg") study_area = gpd.read_file("study_area.gpkg") # 统一坐标系后再裁剪，避免CRS不一致报错 if reserves.crs != study_area.crs: study_area = study_area.to_crs(reserves.crs) # 用空间索引加速，裁剪出研究区内的保护区 clipped = gpd.clip(reserves, study_area) print("裁剪后要素数:", len(clipped)) clipped.to_file("reserves_clipped.gpkg", driver="GPKG")gpd.clip做的是几何裁剪，落在研究区外的部分会被切掉，边界上的保护区会被切成研究区形状。这里的关键参数是 CRS 必须一致，否则 GeoPandas 会直接报错。to_file存成 GeoPackage 而不是 Shapefile，是为了保住中文字段名。如果研究区边界本身也是经纬度，裁剪没问题；但如果后面要算面积，记得先投影。
3.2 投影转换与面积计算
算保护区面积、或者做单位面积统计之前，必须转到等积投影。国内常用 Albers 投影，参数按研究区中央经线调整：
# 定义Albers等积投影，中央经线按研究区选，这里以105°E为例 albers_crs = "+proj=aea +lat_1=25 +lat_2=47 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" # 投影转换 clipped_albers = clipped.to_crs(albers_crs) # 计算每个保护区的面积（平方公里） clipped_albers["area_km2"] = clipped_albers.geometry.area / 1e6 print(clipped_albers[["name", "area_km2"]].head())lat_1和lat_2是标准纬线，一般取研究区南北边界的纬度；lon_0是中央经线。这三个参数决定了投影变形大小，选得越贴合研究区，面积越准。geometry.area返回的是投影平面上的面积，单位是平方米，除以 1e6 得到平方公里。注意：如果几何有自相交或无效环，面积会算错，跑之前用clipped_albers.is_valid.all()检查一下，无效的用buffer(0)修复。
3.3 和栅格数据做分区统计
生态分析里最常见的操作是：把保护区面当分区，统计里面某种栅格（比如 NDVI、土地覆盖、夜间灯光）的均值或面积占比。用rasterstats最省事：
from rasterstats import zonal_stats # 对每个保护区面统计NDVI均值，all_touched=True表示边界像元也算进去 stats = zonal_stats( clipped_albers, "ndvi_2021.tif", stats=["mean", "max", "min"], all_touched=True, nodata=-9999 ) # 把结果挂回GeoDataFrame clipped_albers["ndvi_mean"] = [s["mean"] for s in stats] print(clipped_albers[["name", "ndvi_mean"]].head())all_touched=True是个容易忽略的参数：默认只统计像元中心落在面内的，边界上的保护区会漏掉一圈像元；设成 True 则只要像元碰到面就算，面积小的保护区更适用。nodata要和栅格实际值一致，否则会把填充值当有效值算进去，均值直接失真。如果保护区数量多、栅格又大，这一步会很慢，常见做法是先把栅格按研究区裁小，或者用rasterio.mask批量处理。
4. 避坑与排查：这份数据最容易翻车的五个地方
4.1 中文属性乱码
现象：用 ArcGIS 或某些 Python 环境打开 Shapefile，保护区名称显示成「????」或乱码。原因：Shapefile 的.dbf默认编码是系统本地编码，跨平台读取时对不上，尤其 Windows 中文版是 GBK，Linux 和 Mac 默认 UTF-8。解决：优先用 GeoPackage 格式；如果只有 Shapefile，读取时显式指定编码gpd.read_file("xxx.shp", encoding="gbk")，或者用 QGIS 打开后另存为 UTF-8 的 GeoPackage。
4.2 坐标系缺失或标错
现象：数据能打开，但和底图、其他图层完全对不上，偏移几百米甚至几公里。原因：.prj文件丢失，或者元数据标的是 WGS84 实际却是 CGCS2000，两者在国内有几十米到上百米的差异。解决：先叠加在线影像目视检查，偏移规律如果是整体平移，多半是基准问题；用 QGIS 的「重新投影」工具试 CGCS2000 和 WGS84 哪个对得上。没有.prj时，根据数据来源判断，国内官方数据大概率是 CGCS2000。
4.3 几何无效导致面积和叠加出错
现象：面积算出负数或异常大，zonal_stats报拓扑错误，clip结果缺块。原因：多边形自相交、环方向错误、重复节点，常见于手工数字化或格式转换后的数据。解决：跑之前统一修复，gdf["geometry"] = gdf.geometry.buffer(0)能解决大部分自相交；更彻底的是用shapely.make_valid。修复后重新检查is_valid，再往下做分析。
4.4 保护区范围与最新批复不一致
现象：统计出来的保护区面积和官方公布的对不上，或者某个保护区边界明显是旧版。原因：2021 年是自然保护地整合优化的关键年份，部分保护区的范围、功能区划在这一年有调整，数据可能混用了调整前后的口径。解决：拿几个你熟悉的保护区，对照官方批复文件或最新公告核对边界和面积；如果做时间序列分析，务必确认每一年的数据口径一致，否则趋势分析全是假的。
4.5 叠加分析时要素过多导致内存溢出
现象：全国数据直接和全国栅格做zonal_stats，程序卡死或报 MemoryError。原因：矢量面节点多、栅格分辨率高，分区统计会为每个面生成掩膜，内存占用随要素数和像元数线性增长。解决：分省或分流域批量处理，处理完一个存一个；或者先把矢量转成栅格掩膜（用rasterio.features.rasterize），再用栅格代数做统计，速度快很多但会损失边界精度。
5. 进阶用法：批量处理与结果验证的固定套路
当你不是做一次分析，而是要处理多年份、多区域的数据时，单次脚本就不够用了。我一般会把整个流程封装成一个函数，输入保护区文件、栅格文件、输出路径，中间自动完成坐标系检查、裁剪、投影、分区统计、结果导出。这样换一个研究区只改参数，不用重写代码。
import geopandas as gpd import rasterio from rasterstats import zonal_stats from pathlib import Path def zonal_pipeline(reserve_path, raster_path, out_path, target_crs=None): """保护区矢量与栅格分区统计的通用流程""" gdf = gpd.read_file(reserve_path) # 1. 几何有效性检查与修复 invalid = ~gdf.is_valid if invalid.any(): print(f"修复 {invalid.sum()} 个无效几何") gdf.loc[invalid, "geometry"] = gdf.loc[invalid, "geometry"].buffer(0) # 2. 坐标系处理：没有CRS就报警，有就按需投影 if gdf.crs is None: raise ValueError("数据缺少坐标系，请先确认") if target_crs: gdf = gdf.to_crs(target_crs) # 3. 读取栅格元数据，确认nodata with rasterio.open(raster_path) as src: nodata = src.nodata # 4. 分区统计 stats = zonal_stats( gdf, raster_path, stats=["mean", "max", "min", "count"], all_touched=True, nodata=nodata ) for key in ["mean", "max", "min", "count"]: gdf[f"raster_{key}"] = [s[key] for s in stats] # 5. 导出 gdf.to_file(out_path, driver="GPKG") print(f"完成，输出 {len(gdf)} 个要素到 {out_path}") return gdf # 调用示例 zonal_pipeline( "reserves_clipped.gpkg", "ndvi_2021.tif", "reserves_ndvi_stats.gpkg", target_crs="+proj=aea +lat_1=25 +lat_2=47 +lon_0=105 +datum=WGS84 +units=m +no_defs" )这个函数把前面几章的操作串成了一条线，几个关键点值得说清楚。buffer(0)修复无效几何是通用做法，但会轻微改变边界，如果对精度要求极高，应该用shapely.make_valid保留原始节点。target_crs设成 None 时不做投影，适合只做叠加不涉及面积统计的场景。count统计的是落在面内的有效像元数，可以用来判断某个保护区是不是太小、像元太少导致均值不可靠——如果 count 只有个位数，这个均值基本没有统计意义。
验证结果我有个固定习惯：随机抽 3 到 5 个保护区，用 QGIS 单独打开，把矢量边界叠在栅格上目视检查，再手工框选几个像元算个均值，和脚本结果对一下。这一步看着笨，但能抓出坐标系偏移、nodata 设置错误、边界错位这些脚本本身不会报错的问题。从那以后我每次拿到新的矢量面数据，都强制走一遍「看 CRS → 查几何有效性 → 抽样目视核对」这三步，再开始正式分析。希望帮到你。
本文还有配套的精品资源，点击获取
