Source code for IRData.era5.load

# (C) British Crown Copyright 2017, Met Office
# This code is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
# This code is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# GNU Lesser General Public License for more details.

# The functions in this module provide the main way to load
# ERA5 data.

import os
import os.path
import iris
import iris.time
import datetime
import numpy as np

from .utils import _hourly_get_file_name
from .utils import _translate_for_file_names
from .utils import monolevel_analysis
from .utils import monolevel_forecast

# Need to add coordinate system metadata so they work with cartopy

def _is_in_file(variable,year,month,day,hour,stream='enda'):
    """Is the variable available for this time?
       Or will it have to be interpolated?"""
    if stream=='enda' and hour%3==0:
        return True
    if stream=='oper' and hour%1==0:
        return True
    return False

def _get_previous_field_time(variable,year,month,day,hour,stream='enda'):
    """Get the latest time, before the given time,
                     for which there is saved data"""
    if stream=='enda':
        return {'year':year,'month':month,'day':day,'hour':int(hour/3)*3}
    if stream=='oper':
        return {'year':year,'month':month,'day':day,'hour':int(hour)}
    raise Exception("Unknown stream %s" % stream)

def _get_next_field_time(variable,year,month,day,hour,stream='enda'):
    """Get the earliest time, after the given time,
                     for which there is saved data"""  
    if stream=='enda':
       dr = {'year':year,'month':month,'day':day,'hour':int(hour/3)*3+3}
    elif stream=='oper':
       dr = {'year':year,'month':month,'day':day,'hour':int(hour)+1}
       raise Exception("Unknown stream %s" % stream) 
    if dr['hour']>=24:
        d_next= (['year'],dr['month'],dr['day']) 
                 + datetime.timedelta(days=1) )
        dr = {'year':d_next.year,'month':d_next.month,'day',
    return dr

def _get_slice_at_hour_at_timestep(variable,year,month,day,hour,
    """Get the cube with the data, given that the specified time
       matches a data timestep."""
    if not _is_in_file(variable,year,month,day,hour,stream=stream):
        raise ValueError("Invalid hour - data not in file")
    if not os.path.isfile(file_name):
        raise Exception(("%s for %04d/%02d not available"+
                             " might need era5.fetch") % (variable,
    # This isn't the right error to catch
    except iris.exceptions.ConstraintMismatchError:
       print("Data not available")

    # Enhance the names and metadata for iris/cartopy
    if stream=='enda':
        hslice.dim_coords[0].rename('member') # Consistency with 20CR
    return hslice

[docs]def load(variable,dtime, stream='enda',fc_init=None): """Load requested data from disc, interpolating if necessary. Data must be available in directory $SCRATCH/ERA5, previously retrieved by :func:`fetch`. Args: variable (:obj:`str`): Variable to fetch (e.g. 'prmsl') dtime (:obj:`datetime.datetime`): Date and time to load data for. fc_init (:obj:`int`): If not None; 6 or 18. See below. Returns: :obj:`iris.cube.Cube`: Global field of variable at time. Note that ERA5 data is only output every 3 hours (enda) or hour (oper), so if hour%3!=0, the result may be linearly interpolated in time. Precipitation data in ERA5 is a forecast field: twice a day (at 06:00 and 18:00) forecast data is calculated for the next 18 hours. So at 21:00, for example, there are 2 sets of precipitation available: a 3-hour forecast starting at 18 that day, and a 15-hour forecast starting at 06:00, and there is a discontinuity between the two fields. This function will always load the shortest lead-time forecast available unless fc_init is set (to 6 or 18) in which case it will load the forecast starting at the hour specified. You will only need this if you are making videos, or otherwise need time-continuous forecast fields, in which case you will need to be clever in smoothing over the discontinuity. For analysis fields (everything except prate), this issue does not arise and fc_init is ignored. Raises: StandardError: Data not on disc - see :func:`fetch` | """ if ((variable not in monolevel_analysis) and (variable not in monolevel_forecast)): raise Exception("Unsupported variable %s" % variable) dhour=dtime.hour+dtime.minute/60.0+dtime.second/3600.0 if _is_in_file(variable, dtime.year,dtime.month,,dhour, stream=stream): return(_get_slice_at_hour_at_timestep(variable,dtime.year, dtime.month,, dhour,stream=stream, fc_init=fc_init)) previous_step=_get_previous_field_time(variable,dtime.year,dtime.month,,dhour,stream=stream) next_step=_get_next_field_time(variable,dtime.year,dtime.month,,dhour,stream=stream) dt_current=dtime dt_previous=datetime.datetime(previous_step['year'], previous_step['month'], previous_step['day'], previous_step['hour']) dt_next=datetime.datetime(next_step['year'], next_step['month'], next_step['day'], next_step['hour']) s_previous=_get_slice_at_hour_at_timestep(variable, previous_step['year'], previous_step['month'], previous_step['day'], previous_step['hour'], stream=stream, fc_init=fc_init) s_next=_get_slice_at_hour_at_timestep(variable, next_step['year'], next_step['month'], next_step['day'], next_step['hour'], stream=stream, fc_init=fc_init) # Iris won't merge cubes with different attributes s_previous.attributes=s_next.attributes s_next=iris.cube.CubeList((s_previous,s_next)).merge_cube() s_next=s_next.interpolate([('time',dt_current)],iris.analysis.Linear()) return s_next