Interpolation/regridding/resampling - call for requirements

[RHattersley]

Hello everyone,

Iris already supports some simple forms of interpolation and regridding, but we're keen to build on that.

So, if there is some interpolation/regridding/resampling capability that you would like Iris to have, then please could you let us know. The more specific the better. Whatever requirements we collect over the next week or so will drive our next round of development on this topic.

Kind regards,
Richard

[Phil Elson]

Just as a clarification, you didn't specify spatial regridding (on the horizontal). Was that an intentional omission?
If so, I'd like to highlight that fact in order to get as much diverse information as is possible.

HTH

[RHattersley]

Hi Phil,

Thanks for checking. At this stage I'm interested in pretty much anything which involves resampling in some form or another. If we end up with too much material then we'll just tackle it in several steps.

Richard

[Andy Hartley]

Hi Richard,

Regridding to a new horizontal coordinate system would be very useful. If I'm sharing data with project partners or intersecting with polygons in R (I'd also like to do this in python :-) ), I frequently need to reproject rotated pole to lat/long. If you were to include this feature, I imagine reprojection to any other coordinate system would be equally possible.

I'm an iris newbie, so correct me if this is already included. I had a good search through the documentation, but couldn't find it.

Kind regards,
Andy

[Susan]

Hi Richard,

+1 for the ability to regrid to different co-ordinate systems.

It would be good to be able to regrid from lat, lon to national grid (OS) and vice versa as well as from rotated pole to lat, lon. One constraint on this regridding would be mass conservation.
i.e. Air concentration if regridded from lat, lon to national grid needs to conserve mass.

Susan

[Andrew Dawson]

For horizontal regridding it would be nice if there were a variety of schemes available, as Susan already mentioned. I often need to regrid from high to low horizontal resolution which is best done using an area averaging scheme, rather than linear interpolation which in worst case could simply sub-sample the grid.

[Andy Hartley]

+1 for more aggregation options, when regridding from high to low spatial resolutions. Simple statistics such as mean, sd, median, min, etc I often find are more useful than the various interpolation methods available.

[Prince]

Hi Richard,

It would be ideal to have functions to regrid data on
rectilinear, curvilinear, and unstructured "

source" grids, to "destination" grids of any of these types. Most of these
capacities are developed as part of the Earth System Modeling
Framework-ESMF (http://www.earthsystemmodeling.org/). They also have a
python interface which is quite nice.

For reference, NCAR Command Language (NCL) has
implemented this with some examples here (
http://www.ncl.ucar.edu/Applications/ESMF.shtml).

Cheers,
Prince

[Andrew Dawson]

Most of these capacities are developed as part of the Earth System Modeling
Framework-ESMF (http://www.earthsystemmodeling.org/

ESMF has also been adopted by CDAT/UV-CDAT (a competitor of iris in some sense I suppose) for regridding tasks. The functionality it provides is really very useful. Using ESMF would of course add a new dependency, but it could always be optional I suppose...

[PaulEarnshaw]

I would like to second Andy Hartley's +1 for aggregation methods. These are very useful when going from high to low resolution and are used a lot in the Met Office. It is also the best way of interpolating/aggregating data that is very patchy in nature, e.g. precipitation.

I have attempted to use the Iris aggregation method to do this type of regridding (Note: this requires PR349 or Iris 1.3 when released)

def map_points(oldcoord, newcoord):
    """
    Routine to create mapping from fine resolution to course for area averaging using aggregated_by method
    """
    oldind = np.arange(oldcoord.points.size)
    newind = np.arange(newcoord.points.size)
    if oldcoord.circular:
        mapind = np.around((oldind.astype(float) * newind.size)/oldind.size).astype(int)
    else:
        mapind = np.around((oldind.astype(float) * (newind.size-1))/(oldind.size-1)).astype(int)
    if oldcoord.circular:
        mapind = np.fmod(mapind, newind.size)
    mapcoord = iris.coords.AuxCoord(newcoord.points[mapind], units=newcoord.units, coord_system=newcoord.coord_system)
    mapcoord.rename(newcoord.name())
    return mapcoord

def area_average(cube, aggfunc, samplepoints):
    """
    Calculate crude box area average using aggregated_by method
    - Useful for interpolating from fine to course resolution
    - Does not allow for intersecting grid boxes, or areas of boxes.
   
    aggfunc is one of Iris' aggregator objects
    samplepoints is the list of new grid points cf. iris.analyis.interpolate.linear
    """
    aggcube = cube.copy()
    for (coordname, points) in samplepoints:
       
        if aggcube.coords(coordname):
            aggcube.coord(coordname).rename('old_%s'%coordname)
            oldcoord = aggcube.coord('old_%s'%coordname)
            newcoord = iris.coords.AuxCoord(points, units=oldcoord.units, coord_system=oldcoord.coord_system)
            newcoord.rename(coordname)
            aggcube.add_aux_coord(map_points(oldcoord, newcoord), aggcube.coord_dims(oldcoord))
            # This line requires the fudged iris
            aggcube = aggcube.aggregated_by(coordname, aggfunc)
            aggcube = swap_aux_dim(aggcube, coordname)
   
    return aggcube

[Andrew Munslow]

Hi,

ECMWF grib file output is sometimes in Gaussian Grid (GG) spatial coordinate system.

http://www.ecmwf.int/products/data/technical/gaussian/

I think the requirement would be to convert from Gaussian Grid to regular lat/lon.

Kind regards
Andrew

[RHattersley]

Hi,

When regridding from high to low resolution, I would expect the mean to include an area weighting. Similarly for the standard deviation, variance, and median/percentile. But is there an corresponding modification to min/max (e.g. a minimum area threshold)?

Regards,
Richard

[RHattersley]

Thank you everyone for all the feedback. :-) It's been hugely helpful in setting scope and prioritising.

I do have one further favour to ask though which will help us ensure what we create is actually useful to you...
 1. What performance criteria do you have? e.g. For regridding method X, how much data do you need regridded (what resolution and dimensionality) and in what time frame?
 2. What other non-functional requirements do you have?

Kind regards,
Richard

[LeonH]

Hi Richard

I am a little late in submitting this, but I never thought that the regridding wouldn't support missing data! However, my code needs to put ocean temperature data onto the ocean velocity grid, but when I do it the masked numbers in the array (the continents) are replaced with their numerical values and some interpolation between that and the data inbetween. Not good. No warnings either and nothing in the documentation saying what to do if you have missing data...

Thanks
Leon

[pp-mo]

Hi Leon

Unfortunately Richard is away this week.
However, I think you possibly shouldn't be getting this problem :  Which regridding function are you using, and how are the "missing" points represented?
(If it gets involved, perhaps this discussion should be in a separate thread ...)

Regards
Patrick

[LeonH]

Hi Patrick

I used iris.analysis.interpolate.regrid(original_cube, destination_cube)
It simply interpolates all values that exist in the grid, including creating values where there are missing data in destination_cube by interpolating the missing data value and other surrounding values. So the mask disappears and the returned cube doesn't look like the destination_cube at all!

Leon

[bblay]

Confirmed.
Here's the a test we should add when we fix it:

    def test_regrid_masked(self):
        cs = iris.coord_systems.GeogCS(654321)
       
        src = iris.cube.Cube(np.ma.masked_less(np.arange(25, dtype='f').reshape((5,5)), 10))
        src.add_dim_coord(DimCoord(range(5), "latitude", coord_system=cs), 0)
        src.add_dim_coord(DimCoord(range(5), "longitude", coord_system=cs), 1)
       
        dst = iris.cube.Cube(np.ndarray((2,2)))
        dst.add_dim_coord(DimCoord([1,3], "latitude", coord_system=cs), 0)
        dst.add_dim_coord(DimCoord([1,3], "longitude", coord_system=cs), 1)
       
        r = iris.analysis.interpolate.regrid(src, dst)
        self.assertCMLApproxData(r, ('regrid', 'masked.cml'))

Update: Raised in https://github.com/SciTools/iris/issues/58

[bblay]

Thanks to pp-mo for pointing out the above example uses extrapolation. For data with a hole in the middle:

src = iris.cube.Cube(np.ma.MaskedArray(np.arange(25, dtype='f').reshape((5,5)), mask=np.zeros((5,5))))
src.data.mask[1:-1, 1:-1] = True

[bblay]

Apologies for the fragmented posting. Adding this:

src.data.data[1,1] = 99

highlights the problem better, as the masked 99 is returned in r.

[rsignell]

Richard & Co,

Just wondering how the regridding work is progressing -- which methods are being attacked first, what the timeline might be, how to contribute, etc.

I think I've seen discussion in this thread and others the need for methods that do:
- horizontal gridding  (structured-to-structured, unstructured-to-structured, interpolation to a specified set of points)
- time interpolation
- vertical interpolation (interpolate from dimensionless vertical coordinates to a specific height or pressure vertical coordinate (e.g. slice 5 m below surface in a sigma-layer ocean model).
- combinations of horizontal, vertical and time (interpolation onto an arbitrary (X,Y,Z) slice plane through a 3D mesh of model output,  interpolation onto a 4D (x,y,z,t) ocean glider, AUV, balloon or aircraft track, etc.)

I imagine that some folks in the Iris community have experience with these methods in other packages that might be useful to tap into.

It would be great to have a blog post or update to the roadmap as the path becomes clearer.

Thanks!
Rich

[Bernd Becker]

Hello,

I want to read 4 km resolution pp files from UKV and interpolate/re-grid them to
25 km N512 resolution.

Supposedly this works by use of IDL pp-regrid.

How can I do it using IRIS?

Who has an example to share?

Much obliged,
Bernd.

[MichelleMc]

Hi Andrew, 

I am trying to interpolate my data from high to low resolution using iris.analysis.interpolation.regrid where I have my data in one grid which I want to interpolate to an N96 grid so that I can use it in the UM. In what way do you interpolate and have you ever had any problems/had to work with masked data on this as well?

To give you a better insight into what I am doing, I am trying to create an SST file with which to force my model (HadGEM3 N96) however my data is on a different grid to the model. I also have masked data which when plotted has SSTs around the globe as the same colour and multi-colours around the coastline. I am not sure if you have come across this or have any experience with this but any experience you have had with interpolation in iris would be great

Thanks 
Michelle

[Richard Hattersley]

Hi Michelle,

I can't speak for Andrew, but when going from high to low resolution my first inclination would be to use area-weighted regridding, as described in the user guide. This includes support for masked data via a tolerance parameter.

The strange colours in your plot could be caused by "masked" values actually being represented by a large "fill value". When the colour range extends to cover the fill value (e.g. -9999) the genuine fluctuations in SST are no longer visible. The weird colours around the coast could be caused by interpolating between the fill value and genuine values. Before doing any interpolating/plotting, make sure your fill values have been converted to a real mask, e.g. with masked_data = numpy.ma.masked_equal(data, fill_value).

Regards,
Richard

[Andrew Dawson]

In addition to Richard's comments, you will have to be careful to make sure the mask of your interpolated cube matches the land-sea mask your UM configuration is using. If they don't match you will see errors in the model run.

You can use the mdtol parameter Richard mentioned to make sure you allow some missing data, to avoid masking whole low resolution cells because a few high resolution cells within it have missing values (there tends to be more land in higher resolution models). In all likelihood you wont manage to generate an interpolated filed whose land-mask matches the model land-mask, and you may have to try and create an interpolated field with more ocean points than you might want (high value of mdtol) and then apply the real model land-mask.