从数据到洞察：如何用Python分析这份2023自然保护区数据，发现生态保护热点？

从数据到洞察：Python空间分析揭示2023自然保护区生态价值图谱

自然保护区的边界数据从来不只是地图上的几条线——它们是生态系统的脉搏、物种迁徙的走廊，更是人类与自然对话的桥梁。当我们拿到2023年更新的全国自然保护区矢量数据时，如何用Python将其转化为可量化的生态洞察？本文将带您跨越从数据加载到空间分析的完整链路，通过geopandas、shapely和matplotlib等工具，挖掘保护区分布背后的生态保护优先级逻辑。

1. 数据准备与环境搭建

在开始分析前，我们需要配置专门的Python地理空间分析环境。推荐使用conda创建独立环境，避免库版本冲突：

conda create -n geo-analysis python=3.9 conda activate geo-analysis conda install -c conda-forge geopandas matplotlib contextily folium

关键库的作用说明：

• geopandas：地理空间数据分析的核心库，支持shp文件读取与空间运算
• shapely：处理几何对象的空间关系计算
• contextily：为地图添加底图
• folium：创建交互式地图

提示：自然保护区数据通常包含.shp（几何图形）、.shx（索引）和.dbf（属性表）三个必要文件，需确保文件完整

加载数据时，建议先检查坐标系。我国常用坐标系包括CGCS2000（EPSG:4490）和WGS84（EPSG:4326），转换方法如下：

import geopandas as gpd # 读取数据并检查原始坐标系 reserves = gpd.read_file('nature_reserves.shp') print(f"原始坐标系：{reserves.crs}") # 如有需要转换为WGS84 reserves = reserves.to_crs(epsg=4326)

2. 数据清洗与质量验证

原始保护区的数据质量直接影响分析结果。我们需要系统性地处理以下常见问题：

2.1 缺失值处理策略

2.2 几何有效性验证

保护区边界可能出现自相交、孔洞等拓扑错误，需使用shapely进行修复：

from shapely.validation import make_valid reserves['geometry'] = reserves['geometry'].apply( lambda geom: make_valid(geom) if not geom.is_valid else geom )

2.3 属性一致性检查

特别关注保护区的建立年份和面积单位：

• 年份格式统一（如将"2023年"转换为2023）
• 面积单位统一为平方公里
• 保护类型分类标准化

# 面积单位转换示例 reserves['area_km2'] = reserves['area_ha'] / 100 # 假设原始单位为公顷

3. 空间分析技术实战

3.1 保护区分布密度分析

使用核密度估计（KDE）识别保护区聚集热点：

from scipy.stats import gaussian_kde import numpy as np # 提取所有保护区的几何中心点 centroids = reserves.geometry.centroid coords = np.vstack([centroids.x, centroids.y]) # 计算核密度 kde = gaussian_kde(coords)(coords) reserves['density'] = kde

将密度分析结果可视化：

import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(12, 8)) reserves.plot(column='density', cmap='YlOrRd', legend=True, ax=ax) plt.title('全国自然保护区分布密度热力图') plt.show()

3.2 保护区与生态要素的空间关联

3.2.1 与生物多样性热点区叠加

假设我们有生物多样性优先区矢量数据，可以进行空间连接：

biodiversity = gpd.read_file('biodiversity_hotspots.shp') overlap = gpd.sjoin(reserves, biodiversity, how='inner', op='intersects') print(f"位于生物多样性热点区的保护区占比：{len(overlap)/len(reserves):.1%}")

3.2.2 与人口密度数据的空间关联

通过人口栅格数据提取保护区周边人口指标：

import rasterio from rasterstats import zonal_stats # 假设有人口密度tif文件 stats = zonal_stats( reserves.geometry, 'population.tif', stats=['mean', 'max'] ) reserves['pop_density'] = [x['mean'] for x in stats]

3.3 保护区网络连通性分析

评估保护区之间的生态廊道状况：

from sklearn.neighbors import NearestNeighbors # 计算最近邻保护区距离 coords = np.array([[g.centroid.x, g.centroid.y] for g in reserves.geometry]) nbrs = NearestNeighbors(n_neighbors=2).fit(coords) distances, _ = nbrs.kneighbors(coords) reserves['nearest_dist'] = distances[:,1] # 获取到最近保护区的距离

将结果按生态区进行分组统计：

connectivity = reserves.groupby('eco_region')['nearest_dist'].agg(['mean', 'std']) print(connectivity.sort_values('mean'))

4. 多维数据可视化呈现

4.1 交互式地图制作

使用folium创建可缩放探索的地图：

import folium m = folium.Map(location=[35, 105], zoom_start=5) for _, row in reserves.iterrows(): folium.GeoJson( row['geometry'], tooltip=f"{row['name']}<br>面积：{row['area_km2']:.1f}km²" ).add_to(m) m.save('reserves_map.html')

4.2 保护成效雷达图

评估不同类型保护区的综合表现：

categories = ['森林生态', '野生动物', '野生植物', '地质遗迹', '海洋海岸'] stats = reserves.groupby('type')[['area_km2', 'density', 'nearest_dist']].mean() fig = plt.figure(figsize=(8,8)) ax = fig.add_subplot(111, polar=True) for idx, row in stats.iterrows(): values = row.values.flatten().tolist() values += values[:1] # 闭合雷达图 angles = np.linspace(0, 2*np.pi, len(values), endpoint=False).tolist() angles += angles[:1] ax.plot(angles, values, label=idx) ax.set_xticks(angles[:-1]) ax.set_xticklabels(['面积', '密度', '连通性']) ax.legend() plt.show()

4.3 时间演变动态展示

分析保护区建立的时间分布特征：

timeline = reserves.groupby('establish_year').agg({ 'name': 'count', 'area_km2': 'sum' }).cumsum() fig, ax1 = plt.subplots(figsize=(10,6)) ax1.plot(timeline.index, timeline['name'], 'b-', label='数量') ax2 = ax1.twinx() ax2.plot(timeline.index, timeline['area_km2'], 'r--', label='总面积') ax1.set_xlabel('年份') ax1.set_ylabel('保护区数量', color='b') ax2.set_ylabel('总面积(km²)', color='r') plt.title('自然保护区历年累积增长趋势') plt.show()

5. 生态价值评估模型构建

基于上述分析，我们可以建立简单的保护优先级评分模型：

# 标准化各指标 reserves['score_size'] = (reserves['area_km2'] - reserves['area_km2'].min()) / \ (reserves['area_km2'].max() - reserves['area_km2'].min()) reserves['score_bio'] = reserves['density'] / reserves['density'].max() reserves['score_connect'] = 1 - (reserves['nearest_dist'] / reserves['nearest_dist'].max()) # 加权计算总分（可根据实际调整权重） weights = {'size': 0.4, 'bio': 0.3, 'connect': 0.3} reserves['priority_score'] = ( reserves['score_size'] * weights['size'] + reserves['score_bio'] * weights['bio'] + reserves['score_connect'] * weights['connect'] ) # 查看得分最高的10个保护区 top_reserves = reserves.nlargest(10, 'priority_score')[['name', 'type', 'priority_score']] print(top_reserves)

实际项目中，我们发现面积权重过高会导致小型但生态关键保护区被忽视。后来调整为加入物种丰富度等生物多样性指标后，评估结果更加合理。