以下是Python基础学习第二十九课的完整内容，聚焦Python在环境科学与地球科学中的跨学科应用，涵盖气候变化建模、地理空间数据分析与灾害预警系统，带你探索地球系统的数字化模拟：

 

Python基础学习第二十九课：环境科学与地球科学计算

 

一、课程目标

 

1. 掌握气候时间序列分析（Xarray+StatsModels）

2. 学习地理空间数据处理（GeoPandas+Rasterio）

3. 构建洪水灾害预测模型

4. 开发空气质量实时监测平台

 

二、气候变化建模与时间序列分析

 

1. 全球温度数据集处理（Xarray）

 

(1) 加载CRU TS气候数据

 

import xarray as xr

 

# 打开NetCDF格式的全球温度数据集

ds = xr.open_dataset('cru_ts4.05.1901.2020.tmp.dat.nc')

print(ds) # 查看数据维度（时间、纬度、经度）

 

# 提取中国区域（示例经纬度范围）

china_temp = ds['tmp'].sel(

    latitude=slice(53.5, 18.0), # 北纬18°-53.5°

    longitude=slice(73.5, 135.0) # 东经73.5°-135°

)

 

(2) 计算区域平均温度趋势

 

import numpy as np

 

# 计算时间维度上的线性趋势

def calculate_trend(da):

    """计算xarray DataArray的时间趋势"""

    time = np.arange(len(da.time))

    slope = np.polyfit(time, da.mean(dim=['latitude', 'longitude']), 1)[0]

    return slope * 10 # 转换为每十年变化量(°C/10a)

 

china_trend = calculate_trend(china_temp)

print(f"中国区域温度变化趋势: {china_trend:.3f} °C/10年")

 

2. 极端气候事件检测（Pandas）

 

# 将xarray转换为Pandas DataFrame

df = china_temp.to_dataframe().reset_index()

 

# 计算滑动平均温度

df['rolling_mean'] = df.groupby('latitude')['tmp'].transform(

    lambda x: x.rolling(window=30, center=True).mean()

)

 

# 检测热浪事件（连续5天超过阈值）

threshold = df['rolling_mean'].quantile(0.95)

heatwaves = df[df['tmp'] > threshold].groupby(

    (df['tmp'] <= threshold).cumsum()

).filter(lambda x: len(x) >= 5)

 

print(f"检测到{len(heatwaves)//5}次热浪事件")

 

三、地理空间数据分析

 

1. 洪水风险区划（GeoPandas）

 

(1) 加载数字高程模型(DEM)

 

import rasterio

from rasterio.plot import show

 

# 读取DEM数据

with rasterio.open('dem.tif') as src:

    dem = src.read(1)

    transform = src.transform

 

# 可视化DEM

show(dem, cmap='terrain')

 

(2) 水文分析（Flow Accumulation）

 

from rasterio.features import sieve

from scipy.ndimage import uniform_filter

 

# 地形平滑处理

smoothed_dem = uniform_filter(dem, size=3)

 

# 计算坡度（简化版）

dx, dy = np.gradient(smoothed_dem, transform.a, transform.e)

slope = np.arctan(np.sqrt(dx**2 + dy**2)) * 180 / np.pi

 

# 提取水系（坡度<5°的区域）

river_mask = slope < 5

river_pixels = sieve(river_mask.astype('uint8'), size=100) # 连通区域>100像素

 

# 可视化水系

show(river_pixels, cmap='Blues')

 

2. 洪水淹没模拟（洪水深度计算）

 

# 假设降雨量数据（单位：mm）

rainfall = 200 # 200mm降雨

 

# 计算潜在淹没深度（简化模型）

# 淹没深度 = 降雨量 - 地形储水能力（假设储水能力=坡度*10）

flood_depth = np.where(

    slope < 5,

    np.maximum(rainfall - (slope * 10), 0),

    0

)

 

# 可视化淹没深度

import matplotlib.pyplot as plt

plt.imshow(flood_depth, cmap='Blues', vmin=0, vmax=50)

plt.colorbar(label='淹没深度 (mm)')

plt.title('模拟洪水淹没深度')

plt.show()

 

四、灾害预警系统开发

 

1. 实时降雨数据API接入

 

import requests

import pandas as pd

 

# 从NOAA API获取实时降雨数据

url = "https://www.ncdc.noaa.gov/cdo-web/api/v2/data"

headers = {"token": "YOUR_API_TOKEN"}

params = {

    "datasetid": "GHCND",

    "locationid": "FIPS:US",

    "startdate": "2025-07-27",

    "enddate": "2025-07-27",

    "datatypeid": "PRCP",

    "limit": 1000

}

 

response = requests.get(url, headers=headers, params=params)

data = response.json()['results']

 

# 转换为DataFrame

rain_df = pd.DataFrame(data)

rain_df['date'] = pd.to_datetime(rain_df['date'])

rain_df['prcp'] = rain_df['value'] / 1000 # 转换为mm

 

# 计算每小时降雨量

hourly_rain = rain_df.resample('H', on='date')['prcp'].sum()

 

2. 洪水预警阈值判断

 

# 定义预警阈值（单位：mm/h）

warning_levels = {

    'blue': 10, # 蓝色预警

    'yellow': 25, # 黄色预警

    'orange': 50, # 橙色预警

    'red': 100 # 红色预警

}

 

# 实时判断预警级别

current_rain = hourly_rain.iloc[-1]

warning_level = next(

    (level for level, threshold in warning_levels.items() if current_rain >= threshold),

    'normal'

)

 

print(f"当前降雨量: {current_rain:.1f} mm/h, 预警级别: {warning_level}")

 

五、综合案例：气候变化对农业影响评估

 

1. 系统架构

 

气候数据(CMIP6) → 作物生长模型(DSSAT) → 经济损失评估 → 可视化仪表盘

 

2. 关键代码实现

 

(1) CMIP6数据提取（Xarray）

 

# 打开CMIP6温度预测数据

ds_future = xr.open_dataset('tas_Amon_CESM2_ssp585_r1i1p1f1_gr_202001-204012.nc')

 

# 计算未来30年温度变化

future_trend = calculate_trend(ds_future['tas'].sel(time=slice('2020', '2040')))

historical_trend = calculate_trend(china_temp.sel(time=slice('1980', '2020')))

 

print(f"历史趋势: {historical_trend:.3f} °C/10a, 未来预测: {future_trend:.3f} °C/10a")

 

(2) 作物产量损失模型

 

# 简化模型：温度每升高1°C，小麦减产2%

wheat_yield_loss = (future_trend - historical_trend) * 2

 

# 计算经济损失（假设小麦价格5元/kg，产量1亿吨）

economic_loss = wheat_yield_loss * 0.02 * 1e8 * 5

print(f"预测经济损失: {economic_loss:.2f} 亿元")

 

六、课后练习

 

1. 基础题：

   - 用Xarray分析本地气象站温度数据，计算年际变化趋势。

   - 使用GeoPandas绘制城市洪涝风险分区图。

2. 进阶题：

   - 开发一个基于机器学习的地震震级预测模型（使用USGS地震目录）。

   - 构建实时空气质量预警系统（整合PM2.5、NO2等多污染物数据）。

 

七、常见问题解答

 

1. Q：如何获取高分辨率气候数据？A：

   - 免费源：ERA5（ECMWF）、NASA POWER

   - 商业数据：WorldClim、PRISM

2. Q：地理空间分析的内存优化技巧？A：

   - 使用分块处理（

"rasterio.windows"）

   - 将栅格数据转换为稀疏矩阵格式

3. Q：灾害预警系统的实时性保障？A：

   - 使用消息队列（如Kafka）处理传感器数据流

   - 部署边缘计算节点（如树莓派+LoRa网关）

 

通过本课，你已掌握环境科学与地球科学计算的数字化工具！从气候变化建模到灾害预警系统，具备解决全球可持续发展挑战的技术能力。下一步可探索碳中和路径优化或行星科学数据分析，持续拓展地球计算的前沿边界。 🌍🔍