简介这份资源提供中国首批5个国家公园的边界矢量图SHP图层面向地理信息系统从业者、生态保护研究者、高校师生及规划管理人员可用于自然保护区分析、制图与空间规划。包内共44个文件以shp、shx、dbf、prj为核心配套sbx、sbn、cpg及xml元数据压缩包约1.73MB覆盖大熊猫、东北虎豹、海南热带雨林、三江源和武夷山五个国家公园其中三江源按黄河源、长江源、澜沧江园区分别成层。这些图层可直接在ArcGIS、QGIS等软件中加载与卫星影像、地形图、人口数据叠加评估栖息地状况、人类活动影响与水文变化也可用于制作公众科普分布图。已有2169人学习下载适合需要基础边界数据开展生态研究与空间分析的用户。1. 边界矢量图不是一张图国家公园 shp 图层到底解决什么问题做自然保护地、生态评估、国土空间规划的人几乎都遇到过同一个尴尬手上有遥感影像、有物种分布点、有巡护轨迹但一到“把数据裁到国家公园范围内”这一步就卡住——因为缺一份边界矢量。中国首批 5 个国家公园边界矢量图 shp 图层本质就是把这 5 个片区的管理边界落成 GIS 能直接读、能参与空间运算的矢量面数据。它不是给人看的示意图而是给机器用的“范围定义”。它的价值在于三件事一是做裁剪掩膜把全国尺度的栅格数据裁成园区尺度二是做空间统计算面积、算覆盖比例、算缓冲区三是做叠加分析和道路、村庄、水系、土地利用图层求交。适合谁做生态遥感的研究生、做保护地规划的技术员、做自然教育或科普地图的开发者以及需要给报告配一张准确范围图的从业者。这一章先把“它是什么、能干嘛”讲透后面再动手。2. 拿到 shp 先别急着打开坐标系、拓扑与字段的三重校验很多人拿到一个 shp 压缩包双击打开发现位置飘到几千公里外或者面积算出来是负的第一反应是“数据错了”。血泪经验是八成是坐标系没对上而不是数据本身有问题。这一章讲清楚在动手裁剪之前必须做的三项校验以及为什么它们决定了后面所有分析的可信度。2.1 为什么地理坐标系和投影坐标系必须分清shp 文件本身不“自带”一个绝对坐标它靠同目录下的 .prj 文件声明坐标系。国家公园边界这类数据常见两种状态一种是地理坐标系经纬度单位是度另一种是投影坐标系平面米单位是米。如果你要算面积、算缓冲区距离必须用投影坐标系如果只是做位置展示或和经纬度点位叠加地理坐标系也能用。问题在于很多人把经纬度数据直接丢进“计算几何”工具得到的结果单位是“平方度”这个数字没有任何物理意义。正确做法是先投影。中国范围内常用的是高斯-克吕格投影或 Albers 等积投影。等积投影的好处是面积不变形适合算面积高斯-克吕格适合局部、形变小。我一般会先看 .prj 里是不是有GCS字样有就是地理坐标系需要投影有PCS或Projected字样才是投影坐标系。下面这段用 Python 的 geopandas 做坐标系检查和投影是最小可复现的起点import geopandas as gpd # 读取 shp注意路径不要带中文和空格避免底层驱动报错 gdf gpd.read_file(./data/national_park_boundary.shp) # 打印当前坐标系判断是地理还是投影 print(当前 CRS:, gdf.crs) # 如果是地理坐标系单位度投影到适合中国的等积投影 # EPSG:4490 是 CGCS2000 地理坐标系EPSG:4526 是 CGCS2000 三度带投影示例 if gdf.crs and gdf.crs.is_geographic: gdf_proj gdf.to_crs(epsg4526) print(已投影单位变为米) else: gdf_proj gdf # 投影后再算面积单位才是平方米 gdf_proj[area_m2] gdf_proj.geometry.area print(gdf_proj[[name, area_m2]].head())逻辑说明read_file读取后先看crsis_geographic判断是否为地理坐标系。to_crs做投影转换参数epsg是坐标系编码不同片区可能落在不同投影带实际项目里要按经度范围选带号或者统一用 Albers 等积投影如 EPSG:102025 之类的自定义参数。面积字段必须在投影后计算否则数值不可用。2.2 拓扑检查自相交和缝隙会让面积统计翻车矢量面数据最隐蔽的坑是拓扑错误。一个多边形如果自己和自己交叉自相交或者两个相邻面之间有细缝、有重叠面积统计就会偏。更麻烦的是很多 GIS 软件在显示时看不出来只有做叠加分析时才报错。常见做法是先用is_valid检查几何有效性再用buffer(0)或make_valid修复。注意buffer(0)是老牌修复手段对多数自相交有效但对复杂几何可能产生多部件。修复后要重新检查并对比修复前后的面积差异差异过大说明数据本身有问题不能硬修。# 检查几何有效性 invalid gdf_proj[~gdf_proj.geometry.is_valid] print(无效几何数量:, len(invalid)) # 修复无效几何buffer(0) 是常用手段 gdf_proj[geometry] gdf_proj.geometry.buffer(0) # 修复后再次检查 print(修复后无效数量:, (~gdf_proj.geometry.is_valid).sum()) # 对比修复前后面积差异过大要警惕 gdf_proj[area_fixed] gdf_proj.geometry.area print(gdf_proj[[area_m2, area_fixed]].head())参数说明buffer(0)的 0 表示零距离缓冲作用是触发几何重建不是真的做缓冲。如果修复后仍有无效几何可以换shapely.validation.make_valid。面积对比时如果某个面修复前后差了几个百分点建议回到原始数据核对而不是直接采用修复结果。2.3 字段与属性name 字段为什么不能想当然shp 的属性表字段名有长度限制经典格式是 10 个字符中文名经常被截断或转码成乱码。很多边界数据里园区名称字段可能叫NAME、Name、园区名甚至被截成NAME_1。如果你写脚本时硬编码了字段名换一份数据就报 KeyError。我一般会先打印字段列表确认名称字段和面积字段是否存在再决定用哪个。如果字段缺失可以从文件名或外部对照表补。下面这段是字段探查和重命名的最小操作# 打印所有字段名和类型 print(gdf_proj.dtypes) # 假设名称字段被截断为 NAME重命名为可读的 park_name if NAME in gdf_proj.columns: gdf_proj gdf_proj.rename(columns{NAME: park_name}) # 确认每个园区的名称唯一避免后续分组统计出错 print(gdf_proj[park_name].value_counts())逻辑说明dtypes能看出字段是对象、浮点还是整型。重命名只是让后续代码可读不影响数据。value_counts用来检查是否有重名或空值如果某个园区名称出现多次可能是多部件面需要先融合dissolve再统计。3. 用 shp 做裁剪与统计从全国栅格到园区尺度的完整链路边界矢量图最核心的用途是当“裁剪掩膜”。这一章讲一条完整链路用 shp 裁剪栅格、用 shp 做分区统计、用 shp 生成缓冲区。每一步都给出可抄的代码和参数解释并说明失败时看什么。3.1 用 rasterio 按边界裁剪栅格mask 参数怎么设假设你有一份全国尺度的 NDVI 栅格想裁到某个国家公园范围内。常见做法是用 rasterio 的mask函数传入栅格和几何对象。关键参数是cropTrue它会把输出范围收紧到几何的包围盒减少空值区域all_touchedFalse表示只保留像元中心落在边界内的像元True则会保留与边界相交的所有像元。做面积统计时一般用False做边缘分析时可能用True。import rasterio from rasterio.mask import mask import geopandas as gpd # 读取边界取某一个园区 gdf gpd.read_file(./data/national_park_boundary.shp) one_park gdf[gdf[park_name] 示例园区] # 读取栅格 with rasterio.open(./data/ndvi_2023.tif) as src: # 确保边界和栅格坐标系一致不一致要先投影 if one_park.crs ! src.crs: one_park one_park.to_crs(src.crs) # 裁剪cropTrue 收紧范围all_touchedFalse 只保留中心落入的像元 out_image, out_transform mask( src, one_park.geometry, cropTrue, all_touchedFalse ) out_meta src.meta.copy() # 更新元数据写回新文件 out_meta.update({ height: out_image.shape[1], width: out_image.shape[2], transform: out_transform }) with rasterio.open(./data/ndvi_park.tif, w, **out_meta) as dest: dest.write(out_image)逻辑说明mask返回裁剪后的数组和新的仿射变换。out_meta复制原栅格元数据再更新高、宽和 transform否则写出的文件坐标会错。坐标系不一致时直接裁剪会报错或结果错位所以先to_crs对齐。失败时先看报错是“CRS mismatch”还是“geometry invalid”前者对齐坐标系后者先修几何。3.2 分区统计用 zonal_stats 算园区均值与总和裁剪只是第一步更常见的是分区统计每个园区内 NDVI 均值多少、最大值多少、像元数多少。常用工具是 rasterstats 的zonal_stats它接受矢量、栅格和统计量列表。注意nodata参数要设对否则空值会被当成 0 参与均值结果偏低。from rasterstats import zonal_stats # 统计量均值、最大值、像元计数 stats zonal_stats( one_park, ./data/ndvi_2023.tif, stats[mean, max, count], nodata-9999, # 与栅格实际 nodata 一致 geojson_outFalse ) print(stats)参数说明stats列表里count是有效像元数不是总面积。nodata必须和栅格元数据里的 nodata 一致常见是 -9999 或 0设错会让统计失真。如果结果是 None先检查矢量和栅格坐标系是否一致再检查几何是否落在栅格范围内。3.3 缓冲区分析距离参数和投影单位必须匹配做保护地分析时经常要算边界外 1 公里、5 公里缓冲带。这里最容易翻车的是单位如果数据是地理坐标系buffer(1000)会被解释成 1000 度结果完全错误。所以缓冲区必须在投影坐标系下做且距离单位是米。# 在投影坐标系下做 1 公里缓冲 gdf_proj gdf.to_crs(epsg4526) gdf_proj[buffer_1km] gdf_proj.geometry.buffer(1000) # 缓冲后面积会变大这是正常的 gdf_proj[buffer_area] gdf_proj[buffer_1km].area print(gdf_proj[[park_name, buffer_area]].head())逻辑说明buffer(1000)的 1000 是米前提是 CRS 单位是米。缓冲后几何类型仍是面可以直接参与叠加。如果要做“边界外”缓冲带需要用缓冲面减去原始面得到环状区域。4. 避坑与排查边界 shp 最常见的 5 个翻车现场这一章集中讲踩坑。每一条都按“现象 → 原因 → 解决”写都是实际项目里反复出现的。现象一打开后位置飘到海上或国外。原因.prj 缺失或坐标系声明错误软件按默认坐标系解释。解决先确认 .prj 是否存在不存在就手动指定再用已知点位如园区内一个明显地标反查偏移量判断是坐标系错还是数据本身错。现象二面积算出来是负数或极小值。原因在地理坐标系下直接算面积单位是平方度或者多边形环方向反了。解决先投影到等积投影再算面积环方向问题用make_valid或buffer(0)修复。现象三裁剪结果全是空值。原因矢量和栅格坐标系不一致或者几何范围与栅格不重叠。解决打印两者crs和bounds对比先对齐坐标系再确认栅格覆盖范围是否包含园区。现象四字段名乱码或缺失。原因shp 字段名长度限制和编码问题。解决用 geopandas 读取时指定encodingutf-8或gbk尝试字段名被截断就从外部对照表补全不要硬编码。现象五多个园区融合后面积对不上。原因相邻园区之间有重叠或缝隙直接 dissolve 会重复计算或漏算。解决先做拓扑修复再 dissolve融合后对比各园区面积之和与融合后总面积差异超过阈值就回去查。提示以上五条里坐标系问题占了一半以上。养成“先看 CRS、再算面积、后做叠加”的顺序能省掉大量返工。5. 进阶技巧把边界 shp 变成可复用的空间分析底座走到这里你已经能读、能裁、能统计。最后一章讲一个具体技巧把边界 shp 做成一个可复用的“分析底座”让后续所有分析都基于它而不是每次重新处理。核心思路是统一投影、统一字段、统一几何有效性然后导出为 GeoPackage 或 Parquet减少 shp 的格式限制。具体做法分三步。第一步把 5 个园区统一投影到同一个等积投影保证面积可比。第二步补全字段园区名称、面积、周长、中心点坐标方便后续分组和制图。第三步导出为 GeoPackage因为 shp 有字段名长度和单文件大小限制GeoPackage 没有这些问题。# 统一投影并补全字段 gdf_all gdf.to_crs(epsg4526) gdf_all[area_m2] gdf_all.geometry.area gdf_all[perimeter_m] gdf_all.geometry.length gdf_all[centroid_x] gdf_all.geometry.centroid.x gdf_all[centroid_y] gdf_all.geometry.centroid.y # 导出为 GeoPackage避免 shp 的字段和大小限制 gdf_all.to_file(./data/park_base.gpkg, layerboundary, driverGPKG)参数说明to_crs统一投影area和length在投影坐标系下才有物理意义centroid取几何中心。to_file的driverGPKG指定 GeoPackage 格式layer是图层名。导出后可以用 QGIS 或 geopandas 直接读取后续裁剪、统计都基于这个底座。验证方法导出后重新读取对比面积字段和原始计算是否一致再用一个已知点做空间查询确认点在面内的判断正确。如果要做时间序列分析可以在这个底座上挂多个年份的统计结果形成“园区-年份-指标”的长表。我自己的习惯是每拿到一份新的边界数据先花二十分钟做坐标系、拓扑、字段三项校验再开始任何分析。这二十分钟能省掉后面几小时的返工。边界矢量图看着简单但它是所有空间分析的起点起点歪了后面全歪。希望帮到你。本文还有配套的精品资源点击获取