Python实战：用ERA5数据绘制城市热浪地图（附完整代码）

Python实战：用ERA5数据绘制城市热浪地图（附完整代码）

去年夏天，我在参与一个城市气候适应性规划项目时，遇到了一个棘手的问题：如何量化评估城市热浪的时空分布特征？当时手头有几十年的气象观测站数据，但站点稀疏、空间代表性有限，很难绘制出精细化的热浪风险图。直到我开始尝试使用欧洲中期天气预报中心（ECMWF）的ERA5再分析数据，才真正打开了城市热气候研究的新视野。

ERA5作为目前全球最权威的再分析数据集之一，提供了0.25°×0.25°（约25公里）的高时空分辨率气象数据，覆盖全球范围且时间序列完整。对于城市气候研究者、公共卫生分析师和城市规划师来说，这套数据简直是宝藏——你不再需要依赖零星的气象站点，而是可以直接获取每个网格点的温度、湿度、风速等数十个变量，时间跨度从1950年至今，每小时一次。

但数据到手只是第一步。真正的挑战在于：如何从海量的NetCDF文件中提取出对城市热浪研究有用的信息？如何定义热浪事件？如何将抽象的温度数据转化为直观的空间分布图？更重要的是，如何让这个过程可重复、可自动化，以便应用于不同城市、不同时期的研究？

这篇文章就是我在这个探索过程中的经验总结。我将带你从零开始，用Python构建一套完整的热浪识别与可视化流程。这套方法不仅适用于学术研究，也能为城市规划部门的热岛缓解策略、公共卫生部门的热相关疾病预警提供数据支撑。无论你是刚接触气候数据分析的新手，还是有经验的研究者，都能从中找到实用的代码片段和思路启发。

1. ERA5数据获取与预处理：从云端到本地

1.1 理解ERA5数据结构和访问方式

ERA5数据存储在哥白尼气候数据存储（CDS）平台上，提供两种主要访问方式：手动下载和API自动下载。对于长期研究项目，我强烈推荐使用API方式，虽然初次配置稍显繁琐，但一旦设置完成，后续的数据获取就能完全自动化。

CDS API的核心是cdsapi这个Python库。安装过程很简单：

pip install cdsapi

但安装库只是第一步，真正的关键是获取API密钥。你需要到CDS官网注册账号，然后在用户设置页面找到“API密钥”部分。系统会生成一个包含URL和密钥的文本，你需要将其保存到本地~/.cdsapirc文件中（Windows用户保存在C:\Users\用户名\.cdsapirc）。

注意：API密钥是你的个人凭证，千万不要上传到公开的代码仓库。我通常会在代码中使用环境变量来管理这类敏感信息。

1.2 批量下载城市区域气温数据

假设我们要研究长三角地区2022年的热浪事件，需要下载该区域1970-2022年夏季（6-8月）的2米气温数据。下面是一个完整的下载脚本：

import cdsapi
import calendar
from pathlib import Path
import time

def download_era5_temperature(years, months, area, output_dir):
    """
    批量下载ERA5 2米气温数据
    
    参数:
    years: 年份列表，如[1970, 1971, ..., 2022]
    months: 月份列表，如['06', '07', '08']
    area: 区域范围 [北纬, 西经, 南纬, 东经]
    output_dir: 输出目录路径
    """
    c = cdsapi.Client()
    output_path = Path(output_dir)
    output_path.mkdir(parents=True, exist_ok=True)
    
    for year in years:
        for month in months:
            # 构建文件名
            filename = f"era5_t2m_{year}_{month}_长三角.nc"
            filepath = output_path / filename
            
            # 如果文件已存在，跳过下载
            if filepath.exists():
                print(f"文件已存在: {filename}")
                continue
            
            # 获取该月的天数
            days_in_month = calendar.monthrange(year, int(month))[1]
            days = [f"{day:02d}" for day in range(1, days_in_month + 1)]
            
            # 24小时数据
            hours = [f"{h:02d}:00" for h in range(24)]
            
            print(f"正在下载 {year}年{month}月数据...")
            
            try:
                c.retrieve(
                    'reanalysis-era5-single-levels',
                    {
                        'product_type': 'reanalysis',
                        'variable': '2m_temperature',
                        'year': str(year),
                        'month': month,
                        'day': days,
                        'time': hours,
                        'area': area,  # 长三角区域: [32.5, 118.0, 29.5, 122.0]
                        'format': 'netcdf',
                    },
                    str(filepath)
                )
                print(f"✓ 下载完成: {filename}")
                time.sleep(2)  # 避免请求过于频繁
                
            except Exception as e:
                print(f"✗ 下载失败 {year}-{month}: {e}")
                continue

# 使用示例
if __name__ == "__main__":
    # 定义下载范围
    years = list(range(1970, 2023))  # 1970-2022年
    months = ['06', '07', '08']  # 夏季月份
    
    # 长三角区域 [北纬, 西经, 南纬, 东经]
    yangtze_delta_area = [32.5, 118.0, 29.5, 122.0]
    
    # 输出目录
    output_directory = "./data/era5_raw"
    
    download_era5_temperature(years, months, yangtze_delta_area, output_directory)

这个脚本有几个关键点值得注意：

1. 分月下载：ERA5 API对单次请求的数据量有限制，所以按年月分批下载是稳妥的做法
2. 错误处理：网络请求可能失败，加入try-except结构确保程序不会中途崩溃
3. 速率控制：每次下载后暂停2秒，避免触发API的速率限制
4. 文件检查：下载前检查文件是否已存在，避免重复下载浪费资源

1.3 数据格式转换与质量检查

下载的NetCDF文件需要经过初步处理才能用于分析。我通常会用xarray库来读取和检查数据：

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt

def inspect_era5_file(filepath):
    """检查ERA5数据文件的基本信息"""
    ds = xr.open_dataset(filepath)
    
    print("=== 数据集概览 ===")
    print(ds)
    
    print("\n=== 变量信息 ===")
    for var in ds.data_vars:
        print(f"{var}: {ds[var].attrs.get('long_name', '无描述')}")
        print(f"  单位: {ds[var].attrs.get('units', '未知')}")
        print(f"  形状: {ds[var].shape}")
    
    print("\n=== 坐标信息 ===")
    print(f"时间范围: {ds.time.values[0]} 到 {ds.time.values[-1]}")
    print(f"经度范围: {ds.longitude.values.min():.2f}° 到 {ds.longitude.values.max():.2f}°")
    print(f"纬度范围: {ds.latitude.values.min():.2f}° 到 {ds.latitude.values.max():.2f}°")
    
    # 可视化检查
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    # 时间序列示例
    sample_point = ds.t2m.isel(longitude=10, latitude=10)
    axes[0].plot(sample_point.time, sample_point - 273.15)  # 转换为摄氏度
    axes[0].set_title("单点温度时间序列")
    axes[0].set_ylabel("温度 (°C)")
    axes[0].tick_params(axis='x', rotation=45)
    
    # 空间分布示例（第一个时次）
    temp_map = ds.t2m.isel(time=0) - 273.15
    im = axes[1].imshow(temp_map, cmap='RdBu_r', 
                       extent=[ds.longitude.min(), ds.longitude.max(),
                               ds.latitude.min(), ds.latitude.max()])
    axes[1].set_title("空间分布示例")
    axes[1].set_xlabel("经度")
    axes[1].set_ylabel("纬度")
    plt.colorbar(im, ax=axes[1], label='温度 (°C)')
    
    # 数据统计
    temp_data = ds.t2m.values - 273.15
    axes[2].hist(temp_data.flatten(), bins=50, alpha=0.7)
    axes[2].set_title("温度分布直方图")
    axes[2].set_xlabel("温度 (°C)")
    axes[2].set_ylabel("频数")
    
    plt.tight_layout()
    plt.show()
    
    ds.close()
    return True

# 使用示例
file_path = "./data/era5_raw/era5_t2m_2022_06_长三角.nc"
inspect_era5_file(file_path)

通过这个检查步骤，你可以确认数据是否完整、单位是否正确（ERA5的温度数据默认是开尔文，需要减去273.15转换为摄氏度）、时空范围是否符合预期。

2. 热浪事件的定义与识别算法

2.1 热浪的科学定义与业务标准

在开始编码之前，我们需要明确什么是"热浪"。不同地区、不同研究目的对热浪的定义可能不同，但核心思想是一致的：一段异常高温的持续时期。常见的定义方法包括：

在实际项目中，我通常采用相对阈值法，因为它能更好地反映"相对于当地气候常态的异常情况"。比如，35°C在重庆可能是常态，但在哈尔滨就是极端事件。

2.2 计算气候基准与温度百分位

要使用相对阈值法，首先需要建立气候基准。我推荐使用1971-2000年或1981-2010年作为基准期，这两个时段是WMO推荐的标准气候基准期。

import xarray as xr
import numpy as np
import pandas as pd
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

class HeatwaveDetector:
    """热浪检测器类"""
    
    def __init__(self, data_dir, baseline_start=1971, baseline_end=2000):
        """
        初始化热浪检测器
        
        参数:
        data_dir: ERA5数据目录
        baseline_start: 基准期起始年份
        baseline_end: 基准期结束年份
        """
        self.data_dir = Path(data_dir)
        self.baseline_start = baseline_start
        self.baseline_end = baseline_end
        self.baseline_data = None
        self.thresholds = None
        
    def load_baseline_data(self):
        """加载基准期数据"""
        print("正在加载基准期数据...")
        
        baseline_files = []
        for year in range(self.baseline_start, self.baseline_end + 1):
            for month in [6, 7, 8]:  # 夏季月份
                filename = f"era5_t2m_{year}_{month:02d}_长三角.nc"
                filepath = self.data_dir / filename
                if filepath.exists():
                    baseline_files.append(filepath)
        
        if not baseline_files:
            raise FileNotFoundError("未找到基准期数据文件")
        
        # 使用xarray的open_mfdataset高效读取多个文件
        ds_list = []
        for file in baseline_files:
            ds = xr.open_dataset(file)
            # 提取日最高温
            daily_max = ds['t2m'].resample(time='1D').max() - 273.15
            ds_list.append(daily_max)
        
        # 合并所有基准期数据
        self.baseline_data = xr.concat(ds_list, dim='time')
        print(f"基准期数据加载完成: {len(baseline_files)}个月份文件")
        return self.baseline_data
    
    def calculate_percentile_thresholds(self, percentile=90):
        """
        计算温度百分位阈值
        
        参数:
        percentile: 百分位数，默认90（即第90百分位）
        
        返回:
        每个网格点的阈值数据数组
        """
        if self.baseline_data is None:
            self.load_baseline_data()
        
        print(f"正在计算第{percentile}百分位阈值...")
        
        # 对每个网格点单独计算百分位
        # 使用xarray的apply_ufunc进行并行计算
        def _percentile_func(data, q):
            """计算单个网格点的百分位"""
            return np.nanpercentile(data, q)
        
        self.thresholds = xr.apply_ufunc(
            _percentile_func,
            self.baseline_data,
            input_core_dims=[['time']],
            output_core_dims=[[]],
            vectorize=True,
            kwargs={'q': percentile}
        )
        
        print("阈值计算完成")
        return self.thresholds
    
    def detect_heatwave_days(self, target_year, min_duration=3):
        """
        检测指定年份的热浪日
        
        参数:
        target_year: 目标年份
        min_duration: 最小持续天数，默认3天
        
        返回:
        热浪日标记的xarray数据集
        """
        if self.thresholds is None:
            self.calculate_percentile_thresholds()
        
        print(f"正在检测{target_year}年热浪日...")
        
        # 加载目标年份数据
        target_data = []
        for month in [6, 7, 8]:
            filename = f"era5_t2m_{target_year}_{month:02d}_长三角.nc"
            filepath = self.data_dir / filename
            if not filepath.exists():
                raise FileNotFoundError(f"目标年份数据缺失: {filename}")
            
            ds = xr.open_dataset(filepath)
            daily_max = ds['t2m'].resample(time='1D