尧图网页设计 
简介这份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-8QGIS、ArcGIS Pro、GeoPandas 都能读。如果你后续要做自动化处理我更推荐先把 Shapefile 转成 GeoPackage省掉一堆编码麻烦。2.2 坐标系与投影别在第一步就埋雷拿到任何矢量数据先看.prj文件或元数据里的坐标系。国内这类数据常见两种地理坐标系 WGS84EPSG:4326或 CGCS2000EPSG:4490单位是度投影坐标系常见 Albers 等积投影或 Web MercatorEPSG:3857单位是米。这里有个血泪经验做面积统计必须用等积投影不能用 WGS84 直接算。在经纬度坐标下算面积越往高纬误差越大东北的保护区能给你算出离谱的结果。判断方法很简单用 GeoPandas 读进来打印 CRSimport 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, driverGPKG)gpd.clip做的是几何裁剪落在研究区外的部分会被切掉边界上的保护区会被切成研究区形状。这里的关键参数是 CRS 必须一致否则 GeoPandas 会直接报错。to_file存成 GeoPackage 而不是 Shapefile是为了保住中文字段名。如果研究区边界本身也是经纬度裁剪没问题但如果后面要算面积记得先投影。3.2 投影转换与面积计算算保护区面积、或者做单位面积统计之前必须转到等积投影。国内常用 Albers 投影参数按研究区中央经线调整# 定义Albers等积投影中央经线按研究区选这里以105°E为例 albers_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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_touchedTrue表示边界像元也算进去 stats zonal_stats( clipped_albers, ndvi_2021.tif, stats[mean, max, min], all_touchedTrue, nodata-9999 ) # 把结果挂回GeoDataFrame clipped_albers[ndvi_mean] [s[mean] for s in stats] print(clipped_albers[[name, ndvi_mean]].head())all_touchedTrue是个容易忽略的参数默认只统计像元中心落在面内的边界上的保护区会漏掉一圈像元设成 True 则只要像元碰到面就算面积小的保护区更适用。nodata要和栅格实际值一致否则会把填充值当有效值算进去均值直接失真。如果保护区数量多、栅格又大这一步会很慢常见做法是先把栅格按研究区裁小或者用rasterio.mask批量处理。4. 避坑与排查这份数据最容易翻车的五个地方4.1 中文属性乱码现象用 ArcGIS 或某些 Python 环境打开 Shapefile保护区名称显示成「????」或乱码。原因Shapefile 的.dbf默认编码是系统本地编码跨平台读取时对不上尤其 Windows 中文版是 GBKLinux 和 Mac 默认 UTF-8。解决优先用 GeoPackage 格式如果只有 Shapefile读取时显式指定编码gpd.read_file(xxx.shp, encodinggbk)或者用 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_crsNone): 保护区矢量与栅格分区统计的通用流程 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_touchedTrue, nodatanodata ) for key in [mean, max, min, count]: gdf[fraster_{key}] [s[key] for s in stats] # 5. 导出 gdf.to_file(out_path, driverGPKG) print(f完成输出 {len(gdf)} 个要素到 {out_path}) return gdf # 调用示例 zonal_pipeline( reserves_clipped.gpkg, ndvi_2021.tif, reserves_ndvi_stats.gpkg, target_crsprojaea lat_125 lat_247 lon_0105 datumWGS84 unitsm no_defs )这个函数把前面几章的操作串成了一条线几个关键点值得说清楚。buffer(0)修复无效几何是通用做法但会轻微改变边界如果对精度要求极高应该用shapely.make_valid保留原始节点。target_crs设成 None 时不做投影适合只做叠加不涉及面积统计的场景。count统计的是落在面内的有效像元数可以用来判断某个保护区是不是太小、像元太少导致均值不可靠——如果 count 只有个位数这个均值基本没有统计意义。验证结果我有个固定习惯随机抽 3 到 5 个保护区用 QGIS 单独打开把矢量边界叠在栅格上目视检查再手工框选几个像元算个均值和脚本结果对一下。这一步看着笨但能抓出坐标系偏移、nodata 设置错误、边界错位这些脚本本身不会报错的问题。从那以后我每次拿到新的矢量面数据都强制走一遍「看 CRS → 查几何有效性 → 抽样目视核对」这三步再开始正式分析。希望帮到你。本文还有配套的精品资源点击获取 2021年中国自然保护区矢量面数据：GIS叠加分析与分区统计实战
简介这份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-8QGIS、ArcGIS Pro、GeoPandas 都能读。如果你后续要做自动化处理我更推荐先把 Shapefile 转成 GeoPackage省掉一堆编码麻烦。2.2 坐标系与投影别在第一步就埋雷拿到任何矢量数据先看.prj文件或元数据里的坐标系。国内这类数据常见两种地理坐标系 WGS84EPSG:4326或 CGCS2000EPSG:4490单位是度投影坐标系常见 Albers 等积投影或 Web MercatorEPSG:3857单位是米。这里有个血泪经验做面积统计必须用等积投影不能用 WGS84 直接算。在经纬度坐标下算面积越往高纬误差越大东北的保护区能给你算出离谱的结果。判断方法很简单用 GeoPandas 读进来打印 CRSimport 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, driverGPKG)gpd.clip做的是几何裁剪落在研究区外的部分会被切掉边界上的保护区会被切成研究区形状。这里的关键参数是 CRS 必须一致否则 GeoPandas 会直接报错。to_file存成 GeoPackage 而不是 Shapefile是为了保住中文字段名。如果研究区边界本身也是经纬度裁剪没问题但如果后面要算面积记得先投影。3.2 投影转换与面积计算算保护区面积、或者做单位面积统计之前必须转到等积投影。国内常用 Albers 投影参数按研究区中央经线调整# 定义Albers等积投影中央经线按研究区选这里以105°E为例 albers_crs projaea lat_125 lat_247 lat_00 lon_0105 x_00 y_00 datumWGS84 unitsm 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_touchedTrue表示边界像元也算进去 stats zonal_stats( clipped_albers, ndvi_2021.tif, stats[mean, max, min], all_touchedTrue, nodata-9999 ) # 把结果挂回GeoDataFrame clipped_albers[ndvi_mean] [s[mean] for s in stats] print(clipped_albers[[name, ndvi_mean]].head())all_touchedTrue是个容易忽略的参数默认只统计像元中心落在面内的边界上的保护区会漏掉一圈像元设成 True 则只要像元碰到面就算面积小的保护区更适用。nodata要和栅格实际值一致否则会把填充值当有效值算进去均值直接失真。如果保护区数量多、栅格又大这一步会很慢常见做法是先把栅格按研究区裁小或者用rasterio.mask批量处理。4. 避坑与排查这份数据最容易翻车的五个地方4.1 中文属性乱码现象用 ArcGIS 或某些 Python 环境打开 Shapefile保护区名称显示成「????」或乱码。原因Shapefile 的.dbf默认编码是系统本地编码跨平台读取时对不上尤其 Windows 中文版是 GBKLinux 和 Mac 默认 UTF-8。解决优先用 GeoPackage 格式如果只有 Shapefile读取时显式指定编码gpd.read_file(xxx.shp, encodinggbk)或者用 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_crsNone): 保护区矢量与栅格分区统计的通用流程 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_touchedTrue, nodatanodata ) for key in [mean, max, min, count]: gdf[fraster_{key}] [s[key] for s in stats] # 5. 导出 gdf.to_file(out_path, driverGPKG) print(f完成输出 {len(gdf)} 个要素到 {out_path}) return gdf # 调用示例 zonal_pipeline( reserves_clipped.gpkg, ndvi_2021.tif, reserves_ndvi_stats.gpkg, target_crsprojaea lat_125 lat_247 lon_0105 datumWGS84 unitsm no_defs )这个函数把前面几章的操作串成了一条线几个关键点值得说清楚。buffer(0)修复无效几何是通用做法但会轻微改变边界如果对精度要求极高应该用shapely.make_valid保留原始节点。target_crs设成 None 时不做投影适合只做叠加不涉及面积统计的场景。count统计的是落在面内的有效像元数可以用来判断某个保护区是不是太小、像元太少导致均值不可靠——如果 count 只有个位数这个均值基本没有统计意义。验证结果我有个固定习惯随机抽 3 到 5 个保护区用 QGIS 单独打开把矢量边界叠在栅格上目视检查再手工框选几个像元算个均值和脚本结果对一下。这一步看着笨但能抓出坐标系偏移、nodata 设置错误、边界错位这些脚本本身不会报错的问题。从那以后我每次拿到新的矢量面数据都强制走一遍「看 CRS → 查几何有效性 → 抽样目视核对」这三步再开始正式分析。希望帮到你。本文还有配套的精品资源点击获取