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

WW、forever

已于 2024-10-31 10:14:56 修改

阅读量6.7k

点赞数 6

分类专栏： Python基础处理文章标签： python 开发语言

于 2023-11-04 20:32:59 首次发布

本文链接：https://blog.csdn.net/qq_44246618/article/details/134215773

版权

Python基础处理专栏收录该内容

40 篇文章

订阅专栏

基于Python处理ERA5数据

1 ERA5数据简介
2 数据集处理
参考

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