GRIB: nearest gridpoint

First we ensure the example file is available.

[1]:

import earthkit.data as ekd
from earthkit.geo import nearest_point_haversine, nearest_point_kdtree, GeoKDTree

[2]:

ekd.download_example_file("test.grib")
ds = ekd.from_source("file", "test.grib")

Using the KDTree method

We define a reference point and get the index of the nearest gridpoint from the first field.

[3]:

latlon = ds[0].to_latlon()
lat = latlon["lat"]
lon = latlon["lon"]

p_ref = (51.45, -0.97)
# distance is in m
idx, distance = nearest_point_kdtree(p_ref, (lat, lon))
idx, distance

[3]:

(array([102]), array([218417.94491761]))

With the resulting index we can get the value at the nearest gridpoint.

[4]:

v = ds[0].values[idx]
v

[4]:

array([280.68066406])

The same method works with multiple reference points.

[5]:

p_ref = [
    [51.45, 44.49, 50.73],
    [-0.97, 18.34, -17.1]
]
idx, distance = nearest_point_kdtree(p_ref, (lat, lon))
idx

[5]:

array([102, 144, 116])

[6]:

v = ds[0].values[idx]
v

[6]:

array([280.68066406, 296.4765625 , 285.05761719])

Using a KDTree object

The problem with nearest_point_kdtree() is that the KDTree is rebuilt at each call (since it is not cached at the moment). This can be a costly operation for large grids. To overcome this difficulty we can create a GeoKDTree object and call the nearest point computations directly on it multiple times.

[7]:

tree = GeoKDTree(lat, lon)
idx, distance = tree.nearest_point(p_ref)
idx, distance

[7]:

(array([102, 144, 116]),
 array([218417.94491761, 120066.38079566, 235683.22276184]))

[8]:

v = ds[0].values[idx]
v

[8]:

array([280.68066406, 296.4765625 , 285.05761719])

[9]:

p_ref = (44.49,18.34)
idx, distance = tree.nearest_point(p_ref)
idx, distance

[9]:

(array([144]), array([120066.38079566]))

[10]:

v = ds[0].values[idx]
v

[10]:

array([296.4765625])

Using the Haversine method

The nearest gridpoint can also be determined with the Haversine method, which is a “brute-force” approach since it computes all the distances between the reference point and the other points using the Haversine distance formula, then finds the points with the shortest distance.

The nearest_point_haversine() method can be used exactly in the same way as nearest_point_kdtree().

We define a reference point and get the index of the nearest gridpoint from the first field.

[11]:

latlon = ds[0].to_latlon()
lat = latlon["lat"]
lon = latlon["lon"]

p_ref = (51.45, -0.97)
# distance is in m
idx, distance = nearest_point_haversine(p_ref, (lat, lon))
idx, distance

[11]:

(array([102]), array([218417.94491761]))

With the resulting index we can get the value at the nearest gridpoint.

[12]:

v = ds[0].values[idx]
v

[12]:

array([280.68066406])

This time we use 3 reference points.

[13]:

p_ref = [
    [51.45, 44.49, 50.73],
    [-0.97, 18.34, -17.1]
]
idx, dist = nearest_point_haversine(p_ref, (lat, lon))
idx

[13]:

array([102, 144, 116])

[14]:

v = ds[0].values[idx]
v

[14]:

array([280.68066406, 296.4765625 , 285.05761719])