GRIB: nearest gridpoint

First we ensure the example file is available.

from earthkit.geo.distance import GeoKDTree, nearest_point_haversine, nearest_point_kdtree

import earthkit.data as ekd

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

Using the KDTree method

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

lat, lon = ds[0].geography.latlons()

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

(array([56]), array([175096.60480163]))

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

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

array([285.29272461])

The same method works with multiple reference points.

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

array([56, 73, 53])

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

array([285.29272461, 297.87768555, 284.27905273])

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.

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

(array([56, 73, 53]),
 array([175096.60480163, 142841.04317786, 169632.89688725]))

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

array([285.29272461, 297.87768555, 284.27905273])

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

(array([73]), array([142841.04317786]))

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

array([297.87768555])

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.

latlon = ds[0].geography.latlons()

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

(array([56]), array([175096.60480163]))

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

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

array([285.29272461])

This time we use 3 reference points.

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

array([56, 73, 53])

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

array([285.29272461, 297.87768555, 284.27905273])