简介面向生物多样性与气候变化交叉领域的数据分析学习者这是一套基于Python的完整可视化分析项目涵盖数据采集、清洗、统计与图表生成等环节适合环境科学、生态学专业学生及科研人员上手实践。包体共707个文件约148.22MB其中以Python脚本、Jupyter Notebook分析文档、CSV格式数据集为主辅以HTML交互页面、PNG/JPG可视化结果图DataGather目录存放原始数据ImageResource目录统一汇总生成图表各分析单元按模块分目录组织便于定位与复用。目前已有144人学习浏览。借助这套项目读者能够获取可直接运行的源代码和多年度物种分布、气候因子等真实数据同时得到一套清晰的可视化输出范例有助于快速复现分析流程并在此基础上拓展自己的研究思路与方法。1. 先想清楚再动手生物多样性与气候变化可视化分析到底要分析什么把几百万条物种出现记录和全球逐月气候栅格放在一起做可视化分析最花时间的不是写绘图代码而是先把 CSV 和 NetCDF 两类数据在同一个坐标系里对齐。物种点是表格里的经纬度气候场是网格上的数值矩阵两者精度、单位、时间口径全都不一样一旦对齐画图和做统计反而只要几十行 Python。这类项目通常围绕一条固定的技术路线展开先整理数据集再做地图与曲线可视化最后用统计检验给结论兜底。标题中的源代码与数据集提供的就是这样一条可复现路径。它适合刚迈过 Python 基础门槛、想拿真实数据练手的学习者也适合需要为报告、汇报或课程设计快速出图的人。2. 数据准备把物种出现记录与气候栅格统一成可复现的Python数据集生物多样性侧的数据源最常用的是 GBIF 的物种出现记录气候侧常见的是 CRU TS 这类月值 NetCDF或 WorldClim 的多年平均 GeoTIFF。它们格式不同、坐标意图不同拿到后别急着画图先把数据整理成后续模块都能用的长表。这也是数据集这个概念里最容易被忽略的部分raw 目录保持原始processed 目录放统一口径后的结果。下方工作流按“物种丰富度网格 → 气候年度场 → 栅格值提取”三步推进每一步都产出可继续被下游消费的 DataFrame。2.1 GBIF记录清洗从CSV到物种丰富度网格import pandas as pd gbif pd.read_csv( data/raw/gbif_observations.csv, usecols[scientificName, decimalLongitude, decimalLatitude, year], dtype{decimalLongitude: float32, decimalLatitude: float32}, ) # 缺失坐标和(0,0)假记录必须提前清掉 gbif gbif.dropna(subset[decimalLongitude, decimalLatitude]) gbif gbif.query( (decimalLongitude -180 and decimalLongitude 180) and (decimalLatitude -90 and decimalLatitude 90) ) grid gbif.copy() grid[lon_cell] grid[decimalLongitude].round() grid[lat_cell] grid[decimalLatitude].round() richness ( grid.groupby([scientificName, lon_cell, lat_cell])[year] .size() .reset_index(namerecords) ) species_per_cell ( richness.groupby([lon_cell, lat_cell])[scientificName] .nunique() .reset_index(namerichness) )read_csv 里 usecols 只保留四个列dtype 显式给 float32。GBIF 导出文件动辄几百万行float64 换 float32 可以让 DataFrame 体积缩小不少坐标值域在 [-180,180] 内float32 的有效数字足够。dropna 之后再做范围过滤是为了拦截坐标默认填 0 的脏记录否则图里会在坐标原点附近出现一个异常高值团逐步分析时很难排查。聚合逻辑分两步先按“物种 网格”统计出每个物种在每个 1 度网格里的记录量再用 nunique 统计网格内出现的物种数。用 round() 取整得到 1 度网格如果想做 0.5 度网格把经纬度先乘以 2 再 round最后除以 2 即可。这个思路跟目标检测数据集按类别统计的思路一致只是把类别换成了物种难点从标注框换成了坐标精度。提示如果只想跑通流程可以用随机生成的模拟物种点加一段模拟温升序列替代真实下载效果类似拿 iris 数据集验证代码那样等接口确认无误后再换正式数据。2.2 气候栅格读取NetCDF的变量、单位与时间聚合import xarray as xr ds xr.open_dataset(data/raw/cru_tmp.nc) print(ds) # 取温度变量先切出分析时段 temp ds[tmp].sel(timeslice(1990-01-01, 2020-12-31)) # 月值聚合到年1YE 是 xarray 新推荐的频率写法 annual_temp temp.resample(time1YE).mean(dimtime) # 距平以 1961-1990 为气候基准期 baseline temp.sel(timeslice(1961-01-01, 1990-12-31)).mean(dimtime) anomaly annual_temp - baseline打开 NetCDF 后第一件事永远是 print(ds)确认变量名、单位、时间坐标是否连续三个信息不确认就直接取数后面很容易在单位换算上翻车。CRU 系列的温度变量通常叫 tmp单位是摄氏度降水是 pre单位 mm/month拿到其他产品时以数据集自带 units 属性为准。月值转年值用 resample 的 1YE。旧代码里常见 1Y新版本 xarray 推荐 1YE 以明确表示年末对齐两者实际功能一致但升级依赖时可能看到弃用提示。距平计算选 1961-1990 作为基准期这是气候学里的常用约定你的数据集覆盖范围不足时可以自行调整。后续做全球尺度的可视化分析时距平比绝对温度更适合因为绝对温度本身就有强烈的纬度梯度会把时间变化信号掩盖掉。2.3 坐标对齐把栅格值提取到物种网格上sp species_per_cell.reset_index(dropTrue) sp[pid] sp.index # 用 Dataset 构造成批查询坐标 pts xr.Dataset( { lon: (point, sp[lon_cell]), lat: (point, sp[lat_cell]), } ) # nearest 提取每个网格中心的气候值 matched annual_temp.sel(lonpts.lon, latpts.lat, methodnearest) tmp_df matched.to_dataframe(nametemp_anom).reset_index() merged sp.merge(tmp_df, left_onpid, right_onpoint)sel 一次传入全部网格坐标返回的 DataArray 维度变成 (time, point)随后 to_dataframe 把 MultiIndex 展开再和物种丰富度做 merge。methodnearest 取距网格中心最近的气候格点速度快适合初筛气候场分辨率本身很粗时可以用 interp 做双线性插值代价是计算量明显增加。如果拿到的是 GeoTIFF 等投影坐标先用 rioxarray 的 reproject 转成等距经纬度再走同一套流程。数据形态常见来源Python 读取入口最需要确认的信息CSV / TSVGBIF 导出pandas.read_csv坐标列名、空值与 0,0 记录NetCDFCRU、ERA5 类xarray.open_dataset变量名、单位、时间频率GeoTIFFWorldClim 等rioxarray.open_rasterio投影方式、波段顺序与数值缩放系数这一步得到的 merged 表是后续所有可视化分析的数据基础建议即刻存成 data/processed/merged.parquet避免每次从头读原始文件。降水、日较差等气候特征用同样的提取流程各跑一遍最后合并成一张宽表。3. 可视化分析主体地图、双轴折线和回归带怎么画出气候-物种关系数据对齐完成之后可视化分析的工作量集中在三个层面空间、时间、关系。通常的顺序是先画地图确认分布没有异常再画时间序列看趋势最后把两个变量拉进同一张图找 pattern确认有 pattern 之后才交给统计检验。下面这组图形基本覆盖这类项目九成的出图需求。3.1 全球物种丰富度地图Cartopy 散点把空间分布画出来import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature fig plt.figure(figsize(12, 6)) ax fig.add_subplot(1, 1, 1, projectionccrs.PlateCarree()) ax.add_feature(cfeature.LAND, color#eaeaea) ax.add_feature(cfeature.COASTLINE, linewidth0.4) pc ax.scatter( sp[lon_cell], sp[lat_cell], csp[richness], cmapYlOrRd, s8, alpha0.7, linewidths0, transformccrs.PlateCarree(), rasterizedTrue, ) plt.colorbar(pc, axax, shrink0.7, labelspecies richness)projection 指定画布投影transform 声明数据本身的坐标系这里两者都是 PlateCarree等距经纬度适合全球展示。想要区域效果更好看可以把 projection 换成 Mollweide 或 Orthographic但 transform 保持数据原有坐标系不变。s 控制点大小网格点过万时 s 超过 20 就会糊成一片看不出密度时优先调小 alpha 而不是调大点。rasterizedTrue 在导出 PDF 时会把散点转成栅格文件体积小一个数量级。提示Cartopy 底层依赖 GEOS、Proj 等原生库源码编译容易卡在环境配置上。装不上时先用 conda 建环境或使用预编译 wheel临时替代方案是去掉投影直接用 plt.scatter 画经纬度先跑通业务逻辑再补地图底图。3.2 时间序列双轴图温度异常与物种丰富度同图展示fig, ax1 plt.subplots(figsize(10, 4)) ax1.plot(years, temp_series, color#d62728, lw1.6, labeltem anomaly) ax1.set_xlabel(year) ax1.set_ylabel(temperature anomaly (°C)) ax2 ax1.twinx() ax2.plot(years, richness_series, color#1f77b4, lw1.6, labelrichness) ax2.set_ylabel(species richness per cell) # 年份刻度稀疏化否则横坐标会挤成一排黑色竖线 ax1.set_xticks(range(years[0], years[-1] 1, 10)) ax1.tick_params(axisx, rotation0) fig.tight_layout()twinx 让左右轴共用同一个 x 轴适合把量纲完全不同的两条曲线叠在一起。但这也是可视化分析里最容易误导读者的写法两条轴的起点和刻度范围不同曲线纵向高度并不代表数值可比。更严谨的做法是两个序列各自做 z-score 标准化后画到同一根 y 轴z_t (temp_series - temp_series.mean()) / temp_series.std() z_r (richness_series - richness_series.mean()) / richness_series.std() plt.figure(figsize(10, 4)) plt.plot(years, z_t, color#d62728, labeltemp anomaly (z)) plt.plot(years, z_r, color#1f77b4, labelrichness (z)) plt.legend()横坐标密集的问题在时间跨度长时尤其明显跨度 30 年用 5 年步长跨度 100 年用 10 年步长直接传 range 生成刻度比依赖默认刻度更可控。3.3 关系可视化hexbin 密度图加二次趋势线import numpy as np fig, ax plt.subplots(figsize(8, 5)) hb ax.hexbin( merged[temp_anom], merged[richness], gridsize40, binslog, cmapBlues, ) plt.colorbar(hb, axax, labellog(count)) z np.polyfit(merged[temp_anom], merged[richness], 2) xx np.linspace(merged[temp_anom].min(), merged[temp_anom].max(), 100) ax.plot(xx, np.polyval(z, xx), color#c0392b, lw2) ax.set_xlabel(temperature anomaly (°C)) ax.set_ylabel(richness)关系图最常犯的错误是几万个点直接 scatter最后全变成一个深色聚合块。hexbin 把平面切成六边形格子用颜色表示落点密度binslog 让计数取对数稀疏区域的分布差异才看得出来。二次多项式拟合是这类气候-物种分析的经验选择生物对温度变化的响应经常是驼峰型中间最优、两端下降order2 能还原这个形状再高阶就容易被样本内的噪声带着走。想直接出带置信带的图用 seaborn 一行也能做到import seaborn as sns sns.regplot( xtemp_anom, yrichness, datamerged, order2, scatter_kws{s: 8, alpha: 0.3}, line_kws{color: #c0392b}, )regplot 默认带 95% 置信带适合汇报但数据量很大时它会把所有散点都画出来内存和渲染都会吃力先用 hexbin 看密度、再单独拟合是更节约的选择。3.4 给可视化分析加上交互Plotly 保存成 HTMLimport plotly.express as px fig px.scatter_geo( sp, latlat_cell, lonlon_cell, colorrichness, color_continuous_scaleYlOrRd, projectionnatural earth, ) fig.write_html(richness_map.html)scatter_geo 使用 plotly 自带的底图不需要单独申请地图服务 token比 scatter_mapbox 的接入成本低。交互模式下可以悬浮查看每个网格的坐标与丰富度适合放到网页或交付给不写代码的协作方。离线报告或打印场景里静态 Cartopy 图仍是更可靠的交付格式交互图作为补充。可视化场景推荐图形数据量与参数提示空间分布地图散点网格上万时 s 不超过 8开启 rasterized时间趋势双轴折线或 z-score 单轴刻度步长按跨度设为 5 或 10 年变量关系hexbin 趋势线大样本用 hexbin多项式阶数用 2交互分享plotly 输出 HTML离线展示用 scatter_geo 或退回静态图4. 从图像到结论用统计检验、分箱与降维确认关系强度与形态绘图只负责让 pattern 可见报告里最终要回答的还有三个问题关系有多强、形态是什么、哪些气候维度在起作用。这一步不需要高阶统计量三个常用的工具就能覆盖秩相关做强度检验分箱看非线性形态PCA 加层次聚类做多维气候降维观察。4.1 Spearman 秩相关先确认单调关系是否成立from scipy.stats import spearmanr valid merged.dropna(subset[temp_anom, richness]) rho, p_value spearmanr(valid[temp_anom], valid[richness]) print(fn{len(valid)}, rho{rho:.3f}, p{p_value:.2e})richness 是计数型变量分布偏斜且栅格提取后会带入空间自相关直接用 Pearson 容易被少数极端网格带偏。Spearman 看的是秩相关对极端值不敏感适合做第一道检验。看结果时重点看 n 和 rhon 是有效网格数决定检验功效rho 在 0.1 附近时即使 p 小于 0.001实际解释力也很有限论文和汇报里不要只报 p 值。4.2 分箱再看细节把非线性趋势暴露出来valid[temp_bin] pd.qcut( valid[temp_anom], q10, duplicatesdrop ).astype(str) grp ( valid.groupby(temp_bin, observedFalse)[richness] .agg([mean, std, count]) ) print(grp)qcut 按分位数切分保证每个桶里的网格数量接近适合偏斜的气候变量pd.cut 按数值等宽切极端值会把中间桶挤满。duplicatesdrop 处理的是大量 0 距平导致分位点重复报错的情况。分组后取均值和标准差把桶均值序列画成点线图可以和 4.1 的单调相关互相印证如果均值随桶号先升后降说明关系不是一条直线此时二次拟合的结论比线性相关更可靠。4.3 多维气候的降维观察PCA 投影与层次聚类分群from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from scipy.cluster.hierarchy import linkage, fcluster # 温度、降水、日较差三类气候距平都是同一套提取流程得到的列 clim_cols [temp_anom, prec_anom, diurnal_anom] X StandardScaler().fit_transform(merged[clim_cols]) pc PCA(n_components2).fit_transform(X) merged[PC1], merged[PC2] pc[:, 0], pc[:, 1] print(PCA(n_components2).fit(X).explained_variance_ratio_) Z linkage(X, methodward) merged[cluster] fcluster(Z, t3, criterionmaxclust)StandardScaler 先做标准化是因为温度和降水的单位不同不标准化时 PCA 第一主成分会几乎只反映降水方差失去“综合气候”的意义。explained_variance_ratio_ 打印前两个主成分的累计解释率低于 60% 时说明气候维度信息较分散PC1、PC2 只适合做探索性观察。linkage 的 ward 方法按合并时离差平方和最小化来决定分组风格偏保守得到的 cluster 回填到散点图上可以明显看出“冷湿”与“干热”这类气候组合它们对应的丰富度水平通常也分得很开。分析工具回答的问题使用注意scipy.stats.spearmanr两变量是否存在单调相关n 很大时 p 必显著重点看 rho 量级scipy.stats.pearsonr是否存在线性相关要求近似正态受强异常值影响大pd.qcut groupby关系是否线性、是否偏离分箱数不宜过少重复分位点要去重PCA 层次聚类哪些气候组合与丰富度相关必须标准化再结合载荷解释主成分5. 源代码侧的三个实用技巧模块拆分、并行读栅格与中间结果缓存到这一步分析流程已经完整了。真正决定这套代码能不能复用的是工程组织方式。以下三个技巧不涉及新技术只是把数据读取、可视化和统计解耦让项目在换数据、换年份时不用从头开始。5.1 模块拆分的项目结构project/ ├── data/ │ ├── raw/ # 原始下载数据只读不修改 │ └── processed/ # 清洗对齐后的中间结果 ├── src/ │ ├── preprocess.py # 数据读取与对齐 │ ├── plotting.py # 全部绘图函数 │ └── stats.py # 相关检验与降维分析 └── main.pymain.py 只负责按顺序调用三个模块的函数通常不超过 60 行。换区域、换年份时只需要改 preprocess 里的参数plotting 和 stats 完全不感知数据变化。5.2 懒加载与分块读取import xarray as xr rds xr.open_dataset( data/raw/clim.nc, chunks{time: -1, x: 1024, y: 1024}, ) anom rds[tmp].sel(timeslice(1990-01-01, 2020-12-31))chunks 参数让数据以 Dask 数组的方式懒加载open_dataset 本身不把整个文件读进内存后续的 sel、mean 会在真正需要出结果时才触发计算。块大小没有绝对最优解内存 32GB 的机器x 和 y 维度各设 1024 或 2048 都常见内存紧张时调到 512代价是 IO 次数变多。5.3 中间结果缓存from pathlib import Path import pandas as pd def cache_parquet(key, producer): out Path(data/cache) / f{key}.parquet if out.exists(): return pd.read_parquet(out) df producer() out.parent.mkdir(parentsTrue, exist_okTrue) df.to_parquet(out, indexFalse) return df merged cache_parquet(merge_global_1990_2020_v1, lambda: build_merged(...))缓存 key 的语义化命名建议写成“数据名_空间范围_时间范围_处理版本”例如 merge_global_1990_2020_v1。当重新下载了气候数据、或把年份范围从 1990-2020 改成 2000-2023只需要把这个 key 里的时间部分改掉后面的统计和绘图模块一行都不用动。本文还有配套的精品资源点击获取