从数据到地图用Python复现中国旱区土壤碳分布图附代码与数据当我们在学术论文中看到一张精美的科学图表时常常会好奇这些数据从何而来又是如何被转化为直观可视的图形本文将以中国旱区土壤碳分布研究为例带你用Python完整复现科研论文中的空间分布图从原始数据获取到最终可视化手把手教你掌握地理空间数据分析的核心技能。1. 环境准备与数据获取1.1 搭建Python分析环境地理空间分析需要特定的工具链支持。推荐使用conda创建独立环境避免包冲突conda create -n soil_carbon python3.9 conda activate soil_carbon conda install -c conda-forge geopandas rasterio matplotlib numpy pandas关键库说明geopandas地理空间数据处理的核心库rasterio栅格数据读写与处理matplotlib科学可视化numpy/pandas数值计算与表格处理提示conda-forge通道能确保所有地理空间库的依赖关系正确解析1.2 获取土壤碳数据中国旱区土壤碳数据可从以下公开来源获取HWSDHarmonized 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( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest)3. 空间可视化实现3.1 创建分级色彩图土壤碳分布通常采用分级设色法表示碳含量范围 (kg/m²)颜色代码等级0-2#FFFFCC12-4#C7E9B424-6#7FCDBB36-8#41B6C448-10#2C7FB8510#2534946import 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(axax, linewidth0.5, colorgray) # 绘制土壤碳数据 im ax.imshow(total_carbon, cmapcmap, normnorm, extent[xmin, xmax, ymin, ymax]) # 添加颜色条 divider make_axes_locatable(ax) cax divider.append_axes(right, size5%, pad0.1) cbar fig.colorbar(im, caxcax) cbar.set_label(Soil Carbon Stock (kg/m²)) # 设置标题和注释 ax.set_title(Soil Carbon Distribution in Chinese Drylands, fontsize14) ax.text(0.01, 0.01, Data: HWSD v1.2 | Created with Python, transformax.transAxes, fontsize8) plt.tight_layout() plt.savefig(soil_carbon_distribution.png, dpi300)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, labelSOC) plt.plot(ai_values, sic_values, r-s, labelSIC) plt.xlabel(Aridity Index) plt.ylabel(Carbon Content (kg/m²)) plt.legend() plt.grid(True)4.2 空间自相关分析使用Morans I指数检验土壤碳的空间聚集性from esda.moran import Moran import libpysal as lps # 计算空间权重矩阵 w lps.weights.Queen.from_dataframe(china_gdf) # 计算Morans I moran Moran(china_gdf[total_carbon], w) print(fMorans I: {moran.I:.3f}) print(fP-value: {moran.p_norm:.4f})典型结果解读Morans I 0空间正相关聚集分布Morans 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(windowwindow) # 处理每个分块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_start4) folium.raster_layers.ImageOverlay( imagecarbon_data, bounds[[15, 70], [55, 140]], colormaplambda 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, cmapYlOrBr, vmax10) im2 ax2.imshow(sic_data, cmapBlues, vmax15) # 添加颜色条和标题等
从数据到地图：用Python复现中国旱区土壤碳分布图（附代码与数据）
发布时间：2026/6/3 5:47:03
从数据到地图用Python复现中国旱区土壤碳分布图附代码与数据当我们在学术论文中看到一张精美的科学图表时常常会好奇这些数据从何而来又是如何被转化为直观可视的图形本文将以中国旱区土壤碳分布研究为例带你用Python完整复现科研论文中的空间分布图从原始数据获取到最终可视化手把手教你掌握地理空间数据分析的核心技能。1. 环境准备与数据获取1.1 搭建Python分析环境地理空间分析需要特定的工具链支持。推荐使用conda创建独立环境避免包冲突conda create -n soil_carbon python3.9 conda activate soil_carbon conda install -c conda-forge geopandas rasterio matplotlib numpy pandas关键库说明geopandas地理空间数据处理的核心库rasterio栅格数据读写与处理matplotlib科学可视化numpy/pandas数值计算与表格处理提示conda-forge通道能确保所有地理空间库的依赖关系正确解析1.2 获取土壤碳数据中国旱区土壤碳数据可从以下公开来源获取HWSDHarmonized 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( sourcerasterio.band(src, 1), destinationrasterio.band(dst, 1), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.nearest)3. 空间可视化实现3.1 创建分级色彩图土壤碳分布通常采用分级设色法表示碳含量范围 (kg/m²)颜色代码等级0-2#FFFFCC12-4#C7E9B424-6#7FCDBB36-8#41B6C448-10#2C7FB8510#2534946import 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(axax, linewidth0.5, colorgray) # 绘制土壤碳数据 im ax.imshow(total_carbon, cmapcmap, normnorm, extent[xmin, xmax, ymin, ymax]) # 添加颜色条 divider make_axes_locatable(ax) cax divider.append_axes(right, size5%, pad0.1) cbar fig.colorbar(im, caxcax) cbar.set_label(Soil Carbon Stock (kg/m²)) # 设置标题和注释 ax.set_title(Soil Carbon Distribution in Chinese Drylands, fontsize14) ax.text(0.01, 0.01, Data: HWSD v1.2 | Created with Python, transformax.transAxes, fontsize8) plt.tight_layout() plt.savefig(soil_carbon_distribution.png, dpi300)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, labelSOC) plt.plot(ai_values, sic_values, r-s, labelSIC) plt.xlabel(Aridity Index) plt.ylabel(Carbon Content (kg/m²)) plt.legend() plt.grid(True)4.2 空间自相关分析使用Morans I指数检验土壤碳的空间聚集性from esda.moran import Moran import libpysal as lps # 计算空间权重矩阵 w lps.weights.Queen.from_dataframe(china_gdf) # 计算Morans I moran Moran(china_gdf[total_carbon], w) print(fMorans I: {moran.I:.3f}) print(fP-value: {moran.p_norm:.4f})典型结果解读Morans I 0空间正相关聚集分布Morans 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(windowwindow) # 处理每个分块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_start4) folium.raster_layers.ImageOverlay( imagecarbon_data, bounds[[15, 70], [55, 140]], colormaplambda 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, cmapYlOrBr, vmax10) im2 ax2.imshow(sic_data, cmapBlues, vmax15) # 添加颜色条和标题等
![P10377 [GESP202403 六级] 好斗的牛](http://pic.xiahunao.cn/yaotu/P10377 [GESP202403 六级] 好斗的牛)
