|
| 1 | +""" |
| 2 | +Read data from ECMWF MACC Reanalysis. |
| 3 | +""" |
| 4 | + |
| 5 | +from __future__ import division |
| 6 | +import threading |
| 7 | +import pandas as pd |
| 8 | + |
| 9 | +try: |
| 10 | + import netCDF4 |
| 11 | +except ImportError: |
| 12 | + class netCDF4: |
| 13 | + @staticmethod |
| 14 | + def Dataset(*a, **kw): |
| 15 | + raise ImportError( |
| 16 | + 'Reading ECMWF data requires netCDF4 to be installed.') |
| 17 | + |
| 18 | +try: |
| 19 | + from ecmwfapi import ECMWFDataServer |
| 20 | +except ImportError: |
| 21 | + def ECMWFDataServer(*a, **kw): |
| 22 | + raise ImportError( |
| 23 | + 'To download data from ECMWF requires the API client.\nSee https:/' |
| 24 | + '/confluence.ecmwf.int/display/WEBAPI/Access+ECMWF+Public+Datasets' |
| 25 | + ) |
| 26 | + SERVER = None |
| 27 | +else: |
| 28 | + SERVER = ECMWFDataServer() |
| 29 | + |
| 30 | +#: map of ECMWF MACC parameter keynames and codes used in API |
| 31 | +PARAMS = { |
| 32 | + "tcwv": "137.128", |
| 33 | + "aod550": "207.210", |
| 34 | + 'aod469': '213.210', |
| 35 | + 'aod670': '214.210', |
| 36 | + 'aod865': '215.210', |
| 37 | + "aod1240": "216.210", |
| 38 | +} |
| 39 | + |
| 40 | + |
| 41 | +def _ecmwf(server, startdate, stopdate, params, targetname): |
| 42 | + # see http://apps.ecmwf.int/datasets/data/macc-reanalysis/levtype=sfc/ |
| 43 | + server.retrieve({ |
| 44 | + "class": "mc", |
| 45 | + "dataset": "macc", |
| 46 | + "date": "%s/to/%s" % (startdate, stopdate), |
| 47 | + "expver": "rean", |
| 48 | + "grid": "0.75/0.75", |
| 49 | + "levtype": "sfc", |
| 50 | + "param": params, |
| 51 | + "step": "3/6/9/12/15/18/21/24", |
| 52 | + "stream": "oper", |
| 53 | + "format": "netcdf", |
| 54 | + "time": "00:00:00", |
| 55 | + "type": "fc", |
| 56 | + "target": targetname, |
| 57 | + }) |
| 58 | + |
| 59 | + |
| 60 | +def get_ecmwf_macc(filename, params, startdate, stopdate, lookup_params=True, |
| 61 | + server=SERVER, target=_ecmwf): |
| 62 | + """ |
| 63 | + Download data from ECMWF MACC Reanalysis API. |
| 64 | +
|
| 65 | + Parameters |
| 66 | + ---------- |
| 67 | + filename : str |
| 68 | + full path of file where to save data, ``.nc`` appended if not given |
| 69 | + params : str or sequence of str |
| 70 | + keynames of parameter[s] to download |
| 71 | + startdate : datetime.datetime or datetime.date |
| 72 | + UTC date |
| 73 | + stopdate : datetime.datetime or datetime.date |
| 74 | + UTC date |
| 75 | + lookup_params : bool, default True |
| 76 | + optional flag, if ``False``, then codes are already formatted |
| 77 | + server : ecmwfapi.api.ECMWFDataServer |
| 78 | + optionally provide a server object, default is given |
| 79 | + target : callable |
| 80 | + optional function that calls ``server.retrieve`` to pass to thread |
| 81 | +
|
| 82 | + Returns |
| 83 | + ------- |
| 84 | + t : thread |
| 85 | + a thread object, use it to check status by calling `t.is_alive()` |
| 86 | +
|
| 87 | + Notes |
| 88 | + ----- |
| 89 | + To download data from ECMWF requires the API client. For more information, |
| 90 | + see the `documentation |
| 91 | + <https://confluence.ecmwf.int/display/WEBAPI/Access+ECMWF+Public+Datasets>`_. |
| 92 | +
|
| 93 | + This function returns a daemon thread that runs in the background. Exiting |
| 94 | + Python will kill this thread, however this thread will not block the main |
| 95 | + thread or other threads. This thread will terminate when the file is |
| 96 | + downloaded or if the thread raises an unhandled exception. You may submit |
| 97 | + multiple requests simultaneously to break up large downloads. You can also |
| 98 | + check the status and retrieve downloads online at |
| 99 | + http://apps.ecmwf.int/webmars/joblist/. This is useful if you kill the |
| 100 | + thread. Downloads expire after 24 hours. |
| 101 | +
|
| 102 | + .. warning:: Your request may be queued online for an hour or more before |
| 103 | + it begins to download |
| 104 | +
|
| 105 | + Precipitable water :math:`P_{wat}` is equivalent to the total column of |
| 106 | + water vapor (TCWV), but the units given by ECMWF MACC Reanalysis are kg/m^2 |
| 107 | + at STP (1-atm, 25-C). Divide by ten to convert to centimeters of |
| 108 | + precipitable water: |
| 109 | +
|
| 110 | + .. math:: |
| 111 | + P_{wat} \\left( \\text{cm} \\right) \ |
| 112 | + = TCWV \\left( \\frac{\\text{kg}}{\\text{m}^2} \\right) \ |
| 113 | + \\frac{100 \\frac{\\text{cm}}{\\text{m}}} \ |
| 114 | + {1000 \\frac{\\text{kg}}{\\text{m}^3}} |
| 115 | +
|
| 116 | + The keynames available for the ``params`` argument are given by |
| 117 | + :const:`pvlib.iotools.ecmwf_macc.PARAMS` which maps the keys to codes used |
| 118 | + in the API. The following keynames are available: |
| 119 | +
|
| 120 | + ======= ========================================= |
| 121 | + keyname description |
| 122 | + ======= ========================================= |
| 123 | + tcwv total column water vapor in kg/m^2 at STP |
| 124 | + aod550 aerosol optical depth measured at 550-nm |
| 125 | + aod469 aerosol optical depth measured at 469-nm |
| 126 | + aod670 aerosol optical depth measured at 670-nm |
| 127 | + aod865 aerosol optical depth measured at 865-nm |
| 128 | + aod1240 aerosol optical depth measured at 1240-nm |
| 129 | + ======= ========================================= |
| 130 | +
|
| 131 | + If ``lookup_params`` is ``False`` then ``params`` must contain the codes |
| 132 | + preformatted according to the ECMWF MACC Reanalysis API. This is useful if |
| 133 | + you want to retrieve codes that are not mapped in |
| 134 | + :const:`pvlib.iotools.ecmwf_macc.PARAMS`. |
| 135 | +
|
| 136 | + Specify a custom ``target`` function to modify how the ECMWF API function |
| 137 | + ``server.retrieve`` is called. The ``target`` function must have the |
| 138 | + following signature in which the parameter definitions are similar to |
| 139 | + :func:`pvlib.iotools.get_ecmwf_macc`. :: |
| 140 | +
|
| 141 | +
|
| 142 | + target(server, startdate, stopdate, params, filename) -> None |
| 143 | +
|
| 144 | + Examples |
| 145 | + -------- |
| 146 | + Retrieve the AOD measured at 550-nm and the total column of water vapor for |
| 147 | + November 1, 2012. |
| 148 | +
|
| 149 | + >>> from datetime import date |
| 150 | + >>> from pvlib.iotools import get_ecmwf_macc |
| 151 | + >>> filename = 'aod_tcwv_20121101.nc' # .nc extension added if missing |
| 152 | + >>> params = ('aod550', 'tcwv') |
| 153 | + >>> start = end = date(2012, 11, 1) |
| 154 | + >>> t = get_ecmwf_macc(filename, params, start, end) |
| 155 | + >>> t.is_alive() |
| 156 | + True |
| 157 | +
|
| 158 | + """ |
| 159 | + if not filename.endswith('nc'): |
| 160 | + filename += '.nc' |
| 161 | + if lookup_params: |
| 162 | + try: |
| 163 | + params = '/'.join(PARAMS.get(p) for p in params) |
| 164 | + except TypeError: |
| 165 | + params = PARAMS.get(params) |
| 166 | + startdate = startdate.strftime('%Y-%m-%d') |
| 167 | + stopdate = stopdate.strftime('%Y-%m-%d') |
| 168 | + if not server: |
| 169 | + server = ECMWFDataServer() |
| 170 | + t = threading.Thread(target=target, daemon=True, |
| 171 | + args=(server, startdate, stopdate, params, filename)) |
| 172 | + t.start() |
| 173 | + return t |
| 174 | + |
| 175 | + |
| 176 | +class ECMWF_MACC(object): |
| 177 | + """container for ECMWF MACC reanalysis data""" |
| 178 | + |
| 179 | + TCWV = 'tcwv' # total column water vapor in kg/m^2 at (1-atm,25-degC) |
| 180 | + |
| 181 | + def __init__(self, filename): |
| 182 | + self.data = netCDF4.Dataset(filename) |
| 183 | + # data variables and dimensions |
| 184 | + variables = set(self.data.variables.keys()) |
| 185 | + dimensions = set(self.data.dimensions.keys()) |
| 186 | + self.keys = tuple(variables - dimensions) |
| 187 | + # size of lat/lon dimensions |
| 188 | + self.lat_size = self.data.dimensions['latitude'].size |
| 189 | + self.lon_size = self.data.dimensions['longitude'].size |
| 190 | + # spatial resolution in degrees |
| 191 | + self.delta_lat = -180.0 / (self.lat_size - 1) # from north to south |
| 192 | + self.delta_lon = 360.0 / self.lon_size # from west to east |
| 193 | + # time resolution in hours |
| 194 | + self.time_size = self.data.dimensions['time'].size |
| 195 | + self.start_time = self.data['time'][0] |
| 196 | + self.stop_time = self.data['time'][-1] |
| 197 | + self.time_range = self.stop_time - self.start_time |
| 198 | + self.delta_time = self.time_range / (self.time_size - 1) |
| 199 | + |
| 200 | + def get_nearest_indices(self, latitude, longitude): |
| 201 | + """ |
| 202 | + Get nearest indices to (latitude, longitude). |
| 203 | +
|
| 204 | + Parmaeters |
| 205 | + ---------- |
| 206 | + latitude : float |
| 207 | + Latitude in degrees |
| 208 | + longitude : float |
| 209 | + Longitude in degrees |
| 210 | +
|
| 211 | + Returns |
| 212 | + ------- |
| 213 | + idx_lat : int |
| 214 | + index of nearest latitude |
| 215 | + idx_lon : int |
| 216 | + index of nearest longitude |
| 217 | + """ |
| 218 | + # index of nearest latitude |
| 219 | + idx_lat = int(round((latitude - 90.0) / self.delta_lat)) |
| 220 | + # avoid out of bounds latitudes |
| 221 | + if idx_lat < 0: |
| 222 | + idx_lat = 0 # if latitude == 90, north pole |
| 223 | + elif idx_lat > self.lat_size: |
| 224 | + idx_lat = self.lat_size # if latitude == -90, south pole |
| 225 | + # adjust longitude from -180/180 to 0/360 |
| 226 | + longitude = longitude % 360.0 |
| 227 | + # index of nearest longitude |
| 228 | + idx_lon = int(round(longitude / self.delta_lon)) % self.lon_size |
| 229 | + return idx_lat, idx_lon |
| 230 | + |
| 231 | + def interp_data(self, latitude, longitude, utc_time, param): |
| 232 | + """ |
| 233 | + Interpolate ``param`` values to ``utc_time`` using indices nearest to |
| 234 | + (``latitude, longitude``). |
| 235 | +
|
| 236 | + Parmaeters |
| 237 | + ---------- |
| 238 | + latitude : float |
| 239 | + Latitude in degrees |
| 240 | + longitude : float |
| 241 | + Longitude in degrees |
| 242 | + utc_time : datetime.datetime or datetime.date |
| 243 | + Naive or UTC date or datetime to interpolate |
| 244 | + param : str |
| 245 | + Name of the parameter to interpolate from the data |
| 246 | +
|
| 247 | + Returns |
| 248 | + ------- |
| 249 | + Interpolated ``param`` value at (``utc_time, latitude, longitude``) |
| 250 | +
|
| 251 | + Examples |
| 252 | + -------- |
| 253 | + Use this to get a single value of a parameter in the data at a specific |
| 254 | + time and set of (latitude, longitude) coordinates. |
| 255 | +
|
| 256 | + >>> from datetime import datetime |
| 257 | + >>> from pvlib.iotools import ecmwf_macc |
| 258 | + >>> data = ecmwf_macc.ECMWF_MACC('aod_tcwv_20121101.nc') |
| 259 | + >>> dt = datetime(2012, 11, 1, 11, 33, 1) |
| 260 | + >>> data.interp_data(38.2, -122.1, dt, 'aod550') |
| 261 | + """ |
| 262 | + nctime = self.data['time'] # time |
| 263 | + ilat, ilon = self.get_nearest_indices(latitude, longitude) |
| 264 | + # time index before |
| 265 | + before = netCDF4.date2index(utc_time, nctime, select='before') |
| 266 | + fbefore = self.data[param][before, ilat, ilon] |
| 267 | + fafter = self.data[param][before + 1, ilat, ilon] |
| 268 | + dt_num = netCDF4.date2num(utc_time, nctime.units) |
| 269 | + time_ratio = (dt_num - nctime[before]) / self.delta_time |
| 270 | + return fbefore + (fafter - fbefore) * time_ratio |
| 271 | + |
| 272 | + |
| 273 | +def read_ecmwf_macc(filename, latitude, longitude, utc_time_range=None): |
| 274 | + """ |
| 275 | + Read data from ECMWF MACC reanalysis netCDF4 file. |
| 276 | +
|
| 277 | + Parameters |
| 278 | + ---------- |
| 279 | + filename : string |
| 280 | + full path to netCDF4 data file. |
| 281 | + latitude : float |
| 282 | + latitude in degrees |
| 283 | + longitude : float |
| 284 | + longitude in degrees |
| 285 | + utc_time_range : sequence of datetime.datetime |
| 286 | + pair of start and stop naive or UTC date-times |
| 287 | +
|
| 288 | + Returns |
| 289 | + ------- |
| 290 | + data : pandas.DataFrame |
| 291 | + dataframe for specified range of UTC date-times |
| 292 | + """ |
| 293 | + ecmwf_macc = ECMWF_MACC(filename) |
| 294 | + try: |
| 295 | + ilat, ilon = ecmwf_macc.get_nearest_indices(latitude, longitude) |
| 296 | + nctime = ecmwf_macc.data['time'] |
| 297 | + if utc_time_range: |
| 298 | + start_idx = netCDF4.date2index( |
| 299 | + utc_time_range[0], nctime, select='before') |
| 300 | + stop_idx = netCDF4.date2index( |
| 301 | + utc_time_range[-1], nctime, select='after') |
| 302 | + time_slice = slice(start_idx, stop_idx + 1) |
| 303 | + else: |
| 304 | + time_slice = slice(0, ecmwf_macc.time_size) |
| 305 | + times = netCDF4.num2date(nctime[time_slice], nctime.units) |
| 306 | + df = {k: ecmwf_macc.data[k][time_slice, ilat, ilon] |
| 307 | + for k in ecmwf_macc.keys} |
| 308 | + if ECMWF_MACC.TCWV in df: |
| 309 | + # convert total column water vapor in kg/m^2 at (1-atm, 25-degC) to |
| 310 | + # precipitable water in cm |
| 311 | + df['precipitable_water'] = df[ECMWF_MACC.TCWV] / 10.0 |
| 312 | + finally: |
| 313 | + ecmwf_macc.data.close() |
| 314 | + return pd.DataFrame(df, index=times.astype('datetime64[s]')) |
0 commit comments