Time variable units "day as %Y%m%d.%f" in python iris

Question

I am hoping someone can help. I am running a few climate models (NetCDF files) in python using iris. All was working well until I added my last model which is formatted differently. The units they use for the time variable in the new models is day as %Y%m%d.%f but in the other models it is days since …. This means that when I try to constrain the time variable I get the following error AttributeError: 'numpy.float64' object has no attribute 'year'. I tried adding a year variable using iriscc.add_year(EARTH3, 'time') but that just brings up the error ‘Unit has undefined calendar’.

I’m wondering if you know how I might fix this? Do I need to convert the calendar type? Or is there is there a way around that? Not sure how to do that anyway!

Thank you! Erika

EDIT: here is the full code for my file the model CanESM2 is working, but the model EARTH3 is not - it is the one with the funny time units.

import matplotlib.pyplot as plt
import iris
import iris.coord_categorisation as iriscc
import iris.plot as iplt
import iris.quickplot as qplt
import iris.analysis.cartography
import cf_units
from cf_units import Unit
import datetime
import numpy as np

def main():
    #-------------------------------------------------------------------------

    #bring in all the GCM models we need and give them a name
    CanESM2= '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/GCM_data/tasmin_Amon_CanESM2_historical_r1i1p1_185001-200512.nc'
    EARTH3= '/exports/csce/datastore/geos/users/s0XXXX/Climate_Modelling/GCM_data/tas_Amon_EC-EARTH_historical_r3i1p1_1850-2009.nc'

    #Load exactly one cube from given file
    CanESM2 = iris.load_cube(CanESM2)
    EARTH3 = iris.load_cube(EARTH3)

    print"CanESM2 time"
    print (CanESM2.coord('time'))
    print "EARTH3 time"
    print (EARTH3.coord('time'))

    #fix EARTH3 time units as they differ from all other models
    t_coord=EARTH3.coord('time')
    t_unit = t_coord.attributes['invalid_units']
    timestep, _, t_fmt_str = t_unit.split(' ')
    new_t_unit_str= '{} since 1850-01-01 00:00:00'.format(timestep) 
    new_t_unit = cf_units.Unit(new_t_unit_str, calendar=cf_units.CALENDAR_STANDARD)

    new_datetimes = [datetime.datetime.strptime(str(dt), t_fmt_str) for dt in t_coord.points]
    new_dt_points = [new_t_unit.date2num(new_dt) for new_dt in new_datetimes]
    new_t_coord = iris.coords.DimCoord(new_dt_points, standard_name='time', units=new_t_unit)

    print "EARTH3 new time"
    print (EARTH3.coord('time'))

    #regrid all models to have same latitude and longitude system, all regridded to model with lowest resolution
    CanESM2 = CanESM2.regrid(CanESM2, iris.analysis.Linear())
    EARTH3 =EARTH3.regrid(CanESM2, iris.analysis.Linear())

    #we are only interested in the latitude and longitude relevant to Malawi (has to be slightly larger than country boundary to take into account resolution of GCMs)
    Malawi = iris.Constraint(longitude=lambda v: 32.0 <= v <= 36., latitude=lambda v: -17. <= v <= -8.)   
    CanESM2 =CanESM2.extract(Malawi)
    EARTH3 =EARTH3.extract(Malawi)

    #time constraignt to make all series the same, for ERAINT this is 1990-2008 and for RCMs and GCMs this is 1961-2005
    iris.FUTURE.cell_datetime_objects = True
    t_constraint = iris.Constraint(time=lambda cell: 1961 <= cell.point.year <= 2005)
    CanESM2 =CanESM2.extract(t_constraint)
    EARTH3 =EARTH3.extract(t_constraint)

    #Convert units to match, CORDEX data is in Kelvin but Observed data in Celsius, we would like to show all data in Celsius
    CanESM2.convert_units('Celsius')
    EARTH3.units = Unit('Celsius') #this fixes EARTH3 which has no units defined
    EARTH3=EARTH3-273 #this converts the data manually from Kelvin to Celsius

    #add year data to files
    iriscc.add_year(CanESM2, 'time')
    iriscc.add_year(EARTH3, 'time')

    #We are interested in plotting the data by year, so we need to take a mean of all the data by year
    CanESM2YR=CanESM2.aggregated_by('year', iris.analysis.MEAN)
    EARTH3YR = EARTH3.aggregated_by('year', iris.analysis.MEAN)

    #Returns an array of area weights, with the same dimensions as the cube
    CanESM2YR_grid_areas = iris.analysis.cartography.area_weights(CanESM2YR)
    EARTH3YR_grid_areas = iris.analysis.cartography.area_weights(EARTH3YR)

    #We want to plot the mean for the whole region so we need a mean of all the lats and lons
    CanESM2YR_mean = CanESM2YR.collapsed(['latitude', 'longitude'], iris.analysis.MEAN, weights=CanESM2YR_grid_areas)   
    EARTH3YR_mean = EARTH3YR.collapsed(['latitude', 'longitude'], iris.analysis.MEAN, weights=EARTH3YR_grid_areas) 

    #-------------------------------------------------------------------------
    #PART 4: PLOT LINE GRAPH
    #limit x axis    
    plt.xlim((1961,2005)) 

    #assign the line colours and set x axis to 'year' rather than 'time'
    qplt.plot(CanESM2YR_mean.coord('year'), CanESM2YR_mean, label='CanESM2', lw=1.5, color='blue')
    qplt.plot(EARTH3YR_mean.coord('year'), EARTH3YR_mean, label='EC-EARTH (r3i1p1', lw=1.5, color='magenta')

    #set a title for the y axis
    plt.ylabel('Near-Surface Temperature (degrees Celsius)')

    #create a legend and set its location to under the graph
    plt.legend(loc="upper center", bbox_to_anchor=(0.5,-0.05), fancybox=True, shadow=True, ncol=2)

    #create a title
    plt.title('Tas for Malawi 1961-2005', fontsize=11)   

    #add grid lines
    plt.grid()

    #show the graph in the console
    iplt.show()

if __name__ == '__main__':
    main()

Answer 1

In Iris, unit strings for time coordinates must be specified in the format <time-period> since <epoch>, where <time-period> is a unit of measure of time, such as 'days', or 'years'. This format is specified by udunits2, the library Iris uses to supply valid units and perform unit conversions.

The time coordinate in this case does not have a unit that follows this format, meaning the time coordinate will not have full time coordinate functionality (this partly explains the Attribute Error in the question). To fix this we will need to construct a new time coordinate based on the values and metadata of the existing time coordinate and then replace the cube's existing time coordinate with the new one.

To do this we'll need to:

1. construct a new time unit based on the metadata contained in the existing time unit
2. take the existing time coordinate's point values and format them as datetime objects, using the format string specified in the existing time unit
3. convert the datetime objects from (2.) to an array of floating-point numbers using the new time unit constructed in (1.)
4. create a new time coordinate from the array constructed in (3.) and the new time unit produced in (1.)
5. remove the old time coordinate from the cube and add the new one.

Here's the code to do this...

import datetime
import cf_units
import iris
import numpy as np

t_coord = EARTH3.coord('time')

t_unit = t_coord.attributes['invalid_units']
timestep, _, t_fmt_str = t_unit.split(' ')
new_t_unit_str = '{} since 1850-01-01 00:00:00'.format(timestep)
new_t_unit = cf_units.Unit(new_t_unit_str, calendar=cf_units.CALENDAR_STANDARD)

new_datetimes = [datetime.datetime.strptime(str(dt), t_fmt_str) for dt in t_coord.points]
new_dt_points = [new_t_unit.date2num(new_dt) for new_dt in new_datetimes]
new_t_coord = iris.coords.DimCoord(new_dt_points, standard_name='time', units=new_t_unit)

t_coord_dim = cube.coord_dims('time')
cube.remove_coord('time')
cube.add_dim_coord(new_t_coord, t_coord_dim)

I've made an assumption about the best epoch for your time data. I've also made an assumption about the calendar that best describes your data, but you should be able to replace (when constructing new_t_unit) the standard calendar I've chosen with any other valid cf_units calendar without difficulty.

As a final note, it is effectively impossible to change calendar types. This is because different calendar types include and exclude different days. For example, a 360day calendar has a Feb 30 but no May 31 (as it assumes 12 idealised 30 day long months). If you try and convert from a 360day calendar to a standard calendar, problems you hit include what you do with the data from 29 and 30 Feb, and how you fill the five missing days that don't exist in a 360day calendar. For such reasons it's generally impossible to convert calendars (and Iris doesn't allow such operations).

Hope this helps!

Answer 2

Maybe the answer is not more useful however I write here the function that I made in order to convert the data from %Y%m%d.%f in datetime array.

The function create a perfect datetime array, without missing values, it can be modified to take into account if there are missing times, however a climate model should not have missing data.

def fromEARTHtime2Datetime(dt,timeVecEARTH):
    """
    This function returns the perfect array from the EARTH %Y%m%d.%f time 
    format and convert it to a more useful time, such as the time array
    from the datetime of pyhton, this is WHTOUT any missing data!

    Parameters
    ----------
    dt : string
        This is the time discretization, it can be 1h or 6h, but always it 
        needs to be hours, example dt = '6h'.
        
    timeVecEARTH : array of float
        Vector of the time to be converted. For example the time of the
        EARTH is day as %Y%m%d.%f.
        And only this format can be converted to datetime, for example:
            20490128.0,20490128.25,20490128.5,20490128.75 this will be converted
            in datetime: '2049-01-28 00:00:00', '2049-01-28 60:00:00',
                         '2049-01-28 12:00:00','2049-01-28 18:00:00'

    Returns
    -------
    timeArrNew : datetime
        This is the perfect and WITHOUT any missing data datatime array, 
        for example: DatetimeIndex(['2049-01-28 00:00:00', '2049-01-28 06:00:00',
                                        ...
                                   '2049-02-28 18:00:00', '2049-03-01 00:00:00'],
                                   dtype='datetime64[ns]', length=129, freq='6H')

    """

    dtDay = 24/np.float(dt[:-1])
    partOfDay = np.arange(0,1,1/dtDay)
    hDay = []
    for ip in partOfDay:
        hDay.append('%02.f:00:00' %(24*ip))
    dictHours = dict(zip(partOfDay,hDay))
    
    t0Str = str(timeVecEARTH[0])
    timeAux0 = t0Str.split('.')
    timeAux0 = timeAux0[0][0:4] +'-' + timeAux0[0][4:6] +'-' + timeAux0[0][6:] + ' ' + dictHours[float(timeAux0[1])]
    
    tendStr = str(timeVecEARTH[-1])
    timeAuxEnd = tendStr.split('.')
    timeAuxEnd = timeAuxEnd[0][0:4] +'-' + timeAuxEnd[0][4:6] +'-' + timeAuxEnd[0][6:] + ' ' + dictHours[float(timeAuxEnd[1])]
    
    timeArrNew = pd.date_range(timeAux0,timeAuxEnd, freq=dt)
    
    return timeArrNew