基于Python处理ERA5数据
1 ERA5数据简介
ERA5是ECMWF(欧洲中期天气预报中心)对1950年1月至今全球气候的第五代大气再分析数据集。
- 包含了四个基本变量(日平均温度、降水、比湿度和距离地表2米的气压),这些变量在每日时间尺度上覆盖全球,从而可以对不同地区和时间段进行全面和统一的分析
- 时间分辨率:1940年至今,小时尺度、日尺度、月尺度
- 空间分辨率:0.1°×0.1°(30km)
EAR5数据集的详细介绍及处理可参见另一博客-【数据集】ERA5(欧洲中期天气预报中心)再分析数据介绍及下载。
2 数据集处理
准备工作:xarray库安装
处理ERA5数据的一种常见方法是使用xarray库。
可使用pip list,在cmd控制台查看已安装包(库):
首先,确保已经安装了xarray和netCDF4库,以pip工具(cmd控制台)下载工具箱代码如下:
pip install xarray netCDF4
然后,可以使用xarray的open_dataset()函数加载ERA5数据集:
import xarray as xr
# 加载ERA5数据集
ds = xr.open_dataset('era5_data.nc')
接下来,可以使用xarray的各种功能来处理数据。例如,可以使用sel()函数从数据集中选择特定的经度和纬度:
# 选择经度为-60和纬度为30的数据
ds = ds.sel(longitude=-60, latitude=30)
还可以使用resample()函数对时间进行重新采样:
# 将时间重新采样为每月数据
ds = ds.resample(time='1M').mean()
最后,可以将数据保存到netCDF文件中:
# 将处理后的数据保存到netCDF文件中
ds.to_netcdf('processed_era5_data.nc')
可根据具体需求,使用xarray的其他功能来处理ERA5数据。
2.1 基于小时尺度数据获取日尺度数据-降水+温度
以2020年7月1日-31日24小时降水和温度数据为例,求和得到日尺度降水,求平均得到日尺度温度。
日尺度降水
绘制日尺度降水如下:
Python代码如下:
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import os # 导入 os 模块以处理文件路径
# 设置全局字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'
# 定义时间范围
start_date = pd.Timestamp('2020-07-01')
end_date = pd.Timestamp('2020-07-31')
# print(start_date)
# 替换为 EAR5 数据文件路径
# file_path = 'D:/0 DataBase/7 EAR5/GBA_2020_June/data_0.nc'
file_path = 'D:/0 DataBase/7 EAR5/GBA_2020_July/GBA_2020_July.nc'
# 打开 NetCDF 文件
dataset = nc.Dataset(file_path)
# 提取降水量数据(变量名为 'tp')
precipitation = dataset.variables['tp'][:] # 24小时降水量数据
# 提取纬度和经度
latitude = dataset.variables['latitude'][:]
longitude = dataset.variables['longitude'][:]
# 提取有效时间
time = dataset.variables['valid_time'][:]
time_units = dataset.variables['valid_time'].units # 获取时间单位
# 转换时间为日期
dates = nc.num2date(time, units=time_units)
# 将 cftime 对象转换为字符串
dates_str = [date.strftime('%Y-%m-%d %H:%M:%S') for date in dates]
# 使用 pandas 创建日期索引
dates = pd.to_datetime(dates_str)
# 创建 DataFrame,列为降水量,索引为时间
# 注意:将数据的维度调整为 (时间, 纬度 x 经度)
# 先将降水量数据重塑为 2D 的 (时间, 纬度*经度)
n_time = precipitation.shape[0]
n_lat = precipitation.shape[1]
n_lon = precipitation.shape[2]
# 重新调整数据形状为 (时间, 纬度*经度)
precipitation_reshaped = precipitation.reshape(n_time, n_lat * n_lon)
# 创建 DataFrame
df = pd.DataFrame(data=precipitation_reshaped, index=dates)
# 按日期求和,得到每日总降水量(单位为m,调整为mm)
daily_precipitation = df.resample('D').sum()
daily_precipitation = daily_precipitation * 100
# 自定义导出图片的目录
# output_directory = 'D:/6 Python Codes/GBA_Data_Process/Pre_June_2020/'
output_directory = 'D:/6 Python Codes/GBA_Data_Process/Pre_July_2020/'
os.makedirs(output_directory, exist_ok=True) # 创建目录(如果不存在)
# 循环绘制每一天的总降水量
for i, date in enumerate(daily_precipitation.index):
plt.figure(figsize=(10, 6))
# 获取当前日期的降水量数据
daily_data = daily_precipitation.iloc[i].values.reshape(n_lat, n_lon) # 重新调整为 2D
# 对数据进行逆向转置,使其符合经度(x轴)和纬度(y轴)的绘图要求
daily_data = np.flipud(daily_data)
# 绘制降水量图
plt.imshow(daily_data, cmap='Blues', origin='lower',
extent=[longitude.min(), longitude.max(), latitude.min(), latitude.max()]) # 设置经纬度范围
plt.colorbar(label='Precipitation (mm)')
# 加载并绘制 SHP 文件
GBAshp_file_path = "D:/6 Python Codes/GBA_Data_Process/shpFile/GBA/GBA.shp"
Chinashp_file_path = "D:/6 Python Codes/GBA_Data_Process/shpFile/China/中国.shp"
gdf_GBA = gpd.read_file(GBAshp_file_path)
gdf_China = gpd.read_file(Chinashp_file_path)
# 绘制 SHP 文件
gdf_GBA.boundary.plot(ax=plt.gca(), color='black', linewidth=0.5) # 将边界绘制在图上
gdf_China.boundary.plot(ax=plt.gca(), color='black', linewidth=1) # 将边界绘制在图上
# 设置图形范围与 TIFF 图像边界相同
plt.xlim(longitude.min(), longitude.max())
plt.ylim(latitude.min(), latitude.max())
# 设置标题
plt.title(f'Total Precipitation on {date.strftime("%Y-%m-%d")}', fontsize=16)
plt.xlabel('Longitude', fontsize=14)
plt.ylabel('Latitude', fontsize=14)
# 自定义坐标轴刻度值,后面加上°
def format_ticks(x, pos):
return f'{x:.1f}°'
# 应用自定义格式化函数
plt.gca().xaxis.set_major_formatter(plt.FuncFormatter(format_ticks))
plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(format_ticks))
# 保存图像到自定义目录
plt.savefig(os.path.join(output_directory, f'daily_precipitation_{date.strftime("%Y%m%d")}.png'), dpi=300)
# 显示图像
plt.close() # 关闭当前图形,避免显示后暂停
# 关闭数据集
dataset.close()
日尺度温度
绘制日尺度温度如下:
Python代码如下:
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import os # 导入 os 模块以处理文件路径
# 设置全局字体为 Times New Roman
plt.rcParams['font.family'] = 'Times New Roman'
# 定义时间范围
start_date = pd.Timestamp('2020-07-01')
end_date = pd.Timestamp('2020-07-31')
# 替换为 EAR5 数据文件路径
# file_path = 'D:/0 DataBase/7 EAR5/GBA_2020_June/data_0.nc'
file_path = 'D:/0 DataBase/7 EAR5/GBA_2020_July/GBA_2020_July.nc'
# 打开 NetCDF 文件
dataset = nc.Dataset(file_path)
# 提取 2米处的温度数据(变量名假设为 't2m')
temperature = dataset.variables['t2m'][:] # 2米高度的温度数据
# 提取纬度和经度
latitude = dataset.variables['latitude'][:]
longitude = dataset.variables['longitude'][:]
# 提取有效时间
time = dataset.variables['valid_time'][:]
time_units = dataset.variables['valid_time'].units # 获取时间单位
# 转换时间为日期
dates = nc.num2date(time, units=time_units)
# 将 cftime 对象转换为字符串
dates_str = [date.strftime('%Y-%m-%d %H:%M:%S') for date in dates]
# 使用 pandas 创建日期索引
dates = pd.to_datetime(dates_str)
# 创建 DataFrame,列为温度,索引为时间
n_time = temperature.shape[0]
n_lat = temperature.shape[1]
n_lon = temperature.shape[2]
# 重新调整数据形状为 (时间, 纬度Latitude*经度Longitude)
temperature_reshaped = temperature.reshape(n_time, n_lat * n_lon)
# 创建 DataFrame
df = pd.DataFrame(data=temperature_reshaped, index=dates)
# 按日期求平均,得到每日平均温度
daily_temperature = df.resample('D').mean()
# 将温度从开尔文转换为摄氏度
daily_temperature = daily_temperature - 273.15
# 自定义导出图片的目录
output_directory = 'D:/6 Python Codes/GBA_Data_Process/Tem_July_2020/'
os.makedirs(output_directory, exist_ok=True) # 创建目录(如果不存在)
# 循环绘制每一天的平均温度
for i, date in enumerate(daily_temperature.index):
plt.figure(figsize=(10, 6))
# 获取当前日期的温度数据,按照经度、纬度重新调整为 2D
daily_data = daily_temperature.iloc[i].values.reshape(n_lat, n_lon) # 先按经度再按纬度
# 对数据进行逆向转置,使其符合经度(x轴)和纬度(y轴)的绘图要求
daily_data = np.flipud(daily_data)
# 绘制温度图
plt.imshow(daily_data, cmap='hot_r', origin='lower',
extent=[longitude.min(), longitude.max(), latitude.min(), latitude.max()]) # 设置经纬度范围
plt.colorbar(label='Temperature (°C)')
# 加载并绘制 SHP 文件
GBAshp_file_path = "D:/6 Python Codes/GBA_Data_Process/shpFile/GBA/GBA.shp"
Chinashp_file_path = "D:/6 Python Codes/GBA_Data_Process/shpFile/China/中国.shp"
gdf_GBA = gpd.read_file(GBAshp_file_path)
gdf_China = gpd.read_file(Chinashp_file_path)
# 绘制 SHP 文件
gdf_GBA.boundary.plot(ax=plt.gca(), color='black', linewidth=0.5) # 将边界绘制在图上
gdf_China.boundary.plot(ax=plt.gca(), color='black', linewidth=1) # 将边界绘制在图上
# 设置图形范围与 TIFF 图像边界相同
plt.xlim(longitude.min(), longitude.max())
plt.ylim(latitude.min(), latitude.max())
# 设置标题
plt.title(f'Average Temperature on {date.strftime("%Y-%m-%d")}', fontsize=16)
plt.xlabel('Longitude') # 设置 x 轴为经度
plt.ylabel('Latitude') # 设置 y 轴为纬度
# 自定义坐标轴刻度值,后面加上°
def format_ticks(x, pos):
return f'{x:.1f}°'
# 应用自定义格式化函数
plt.gca().xaxis.set_major_formatter(plt.FuncFormatter(format_ticks))
plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(format_ticks))
# 保存图像到自定义目录
plt.savefig(os.path.join(output_directory, f'daily_temperature_{date.strftime("%Y%m%d")}.png'), dpi=300)
# 关闭当前图形
plt.close()
# 关闭数据集
dataset.close()