我在特定年份每天都有来自ERA5的每小时数据。我想将这些数据从每小时转换为每日。我知道做这件事的漫长而艰难的方法,但我需要一些容易做到的事情。
哥白尼在这里有一个代码https://confluence.ecmwf.int/display/CKB/ERA5%3A+How+to+calculate+daily+total+precipitation,如果数据集只转换了一天,它可以正常工作,但是当转换为全年时,我遇到了问题。
链接下载ERA5数据集,该数据集可在https://cds.climate.copernicus.eu/cdsapp#!/home获得
https://confluence.ecmwf.int/display/CKB/How+to+download+ERA5
This script downloads the houly data for only 2 days (1st and 2nd of January 2017):#!/usr/bin/env python
"""
Save as get-tp.py, then run "python get-tp.py".
Input file : None
Output file: tp_20170101-20170102.nc
"""
import cdsapi
c = cdsapi.Client()
r = c.retrieve(
'reanalysis-era5-single-levels', {
'variable' : 'total_precipitation',
'product_type': 'reanalysis',
'year' : '2017',
'month' : '01',
'day' : ['01', '02'],
'time' : [
'00:00','01:00','02:00',
'03:00','04:00','05:00',
'06:00','07:00','08:00',
'09:00','10:00','11:00',
'12:00','13:00','14:00',
'15:00','16:00','17:00',
'18:00','19:00','20:00',
'21:00','22:00','23:00'
],
'format' : 'netcdf'
})
r.download('tp_20170101-20170102.nc')
## Add multiple days and multiple months to donload more data
Below script will create a netCDF file for only one day
#!/usr/bin/env python
"""
Save as file calculate-daily-tp.py and run "python calculate-daily-tp.py".
Input file : tp_20170101-20170102.nc
Output file: daily-tp_20170101.nc
"""
import time, sys
from datetime import datetime, timedelta
from netCDF4 import Dataset, date2num, num2date
import numpy as np
day = 20170101
d = datetime.strptime(str(day), '%Y%m%d')
f_in = 'tp_%d-%s.nc' % (day, (d + timedelta(days = 1)).strftime('%Y%m%d'))
f_out = 'daily-tp_%d.nc' % day
time_needed = []
for i in range(1, 25):
time_needed.append(d + timedelta(hours = i))
with Dataset(f_in) as ds_src:
var_time = ds_src.variables['time']
time_avail = num2date(var_time[:], var_time.units,
calendar = var_time.calendar)
indices = []
for tm in time_needed:
a = np.where(time_avail == tm)[0]
if len(a) == 0:
sys.stderr.write('Error: precipitation data is missing/incomplete - %s!\n'
% tm.strftime('%Y%m%d %H:%M:%S'))
sys.exit(200)
else:
print('Found %s' % tm.strftime('%Y%m%d %H:%M:%S'))
indices.append(a[0])
var_tp = ds_src.variables['tp']
tp_values_set = False
for idx in indices:
if not tp_values_set:
data = var_tp[idx, :, :]
tp_values_set = True
else:
data += var_tp[idx, :, :]
with Dataset(f_out, mode = 'w', format = 'NETCDF3_64BIT_OFFSET') as ds_dest:
# Dimensions
for name in ['latitude', 'longitude']:
dim_src = ds_src.dimensions[name]
ds_dest.createDimension(name, dim_src.size)
var_src = ds_src.variables[name]
var_dest = ds_dest.createVariable(name, var_src.datatype, (name,))
var_dest[:] = var_src[:]
var_dest.setncattr('units', var_src.units)
var_dest.setncattr('long_name', var_src.long_name)
ds_dest.createDimension('time', None)
var = ds_dest.createVariable('time', np.int32, ('time',))
time_units = 'hours since 1900-01-01 00:00:00'
time_cal = 'gregorian'
var[:] = date2num([d], units = time_units, calendar = time_cal)
var.setncattr('units', time_units)
var.setncattr('long_name', 'time')
var.setncattr('calendar', time_cal)
# Variables
var = ds_dest.createVariable(var_tp.name, np.double, var_tp.dimensions)
var[0, :, :] = data
var.setncattr('units', var_tp.units)
var.setncattr('long_name', var_tp.long_name)
# Attributes
ds_dest.setncattr('Conventions', 'CF-1.6')
ds_dest.setncattr('history', '%s %s'
% (datetime.now().strftime('%Y-%m-%d %H:%M:%S'),
' '.join(time.tzname)))
print('Done! Daily total precipitation saved in %s' % f_out)
结果应该是全年计算变量的每日值(例如降水等)。
示例:假设我的全年降水量数据为每天1毫米/小时,全年将有2928个值。
我想要的是全年24毫米/天,非闰年只有365个值。
示例输入数据集:数据的子集可以从这里下载(2017年1月1日和2日)https://www.dropbox.com/sh/0vdfn20p355st3i/AABKYO4do_raGHC34VnsXGPqa?dl=0。在此之后使用第二个脚本来检查代码。 {全年代码> 10 GB因此无法上传
提前致谢
xarray resample只是你的工具。它将netCDF数据从一个时间分辨率(例如每小时)转换为另一个(例如每天)一行。使用您的示例数据文件,我们可以使用以下代码创建daily-means:
import xarray as xr
ds = xr.open_dataset('./tp_20170101-20170102.nc')
tp = ds['tp'] # dimensions [time: 48, latitude: 721, longitude: 1440]
tp_daily = tp.resample(time='D').mean(dim='time') # dimensions (time: 2, latitude: 721, longitude: 1440)
你会看到resample
命令接受一个时间码,在这种情况下'D'
意味着每天,然后我们指定我们想用.mean(dim='time')
使用那天的每小时数据计算每一天的平均值。
相反,例如,如果你想计算每日最大值而不是每日均值,你可以用.mean(dim='time')
替换.max(dim='time')
。您还可以从小时到月(MS
或月开始),年度(AS
或年度开始)等等。时间频率代码可以在Pandas docs中找到。