ARTICLE · INTELLIGENCE

战地情报 · 详情页

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

自然保护区空间分布shp数据处理全攻略:下载、坐标系与叠加分析

自然保护区空间分布shp数据处理全攻略:下载、坐标系与叠加分析 简介面向GIS制图与空间分析的实用数据资源提供全国国家级和地方级自然保护区空间分布矢量数据shp格式数据按保护等级分别组织适合开展保护区分布制图、区域生态格局研究、规划前期底图叠加等场景也可作为ArcGIS、QGIS等平台的教学实操素材。压缩包共16个文件完整包含shp几何数据、dbf属性表、prj投影参数、shx索引以及sbn/sbx空间索引、cpg编码参数与xml元数据等配套文件整体约7.95MB解压后可直接加载使用。已有948人学习浏览说明在生态、规划与GIS教学领域具备较高参考价值。数据内容结构规范既可用于科研论文配图与专题地图制作也可用于国家级/地方级保护区关联属性查询、面积统计等基础空间分析对GIS初级学习者亦能起到熟悉矢量数据组织方式的作用对中高级使用者则可快速获取制图底图节省数据整理时间。1. 中国自然保护区空间分布数据shp一张能救场也最容易翻车的底图“中国自然保护区空间分布数据shp”这个标题听起来像是一个现成的下载包实际上却是GIS圈子里流通最散、版本最杂的一份底图。做过环评、生态评估或科研项目的人大概都经历过明明只想要一张全国保护区边界用来核项目选址是否压占结果从下载站翻到网盘折腾一两个小时要么坐标系没标注要么字段全是拼音要么边界和官方名录对不上。这篇文章按真实工作流来讲清楚这份数据长什么样能从哪些渠道拿到拿到之后怎么体检、怎么标准化、叠加分析时怎么用以及我踩过的那些坑。适合GIS从业者、生态方向的科研人员以及规划/环评工作者直接照着操作。2. 先把shp拆开看6个配套文件、坐标系与字段表2.1 shp不是单个文件丢失prj的shp等于废数据很多人第一次下载“自然保护区shp”时解压出来发现里面没有单个.shp文件而是一堆名字相同的文件第一反应是“下错了”。其实shp格式生来就是多文件结构一个完整的矢量面数据至少要凑齐下面这张表里的几个文件扩展名作用缺失或损坏的后果.shp存储几何要素本身点/线/面坐标都在这里面整个文件打开即报错.shx几何要素的位置索引用来快速定位部分软件能硬读但速度极慢且易闪退.dbf属性表保护区名称、级别、类型、面积都在这图形能显示但要素完全没法查属性.prj坐标系描述告诉软件这套坐标是什么参考最隐蔽的坑图形能打开但坐标系显示unknown.cpg声明dbf里中文属性用的字符编码没有它中文属性表读出来就是乱码.sbn/.sbx等空间索引由ArcGIS生成缺失不影响使用下次保存会自动生成我一般判断一份共享数据“是否完整”就直接看有没有.prj和.cpg。shp本身只是个几何外壳真正决定它能落到正确位置的是.prj决定属性表能不能读的是.cpg和dbf的编码声明。下载自然保护区数据时如果压缩包里只有.shp/.shx/.dbf三个文件我会直接把它归入“需要花半小时抢救”的那一类。另外要注意shp文件的几何类型必须唯一一个面shp里不能混着点和线同一个文件里不允许出现线要素和面要素共存的情况这也是后面用叠加分析时经常遇到报错“Geometry type mismatch”的根源。2.2 坐标系与字段决定叠加结果可信度的两个变量中国自然保护区数据的坐标系来源非常杂。常见的有三种WGS84经纬度坐标很多国际平台和老ArcGIS数据用这个CGCS2000坐标近年国内新矢量化数据的标准选择西安80、北京54这类历史坐标系存量数据里仍然大量存在。如果一份shp的.prj文件声明了坐标参考软件可以自动做动态投影所以表面上叠加不会错但真正计算面积和距离时度量结果会按照“数据源本身的坐标系”来算。也就是说一个WGS84地理坐标系的面shp和一个CGCS2000投影坐标系的面shp叠加显示没问题输出面积却可能差得非常远。字段层面自然保护区shp通常会有名称、级别、类型、面积、批复时间这几类信息。但来源不同字段名差异很大。正规数据库出来的数据字段是英文/拼音组合比如NAME、PROVINCE、LEVEL、TYPE、AREA二传手加工过的数据可能直接是中文字段名比如“保护区名称”“级别”“类型”。还会遇到同一份数据里级别字段写的是“National”“Provincial”“Local”的情况这是WDPA风格的属性。做标准化之前一定先把字段词典摸清楚不然统计结果根本对不上。2.3 常见数据源对比保护区shp都是从哪些渠道流出来的跑过一遍“中国自然保护区shp怎么下载”的人都知道这个数据没有官方一键发布的入口。我整理几个实际可用的方向。数据渠道质量特点适合谁WDPA全球保护区库属性规范有IUCN分类和designation字段全球统一框架但中国区域的边界更新频率不稳定且有使用授权条款做跨国对比、全球生态研究的人国内科研数据共享平台由科研机构整理发布坐标系统和字段较规整但部分数据集需要按科研用途申请审批高校和科研院所项目数据论文附件以公开数据集论文的形式发布有明确的时间、版本、制作方法说明边界可信度中等需要可引用依据的论文和课题官方名录/公告转矢量只有名称、面积、批文时间没有现成边界需要自己按边界描述矢量化对最新审批状态要求极高的人各类GIS资源站/网盘获取最快但prj缺失、边界过时、层级缺失是常态临时出图、内部自用从最终项目验收角度我建议正式报告能用科研数据平台或WDPA内部初判和快速出图用资源站数据一旦涉及项目边界对接务必以官方最新名录和审批文件为准。顺带说一句国界和行政区划边界这类敏感数据不要从非官方渠道拿保护区边界虽然不是国界但涉及具体坐标位置来源尽量规范至少能在元数据里说清出处。3. 自然保护区shp文件怎么找数据源筛选下载后先做体检3.1 四个判断标准数据能不能用下载前先看这四点面对一份标注着“全国自然保护区shp”的下载包别急着解压先按四个标准过滤一遍。第一看有没有面状几何。自然保护区空间分布的本质是边界多边形如果解压出来是点文件那就只能知道保护区中心在哪没法做压占分析。第二看属性表字段是否超过5个。正常的分级数据至少该有保护区名称、级别国家级/省级/市县、类型、面积、所在省份只有三个字段以内的多半是精简加工版做不了分类统计。第三看是否有时间标识。自然保护区体系一直在调整范围变更、晋升、整合经常发生一套没有年份的数据你无法判断它是否还符合现状。第四看.prj是否存在且内容合理。完全没有坐标系声明的直接放弃有.prj但坐标为CGCS2000的是近年规范化的数据可以优先用。这四条是硬标准任何一条不满足都会在后面某个环节爆雷。最省时间的方法就是下载前先看元数据说明没有元数据说明的文件哪怕名字起得再正统也谨慎对待。3.2 下载后先体检用geopandas批量检查shp能不能读拿到一批shp文件后我习惯先用geopandas跑一遍批量体检不打开ArcGIS/QGIS也能快速知道每份文件的几何类型、要素数量、空几何情况。下面这个脚本是整理数据时的常驻工具。import geopandas as gpd from pathlib import Path data_dir Path(./natura_data) # 遍历目录下所有shp逐个尝试读取并打印基础信息 for shp in data_dir.glob(*.shp): try: # 先用gbk编码读多数国内矢量数据是gbk编码 gdf gpd.read_file(shp, encodinggbk) print( shp.name, 几何类型:, gdf.geometry.geom_type.unique(), 要素数:, len(gdf), 空几何:, int(gdf.geometry.is_empty.sum()) ) except Exception: try: # gbk读不了再试utf-8避免漏掉utf-8声明的文件 gdf gpd.read_file(shp, encodingutf-8) print( shp.name, (utf-8), 几何类型:, gdf.geometry.geom_type.unique(), 要素数:, len(gdf), 空几何:, int(gdf.geometry.is_empty.sum()) ) except Exception as e: print(shp.name, 读取失败原因:, e)这段代码的逻辑很简单对目录下每个.shp先按GBK编码读属性表读不了再按UTF-8重试最后把几何类型、要素条数和空几何数量打出来。参数说明里最值得注意的是encoding参数和is_empty属性。国内老数据dbf编码大量是CP936也就是GBK延伸字符集用默认UTF-8去读会出现乱码或者直接抛异常。is_empty则是判断有没有空几何正常的面文件空几何应该是0出现非0值说明原数据在编辑时产生了坏要素。读不出来不代表文件报废先看报错信息。如果报“No features found”则多半是shp文件本身已损坏报“Invalid field type”则是dbf字段类型定义异常报编码解码错误则只需要换个encoding再试。体检脚本跑完后把能读出来的文件统一放进一个工作目录剩下的事就是坐标标准化。3.3 统一坐标系把WGS84和CGCS2000混用的数据拉到同一张底图上体检通过之后下一步是把坐标系统一。全国尺度做叠加分析我一般统到WGS84或CGCS2000地理坐标系这样任何下游软件都不会认错。用geopandas做这件事很简单。import geopandas as gpd # 读取一份原始数据注意此处仍要指定正确编码 gdf gpd.read_file(./natura_data/protect_2020.shp, encodinggbk) # 先看原始坐标系 print(原始坐标系:, gdf.crs) # 统一到CGCS2000地理坐标EPSG:4490 gdf_4490 gdf.to_crs(epsg4490) gdf_4490.to_file(./natura_std/地区级保护区_4490.shp, encodingutf-8) # 单独转一份WGS84供谷歌地球/GMT等工具使用 gdf_4326 gdf.to_crs(epsg4326) gdf_4326.to_file(./natura_std/地区级保护区_4326.shp, encodingutf-8)这里有个核心逻辑to_crs做的是重投影而不是简单的“给坐标套个坐标系”。如果原始shp缺失.prj文件gdf.crs会显示None这时候直接to_crs会立刻报错因为软件无从得知原始坐标到底属于哪个参考系统。遇到这种情况只能根据数据源说明手动指定crs比如gdf gdf.set_crs(epsg4326)但这一步非常危险如果原始数据是西安80或北京54你硬标成WGS84所有要素会偏移几百米。所以缺失prj的文件我基本直接放弃不值得为它冒位置偏移的风险。写入时encoding一定要显式设置。to_file里不写encoding默认按本机区域设置写dbf中文Windows环境下会写成GBK跨平台交换时又乱码。我都习惯统一写utf-8代价是ArcMap老版本读会乱码但现在主流GIS都兼容UTF-8利大于弊。4. 把保护区shp用起来叠加统计、格式转换与跨软件输出4.1 项目压占查询ArcGIS里的叠加统计标准操作拿到标准化的保护区shp后最多使用的场景就是“我的项目地块有没有压占自然保护区”。在ArcGIS里做这个操作的常规路径是identity或intersect把保护区属性附加到项目地块上再分组统计。打开ArcMap把项目面shp和保护区面shp都加进来先用“分析工具→叠加→相交”输入要素选项目地块叠加要素选保护区输出要素里就会包含“项目地块自身字段保护区字段”。关键一步是设置“JoinAttributes”默认是ALL会把保护区全部属性带过来包括名称和级别。如果只想查哪些地块在保护区内输出后按保护区名称字段做“符号系统→唯一值”可视化结果一目了然。这里要注意一个问题叠加分析前两个图层必须统一坐标系。ArcMap界面上看它们叠得挺好是因为软件做了动态投影但算面积、算长度时很可能会按数据源坐标直接算导致结果偏得离谱。所以在做相交之前我会先把保护区shp和项目shp都批处理到同一个EPSG。用ArcToolbox里的“投影”功能批量处理避免手动逐个操作。4.2 ArcGIS shp转kml和json跨软件交付时怎么保住属性保护区数据经常要交给非GIS背景的同事最常见的要求是“帮我导成kml我放到地图软件里看”。ArcGIS里shp转kml的正式做法是“转换工具→转为KML→图层转KML”输出时需要设置一个参数就是每个要素的名称字段。图层转KML的参数设置示例 图层: 保护区面 输出文件: 保护区.kml 图层输出比例: 1 将要素限制为: (默认全部) 每个要素一个KML Feature: 勾选 要素名称字段: 保护区名称这里有个非常容易踩的细节KML格式内部强制WGS84坐标而且一个kml只能承载单一几何类型。如果你的保护区shp里有MultiPolygon和Polygon混合生成时会出现部分要素丢失。更麻烦的是属性表里的中文字段名在转KML时会被折叠到HTML描述块里下游软件搜索时搜不到。我的做法是转KML之前先把保护区名称重命名成英文或拼音字段再把它作为名称字段传入出了图再配一张字段对照表。用命令行也能完成KML转换GDAL没装全时会比ArcGIS更可控。ogr2ogr -f KML -dsco NameFieldname 保护区.kml 保护区.shp这个命令里-f KML指定输出格式-dsco NameField负责告诉转换器哪个字段当作KML里要素显示名。用这个方式能直接过滤掉中文字段名带来的显示问题。同理shp转GeoJSON也常用ogr2ogr处理ogr2ogr -f GeoJSON -lco RFC7946YES -lco WRITE_BOMYES 保护区.json 保护区.shpGeoJSON格式永远使用UTF-8编码WRITE_BOMYES是为了让Windows上部分浏览器和Excel能正确识别中文。转换后建议用文本编辑器打开确认前几百字节看看字段名是否正常。4.3 不装ArcGIS也能干活QGIS与命令行工具的组合拳ArcGIS正版授权不是每个团队都有这时候QGIS加GDAL命令行是完全够用的。QGIS把shp拖进去右键图层→另存为可以选GeoJSON、KML、GeoPackage等格式编码可以在“字符编码”下拉框里选UTF-8。整个过程图形化适合偶尔做一次数据处理的同事。如果项目用地还是dwg格式先别急着往GIS里放。AutoCAD的dwg一般是局部坐标系或无坐标系状态直接转shp后要素位置无法与保护区数据对齐。正确做法是在CAD里先把图纸在世界坐标系下另存为带坐标的文件或者确认好中央经线和投影带再用ArcGIS或QGIS的dwg导入工具转成shp。dwg转shp后多段线会变成密集折点线和面要素被打散后面做拓扑求交会慢很多所以务必在CAD端清理冗余节点后再转。还有一类需求是json转shp。网上那些json转shp的在线工具适合临时转小文件但坐标系、字段类型经常被猜错正式数据我还是习惯用ogr2ogr命令就是把输出格式换成ESRI Shapefileogr2ogr -f ESRI Shapefile 保护区.shp 保护区.json -lco ENCODINGUTF-8顺带提一个常见组合把自然保护区shp和流域边界shp叠加分析时两套数据经常来自不同平台坐标系不统一是常态。拿塔里木河流域、淮河这类流域边界数据来叠一定先各自用ogr2ogr转成同一EPSG再进叠加流程。5. 自然保护区数据避坑指南我踩过且还在踩的5个坑5.1 带属性的shp读到GIS里变成了一张白纸现象同一个面shp在QGIS里打开能显示在ArcMap里打开后却是空白缩放全图也无济于事。最气人的是属性表打开后有几百行记录。原因shp文件本身的几何部分和属性部分是分开存储的。能读到属性表说明.dbf是好的显示空白则说明.shp几何数据有问题最常见的是.shp文件缺失或损坏ArcMap读不到几何信息就直接显示空图层。另一种可能是几何类型与包含数据不符比如文件头声明了Polygon但实际存的是Polyline数据。解决先检查主文件的大小.shp文件只有几KB而.dbf有上百KB这基本就是几何数据丢了。老办法是重新下载原始压缩包。如果几何数据还在但文件结构错误可以用ESRI的Shapechk工具跑一遍修复命令行执行shapechk 文件路径.shp它会自动重建损坏的索引和文件头。跑完后重新加载一遍绝大多数白纸问题都能解决。5.2 保护区面积算出来和公告值差两成现象用属性表里的面积字段做统计得到的结果和官方公告的保护区面积差得非常远有的甚至差出一倍。原因属性表里的area字段是数据生产方在某种坐标系下算的面积可能是墨卡托投影的直接平面面积。Web墨卡托在低纬度还凑合到中高纬度面积被放大了很多拿它做全国统计必然偏大。另一个常见原因是数据源本身坐标是WGS84经纬度直接用经纬度当平面坐标算面积单位还会变成平方米但是按“度”的度量来算结果完全不能用。解决面积计算必须用等积投影坐标系。全国尺度我一般用Albers等积圆锥投影设置中央经线105°E标准纬线25°N和47°N投影参数如下import geopandas as gpd gdf gpd.read_file(./natura_std/保护区.shp, encodingutf-8) # Albers等积圆锥投影适用于全国尺度面积计算 albers projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs gdf_proj gdf.to_crs(albers) # 重新计算面积单位转为平方公里 gdf[area_km2_calc] gdf_proj.geometry.area / 1e6 print(gdf[area_km2_calc].sum())注意这个计算结果才有参考价值。对比时还要考虑官方面积一般是按批复文件里的坐标边界计算手工矢量化时边界平滑程度不同会有微差但不应差到20%以上。如果算完仍然差很多就要怀疑边界数据本身是否完整。5.3 dbf乱码中文属性表秒变拼音火星文现象打开属性表保护区名称显示成一堆拼音或完全不可读的符号有的字段名直接是“鍖哄潙”这类乱码。原因dbf文件里保存中文属性有两种编码路线GBK和UTF-8。ArcGIS老版本在中文系统里写出来的都是GBKQGIS默认按UTF-8读读GBK就乱。反过来新数据用UTF-8写ArcGIS老版本按系统ANSI去读出来的也是乱码。核心是对不上编码就乱。解决读文件时显式指定编码。QGIS图层右键“图层属性→数据源→字符编码”改成GBK或UTF-8即可。用geopandas读的时候gpd.read_file里加上encodinggbk或encodingutf-8。处理完数据统一重新导出时我用to_file写好encodingutf-8并在dbf同目录保留一个.cpg文件声明编码这样换软件也不容易踩坑。5.4 边界和官方公告对不上数据版本问题现象项目核查时发现某些保护区边界比公众认知范围大出一圈或者该保护区的名称已经在新的公告里变更数据里还是老名字。原因这是数据版本滞后造成的。自然保护区体系的调整一直在推进一批又一批的保护区改变了范围、晋升了级别或者更名。下载的shp制作时间早于最新公告就必然出现偏差。解决使用数据前在元数据里记下下载时间和数据版本。报告里引用边界图时建议在数据说明里写清“边界来源自某平台年某版本与最新名录的差异以主管部门公告为准”。如果项目涉及审批必须用官方最新的名录核对名称、级别和批复面积不能直接拿网上shp做结论。5.5 叠加分析后出现一堆破碎小多边形现象项目面shp和保护区面shp做intersect后输出要素数量暴涨出现大量几平方米甚至几乎零面积的小碎片统计结果出现荒谬值。原因两个数据源在相同边界处节点坐标不一致导致求交时产生狭长或微小多边形。这在两套不同来源的边界数据之间几乎必然发生一个用WGS84采集一个用CGCS2000矢量化即便重投影后位置看起来一致节点不重合就是问题。解决尽量不做“精确相交”而是改用Spatial Join或按中心点判断归属。如果非要求交在ArcGIS里把“分析工具→叠加→相交”属性里的“容差”调到合理值比如0.00001度对应大约一米但不要调过头。更严谨的做法是在相交前对两套数据先做一步“集成”用ArcGIS的“修复几何”和“整合”工具统一节点再跑相交。跑完后过滤掉面积小于阈值的小碎面。6. 进阶验证自己的保护区库能不能拿出去用6.1 属性自检按级别和类型分组统计标准化后的保护区库不能直接用先跑一遍属性逻辑校验。全国保护区按级别分为国家级、省级、地市级和区县级属性字段里应该有对应层级。用geopandas做分组统计几秒钟就能看出数据是否完整。import geopandas as gpd gdf gpd.read_file(./natura_std/保护区_4326.shp, encodingutf-8) # 如果字段名不是level改成实际字段名比如dengji或protected_level print(gdf.groupby(level)[name].count()) print(gdf.groupby(type)[name].count()) # 检查名称是否有重复 print(重复名称数量:, gdf[name].duplicated().sum())这段代码的逻辑是统计每个级别下的保护区数量和每个类型下的数量再检查名称重复情况。正常数据里同名保护区会有跨省重名但同名同时同位就要警惕。如果分组结果里某一级别的数量明显偏少多半是属性字段被错读了或者是分级字段本身映射不正确需要回到原始属性表确认。6.2 空间自检无效几何与面积奇异值属性正常不代表几何一定正常空间自检专门揪“能读但经不起分析”的坏要素。重点看两类问题一类是无效几何比如自相交、环方向错误另一类是面积奇异值比如面积字段显示几万平方公里但实际只能容纳个村庄。import geopandas as gpd gdf gpd.read_file(./natura_std/保护区_4326.shp, encodingutf-8) # 有效几何检查无效几何会破坏后续叠加分析 print(无效几何数量:, int(~gdf.geometry.is_valid.sum())) # 面积粗算统一到米制等积投影再算 albers projaea lat_125 lat_247 lat_00 lon_0105 datumWGS84 unitsm no_defs gdf_proj gdf.to_crs(albers) area_km2 gdf_proj.geometry.area / 1e6 # 打印面积最大的前5条记录检查是否符合常识 max_i area_km2.nlargest(5).index print(gdf.loc[max_i, [name, level]].assign(area_km2area_km2.loc[max_i]))is_valid是shapely基于GEOS库做的拓扑有效性检查自相交、线条闭合异常都会被标出来。面积奇异值则要靠Albers等积投影重算来对比。如果最大面积比官方目录里最大的保护区还大出很多说明边界数据里可能混入了跨多个保护区的多边形这时要直接回到原始数据源排查。6.3 交叉验证与名录、行政区界互相印证最后一步是拿数据库和外部权威参考叠加全国行政区划shp和官方保护区名录。把保护区图层叠在行政区界上看是否有明显跨省界的大面元素。并不是每个保护区都不能跨省交界处的保护区跨省是正常的但大面积横穿多个省份就需要怀疑边界是否被错误拼接。具体做法是在QGIS里加两层保护区shp设置半透明填充行政区shp用深色细线符号缩放观察。再抽查几个典型保护区把shp里的名称、级别和公开名录逐条对照。我一般抽三个类型国家级名录里的老牌保护区容易查到近年新晋升的保护区检验数据时效性身边熟悉的小保护区检验空间精度。三个都过基本可以认为这套数据能用于常规分析。说实话我自己在这套流程上也走过弯路前两年拿一份没有prj的保护区数据硬叠加项目汇报时被专家质疑空间位置有误场面相当尴尬。现在拿到任何一份自然保护区的shp我都不急着打开看边界而是先跑一遍体检脚本确认坐标系、编码、几何类型都正常再往下走。养成这个习惯之后因为数据质量半路翻车的次数少了很多。希望帮到你。本文还有配套的精品资源点击获取
RELATED READING

延伸阅读

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