#!/usr/bin/env python3
"""
Convenient wrappers for python APIs used to download climatological
datasets.
"""
import calendar
import numpy as np
__all__ = [
'cmip',
'era',
'merra',
'ncar',
'satellite',
]
[docs]def era(
params,
stream,
levtype,
daterange=None,
yearrange=None,
monthrange=None,
# dayrange=None, # not yet used
years=None,
months=None,
format='netcdf',
forecast=False,
step=12,
levrange=None,
levs=None,
grid=None,
hours=(0, 6, 12, 18),
hour=None,
res=None,
box=None,
filename='era.nc',
):
"""
Retrieves ERA reanalysis data using the provided API. User must have, in
home directory, a file named ``.ecmwfapirc``. See API documentation, but
should look like:
::
{
"url" : "https://api.ecmwf.int/v1",
"key" : "abcdefghijklmnopqrstuvwxyz",
"email" : "email@gmail.com"
}
with the key found on your user/profile page on the ECMWF website.
Parameters
----------
params : str or list of str
Variable name. Gets translated to MARS id name by dictionary below.
Add to this from the `online GRIB table <https://rda.ucar.edu/datasets/ds627.0/docs/era_interim_grib_table.html>`_.
Pay attention to *available groups*. If not available for the group
you selected (e.g. pressure levs, moda), get ``ERROR 6 (MARS_EXPECTED_FIELDS)``.
For rates of change of *parameterized* processes (i.e. diabatic) see
`this link <https://confluence.ecmwf.int/pages/viewpage.action?pageId=57448466>`_.
stream : {'oper', 'moda', 'mofm', 'mdfa', 'mnth'}
The data stream.
levtype : {'ml', 'pl', 'sfc', 'pt', 'pv'}
Level type: model, pressure, surface, potential temperature, and 2PVU
surface, respectively.
levrange : float or (float, float), optional
Individual level or range of levels.
levs : float or ndarray, optional
Individual level or list of levels.
yearrange : int or (int, int)
Individual year or range of years.
years : int or ndarray, optional
Individual year or list of years.
monthrange : int or (int, int), optional
Individual month or range of months.
months : int or ndarray, optional
Individual month or list of months.
daterange : (datetime.datetime, datetime.datetime), optional
Range of dates.
hours : {0, 6, 12, 18} or list thereof, optional
Hour(s) (UTC) of observation.
forecast : bool, optional
Whether we want forecast `'fc'` or analysis `'an'` data. Note that
some data is only available in `'fc'` mode, e.g. diabatic heating.
grid : str, optional
The grid type. Default is ``N32`` which returns data on 64 latitudes.
res : float, optional
Alternative to `grid` that specifies the desired output grid resolution
in degrees. ERA-Interim has a few valid preset resolutions and will
choose the resolution that most closely matches the input.
box : str or length-4 list of float, optional
String name for particular region, e.g. ``'europe'``, or the west,
south, east, and north boundaries, respectively.
format : {'grib1', 'grib2', 'netcdf'}, optional
Output format.
filename : str, optional
Name of file output.
Notes
-----
Some fields (seems true for most model fields) are not archived as
monthly means for some reason! Have no idea why because it would need
almost zero storage requirements.
""" # noqa: E501
# Data stream
import ecmwfapi as ecmwf # only do so inside function
# Variable id conversion (see:
# https://rda.ucar.edu/datasets/ds627.0/docs/era_interim_grib_table.html)
if isinstance(params, str) or not np.iterable(params):
params = (params,)
params = [
{
'tdt': '110.162',
# model-level surface pressure; use this whenever getting tdt and stuff
# (requires lev=1)
'msp': '152.128',
'sp': '134.128', # surface pressure
't2m': '167.128', # 2m temp
'd2m': '168.128', # 2m dew point
'sst': '34.128', # sst
'msl': '151.128', # sea level pressure
'slp': '151.128', # same
'z': '129.128', # geopotential
't': '130.128', # temp
'u': '131.128', # u wind
'v': '132.128', # v wind
'w': '135.128', # w wind
'q': '133.128', # specific humidity
'r': '157.128', # relative humidity
'vort': '138.128', # relative vorticity
'vo': '138.128', # same
'zeta': '138.128', # same
'pt': '3.128', # potential temp (available on 2pvu surf)
'theta': '3.128', # same
'p': '54.128', # pressure (availble on pt, 2pvu surfaces)
'pres': '54.128', # same
# potential vorticity (available on p, pt surfaces)
'pv': '60.128',
'precip': '228.128',
}.get(p)
for p in params
] # returns generator object for each param
if None in params:
raise ValueError(
'MARS id for variable is unknown (might need to be added to this script).'
)
params = '/'.join(params)
# Time selection as various RANGES or LISTS
# Priority; just use daterange as datetime or date objects
if daterange is not None:
if not np.iterable(daterange):
daterange = (daterange,) # want a SINGLE DAY
# options for monthly or daily data
if stream != 'oper':
y0, m0, y1, m1 = (
daterange[0].year,
daterange[0].month,
daterange[1].year,
daterange[1].month,
)
N = max(y1 - y0 - 1, 0) * 12 + (13 - m0) + m1 # number of months in range
dates = '/'.join(
'%04d%02d00' % (y0 + (m0 + n - 1) // 12, (m0 + n - 1) % 12 + 1)
for n in range(N)
)
else:
# MARS will get calendar days in range
dates = '/to/'.join(d.strftime('%Y%m%d') for d in daterange)
# Alternative; list the years/months desired, and if synoptic, get all
# calendar days within
else:
# First, years
if years is not None:
if not np.iterable(years):
years = (years,) # single month
elif yearrange is not None:
if not np.iterable(yearrange):
years = (yearrange,)
else:
years = tuple(range(yearrange[0], yearrange[1] + 1))
else:
raise ValueError('You must use "years" or "yearrange" kwargs.')
# Next, months (this way, can just download JJA data, for example)
if months is not None:
if not np.iterable(months):
months = (months,) # single month
elif monthrange is not None:
if not np.iterable(monthrange):
months = (monthrange, monthrange)
else:
months = tuple(range(monthrange[0], monthrange[1] + 1))
else:
months = tuple(range(1, 13))
# And get dates; options for monthly means and daily stuff
if stream != 'oper':
dates = '/'.join(
'/'.join('%04d%02d00' % (y, m) for m in months) for y in years
)
else:
dates = '/'.join(
'/'.join(
'/'.join(
'%04d%02d%02d' % (y, m, i + 1)
for i in range(calendar.monthrange(y, m)[1])
)
for m in months
)
for y in years
)
# Level selection as RANGE or LIST
# Update this list if you modify script for ERA5, etc.
levchoices = {
'sfc': None,
'pv': None,
'ml': np.arange(1, 137 + 1),
'pl': np.array(
[
1, 2, 3, 5, 7, 10, 20, 30, 50, 70,
100, 125, 150, 175, 200, 225, 250, 300, 350, 400, 450, 500, 550,
600, 650, 700, 750, 775, 800, 825, 850, 875, 900, 925, 950, 975, 1000,
]
),
'pt': np.array(
[265, 270, 285, 300, 315, 330, 350, 370, 395, 430, 475, 530, 600, 700, 850]
),
}.get(levtype, [])
if levchoices == []:
raise ValueError('Invalid level type. Choose from "pl", "pt", "pv", "sfc".')
# More options
# WARNING: Some vars are technically on "levels" like hybrid level surface
# pressure, but we still need 60.
if levtype not in ('sfc', 'pv'): # these have multiple options
# Require input
if levs is None and levrange is None:
raise ValueError(
'Must specify list of levels to "levs" kwarg, range of levels '
'to "levrange" kwarg, or single level to either one.'
)
# Convert levels to mars request
if levs is not None:
if not np.iterable(levs):
levs = (levs,)
else:
if not np.iterable(levrange):
levs = (levrange,)
else:
levs = levchoices[
(levchoices >= levrange[0]) & (levchoices <= levrange[1])
].flat
levs = '/'.join(str(l) for l in levs)
# Other parameters
# Resolution
# same in latitude/longitude required, for now
if res is not None:
grid = '%.5f/%.5f' % (res, res)
elif grid is None:
grid = 'N32'
# Area - can be specified as pre-defined region (e.g. string 'europe') OR
# n/s/w/e boundary
if box is not None and not isinstance(box, str):
box = '/'.join(str(b) for b in (box[3], box[0], box[2], box[1]))
# Hour conversion
if not np.iterable(hours):
hours = (hours,)
hours = '/'.join(str(h).zfill(2) for h in hours) # zfill padds 0s on left
# Forecast type
# TODO: Fix the 12 hour thing. Works for some params (e.g. diabatic
# heating, has 3, 6, 9, 12) but others have 0, 6, 12, 18.
if forecast:
dtype, step = 'fc', str(step)
else:
dtype, step = 'an', '0'
# Server instructions
# Not really sure what happens in some situations: list so far:
# 1. Evidently if you provide with variable string-name instead of numeric ID,
# MARS will search for correct one; if there is name ambiguity/conflict will
# throw error.
# 2. On GUI framework, ECMWF only offers a few resolution options, but program
# seems to run when requesting custom resolutions like 5deg/5deg
# Can also spit raw output into GRIB; apparently ERA-Interim uses
# bilinear interpolation to make grid of point obs, which makes sense,
# because their reanalysis model just picks out point observations
# from spherical harmonics; so maybe grid cell concept is dumb? Maybe need
# to focus on just using cosine weightings, forget about rest?
request = {
# Important ones
'class': 'ei', # ecmwf classifiction; choose ERA-Interim
'expver': '1',
'dataset': 'interim', # thought we already did that; *shrug*
'type': dtype, # type of field; analysis 'an' or forecast 'fc'
'resol': 'av', # prevents truncation before transformation to geo grid
'gaussian': 'reduced',
'format': format,
'step': step, # NOTE: ignored for non-forecast type
# 'grid': res,
'grid': grid, # 64 latitudes, i.e. T42 truncation
'stream': stream, # product monthly, raw, etc.
'date': dates,
'time': hours,
'levtype': levtype,
'param': params,
'target': filename, # save location
}
maxlen = max(map(len, request.keys()))
if levs is not None:
request.update(levelist=levs)
if box is not None:
request.update(area=box)
if stream == 'oper': # TODO: change?
request.update(hour=hour)
print(
'REQUEST\n'
+ '\n'.join(
f'"{key}": ' + ' ' * (maxlen - len(key)) + f'{value}'
for key, value in request.items()
)
)
# Retrieve DATA with settings
server = ecmwf.ECMWFDataServer()
server.retrieve(request)
return
[docs]def merra():
"""
Download MERRA data. Is this possible?
Warning
-------
Not yet implemented.
"""
raise NotImplementedError
[docs]def ncar():
"""
Download NCAR CFSL data. Is this possible?
Warning
-------
Not yet implemented.
"""
raise NotImplementedError
[docs]def satellite():
"""
Download satellite data.
Warning
-------
Not yet implemented.
"""
raise NotImplementedError
[docs]def cmip():
"""
Download CMIP5 model data.
Warning
-------
Not yet implemented.
"""
raise NotImplementedError