土壤重金属数据实战分析Python自动化处理与空间可视化指南当面对数百个土壤采样点的重金属数据时传统的手工统计和Excel操作不仅效率低下还容易出错。本文将带你用Python生态中的Pandas、GeoPandas和Matplotlib工具链实现从原始数据到专业分析报告的全流程自动化处理。1. 数据准备与初步探索拿到土壤重金属数据的第一件事是检查数据质量和结构。典型的土壤重金属数据集通常包含采样点ID、经纬度坐标、采样深度以及各种元素Cr、Pb、Cu等的浓度值。假设我们有一个CSV格式的原始数据文件soil_samples.csvimport pandas as pd # 读取数据并显示前5行 df pd.read_csv(soil_samples.csv) print(df.head()) # 检查数据基本信息 print(df.info()) # 描述性统计 stats df.describe(percentiles[.25, .5, .75, .95]) print(stats.loc[[mean, std, min, 50%, max]])常见的数据质量问题包括缺失值特别是偏远地区的采样点异常值可能是录入错误或真实污染单位不统一ppm与mg/kg混用坐标系统不一致WGS84与GCJ02混用数据清洗示例# 处理缺失值 - 按元素中位数填充 elements [Cr, Pb, Cu, Zn, As, Hg, Cd] for elem in elements: df[elem].fillna(df[elem].median(), inplaceTrue) # 移除极端异常值超过99百分位 for elem in elements: upper_limit df[elem].quantile(0.99) df df[df[elem] upper_limit]2. 重金属统计分析技术2.1 基础统计量计算除了常规的平均值、标准差外土壤重金属分析需要特别关注超标倍数相对于区域背景值的比率变异系数反映元素空间分布均匀性偏态系数判断污染源分布特征# 定义区域背景值单位mg/kg background { Cr: 61.0, Pb: 26.0, Cu: 22.6, Zn: 74.2, As: 11.2, Hg: 0.065, Cd: 0.097 } # 计算各元素超标倍数 exceedance {} for elem in elements: mean_val df[elem].mean() exceedance[elem] mean_val / background[elem] # 转换为DataFrame便于查看 exceed_df pd.DataFrame.from_dict(exceedance, orientindex, columns[Exceedance]) print(exceed_df.sort_values(Exceedance, ascendingFalse))2.2 元素相关性分析重金属元素之间的相关性可以揭示共同污染源。使用Pandas计算相关系数矩阵corr_matrix df[elements].corr() # 可视化热力图 import seaborn as sns import matplotlib.pyplot as plt plt.figure(figsize(10,8)) sns.heatmap(corr_matrix, annotTrue, cmapcoolwarm, center0) plt.title(重金属元素相关性热力图) plt.show()典型的相关模式工业污染特征Pb、Zn、Cd高度相关农业污染特征As、Hg与其它元素相关性低自然背景特征Cr与其它元素相关性弱3. 空间分布可视化3.1 地理数据处理使用GeoPandas将普通数据框转换为地理数据框import geopandas as gpd from shapely.geometry import Point # 创建几何对象 geometry [Point(xy) for xy in zip(df[经度], df[纬度])] gdf gpd.GeoDataFrame(df, geometrygeometry, crsEPSG:4326) # 转换为适合分析的投影如Web墨卡托 gdf gdf.to_crs(epsg3857)3.2 制作元素分布图绘制Cd元素的空间分布散点图fig, ax plt.subplots(figsize(12,10)) # 背景地图需安装contextily import contextily as ctx gdf.plot(axax, alpha0.5) ctx.add_basemap(ax, sourcectx.providers.Stamen.Terrain) # 按Cd浓度分级着色 scatter gdf.plot(columnCd, cmapviridis, legendTrue, markersize50, axax, schemequantiles, edgecolorwhite, linewidth0.3) ax.set_title(土壤Cd含量空间分布, fontsize16) ax.set_axis_off() plt.show()3.3 空间插值可选虽然GeoPandas不直接支持空间插值但可以结合scipy进行简单插值from scipy.interpolate import griddata import numpy as np # 准备插值网格 x gdf.geometry.x y gdf.geometry.y z gdf[Cd] # 创建网格 xi np.linspace(x.min(), x.max(), 100) yi np.linspace(y.min(), y.max(), 100) xi, yi np.meshgrid(xi, yi) # 执行插值 zi griddata((x, y), z, (xi, yi), methodcubic) # 绘制等高线图 plt.figure(figsize(12,10)) contour plt.contourf(xi, yi, zi, levels15, cmapReds) plt.colorbar(contour) plt.scatter(x, y, ck, s5) plt.title(Cd含量空间插值结果) plt.show()4. 高级分析与报告生成4.1 污染评价指数计算常用的土壤污染评价方法包括单因子污染指数$P_i C_i / S_i$内梅罗综合指数$P_{综合} \sqrt{(P_{avg}^2 P_{max}^2)/2}$# 定义土壤环境质量标准二级标准pH7.5 standard { Cr: 250, Pb: 350, Cu: 100, Zn: 300, As: 25, Hg: 1.0, Cd: 0.6 } # 计算单因子指数 for elem in elements: df[f{elem}_PI] df[elem] / standard[elem] # 计算综合指数 df[P_avg] df[[f{e}_PI for e in elements]].mean(axis1) df[P_max] df[[f{e}_PI for e in elements]].max(axis1) df[P_nemerow] np.sqrt((df[P_avg]**2 df[P_max]**2)/2) # 污染等级划分 bins [0, 0.7, 1.0, 2.0, 3.0, float(inf)] labels [安全, 警戒线, 轻度污染, 中度污染, 重度污染] df[Pollution_Level] pd.cut(df[P_nemerow], binsbins, labelslabels)4.2 自动化报告生成使用Jupyter Notebook结合Python代码可以生成交互式分析报告。以下是一个静态报告生成示例from fpdf import FPDF class PDF(FPDF): def header(self): self.set_font(Arial, B, 15) self.cell(0, 10, 土壤重金属分析报告, 0, 1, C) def chapter_title(self, title): self.set_font(Arial, B, 12) self.cell(0, 10, title, 0, 1, L) self.ln(4) def chapter_body(self, body): self.set_font(Arial, , 12) self.multi_cell(0, 10, body) self.ln() # 创建PDF pdf PDF() pdf.add_page() # 添加内容 pdf.chapter_title(1. 分析概述) pdf.chapter_body(f本次分析共处理{len(df)}个采样点覆盖{len(df[地区].unique())}个地区。 主要分析了Cr、Pb、Cu、Zn、As、Hg、Cd七种重金属元素的含量分布特征。) pdf.chapter_title(2. 主要发现) pdf.chapter_body(f- Cd的平均超标倍数最高达到{exceedance[Cd]:.1f}倍\n f- {df[Pollution_Level].value_counts().idxmax()}区域占比最大) # 保存文件 pdf.output(soil_analysis_report.pdf)4.3 交互式可视化可选使用Plotly创建交互式地图import plotly.express as px fig px.scatter_mapbox(df, lat纬度, lon经度, colorCd, sizeCd, hover_name采样点ID, hover_dataelements, color_continuous_scalepx.colors.sequential.Viridis, zoom10) fig.update_layout(mapbox_styleopen-street-map) fig.update_layout(margin{r:0,t:0,l:0,b:0}) fig.show()这套分析方法已经成功应用于多个农业土壤调查项目。在实际项目中建议将上述代码封装成函数和类并通过配置文件管理不同地区的背景值和标准值。对于超大规模数据集10万采样点可以考虑使用Dask替代Pandas进行分布式计算。