ERA5数据获取和处理

最新推荐文章于 2024-05-11 22:09:02 发布

取间

最新推荐文章于 2024-05-11 22:09:02 发布

阅读量5.3k

点赞数 6

文章标签：绘图

本文链接：https://blog.csdn.net/qq_40497403/article/details/112509547

版权

ERA5数据的获取地址：https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview
这时其中ERA5的其中一个产品，我们也可以在主页https://cds.climate.copernicus.eu/cdsapp#!/home上进行搜索，然后选择自己想要的产品。
ERA5数据的获取步骤：
我们通过python代码进行批量处理，第一步先打开Internet Download Manager软件，第二步在python代码上进行相应修改，主要是修改时间，想获取的变量类型，数据的经纬度范围等。第三步：运行python代码即可，这时候就会在Internet Download Manager上下载相应的文件，这样会比在python上直接下载快一些。下载所用的代码如下：


import cdsapi
import calendar
from subprocess import call

def idmDownloader(task_url, folder_path, file_name):
    """
    IDM下载器
    :param task_url: 下载任务地址
    :param folder_path: 存放文件夹
    :param file_name: 文件名
    :return:
    """
    # IDM安装目录
    idm_engine = "C:/Program Files (x86)/Internet Download Manager/IDMan.exe"
    # 将任务添加至队列
    call([idm_engine, '/d', task_url, '/p', folder_path, '/f', file_name, '/a'])
    # 开始任务队列
    call([idm_engine, '/s'])

if __name__ == '__main__':
    c = cdsapi.Client(url="https://cds.climate.copernicus.eu/api/v2",
                  key = "73082:585e6691-6763-4ab3-bcd1-c5347276f09f")
    # 数据信息字典
    dic = {
        'product_type': 'reanalysis',  # 产品类型
        'format': 'netcdf',  # 数据格式
        'variable': [
            '10m_u_component_of_wind', '10m_v_component_of_wind',
        ],  #下载的变量
        'year': '',  # 年，设为空
        'month': '',  # 月，设为空
        'day': [],  # 日，设为空
        'time': [
            '00:00', '06:00', '12:00',
            '18:00',
        ],
        'area': [
            30, -150, -70,
            -130,
        ],   #下载数据的区域
    }

    # 通过循环批量下载1979年到2020年所有月份数据
    for y in range(2016,2017):  # 遍历年
        for m in range(1,2):  # 遍历月
            day_num = calendar.monthrange(y, m)[1]  # 根据年月，获取当月日数
            # 将年、月、日更新至字典中
            dic['year'] = str(y)
            dic['month'] = str(m).zfill(2)
            dic['day'] = [str(d).zfill(2) for d in range(1, day_num + 1)]

            r = c.retrieve('reanalysis-era5-single-levels', dic, )  # 文件下载器
            url = r.location  # 获取文件下载地址
            path = 'E:/ERA5data'  # 存放文件夹
            filename = str(y) + str(m).zfill(2) + '.nc'  # 文件名
            idmDownloader(url, path, filename)  # 执行这个函数，添加进IDM中下载

其中与服务器建立连接所用的url和key是在主页（https://cds.climate.copernicus.eu/cdsapp#!/home）中点击“Climate Data Store API”可以获取。

ERA5数据的解析

from netCDF4 import Dataset
import numpy as np
import os

path = "D:/ERAinterim/2020.nc"
dst = Dataset(path, mode='r', format="netCDF4")

# 查看nc文件中包含了什么
print(dst)
print('---------------------------------------------------------')
# 查看nc文件有哪些变量
print(dst.variables.keys())
print('--------------------------------------------------------')
# 查看nc文件中变量的属性名称
print(dst.variables.keys())
for i in dst.variables.keys():
    print('%s: %s' % (i, dst.variables[i].ncattrs()))
print('--------------------------------------------------------')
# 查看nc文件中变量的属性的具体信息
print(dst.variables.keys())
for i in dst.variables.keys():
    print(dst.variables[i])
    print('\n')
print('-------------------------------------------------------')

# 获取维度的sizes信息，获得索引范围
for i in dst.dimensions.keys():
    print('%s_sizes: %s' % (i, dst.dimensions[i].size))
print('------------------------------------------------------')
# 循环提取数据
# for i in dst.variables.keys():
#     if i not in dst.dimensions.keys():

通过这个代码，我们可以知道我们下载的nc文件的结构，变量有10m u-component of wind（u10）, 10m v-component of wind（v10）, 10m wind speed（si10）, Mean direction of total swell(mdts), Mean direction of wind waves(mdww), Mean period of total swell(mpts), Mean period of wind waves(mpww), Mean wave direction(mwd), Mean wave period(mwp), Significant height of combined wind waves and swell(swh), Significant height of total swell(shts), Significant height of wind waves(shww)
并且每个变量都有四个维度时间（492）expver（2） 纬度（261）*经度（325）

接下来我们将数据写入到csv文件中方便我们浏览

import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
import csv

plt.rcParams['font.sans-serif'] = ['SimHei']        # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False            # 用来正常显示负号

def to_csv(source_file):
    f = nc.Dataset(source_file)
    all_vars = f.variables.keys()   #获取所有变量名称
    print(all_vars)  #长度为18
 
        # 读取nc数据
    lon=f['longitude'][:]  
    lat=f['latitude'][:]
    time=f['time'][:]  
    expver=f['expver'][:]#上面三个变量都是一维变量,都是维度
    u10=f['u10'][:]
    v10=f['v10'][:]
    si10=f['si10'][:]
    mdts=f['mdts'][:]
    mdww=f['mdww'][:]
    mpts=f['mpts'][:]
    mpww=f['mpww'][:]
    mwd=f['mwd'][:]
    mwp=f['mwp'][:]
    swh=f['swh'][:]
    shts=f['shts'][:]
    shww=f['shww'][:]

        lon=np.array(lon)
    lat=np.array(lat)
    time=np.array(time)
    expver=np.array(expver)
    u10=np.array(u10)
    v10=np.array(v10)
    si10=np.array(si10)
    mdts=np.array(mdts)
    mdww=np.array(mdww)
    mpts=np.array(mpts)
    mpww=np.array(mpww)
    mwd=np.array(mwd)
    mwp=np.array(mwp)
    swh=np.array(swh)
    shts=np.array(shts)
    shww=np.array(shww)
 
        print('-------------------------------------------------------------------')
    # 文件名不含扩展名
    source_file = source_file.split('.')
    file_name = source_file[0]
    print(file_name)   

            # 创建csv目标文件
    try:
        # 打开目标csv文件
        with open(file_name+'.csv', 'a', newline='') as targetFile:
            # 创建写入流
            writer = csv.writer(targetFile)
            # 写入表头
            writer.writerow(('lon', 'lat','expver','time', 'u10', 'v10','si10','mdts','mdww','mpts','mpww','mwd','mwp','swh','shts','shww'))
            # 写入数据
            for i in range(len(lon)):
                for j in range(len(lat)):
                    for m in range(len(expver)):
                       for k in range(len(time)):           
                                     writer.writerow((lon[i], lat[j],expver[m],time[k], u10[k][m][j][i], v10[k][m][j][i],u10[k][m][j][i],si10[k][m][j][i],mdts[k][m][j][i],mdww[k][m][j][i],mpts[k][m][j][i],mpww[k][m][j][i],mwd[k][m][j][i],mwp[k][m][j][i],swh[k][m][j][i],shts[k][m][j][i],shww[k][m][j][i]))#dates.num2date(time[i], tz=None)
        targetFile.close()              # 关闭文件
        print('Get'+file_name+'.csv Success!')
    except Exception as e:              # 抛异常
        print('Error :'+str(e))

if '__name__ ==__main__':
    print("start transformation...")
    to_csv('D:/ERAinterim/2020.nc')
    print('Transformed successfully')

我们此时读取csv文件可以发现，里面的数据有很多是缺失的，因为读取区域内有一部分是陆地，所以我们要删掉这部分数据。

import pandas as pd

data = pd.read_csv('D:/ERAinterim/regin/2020.csv')
data=data[~data['mwp'].isin([-32767.000000])]  #通过~取反，选取不包含数字1的行
data=data[~data['swh'].isin([-32767.000000])]  #我们将swh和mwp中缺失数据的陆地数据删除
data=data.reset_index()
print(data)
data.to_csv("D:/ERAinterim/regin/2020.csv")

接下来我们通过文件中的u10和v10计算出风向，并将这个风向当做新的一列加到2020.csv这个文件中去。

import pandas as pd
import math

def fengxiang (u,v):  #这个函数是利用u，v分量来计算风的风向
    if u>0:
        return math.degrees(math.atan2(v,u))
    elif u<0 and v>=0:
        return math.degrees(math.atan2(v,u))+180
    elif u<0 and v<0:
        return math.degrees(math.atan2(v,u))-180
    elif v>0 and u==0:
        return 90
    elif v<0 and u==0:
        return -90
    else:
        return 0

data = pd.read_csv('D:/ERAinterim/2020.csv')
#此时我们要通过u10和v10这两列建立新的一列表示风向
u10=data['u10']
v10=data['v10']
WindDir=[]
for i in range(len(u10)):
    WindDir.append(fengxiang(u10[i],v10[i]))
data['WindDir']=WindDir
data.to_csv('D:/ERAinterim/2020.csv',index=0)

接着是求分析数据文件

import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
import csv
import pandas as pd
import math

#分析西北太平洋的盛行风速风向，大风频率，极值风速
#（2）缝隙风浪，涌浪，混合浪的盛行波高波向，波长，波周期
#盛行风向也称最多风向，是一个地区在某一段时间内出现频率最多的风向。
#通常按日、月、季和年的时段运用统计学的办法求出相应时段的盛行风向。

#我们首先要知道风速为si10,风向为u10，v10. 
data = pd.read_csv('D:/ERAinterim/2020.csv')
MaxWindSpeed=max(data['si10'])  #得到极值风速
ShengxingWindSpeed=data['si10'].value_counts().index[1]  #此时得到除了缺失值以外最多的数
ShengxingWindDir=data['WindDir'].value_counts().index[1]

#我们认为风速大于10的为大风，并求出大风频率
WindSpeed=data['si10'].value_counts().drop([-32767.000000])  #我们不考虑缺失值
zong=sum(WindSpeed.values)
Windnum=list(WindSpeed.values)
WindValue=np.array(WindSpeed.index)
Windpin=np.array(Windnum/zong)
print(Windpin[WindValue>10])  #这里我们得到风速大于10的风的频率

#接下来求风浪，涌浪，混合朗的盛行波高，波向，波周期
#风浪的波高为 shww ,波向为 mdww， 波周期为 mpww
#涌浪的波高为 shts ，波向为mdts，波周期为mpts
#混合浪的波高为 swh ，波向为mwd，波周期为mwp
WindWaveH=data['shww'].value_counts().index[1]
WindWaveD=data['mdww'].value_counts().index[1]
WindWaveP=data['mpww'].value_counts().index[1]

SurWaveH=data['shts'].value_counts().index[1]
SurWaveD=data['mdts'].value_counts().index[1]
SurWaveP=data['mpts'].value_counts().index[1]

MixWaveH=data['swh'].value_counts().index[1]
MixWaveD=data['mwd'].value_counts().index[1]
MixWaveP=data['mwp'].value_counts().index[1]

接下来的问题就是怎么处理40年的数据文件，我们将最近40年的数据每5年放在一个文件中，分别是1.nc,2.nc…直到8.nc。我们将其先转换成1.csv等等。转换完成后再得到各个文件的风向列。
因为每一个文件都有5年的数据，所以我们以每5年一个单位获得盛行风速风向，大风频率，极值风速；风浪，涌浪，混合浪的盛行波高波向，波长，波周期并将这些数据存放到列表中

import matplotlib.pyplot as plt
import numpy as np
import netCDF4 as nc
import csv
import pandas as pd
import math

def getWind(num):
    data = pd.read_csv("D:/ERAinterim/%d.csv"%(num))
    MaxWindSpeed=max(data['si10'])  #得到极值风速
    ShengxingWindSpeed=data['si10'].value_counts().index[1]  #得到盛行风速
    ShengxingWindDir=data['WindDir'].value_counts().index[1] #得到盛行风向
    #我们认为风速大于10的为大风，并求出大风频率
    WindSpeed=data['si10'].value_counts()  
    zong=sum(WindSpeed.values)
    Windnum=list(WindSpeed.values)
    Windpin=np.array(Windnum/zong)
    BigWindpin=max(Windpin)
    return (MaxWindSpeed,ShengxingWindSpeed,ShengxingWindDir,BigWindpin)

def getWave(num):
    data = pd.read_csv("D:/ERAinterim/%d.csv"%(num))
    WindWaveH=data['shww'].value_counts().index[1]
    WindWaveD=data['mdww'].value_counts().index[1]
    WindWaveP=data['mpww'].value_counts().index[1]
    SurWaveH=data['shts'].value_counts().index[1]
    SurWaveD=data['mdts'].value_counts().index[1]
    SurWaveP=data['mpts'].value_counts().index[1]
    MixWaveH=data['swh'].value_counts().index[1]
    MixWaveD=data['mwd'].value_counts().index[1]
    MixWaveP=data['mwp'].value_counts().index[1]
    return (WindWaveH,WindWaveD,WindWaveP,SurWaveH,SurWaveD,SurWaveP,MixWaveH,MixWaveD,MixWaveP)

MaxWind=[] #将8个文件中的最大风速放到这个列表中
ShengWindS=[]
ShengWindD=[]
BigWindP=[]
for i in range(8):
    value1,value2,value3,value4=getWind(i+1)
    MaxWind.append(value1)
    ShengWindS.append(value2)
    ShengWindD.append(value3)
    BigWindP.append(value4)

print(MaxWind)
print(BigWindP) 

WWH=[]
WWD=[]
WWP=[]
SWH=[]
SWD=[]
SWP=[]
MWH=[]
MWD=[]
MWP=[]

for i in range(8):
    value1,value2,value3,value4,value5,value6,value7,value8,value9=getWave(i+1)
    WWH.append(value1)
    WWD.append(value2)
    WWP.append(value3)
    SWH.append(value4)
    SWD.append(value5)
    SWP.append(value6)
    MWH.append(value7)
    MWD.append(value8)
    MWP.append(value9)

print(WWH)

取间

关注

6
点赞
踩
78

收藏

觉得还不错? 一键收藏
5
评论
ERA5数据获取和处理

ERA5数据的获取地址：https://www.ecmwf.int/en/forecasts/datasets/browse-reanalysis-datasetsERA5数据的解析from netCDF4 import Datasetimport numpy as npimport ospath = "D:/ERAinterim/adaptor.mars.internal-1610440615.8431206-6767-1-48905033-3384-497a-948b-dc8b115cf.
复制链接

扫一扫