This knowledge base article shows you how to calculate daily total precipitation using ERA-Interim data. If you just want monthly means, then you can simply download it from http://apps.ecmwf.int/datasets/data/interim-mdfa/levtype=sfc/.
Before you continue, make sure you read through knowledge base articles listed below:
You are also supposed to know how to work with Python under Linux, in particular, how to install packages using pip. You are recommended to use the latest release of packages listed here:
#!/usr/bin/env python
"""
Save as get-tp.py, then run "python get-tp.py".
Input file : None
Output file: tp_20180101.nc
"""
from ecmwfapi import ECMWFDataServer
server = ECMWFDataServer()
server.retrieve({
"class" : "ei",
"dataset": "interim",
"date" : "2018-01-01",
"expver" : "1",
"grid" : "0.75/0.75",
"levtype": "sfc",
"param" : "228.128",
"step" : "12",
"stream" : "oper",
"time" : "00:00:00/12:00:00",
"type" : "fc",
"format" : "netcdf",
"target" : "tp_20180101.nc",
}) |
Run a second script to calculate daily total precipitation. All it does is to add up the two values for 00 - 12 UTC and 12 - 24 UTC for a given day.
#!/usr/bin/env python
"""
Save as file calculate-daily-tp.py and run "python calculate-daily-tp.py".
Input file : tp_20180101.nc
Output file: daily-tp_20180101.nc
"""
import time, sys
from datetime import datetime, timedelta
from netCDF4 import Dataset, date2num, num2date
import numpy as np
day = 20180101
f_in = 'tp_%d.nc' % day
f_out = 'daily-tp_%d.nc' % day
d = datetime.strptime(str(day), '%Y%m%d')
time_needed = [d + timedelta(hours = 12), d + timedelta(days = 1)]
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)
var = ds_dest.createVariable(var_tp.name, np.double, var_tp.dimensions)
var[0, :, :] = data
#var.units = units
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)
|
| For simplicity, data in the output NetCDF file of the above script is unpacked. You may want to pack the data to save some disk spaces. Refer to https://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html#bp_Packed-Data-Values for detailed information. |
Related articles appear here based on the labels you select. Click to edit the macro and add or change labels.
|