interpolate cube data to a given (list of) point(s)

Question

Given a cube with lat/lon coordinates

surface altitude [m] / (m)          (latitude: 21600; longitude: 43200)
     Dimension coordinates:
          latitude                           x                 -
          longitude                          -                 x
     Attributes:
          Conventions: CF-1.4
          description: GTOPO30 surface elevation dataset
          history: Mon Aug 19 14:04:58 2013: ncrename -v surface altitude,surface_altitude...
          institution: Institute of Environmental Physics, University of Bremen,     Germany.
          label: surface altitude [m]
          least_significant_digit: 4
          source: http://eros.usgs.gov/#/Find_Data/Products_and_Data_Available/gtopo30_i...
          title: Global surface elevation from the GTOPO30

and given a list of lat/lon coordinates, how can I interpolate the cube's data to the given points? (For now, nearest-neighbor would be fine, but for coarser cubes, splines would be nice)

PS: For nearest-neighbor, this interpolation should not depend on loading the whole cube into memory.

score 2 · Answer 1 · answered Aug 19 '13 at 15:35

The nearest neighbour code in Iris, sadly, currently loads the data to identify the indexes needed. I've submitted a trivial pull request (made complex by testing) to fix this (https://github.com/SciTools/iris/pull/707) which you could try using to work with this sized dataset.

I'm going to work with a cube from the sample data:

import iris
cube = iris.load_cube(iris.sample_data_path('air_temp.pp'))

And I can check whether there is data loaded with the following function:

def cube_data_is_loaded(cube):
    # A None data manger means the data is loaded...
    return cube._data_manager is None

So:

>>> print cube_data_is_loaded(cube)
False

Basically, the interface (http://scitools.org.uk/iris/docs/latest/iris/iris/analysis/interpolate.html#iris.analysis.interpolate.extract_nearest_neighbour) for nearest neighbour allows you to do a point extraction:

from iris.analysis.interpolate import extract_nearest_neighbour
smaller_cube = extract_nearest_neighbour(cube,
                        (('longitude', -180), ('latitude', 1.5)))

>>> print smaller_cube
air_temperature / (K)               (scalar cube)
     Scalar coordinates:
          forecast_period: 6477 hours
          forecast_reference_time: 1998-03-01 03:00:00
          latitude: 2.50002 degrees
          longitude: 180.0 degrees
          pressure: 1000.0 hPa
          time: 1996-12-01 00:00:00
     Attributes:
          STASH: m01s16i203
          source: Data from Met Office Unified Model
     Cell methods:
          mean: time

Notice how the extraction has actually picked the nearest latitude value to my requested point. One thing that is really important is to note that this function really doesn't handle wrapping if your longitude coordinate is not circular:

cube.coord('longitude').circular = False
smaller_cube = extract_nearest_neighbour(cube,
                   (('longitude', -180), ('latitude', 1.5)))
cube.coord('longitude').circular = True
>>> print smaller_cube
air_temperature / (K)               (scalar cube)
     Scalar coordinates:
          forecast_period: 6477 hours
          forecast_reference_time: 1998-03-01 03:00:00
          latitude: 2.50002 degrees
          longitude: 0.0 degrees
          pressure: 1000.0 hPa
          time: 1996-12-01 00:00:00, bound=(1994-12-01 00:00:00, 1998-12-01 00:00:00)
     Attributes:
          STASH: m01s16i203
          source: Data from Met Office Unified Model
     Cell methods:
          mean: time

Notice how the longitude range in the original cube (0-360) now means that the nearest value to -180 is actually 0.

There is also a function to do a trajectory extraction (http://scitools.org.uk/iris/docs/latest/iris/iris/analysis/trajectory.html?highlight=trajectory#iris.analysis.trajectory.interpolate):

smaller_traj = interpolate(cube,
                  (('longitude', [-180, -180]), ('latitude', [1.5, 3.5])),
                  'nearest')
>>> print smaller_traj
air_temperature / (K)               (*ANONYMOUS*: 2)
     Auxiliary coordinates:
          latitude                              x
          longitude                             x
     Scalar coordinates:
          forecast_period: 6477 hours
          forecast_reference_time: 1998-03-01 03:00:00
          pressure: 1000.0 hPa
          time: 1996-12-01 00:00:00, bound=(1994-12-01 00:00:00, 1998-12-01 00:00:00)
     Attributes:
          STASH: m01s16i203
          source: Data from Met Office Unified Model
     Cell methods:
          mean: time

Finally, notice how the original cube's data has not been loaded throughout (using my branch), and indeed, the data from smaller_cube has also been deferred:

>>> print cube_data_is_loaded(cube)
False
>>> print cube_data_is_loaded(smaller_cube)
False

For the trajectory, in general, deferred loading is not possible, but it is worth noting that when using NetCDF the indices are being passed right through to the underlying NetCDF library so that at no point is the entire array in memory.

HTH!

P.S. I'm not aware of any spline interpolation algorithms which work with Iris yet, though there is a similar interface to do linear interpolation, should that be of interest.

Thanks, @pelson, for this detailed answer! However, I don't understand your comment about trajectories, netCDF, and indices. When I try the `trajectory.interpolate` method on my cube to extract four points using `nearest` then the 4G RAM of my desktop aren't enough and the system starts swapping (or do I need your newest PR to take advantage of this also in the trajectory interface?) — andreas-h, Aug 20 '13 at 08:55
Yes. The new PR is absolutely necessary to get the interpolation without loading all of the data. It was a trivial but necessary change for this to work as you would like. — pelson, Aug 20 '13 at 09:47

interpolate cube data to a given (list of) point(s)

1 Answers1