ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

自然保护区矢量边界数据处理:从数据清洗到可视化全流程指南

自然保护区矢量边界数据处理:从数据清洗到可视化全流程指南

简介:一份面向生物多样性、生态保护评估及地理信息分析人员的中文自然保护区空间数据集,整合全国自然保护区名录与矢量边界,属性涵盖保护级别、保护类型、建立年份、保护面积、所属部门、所在城市等关键字段,可用于区域生态制图、保护区空间统计与科研教学。资源为rar压缩包,共7个文件、约5.46MB,包含全国保护区数据xlsx名录表、shp矢量边界及配套的prj投影、dbf属性、shx索引、xml元数据和kmz地标文件,各类型相互关联,便于在ArcGIS、QGIS等平台直接加载使用。目前已有1426人学习下载。读者可据此快速构建全国保护区基础底图,结合属性表进行分级分类展示、面积计算和时空分布分析,也可筛选指定省份或保护级别展开专题研究,省去自行爬取整理数据的环节。

1. 自然保护区矢量边界:GIS从业者最绕不开的底图数据

做国土空间规划、生物多样性评估或环境影响评价时,很少有人没被“中国自然保护区名录和边界数据”折磨过。这个标题听起来像是一份简单的公开资料,实际落地时却牵涉数据源筛选、坐标系转换、边界精度校验、名录与边界对不齐等一系列问题。我这些年经手过不少项目,从某跨区域生态廊道规划到某图像处理Demo附带的地理底图叠加,几乎每一次都需要把这套数据从原始零散状态整理成可用的矢量图层。

这套数据的核心价值在于:它既是国家级生态保护红线的空间表达,又是环评、矿权退出、生态修复项目最常引用的权威空间范围。但网上能找到的版本五花八门,有的来自某平台手工整理,有的出自某实验室抓取加工,字段缺漏严重,边界偏移几公里都不罕见。后面我会把整理路数、坐标系处理、字段清洗和可视化落地完整讲清楚,让新手能照做跑通,熟手也能对边界坑和精度判断有个参照。

2. 数据从哪里来:名录与矢量边界的数据源选型

2.1 名录数据的官方渠道与抓取策略

中国自然保护区名录的权威源头是主管部门定期发布的批复文件和统计年报,但公开端口往往只提供网页表格或PDF文档,直接当数据用需要大量清洗。常见做法是我会先锁定国家级自然保护区名录,以批复文件为准建立基础台账,再补充省级、市县级名录作为扩展层。

具体操作上,名录字段至少要包含:保护区名称、所在省份、级别(国家级/省级/市县级)、类型(森林生态、草原草甸、荒漠生态、内陆湿地、海洋海岸、野生动物、野生植物、地质遗迹、古生物遗迹)、主管部门、批复文号和批复时间。单纯从网页复制粘贴得到的版本经常缺字段,尤其批复文号和历史调整记录容易丢,而这恰恰是判断数据时效性的关键。

抓取方面,我习惯用Python的requests加BeautifulSoup做静态页面解析,遇到动态加载的表格就换成Selenium。请求头要伪装成正常浏览器,频率控制在每秒一两个请求以内,否则很容易被限流。历史名录调整信息比当前版本更难拿,建议优先找标注了“调整”“扩建”“晋升”字样的页面,逐条补齐变更记录。

2.2 矢量边界数据的公开来源与取舍逻辑

矢量边界数据比名录复杂得多。正规渠道下载的自然保护区边界shp文件往往只覆盖国家级或部分省级保护区,而且年份较旧,最近批复的新设或扩界区域不一定收录。公开数据集里能找到的边界数据来自多个渠道,质量参差不齐,有些干脆是某高校科研团队解译的成果,精度和现势性都要打折。

我的选型顺序是:先找官方发布的基础地理数据或生态保护红线数据,再看科研机构共享数据,最后才考虑第三方整理的合集。判断数据可用性有几个硬指标:是否带投影信息、边界是否经过拓扑检查、属性表里是否有与名录匹配的统一编码、时间戳是哪一年。

拿到手第一步就是核对边界数量和名录数量是否对得上。经常出现的状况是名录里有某个保护区,边界图层里却找不到;或者反过来,边界图层有坐标范围但名录里没这个名字。这时我会用“名录中的管理面积”对照“矢量边界计算面积”,差太多就说明数据源不匹配或者边界不完整。面积比对可以在QGIS里直接用字段计算器跑,不需要写复杂代码。

3. 数据清洗与统一编码:把名录和边界对齐的预处理流程

3.1 坐标系与投影的统一规范

各地拿到的保护区矢量数据投影五花八门,有的是WGS84经纬度,有的是CGCS2000高斯投影,还有个别老数据用北京54或西安80。统一坐标系这一步不能省,我通常先全部转成CGCS2000地理坐标系,再按项目要求做投影。如果项目跨多个省份,用Albers等面积投影最稳,面积统计不会变形太大。

转换工具用QGIS或GDAL都行,命令行写法比较适合批量处理。以GDAL为例:

ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:4490 保护区_CGCS2000.shp 原始边界.shp

ogr2ogr -t_srs EPSG:4490 -nlt PROMOTE_TO_MULTI 最终边界.shp 保护区_CGCS2000.shp

逻辑说明:第一条命令把WGS84经纬度坐标转成CGCS2000地理坐标系,EPSG:4490是它对应的地理坐标代码;第二条强制把要素类型转成MULTIPOLYGON,避免后续叠加分析时因为单多边形与多多边形混用而报错。参数里最关键的是-nlt PROMOTE_TO_MULTI,很多保护区边界被切分成多个面要素,不统一类型后面做相交、裁剪都会出问题。

转换完成后要在QGIS里用“检查几何”工具跑一遍,重点看有没有碎片多边形和自相交。碎边界一般是原始数据切割时产生的,自相交常见于手工勾绘的版本,这两个问题不解决,后续计算面积和出图都会翻车。

3.2 字段对齐与统一编码规则

名录和边界要能关联起来,必须有一套共同的编码规则。常见做法是用行政代码加保护区类别的组合编码,比如省份代码(2位)+保护区区划代码(4位)+类型代码(2位),生成一个唯一ID。手工整理的数据往往没有这个ID,需要自己拼接。

我一般先在Excel里给名录表加一列“统一编码”,然后在QGIS里用属性连接把名录信息关联到边界要素上。具体步骤是:右键边界图层→属性连接→选择名录表的编码字段做连接键。连接完成后检查一遍未匹配的记录数,未匹配的要逐个看是名称写法不同还是空间范围确实缺失。

字段清洗时注意几个细节:保护区名称可能有全称、简称、括号备注等不同写法,建议做一次模糊匹配;类型字段不同数据源用词不一致,有的写“森林生态”,有的写“森林生态系统类型”,要统一映射成标准分类;批复面积和矢量计算面积放在同一张表里方便后续对比。

3.3 面积计算与疑似错位边界的快速排查

清洗完字段后,我用QGIS的字段计算器给每个要素计算投影面积,代码是:

area($geometry) / 10000

结果单位是公顷。行政批复文件里管理面积通常以公顷为单位,可以直接对比。面积相差超过20%的要素需要重点排查,常见原因有三种:边界图层只覆盖了核心区而非全部管理区;坐标系搞错导致面积虚增或缩小;数据源用了某旧版本边界而名录是扩界后版本。

错位边界更隐蔽,边界图形看起来正常,但和遥感影像或地形对不上。排查方法是用QGIS加载在线影像底图,把边界叠加到遥感影像上检查是否明显偏离自然地理边界。保护区的边界一般会和山脊线、河流、道路等明显地理要素吻合,如果边界划到农田中间或者穿越城镇建成区,多半数据有问题。

4. 从数据到产出:R与Python绘制保护区分布图的两种路线

4.1 用R的sf和ggplot2快速出图

R里处理矢量数据最顺手的是sf包,绘图用ggplot2加geom_sf。适合快速探索式制图,不用写太多代码就能看到整体分布。典型流程是读入shp文件,关联名录表,然后按级别或类型着色。

library(sf)

library(ggplot2)

library(dplyr)

reserve_sf <- st_read("保护区_清洗后.shp", options = "ENCODING=UTF-8")

reserve_sf <- left_join(reserve_sf, reserve_info, by = "统一编码")

ggplot() +

geom_sf(data = reserve_sf, aes(fill = 级别), linewidth = 0.2) +

scale_fill_manual(values = c("国家级" = "#d73027", "省级" = "#fc8d59", "市县级" = "#fee08b")) +

theme_minimal()

这里linewidth控制边界线粗细,数值调太小省级保护区边界会看不清,调太大会盖住填充色。我一般用0.2到0.3之间。options = "ENCODING=UTF-8"是无数次翻车换来的参数,shp文件的dbf属性表编码经常是GBK或UTF-8,不指定编码中文标签会乱码。

4.2 用Python的GeoPandas做批量制图与空间分析

项目需要批量出图或做空间叠加时,我偏向用Python,GeoPandas配合Matplotlib能实现自动化流程。读数据、合并、制图、导出一步到位,改参数重新跑也方便。

import geopandas as gpd

import matplotlib.pyplot as plt

reserve_gdf = gpd.read_file("保护区_清洗后.shp", encoding="utf-8")

# 确保坐标系是投影坐标系,否则面积计算不准确

if reserve_gdf.crs.is_geographic:

reserve_gdf = reserve_gdf.to_crs("EPSG:4490")

# 按级别分组统计数量

stats = reserve_gdf.groupby("级别").size()

fig, ax = plt.subplots(1, 1, figsize=(10, 8))

reserve_gdf.plot(ax=ax, column="级别", cmap="YlOrRd", legend=True)

ctx.add_basemap(ax, source=ctx.providers.CartoDB.Positron)

plt.savefig("保护区分布图.png", dpi=300, bbox_inches="tight")

GeoPandas里最容易踩的坑是坐标系判断不够严谨。有些数据虽然声明了地理坐标系但实际上坐标值已经是投影后的米制单位,直接用to_crs转会得到离谱结果。我的经验是先看一眼坐标范围,如果经纬度应该在80到135、18到54左右,但如果出现几百万的数字,说明数据本身就是投影坐标,只需要定义投影而不是转换投影。这个判断步骤看似简单,能省掉后面一大堆无效重跑。

4.3 分类出图与专题图配色

出专题图时配色需要兼顾可读性和打印效果。我的习惯是:国家级用深红色系、省级用橙色系、市县级用浅黄色系;按保护区类型着色时用色系区分大类,比如森林生态用绿色系、内陆湿地用蓝色系、荒漠生态用棕黄色系。

图例和比例尺放在图上要避开数据密集区域,可以用matplotlib调整legend_loc参数。QGIS里更灵活,直接在图层面板拖图例位置。出图分辨率至少300dpi,否则论文或报告压缩后边界线会发虚。输出PDF格式能保证矢量清晰,插到文档里缩放到多大都不失真。

5. 数组边界细节校验:5个常见问题与处理办法

5.1 边界数据与行政区划边界不吻合

现象:保护区边界局部跨越省界或县界,但行政边界图层里对应区域不属于相邻行政区。

原因:数据源用了不同年份的行政边界底图,而保护区边界是依据更早的调查数据勾绘的。还有一种情况是原始数据本身是示意图纸数字化翻版,边界做了平滑简化。

解决:先核对边界图层的时间戳和来源说明,确认是旧版还是错版。如果是历史版本且保护区后续做过范围和功能区调整,应查找最新批复文件确认调整范围。找不到替代数据时,按相邻行政区最新的境界线裁掉越界部分,并在属性表备注里标记“边界已修正”。

5.2 名录里有但边界图层没有的保护区

现象:名录表收录了某国家级保护区,但矢量边界图层里完全没有对应要素。

原因:数据源之间收录标准不一致,有的只收录有公开发布边界的数据,新批建或扩界的保护区边界尚未更新入库。

解决:单独建一个“名录中有边界缺”清单,记录保护区名称、级别、类型和管理面积。如果项目必须要这块边界,用管理面积加批复文件附图做手头最可靠的近似边界。手勾边界是下策,但比完全缺失要好。勾完后要标记数据来源为“人工参考批复文件数字化”,不可以混入正式边界数据一起分析。

5.3 属性表中文乱码或字段名不可读

现象:QGIS打开shp文件后,保护区名称显示为乱码,或者字段名变成拼音缩写,看不懂含义。

原因:dbf属性表编码和GIS软件默认编码不一致,常见情况是UTF-8文件被按GBK读取或反过来。字段名是早期系统限定的10字节缩写,原始含义在生产时已经丢失。

解决:QGIS里在图层属性→数据源→编码方式手动切换“UTF-8”和“GBK”一般能解决乱码。字段名丢失只能靠字段值内容猜,比如value里都是数字和字符组合的一般是编码,汉字的是名称和类型。拿到新数据第一件事就是检查这两个问题,别等分析做完才发现。

5.4 同一保护区出现多个重叠面要素

现象:筛选某个保护区名称,返回多个要素,而且这些要素空间范围有重叠。

原因:原始数据在多次更新编辑过程中没有做好拓扑维护,新增边界直接叠加在旧边界上;也有的是同一保护区按功能区拆成了多个部分,但属性表里没有功能区字段。

解决:用QGIS的“按位置合并”功能先确认是否空间重叠。如果有功能区字段,保留字段并按功能区属性区分;如果没有且确实是重复边界,保留面积较大的要素,删除其他。这个操作执行前要备份原始图层,否则误删了某块核心区没后悔药。

5.5 边界数据在ArcGIS里正常但在QGIS里显示异常

现象:同一个shp文件,ArcMap打开显示正常,QGIS打开后要素跑到大洋中间或者图形完全变形。

原因:shp文件缺少投影信息文件(.prj),ArcGIS会按默认坐标系猜测,QGIS也猜但猜测逻辑不同。另一个常见原因是数据含多部分几何,两个软件对multipolygon的解释有差异。

解决:先用文本编辑器打开原始数据的.prj文件,确认其中WKT字符串描述的是哪个坐标系。缺失的话尝试从数据来源说明或数据值范围反推。实在找不到原坐标系,按坐标值量级做一个合理推断:6位数坐标大概率是高斯投影,7位数带带号可能是CGCS2000分带投影。定错坐标系比不定义还麻烦,务必谨慎。

6. 进阶用法:从静态数据到带时间维度的保护地变化分析

做项目到一定程度,单期边界数据的价值就不够看了。领导或客户常会问“这个保护区近几年范围有没有变化”“扩界面积有多大”“新增区域占了什么地类”。这需要把历史各版本的名录和边界数据整理成时间序列。我在处理跨年份数据时,会维护一个结构统一的数据库,核心有三张表:名录表、边界版本表、变更记录表。名录表是主表,存保护区属性和唯一编码;边界版本表每条记录存一个保护区的某一版本边界及对应年份;变更记录表存批次、扩界、降级等类型及前后面积对比。

做变化分析时先把不同年份的边界做空间叠加,用GeoPandas跑intersection和difference就能得到新增、退出、保持不变三类区域:

gdf_2023 = gpd.read_file("保护区_2023.shp")

gdf_2015 = gpd.read_file("保护区_2015.shp")

新增区域 = gpd.overlay(gdf_2023, gdf_2015, how="difference")

新增区域.to_file("扩界区域.shp", encoding="utf-8")

叠加结果面积占比能说明扩界强度,再把新增区域和土地利用数据做intersection就能看出新增范围主要是什么地类。这块分析做扎实,报告里一张图可以顶三页文字。

变化分析最常遇到的坑是不同年份的坐标系不一致,叠加前必须统一到同一投影坐标系。另一坑是同一保护区改名了,比如“A自然保护区”改成“B国家级自然保护区”,空间范围没有变化,但按名称关联会漏掉。我会在变更记录表里建立名称变更链,关联键始终用统一编码而非名称字段。

我的习惯性动作是用SQLite存这个库,名录表和变更表用SQL管理数据完整性,边界文件单独存文件夹,在SQLite里只存文件路径和版本号。这样做的好处是空间数据大文件不拖慢数据库备份,而名录表里保存的是轻量信息,查询速度更快。最近也在尝试把这个流程包装成可复用脚本,输入一个年份就能自动输出该年份的全量保护地名录及边界图层,省得每次手工操作。

说回最根本的经验:所有边界数据拿回来,第一件事就是统一坐标系、统一编码、备份原始文件。这三步做完,后面写代码、出图、做分析基本顺畅。希望这篇文章里提到的处理流程和数据坑能帮你少走几趟弯路,在下次接到保护区数据需求时,直接照这套流程做就能落地。

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

RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进