下图表征2022年夏季（6月-8月）长三角地区（YRD）每天的 2 米气温异常值（°C）与受影响面积之间的关系。

图形出自论文：J2024-Impacts of urban expansion on air temperature and humidity during 2022 mega-heatwave over the Yangtze River Delta, China-Science of The Total Environment

图的轴说明：

• X轴（横轴）：时间（从6月1日到8月31日）。
• Y轴（纵轴）：气温异常（Temperature Anomaly, 单位为°C），即某地当天的气温与历史气候平均值之差。
• 颜色（色块）：代表该温度异常值对应的受影响面积（Affected Area），单位为平方公里（km²）。颜色越深，面积越大。
• 右侧色标：颜色代表不同的受影响面积（×10⁴ km²）。

本博客复现上述图形，详细介绍ERA5数据的下载、热浪事件的识别等多个步骤，并提供相应Python代码。

Step1：下载多年夏季气温/湿度数据-ERA5

ERA5数据下载步骤

ERA5数据的详细介绍可参考另一博客-【数据集】ERA5（欧洲中期天气预报中心）再分析数据介绍及下载。

• 数据集：ERA5 hourly data on single levels from 1959 to present
• 变量：2m air temperature（2m_temperature）

ERA5-Climate Data Store-查看个人API等进行配置。

API的配置教程-CDSAPI setup

基于Python下载ERA5数据

下载1970年至2024年6、7、8月的气温（2m_temperature）和湿度（2m_specific_humidity）数据，范围包括每个月的所有天数。

Areaf范围说明：北边界 (North), 西边界 (West), 南边界 (South), 东边界 (East)

完整Python代码如下：

import cdsapi
import calendar

c = cdsapi.Client()

years = list(range(1970, 2025))
months = ['06', '07', '08']

for year in years:
    for month in months:
        # 获取该月所有日期
        days_in_month = calendar.monthrange(year, int(month))[1]
        days = [f"{day:02d}" for day in range(1, days_in_month + 1)]
        hours = [f"{h:02d}:00" for h in range(24)]

        # 下载 2m 气温
        c.retrieve(
            'reanalysis-era5-single-levels',
            {
                'product_type': 'reanalysis',
                'variable': '2m_temperature',
                'year': str(year),
                'month': [month],
                'day': days,
                'time': hours,
                'area': [30.5, 106.0, 29.0, 107.5],  # 重庆主城区，N, W, S, E
                'format': 'netcdf',
            },
            f'chongqing_{year}_{month}_temperature.nc'
        )

基于Panoply查看.nc数据结构，如下：

读取气温数据+计算日均/月均气温并绘图（Python代码）

首先，对每天取24小时平均，得到日平均气温；此外，计算各月的月均温。

绘制的月均气温如下：（此案例中选取的区域较小，为7×7方块）

完整代码如下：

import os
import numpy as np
import pandas as pd
import xarray as xr
import rasterio
from rasterio.transform import from_origin
import geopandas as gpd
import matplotlib.pyplot as plt
from tqdm import tqdm

# 设置全局路径
basic_dir = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_nc_1970-2024'
data_dir = os.path.join(basic_dir, 'Ta_nc_1970-2024')

output_dir = os.path.join(basic_dir, 'Ta_1970-2024')
os.makedirs(output_dir, exist_ok=True)

figure_dir = os.path.join(basic_dir, 'Figures')
os.makedirs(figure_dir, exist_ok=True)

# 设置重庆边界路径
shp_path = r'D:\0 DataBase\0 Chongqin Database\1 Boundary\CQcity.shp'
gdf = gpd.read_file(shp_path)

# 设置字体
plt.rcParams['font.family'] = 'Times New Roman'


# 保存 tif 的函数
def save_tif(data_array, out_path, lon, lat):
    transform = from_origin(lon.min(), lat.max(), abs(lon[1] - lon[0]), abs(lat[1] - lat[0]))
    with rasterio.open(
            out_path,
            'w',
            driver='GTiff',
            height=data_array.shape[1],
            width=data_array.shape[2],
            count=data_array.shape[0],
            dtype='float32',
            crs='EPSG:4326',
            transform=transform,
    ) as dst:
        for i in range(data_array.shape[0]):
            dst.write(data_array[i, :, :].values.astype('float32'), i + 1)


# 绘图函数
def plot_map(data, title, save_path, lon, lat, cmap='hot_r', vmin=None, vmax=None):
    fig, ax = plt.subplots(figsize=(10, 8))
    im = ax.imshow(
        np.flipud(data),
        cmap=cmap,
        extent=[lon.min(), lon.max(), lat.min(), lat.max()],
        vmin=vmin,
        vmax=vmax,
        origin='lower'
    )
    gdf.boundary.plot(ax=ax, edgecolor='black', linewidth=1)

    cbar = plt.colorbar(im, ax=ax, shrink=0.95, fraction=0.036, pad=0.04)
    cbar.set_label('Temperature (°C)', fontsize=12)
    cbar.ax.tick_params(labelsize=12)

    ax.set_title(title, fontsize=16)
    ax.set_xlabel('Longitude(°)', fontsize=14)
    ax.set_ylabel('Latitude(°)', fontsize=14)
    ax.tick_params(axis='both', labelsize=12)

    plt.savefig(save_path, dpi=300)
    plt.close()


# 主处理函数
def process_temperature_file(year, month):
    month_str = f"{int(month):02d}"
    filename = f'chongqing_{year}_{month_str}_temperature.nc'
    filepath = os.path.join(data_dir, filename)

    if not os.path.exists(filepath):
        print(f"文件不存在：{filename}")
        return

    # 输出文件路径
    tmax_path = os.path.join(output_dir, f'CQcity_Tmax_{year}{month_str}.tif')
    tmin_path = os.path.join(output_dir, f'CQcity_Tmin_{year}{month_str}.tif')
    tmean_path = os.path.join(output_dir, f'CQcity_Tmean_{year}{month_str}_mean.tif')

    # 如果所有文件都已存在，跳过
    if os.path.exists(tmax_path) and os.path.exists(tmin_path) and os.path.exists(tmean_path):
        print(f"{year}-{month_str} 已处理，跳过。")
        return

    ds = xr.open_dataset(filepath)

    # 正确选择变量
    t2m = ds['t2m'] - 273.15  # 转为摄氏度
    time = pd.to_datetime(ds['valid_time'].values)
    lat = ds['latitude'].values
    lon = ds['longitude'].values

    days = pd.to_datetime(np.unique(time.date))

    # 使用 groupby 将 t2m 按日期分组
    grouped = t2m.groupby('valid_time.date')

    tmax_list = []
    tmin_list = []

    for date, day_data in grouped:
        tmax = day_data.max(dim='valid_time')
        tmin = day_data.min(dim='valid_time')
        tmax_list.append(tmax)
        tmin_list.append(tmin)

    # 合并为 DataArray
    tmax_all = xr.concat(tmax_list, dim='time')
    tmin_all = xr.concat(tmin_list, dim='time')
    tmax_all['time'] = pd.to_datetime([pd.Timestamp(str(d)) for d in grouped.groups.keys()])
    tmin_all['time'] = pd.to_datetime([pd.Timestamp(str(d)) for d in grouped.groups.keys()])

    # 合并为 3D 数据
    tmax_all = xr.concat(tmax_all, dim='time')
    tmin_all = xr.concat(tmin_all, dim='time')
    tmean = t2m.mean(dim='valid_time')

    print("Tmax shape:", tmax_all.shape)
    print("Tmin shape:", tmin_all.shape)

    # 保存 tif
    save_tif(tmax_all, tmax_path, lon, lat)
    save_tif(tmin_all, tmin_path, lon, lat)

    with rasterio.open(
            tmean_path,
            'w',
            driver='GTiff',
            height=tmean.shape[0],
            width=tmean.shape[1],
            count=1,
            dtype='float32',
            crs='EPSG:4326',
            transform=from_origin(lon.min(), lat.max(), abs(lon[1] - lon[0]), abs(lat[1] - lat[0])),
    ) as dst:
        dst.write(tmean.values.astype('float32'), 1)

    # 绘图（每日）
    for i, day in enumerate(days):
        day_str = day.strftime('%Y-%m-%d')
        # plot_map(tmax_all[i].values, f'Daily Tmax on {day_str}', os.path.join(figure_dir, f'Tmax_{day_str}.png'), lon, lat)
        # plot_map(tmin_all[i].values, f'Daily Tmin on {day_str}', os.path.join(figure_dir, f'Tmin_{day_str}.png'), lon, lat)

    # 绘制月均温
    plot_map(tmean.values, f'Monthly Mean Temperature {year}-{month_str}',
             os.path.join(figure_dir, f'Tmean_{year}-{month_str}.png'), lon, lat)

    ds.close()
    print(f"{year}-{month_str} 处理完成 ✅")


# ---------------------------
# 设置起始和结束年份与月份
start_year, start_month = 1970, 6
end_year, end_month = 1976, 8

# 自动调用处理函数
for year in range(start_year, end_year + 1):
    for month in range(6, 9):  # 6, 7, 8月
        if year == start_year and month < start_month:
            continue
        if year == end_year and month > end_month:
            break
        process_temperature_file(year, month)

Step2：计算每日气温异常（Temperature Anomaly）

计算原理

气候基准期的多年平均气温（例如1970–2024年），用于计算“气温异常”：

Temperature Anomaly(x,y,t) = T2022(x,y,t) − Tclimate(x,y,t)

两个版本的气温异常（Temperature Anomaly）分析：

方法编号	方法描述
1️⃣	逐日序号法：对每年同一个日序（如6月15日）在1970–2024年间的值求平均，作为该日气候均值
2️⃣	滑动窗口法：对该日及其前后各7天的共15天数据求平均（跨月、跨年闭环处理），作为该日气候均值

完整的Python代码如下：

import os
import numpy as np
import rasterio
import datetime
import pandas as pd
from glob import glob
from tqdm import tqdm
import matplotlib.pyplot as plt
import calendar

"""
代码说明：
1、根据计算的日最高/日最低气温，计算日均气温（最大值/最小值）
方法1：按年序号平均
方法2：15日滑动平均（闭环）
2、计算气温异常：可指定时间段
"""

# 设置字体
plt.rcParams['font.family'] = 'Times New Roman'

# 加载 tif 为 numpy 数组
def load_tif(path):
    with rasterio.open(path) as src:
        data = src.read()
        profile = src.profile
    return data, profile

# 保存 tif
def save_tif(path, data, profile):
    profile.update(count=data.shape[0], dtype='float32')
    with rasterio.open(path, 'w', **profile) as dst:
        dst.write(data.astype('float32'))


### 方法一：按年序号平均
def compute_temperature_anomaly_dayofyear(tif_dir, output_dir, year_target=2022):
    os.makedirs(output_dir, exist_ok=True)

    # 构建年份和月日映射
    years = list(range(1970, 2025))
    date_range = pd.date_range(f'{year_target}-06-01', f'{year_target}-08-31', freq='D')

    for date in tqdm(date_range):
        day_index = f"{date.month:02d}{date.day:02d}"  # MMDD
        all_years_data = []
        profile = None

        for year in years:
            tif_path = os.path.join(tif_dir, f'CQcity_Tmax_{year}{date.month:02d}.tif')
            if not os.path.exists(tif_path):
                continue
            data, profile = load_tif(tif_path)
            band_index = date.day - 1
            if band_index < data.shape[0]:
                all_years_data.append(data[band_index])

        if len(all_years_data) == 0:
            continue

        all_years_data = np.stack(all_years_data, axis=0)
        clim_mean = np.nanmean(all_years_data, axis=0)

        # 提取目标年当天数据
        current_path = os.path.join(tif_dir, f'CQcity_Tmax_{year_target}{date.month:02d}.tif')
        current_data, _ = load_tif(current_path)
        today_data = current_data[date.day - 1]

        anomaly = today_data - clim_mean

        out_path = os.path.join(output_dir, f'CQcity_Tmax_anomaly_{date.strftime("%Y%m%d")}.tif')
        save_tif(out_path, anomaly[np.newaxis, :, :], profile)


### 方法二：15日滑动平均（闭环）
def compute_temperature_anomaly_rolling15(tif_dir, output_dir, year_target=2022):
    os.makedirs(output_dir, exist_ok=True)

    years = list(range(1970, 2025))
    date_range = pd.date_range(f'{year_target}-06-01', f'{year_target}-08-31', freq='D')
    all_dates = pd.date_range('06-01', '08-31', freq='D')

    for date in tqdm(date_range):
        # 构建 15 天窗口
        day_idx = (date - pd.Timestamp(f'{year_target}-06-01')).days
        window_indices = [(day_idx + i) % len(all_dates) for i in range(-7, 8)]
        window_dates = [all_dates[i] for i in window_indices]

        all_window_data = []
        profile = None

        # 每年每个窗口日
        for year in years:
            for d in window_dates:
                month = d.month
                day = d.day
                tif_path = os.path.join(tif_dir, f'CQcity_Tmax_{year}{month:02d}.tif')
                if not os.path.exists(tif_path):
                    continue
                data, profile = load_tif(tif_path)
                if day - 1 < data.shape[0]:
                    all_window_data.append(data[day - 1])

        if len(all_window_data) == 0:
            continue

        all_window_data = np.stack(all_window_data, axis=0)
        clim_mean = np.nanmean(all_window_data, axis=0)

        # 目标年当天数据
        current_path = os.path.join(tif_dir, f'CQcity_Tmax_{year_target}{date.month:02d}.tif')
        current_data, _ = load_tif(current_path)
        today_data = current_data[date.day - 1]

        anomaly = today_data - clim_mean

        out_path = os.path.join(output_dir, f'CQcity_Tmax_anomaly_{date.strftime("%Y%m%d")}.tif')
        save_tif(out_path, anomaly[np.newaxis, :, :], profile)


### 主函数
# --------------------------------------------------------
# Tmax 输入输出路径
tmax_dir = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_1970-2024\Tmax'
tmax_out_dayofyear = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_Anomaly\Tmax\Method1_DayOfYear'
tmax_out_rolling = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_Anomaly\Tmax\Method2_Rolling15'

# Tmin 同理
tmin_dir = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_1970-2024\Tmin'
tmin_out_dayofyear = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_Anomaly\Tmin\Method1_DayOfYear'
tmin_out_rolling = r'D:\0 DataBase\7 EAR5\CQ_Summer\Ta_Anomaly\Tmin\Method2_Rolling15'

# 调用方法一
compute_temperature_anomaly_dayofyear(tmax_dir, tmax_out_dayofyear, year_target=2022)
compute_temperature_anomaly_dayofyear(tmin_dir, tmin_out_dayofyear, year_target=2022)

# 调用方法二
compute_temperature_anomaly_rolling15(tmax_dir, tmax_out_rolling, year_target=2022)
compute_temperature_anomaly_rolling15(tmin_dir, tmin_out_rolling, year_target=2022)

Step3：识别热浪事件

热浪计算原理

常见热浪识别准则（任选其一）：

• 百分位法：某地日最高温连续 ≥ 3天超过90/95百分位。
• 绝对阈值法：如重庆主城区连续3天 Tmax ≥ 35°C。

可以先标记出热浪日，再聚合其空间影响范围。

完整的Python代码如下：

Step4：计算每一天的“受影响面积”

• 对每一天的温度异常字段，使用阈值（如1°C, 2°C…）筛选出“异常升温”的格点。
• 每个格点面积可根据经纬度估算（例如ERA5分辨率为0.25°≈25×25 km²）。
• 统计各温度异常等级下的面积总和，绘制等值线。

完整的Python代码如下：

参考