how to extract ECMWF ERA-5 data along animal track points?

Question

I want to extract / interpolate an air temperature value for each track location (in space [xy], time [t], and altitude) from ERA5 hourly pressure level data.

I retrieved the ERA5 dataset using CDS toolbox like below. But, I could not figure out how to extract the values for each point. I tried to use the tool: 'ct.observation.interp_from_grid()' within CDS toolbox without success.

import cdstoolbox as ct

# Initialise the application
@ct.application(title='my trial to retrieve and annotate movement data')

# Define a download output for the application
@ct.output.download()

# Define application function
def application():
    """Define a function that extracts hourly Air Temperature in 2018 for track points and provides a download link.

    # Retrieve hourly air temperature
    data = ct.catalogue.retrieve(
        'reanalysis-era5-pressure-levels',
        {
            'variable': 'temperature',
            'product_type': 'reanalysis',
            'pressure_level': [
            '900', '925', '950',
        ],
            'year': 2018,
            'month': '05',
            'day': '01',
            'time': [
                '00:00', '01:00', '02:00', '03:00',
                '04:00', '05:00', '06:00', '07:00',
                '08:00', '09:00', '10:00', '11:00',
                '12:00', '13:00', '14:00', '15:00',
                '16:00', '17:00', '18:00', '19:00',
                '20:00', '21:00', '22:00', '23:00',
            ],
        'area': [
            48, 111, 47,
            112,
        ],
            }
    )
    # Interpolate data for track points
    indexers = {'lon': [146.29, 147.10], 'lat': [-6.689, -7.644], 'plev':['891', '653'], 'time':['2019-03-23 18:52:29', '2019-03-23 21:52:30']}
    points = ct.observation.interp_from_grid(data, method='nearest', drop='False', wrap_lon='False', **indexers)    
    
    print(points)
    return points

Alternatively, I can download the ERA5 data first and then use the extract function of the raster package in R. However, I do not want to download a huge amount of datasets (probably hundreds of GB, even TB) on my computer since my track points cover large spatial and temporal scales.

here is a dummy track points just for demo.

structure(list(Latitude = c(-6.689718, -7.644683, -8.31021, -9.177921, 
-9.493564), Longitude = c(146.297638, 147.107101, 148.211472, 
148.670151, 149.00795), timestamp = c("2019-03-23 15:52:14", 
"2019-03-23 18:52:29", "2019-03-23 21:52:30", "2019-03-24 00:52:29", 
"2019-03-24 03:52:15"), altitude_hPa = c(891, 653, 521, 910, 
711)), class = "data.frame", row.names = c(NA, -5L)) 

I would be so grateful for any suggestions or other ways of doing this.

thanks in advance,

Bat

Answer 1

Hy Bat

I was not aware of the cdstoolbox so far but have had a bit of a deeper look given your demo request (used cdstoolbox-remote; very convenient!). I got the problem tracked down to the interp_from_grids method which contains the following line(s) of code:

if 'time' in indexers:
   indexers['time'] = indexers['time'].astype('float64')

If indexers contains "time" the method tries to convert it into float64 - which does not work with the str list in your demo. To get around this, I converted the "time" array into a numpy.datetime64 object. Something like:

numpy.array(['2019-03-23 18:52:29', '2019-03-23 21:52:30'], dtype = 'datetime64')

This resolves the "AttributeError: 'list' object has no attribute 'astype'" error (as it is now convertable into float64), however, it is not JSON serializable (new Error: "AttributeError: 'list' object has no attribute 'astype'").

At this point I am a bit lost - does time-interpolation even work? The method (did not have time to dig even deeper) seems to handle "time" somehow, however, I can not find an example on the cdstoolbox website. Is time-point-interpolation even possible?

All best, R