从数据到地图：用Python复现中国旱区土壤碳分布图（附代码与数据）
从数据到地图：用Python复现中国旱区土壤碳分布图（附代码与数据）
当我们在学术论文中看到一张精美的科学图表时，常常会好奇：这些数据从何而来？又是如何被转化为直观可视的图形？本文将以中国旱区土壤碳分布研究为例，带你用Python完整复现科研论文中的空间分布图，从原始数据获取到最终可视化，手把手教你掌握地理空间数据分析的核心技能。
1. 环境准备与数据获取
1.1 搭建Python分析环境
地理空间分析需要特定的工具链支持。推荐使用conda创建独立环境，避免包冲突：
conda create -n soil_carbon python=3.9 conda activate soil_carbon conda install -c conda-forge geopandas rasterio matplotlib numpy pandas关键库说明：
- geopandas：地理空间数据处理的核心库
- rasterio：栅格数据读写与处理
- matplotlib：科学可视化
- numpy/pandas：数值计算与表格处理
提示：conda-forge通道能确保所有地理空间库的依赖关系正确解析
1.2 获取土壤碳数据
中国旱区土壤碳数据可从以下公开来源获取：
HWSD（Harmonized World Soil Database）：
- 全球1km分辨率土壤数据库
- 包含有机碳(SOC)和无机碳(SIC)含量数据
- 下载地址：FAO官网土壤门户
WorldClim气候数据：
- 提供降水、温度等气候指标
- 可用于验证干旱梯度变化
- 下载地址：WorldClim官网
中国行政区划数据：
- 从国家基础地理信息中心获取省界矢量数据
- 用于裁剪和可视化
import geopandas as gpd from rasterio.warp import calculate_default_transform, reproject # 示例代码：读取HWSD数据 hwsd_path = "path/to/HWSD_China.tif" with rasterio.open(hwsd_path) as src: china_carbon = src.read(1) # 读取土壤碳含量波段2. 数据处理与空间分析
2.1 定义中国旱区范围
根据联合国环境规划署定义，旱区包括：
- 干旱指数(AI) < 0.65的区域
- 涵盖中国北方大部分地区
# 计算干旱指数AI = P/PET def calculate_ai(precip, pet): return precip / pet # 从WorldClim数据计算AI precip = rasterio.open("wc2.1_30s_prec.tif") # 年降水量 pet = rasterio.open("global_et0.tif") # 潜在蒸散发 ai = calculate_ai(precip.read(1), pet.read(1)) # 创建旱区掩膜 dryland_mask = ai < 0.652.2 土壤碳数据预处理
原始数据通常需要以下处理步骤：
- 投影转换：
- 统一所有数据到相同坐标系（如WGS84）
- 重采样：
- 将不同分辨率数据统一到相同网格
- 异常值处理：
- 过滤明显错误的土壤碳值
# 投影转换示例 def reproject_raster(src_file, dst_crs): transform, width, height = calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds) kwargs = src.meta.copy() kwargs.update({ 'crs': dst_crs, 'transform': transform, 'width': width, 'height': height }) with rasterio.open('reprojected.tif', 'w', **kwargs) as dst: reproject( source=rasterio.band(src, 1), destination=rasterio.band(dst, 1), src_transform=src.transform, src_crs=src.crs, dst_transform=transform, dst_crs=dst_crs, resampling=Resampling.nearest)3. 空间可视化实现
3.1 创建分级色彩图
土壤碳分布通常采用分级设色法表示：
| 碳含量范围 (kg/m²) | 颜色代码 | 等级 |
|---|---|---|
| 0-2 | #FFFFCC | 1 |
| 2-4 | #C7E9B4 | 2 |
| 4-6 | #7FCDBB | 3 |
| 6-8 | #41B6C4 | 4 |
| 8-10 | #2C7FB8 | 5 |
| >10 | #253494 | 6 |
import matplotlib.colors as colors # 自定义颜色映射 cmap = colors.ListedColormap([ '#FFFFCC', '#C7E9B4', '#7FCDBB', '#41B6C4', '#2C7FB8', '#253494']) bounds = [0, 2, 4, 6, 8, 10, 15] norm = colors.BoundaryNorm(bounds, cmap.N)3.2 绘制专题地图
完整地图应包含以下要素：
- 主图：土壤碳空间分布
- 图例：颜色与数值对应关系
- 比例尺和指北针
- 数据来源说明
import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable fig, ax = plt.subplots(figsize=(12, 8)) # 绘制底图 china = gpd.read_file('china_province.shp') china.boundary.plot(ax=ax, linewidth=0.5, color='gray') # 绘制土壤碳数据 im = ax.imshow(total_carbon, cmap=cmap, norm=norm, extent=[xmin, xmax, ymin, ymax]) # 添加颜色条 divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.1) cbar = fig.colorbar(im, cax=cax) cbar.set_label('Soil Carbon Stock (kg/m²)') # 设置标题和注释 ax.set_title('Soil Carbon Distribution in Chinese Drylands', fontsize=14) ax.text(0.01, 0.01, 'Data: HWSD v1.2 | Created with Python', transform=ax.transAxes, fontsize=8) plt.tight_layout() plt.savefig('soil_carbon_distribution.png', dpi=300)4. 进阶分析与验证
4.1 沿干旱梯度的碳变化分析
为验证论文中"有机碳下降、无机碳增加"的结论，我们可以提取沿干旱梯度的变化曲线：
# 创建采样点沿干旱梯度 transect_points = [ (85.0, 44.5), # 新疆北部 (100.3, 40.2), # 河西走廊 (112.5, 41.5), # 内蒙古中部 (116.4, 39.9) # 华北平原北部 ] # 提取各点碳含量 soc_values = [extract_value(soc_data, point) for point in transect_points] sic_values = [extract_value(sic_data, point) for point in transect_points] ai_values = [extract_value(ai_data, point) for point in transect_points] # 绘制变化曲线 plt.figure(figsize=(10, 5)) plt.plot(ai_values, soc_values, 'b-o', label='SOC') plt.plot(ai_values, sic_values, 'r-s', label='SIC') plt.xlabel('Aridity Index') plt.ylabel('Carbon Content (kg/m²)') plt.legend() plt.grid(True)4.2 空间自相关分析
使用Moran's I指数检验土壤碳的空间聚集性：
from esda.moran import Moran import libpysal as lps # 计算空间权重矩阵 w = lps.weights.Queen.from_dataframe(china_gdf) # 计算Moran's I moran = Moran(china_gdf['total_carbon'], w) print(f"Moran's I: {moran.I:.3f}") print(f"P-value: {moran.p_norm:.4f}")典型结果解读：
- Moran's I > 0：空间正相关（聚集分布）
- Moran's I < 0：空间负相关（分散分布）
- P-value < 0.05：统计显著
5. 完整代码架构
推荐的项目目录结构：
soil_carbon_map/ ├── data/ # 原始数据 │ ├── HWSD_China.tif │ ├── china_province.shp │ └── worldclim/ ├── notebooks/ # Jupyter笔记本 │ └── analysis.ipynb ├── scripts/ # Python脚本 │ ├── data_processing.py │ └── visualization.py ├── output/ # 结果输出 │ ├── figures/ │ └── processed_data/ └── requirements.txt # 依赖列表核心处理流程：
数据预处理脚本(
data_processing.py)：- 数据下载与格式转换
- 坐标系统一
- 掩膜裁剪
分析脚本(
analysis.py)：- 空间统计计算
- 干旱梯度分析
- 相关性检验
可视化脚本(
visualization.py)：- 专题地图生成
- 变化曲线绘制
- 交互式地图输出
提示：使用Jupyter Notebook进行探索性分析，成熟后迁移到.py脚本实现自动化
6. 常见问题与优化技巧
6.1 性能优化策略
处理全国范围高分辨率数据时，可能遇到内存问题：
- 分块处理：使用rasterio的窗口读取功能
with rasterio.open('large.tif') as src: for window in src.block_windows(): data = src.read(window=window) # 处理每个分块- Dask并行：适用于大规模计算
import dask.array as da carbon_dask = da.from_array(carbon_data, chunks=(1000, 1000)) result = (carbon_dask * factor).compute()6.2 可视化增强技巧
- 添加地形阴影：
- 使用SRTM高程数据增强立体感
- 动态可视化：
- 用folium创建交互式地图
import folium m = folium.Map(location=[35, 105], zoom_start=4) folium.raster_layers.ImageOverlay( image=carbon_data, bounds=[[15, 70], [55, 140]], colormap=lambda x: (1, 0, 0, x) ).add_to(m)- 多图组合：
- 使用subplots同时显示SOC和SIC分布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6)) im1 = ax1.imshow(soc_data, cmap='YlOrBr', vmax=10) im2 = ax2.imshow(sic_data, cmap='Blues', vmax=15) # 添加颜色条和标题等