【数据集处理】基于Python处理EAR5数据

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()

参考

### ERA5 数据处理成日值的方法 对于ERA5数据处理,尤其是将其转换为日值的过程,可以采用多种方法和工具实现。这些工具包括Python、MATLAB 和 NCL。 #### 使用 Python 进行 ERA5 数据的日值计算 在Python中,`netCDF4`库能够方便地读取 NetCDF 文件中的ERA5数据[^1]。为了获取每日平均温度或其他变量的日值,通常会执行如下操作: ```python from netCDF4 import Dataset import numpy as np def daily_mean_temperature(nc_file_path, variable_name='t2m'): dataset = Dataset(nc_file_path) time_var = dataset.variables['time'] temp_data = dataset.variables[variable_name][:] # 假设每小时有记录,则一天内共有24个时间点 days_count = int(len(time_var)/24) reshaped_temp = temp_data.reshape(days_count, 24, *temp_data.shape[1:]) mean_daily_temps = np.mean(reshaped_temp, axis=1) return mean_daily_temps ``` 这段代码展示了如何通过重塑数组结构来获得每天的平均气温[^3]。 #### 利用 MATLAB 实现相同目标 同样,在MATLAB环境下,可以通过调用特定函数如 `ncread` 来加载NetCDF格式下的ERA5资料,并进一步加工得到所需的时间尺度上的统计数据[^2]。 ```matlab % 加载nc文件 filename = 'era5_example.nc'; varname = 'temperature'; % 或者其他感兴趣的物理量名称 data = ncread(filename,varname); % 计算日均温 (假设数据按小时存储) dailyMeans = mean(reshape(data,[],24,size(data,end)),2); ``` 此脚本片段实现了从原始逐时观测到日度统计指标转变的功能[^4]。 #### 应用 NCL 完成任务 NCL(NCAR Command Language)作为一种专用于科学数据分析的语言,提供了简洁而强大的功能去处理像ERA5这样的大型气候模型输出产品。下面给出了一种可能的方式来进行日汇总运算: ```ncl load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl" load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" begin file = addfile("your_era5_dataset.nc","r") var = file->variable_of_interest@values ; 将一维序列重新排列为二维矩阵,其中每一列代表一日内的不同时刻 dim0 = dimsizes(var)(0) / 24 reshape_var = conform_dims((/dim0,24,dimsizes(var)(1:)/),var,False) avg_dayly_values = dim_avg_wgt(reshape_var,True,(/True,False,False/)) end ``` 上述三种方式分别适用于不同的编程环境和个人偏好;然而无论选用哪种途径,核心思路都是相似的——即先按照日期分组再求得各组内部数值特征描述符(比如这里所说的“平均”)。值得注意的是实际应用过程中还需考虑更多细节因素,例如缺失值处理、质量控制等。
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

WW、forever

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值