Source code for IRData.cera20c.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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#

#Load CERA20C data from previously downloaded files.

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

# Eliminate incomprehensible warning message
iris.FUTURE.netcdf_promote='True'

# Need to add coordinate system metadata so they work with cartopy
coord_s=iris.coord_systems.GeogCS(iris.fileformats.pp.EARTH_RADIUS)

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

def _get_previous_field_time(variable,year,month,day,hour):
    """Get the latest time, before the given time,
                     for which there is saved data"""
    return {'year':year,'month':month,'day':day,'hour':int(hour/3)*3}

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

def _get_slice_at_hour_at_timestep(variable,year,month,day,hour,
                                   fc_init=None,deaccumulate=True):
    
    # 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):
        raise ValueError("Invalid hour - data not in file")

    # Precipitation is accumulated over the forecast, convert to
    # Accumulation over last (3-hour timestep)
    # Still in weird units, be careful.
    if (variable=='prate' and deaccumulate and
       (hour!=21 or (fc_init is not None and fc_init=='last'))):
        r1=_get_slice_at_hour_at_timestep(variable,year,month,day,hour,
                                           fc_init=fc_init,deaccumulate=False)
        # Subtract the values from 3 hours ago
        lt=datetime.datetime(year,month,day,int(hour))-datetime.timedelta(hours=3)
        r2=_get_slice_at_hour_at_timestep(variable,lt.year,lt.month,lt.day,lt.hour,
                                           fc_init=fc_init,deaccumulate=False)
        r1.data=(r1.data-r2.data) # to m/(3-hour period)
        # convert to kg m**-2 s**-1 (same as 20CR)
        # Decided against doing this.
        #r1.data=r1.data/10.8    # to kg m**-2 s**-1 - same as 20CR
        #r1.units='kg m**-2 s**-1'
        return r1

    # Not precipitation - just get the data for this timestep
    file_name=_hourly_get_file_name(variable,year,month,day,hour,
                                    type=type,fc_init=fc_init)
    if not os.path.isfile(file_name):
        raise StandardError(("%s for %04d/%02d not available"+
                             " might need cera20c.fetch") % (variable,
                                                             year,month))
    time_constraint=iris.Constraint(time=iris.time.PartialDateTime(
                                   year=year,
                                   month=month,
                                   day=day,
                                   hour=hour))
    try:
        with iris.FUTURE.context(cell_datetime_objects=True):
            hslice=iris.load_cube(file_name,
                                  time_constraint)
    # This isn't the right error to catch
    except iris.exceptions.ConstraintMismatchError:
       print("Data not available")

    # Enhance the names and metadata for iris/cartopy
    hslice.coord('latitude').coord_system=coord_s
    hslice.coord('longitude').coord_system=coord_s
    hslice.dim_coords[0].rename('member') # Consistency with 20CR
    return hslice

[docs]def load(variable,dtime,fc_init=None): """Load requested data from disc, interpolating if necessary. Data must be available in directory $SCRATCH/CERA-20C, previously retrieved by :func:`fetch`. Args: variable (:obj:`str`): Variable to fetch (e.g. 'prmsl') dtime (:obj:`datetime.datetime`): Date-time to get data for. fc_init (:obj:`str`): See below Returns: :obj:`iris.cube.Cube`: Global field of variable at time. Note that CERA-20C data is only output every 3 hours, so if hour%3!=0, the result will be linearly interpolated in time. Precipitation data in CERA is a forecast field: once a day (at 18:00) 3-hourly forecast data is calculated for the next 27 hours. So at 21:00, there are 2 sets of precipitation available: a 3-hour forecast starting at 18 that day, and a 27-hour forecast starting at 18:00 the day before; and there is a discontinuity in the fields at that time. This function will always load the shortest lead-time forecast available unless fc_init is set to 'last'. 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 StandardError("Unsupported variable %s" % variable) dhour=dtime.hour+dtime.minute/60.0+dtime.second/3600.0 if _is_in_file(variable,dtime.year,dtime.month, dtime.day,dhour): return(_get_slice_at_hour_at_timestep(variable,dtime.year, dtime.month,dtime.day, dhour, fc_init=fc_init)) previous_step=_get_previous_field_time(variable,dtime.year,dtime.month, dtime.day,dhour) next_step=_get_next_field_time(variable,dtime.year,dtime.month, dtime.day,dhour) 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'], 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'], 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