Calculate a value at specific latitude and longitude via interpolation

[mcrimaldi]

I am quite new to using iris so I apologize for the newbie question.
I'm trying to interpolate GRIB temperature data at a certain point with known latitude and longitude.

For example, I've downloaded a GRIB file for air temperature over 2 days, hourly, converted in celsius degrees:

air_temperature / (celsius)         (time: 48; latitude: 21; longitude: 31)
    Dimension coordinates:
        time                             x             -              -
        latitude                         -             x              -
        longitude                        -             -              x
    Auxiliary coordinates:
        forecast_period                  x             -              -
    Scalar coordinates:
        height                      2 m
        originating_centre          European Centre for Medium Range Weather Forecasts

Now, I made a daily mean of temperature, so the new cube has 2 values of temperature:

iris.coord_categorisation.add_day_of_year(temperature, "time")
daily_Tmean = temperature.aggregated_by(["day_of_year"], iris.analysis.MEAN)

The new temperature values are calculated over the entire area (latitude and longitude) of original GRIB file:

temperature_iter = daily_Tmean.slices(("latitude", "longitude"))
ani = animate.animate(temperature_iter, qplt.contourf)
plt.show()

Now, I'm trying to calculate a new temperature value at given latitude and longitude interpolating the temperatures already calculated as daily mean. What I am trying to achieve is something like "resizing" the pixel resolution.
I did this before in MatLab and R applying Delanauy triangulation (https://it.mathworks.com/help/matlab/ref/delaunay.html)

With iris I tried this:

ws_coords = [("latitude", 40.94), ("longitude", 14.02)]
t_atCoords = daily_Tmean.interpolate(ws_coords, iris.analysis.Nearest())

But when I print the calculated temperature values I have this (no Delanauy triangulation so far):

print(t_atCoords.data)
[-- --]

I am quite sure there is something I did not understand or doing wrong, Any help or suggestions would be very much appreciated. Many thanks.

[TonyP-BAS]

As Laura Dreyer has already mentioned in response to this question in the google group, the answer is that your GRIB data are masked to land points only - that is why your contourf plot has the white regions. I'm guessing that your data are ERA5-Land and that you're looking at a point on the coast of Italy.

I have some daily ERA5-Land global data, and if I load a day of that I get a cube I'm calling t2m, which looks like this:

In [24]: print(t2m)
2 metre temperature / (K) (latitude: 1801; longitude: 3600)
Dimension coordinates:
latitude x -
longitude - x

I can find the indices of the point closest to your point of interest like this (there may be better ways):

In [30]: x = t2m.coord(axis='x')

In [31]: np.where(np.abs(x.points - 14.02) == np.abs(x.points - 14.02).min())
Out[31]: (array([140]),)

In [32]: y = t2m.coord(axis='y')

In [33]: np.where(np.abs(y.points - 40.94) == np.abs(y.points - 40.94).min())
Out[33]: (array([491]),)

Which tells me that the closest point is at x index 140 and y index 491. I can then print the surrounding few points from the data array:

In [34]: t2m.data[489:494, 138:143]
Out[34]:
masked_array(
data=[[--, --, 280.2789001464844, 279.6180114746094, 278.9607238769531],
[--, --, 281.61517333984375, 280.8266296386719,
279.7364807128906],
[--, --, --, 282.213134765625, 280.9557800292969],
[--, --, --, --, --],
[--, --, --, --, --]],
mask=[[ True, True, False, False, False],
[ True, True, False, False, False],
[ True, True, True, False, False],
[ True, True, True, True, True],
[ True, True, True, True, True]],
fill_value=1e+20,
dtype=float32)

Which does indeed show that that point (the one in the middle of the extracted region) is indeed masked but has non-masked neighbours. So there is no issue with your code, but you do need to choose a point for which the field has data.

In passing, I don't think that the interpolation will perform a triangulation in any case: I think it will just select the point with the closest coordinate values to your specified point (sort of like I did to find the indices above), essentially treating the field's latitude and longitude coordinates as Cartesian coordinates (with appropriate handling of longitude wrap-around).

I hope this is useful.

[trexfeathers]

How would I go about calculating the temperature datum at my point using the data from the neighboring points? (the datum in the center of the neighboring pixels).

I thought the below was evidence that everything was working properly:

If I change the coordinates to some that fall within the GRIB's coverage, as you pointed out, the calculation happens and the interpolation returns me non-masked values.

Are you still experiencing some problems even when you stay in the un-masked area? Or are you saying that you want to be able to use interpolation in the masked area (I don't think that would be possible)?

Therefore, considering the GRIB as an array of points and not as a raster (the measured points are in the center of the pixels)...

Interpolation should work fine on point data, as far as I'm aware.

[mcrimaldi]

Sorry for the extreme delay I completely missed the notification.

I've managed the interpolation of data and It looks like it's working.

Just a little thing that I don't understand is why when I try to plot the results, the map looks strange:

%matplotlib notebook

temperature_iter_ws_samplepoints = t_atCoords2.slices(("latitude", "longitude"))

ani = animate.animate(temperature_iter_ws_samplepoints, qplt.contourf)

plt.show()

Checking the cube data, I have:

print(t_atCoords2)

air_temperature / (celsius) (time: 2; -- : 18; -- : 18)
Dimension coordinates:
time x - -
Auxiliary coordinates:
day_of_year x - -
forecast_period x - -
latitude - x -
longitude - - x
Scalar coordinates:
height 2 m
originating_centre European Centre for Medium Range Weather Forecasts
Cell methods:
mean day_of_year

And It looks like that the coordinates of my points are stored as "Auxiliary coordinates" instead of Dimension coordinates (the 18 masked (?) ones).
I am assuming that the strange plot depends on this, so, how can I plot the new calculated (interpolated) temperatures the right way?

[trexfeathers]

-- on a Cube summary simply states that the dimension isn't 'named' - because it doesn't have an associated Dimension Coordinate - it is unrelated to masking. Iris is fine to pick up Auxiliary Coordinates when plotting, which is how your plot has been created.

I expect that the whitespace is around masked data points. If you want to confirm this you will need:

for i in range(2):
    print(t_atCoords2.data[i, ...])

I have lost track of why these points might be masked. How did you generate t_atCoords2?

[mcrimaldi]

t_atCoords2 are interpolated temperature values at some sample points I have.
More in detail:

These are the coordinates of my sample points (lat and lon). The values of coordinates (ws_latitude and ws_longitudes) are numpy arrays
ws_samplepoints = [("latitude", ws_latitudes), ("longitude", ws_longitudes)]

then I calculate the interpolation of temperature at these points:
t_atCoords2 = daily_Tmean.interpolate(ws_samplepoints, iris.analysis.Nearest(extrapolation_mode="extrapolate"))

print(t_atCoords2)
air_temperature / (celsius)         (time: 3; -- : 18; -- : 18)
    Dimension coordinates:
        time                             x       -        -
    Auxiliary coordinates:
        day_of_year                      x       -        -
        forecast_period                  x       -        -
        latitude                         -       x        -
        longitude                        -       -        x
    Scalar coordinates:
        height                      2 m
        originating_centre          European Centre for Medium Range Weather Forecasts
    Cell methods:
        mean                        day_of_year

"time" coordinate is 3 because the temperature has been aggregated daily (3 days, the temperature at each day is equal to the daily mean)

With these t_atCoords2 values I have the plot above.

for i in range(2):
    print(t_atCoords2.data[i, ...])

returns:

[[-- 29.231072998046898 28.998223876953148 -- 28.677341715494816 --
  29.341455078125023 26.975325520833355 -- 28.677341715494816 --
  29.341455078125023 -- -- 28.998223876953148 28.070245361328148 --
  29.096032714843773]
 [28.260237630208355 28.053999837239605 28.87707926432294 --
  28.396234130859398 27.648411051432316 28.567478434244816
  31.64790242513023 27.992242431640648 28.396234130859398 --
  28.567478434244816 28.75251668294273 28.694563802083355
  28.87707926432294 29.975315348307316 28.75251668294273
  28.99664713541669]
 [28.260237630208355 28.053999837239605 28.87707926432294 --
  28.396234130859398 27.648411051432316 28.567478434244816
  31.64790242513023 27.992242431640648 28.396234130859398 --
  28.567478434244816 28.75251668294273 28.694563802083355
  28.87707926432294 29.975315348307316 28.75251668294273
  28.99664713541669]

...continues

By the way, I am going to open another question here because I am trying something else with some cubes math.
In the meanwhile, thanks for helping me.

[trexfeathers]

That confirms that the blank spaces in the map are indeed coming from masked data. As for why it would appear in that pattern: I am less sure. Doesn't seem to match up with your source data (based on the original plot). Having not worked closely with interpolation I can't say why; perhaps @stephenworsley can comment.

[trexfeathers]

Please re-open if you would like further help with this 🙂