Potential nearest neighbour interpolation error

[Duncan Watson-Parris]

I've been updating our code to use the new interpolation interface available through cube.interpolate() but I've hit an interesting issue and was wondering if one of the devs could comment.

We have a few unit tests which checks that interpolation of a given point on a cube returns the correct value, and one of those checks the values when the points near on the bounds of one (or more) of the cube's coordinates. It appears the behaviour has changed in the new interface (compared to the nearest_neighbour_data_value) and I'm not sure it seems correct.

Running the following snippet:

import numpy as np
from iris.cube import Cube
from iris.coords import DimCoord
from iris.analysis import Nearest
from iris.analysis.interpolate import nearest_neighbour_data_value

latitude = DimCoord(np.arange(-10, 11, 5), var_name='lat', standard_name='latitude', units='degrees')
longitude = DimCoord(np.arange(-5, 6, 5), var_name='lon', standard_name='longitude', units='degrees')
data = np.reshape(np.arange(15)+1.0,(5,3))
cube = Cube(data, dim_coords_and_dims=[(latitude, 0), (longitude, 1)], var_name='dummy')

sample_points = [[['latitude',2.51],['longitude',2.5]],[['latitude',-2.51],['longitude',2.5]],
                 [['latitude',2.5],['longitude',-2.51]], [['latitude',-2.51],['longitude',-2.51]]]

for point in sample_points:
    print 'New interface: {}'.format(cube.interpolate(point,scheme=Nearest()).data)
    print 'Original: {}'.format(nearest_neighbour_data_value(cube, point))

Should give the following output for iris 1.8.0:

New interface: 8.0
Original: 11.0
New interface: 8.0
Original: 5.0
New interface: 8.0
Original: 7.0
New interface: 8.0
Original: 4.0

As you can see the new interface isn't respecting the bounds of the coordinates (which should be at 2.5, and -2.5). Interestingly, if I set the points to be at 3.0 rather than 2.51 the new interface agrees, so perhaps this is a rounding issue?

[Phil Elson]

Hi Duncan,

> perhaps this is a rounding issue?

This was the clue I needed. I was able to take your code (which made finding the problem so, so much easier, thank you!) and get the expected values by ensuring that the coordinates were floating point, not integers. I suspect the old interface didn't care about the dtype, whereas the new algorithm honours the dtype correctly.

Simply changing the arange command to include a couple of decimal places does the job. See https://gist.github.com/pelson/8d6bf77805aadbda52a1.

HTH,

Phil

--
You received this message because you are subscribed to the Google Groups "Iris" group.
To unsubscribe from this group and stop receiving emails from it, send an email to scitools-iri...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.