How can we extract the maximum value from the neighborhood value of a specified coordinate?

Viewed 118

How can I extract the maximum value from four points in the neighborhood of a specified coordinate?

import xarray as xr 
import numpy as np

lat = [0, 10, 20]
lon = [50, 60, 70, 80]

#sample data
test_data = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12]])

#to xarray
data_xarray = xr.DataArray(test_data, dims=("lat","lon"), coords={"lat":lat, "lon":lon})

#<xarray.DataArray (lat: 3, lon: 4)>
#array([[ 1,  2,  3,  4],
#       [ 5,  6,  7,  8],
#       [ 9, 10, 11, 12]])
#Coordinates:
#  * lat      (lat) int64 0 10 20
#  * lon      (lon) int64 50 60 70 80

data_xarray.plot()

What I want to implement When 5.5 and 52 are specified for lat and lon respectively, extract 10, the maximum value of the four surrounding points.

1 Answers

Labeled indexing using sel supports nearest neighbour lookup. You can use that to look up the four values of interest, reconcatenate them and then compute the max:

lat_search = 5.5
lon_search = 52

# Select the four nearest values
llat_llon = data_xarray.sel(lat=lat_search, lon=lon_search, method="pad")
ulat_ulon = data_xarray.sel(lat=lat_search, lon=lon_search, method="backfill")
ulat_llon = data_xarray.sel(lat=lat_search, method="backfill").sel(
    lon=lon_search, method="pad"
)
llat_ulon = data_xarray.sel(lat=lat_search, method="pad").sel(
    lon=lon_search, method="backfill"
)

# Combine the four values providing them in the correct order
ds_grid = [[llat_llon, ulat_llon], [llat_ulon, ulat_ulon]]
neighbours = xr.combine_nested(ds_grid, concat_dim=("lon", "lat"))

# Alternatively, combine them automatically
neighbours = xr.combine_by_coords(
    [
        x.to_dataset(name="foo").expand_dims(["lat", "lon"])
        for x in [llat_llon, ulat_llon, llat_ulon, ulat_ulon]
    ]
)

# Compute the maximum value
neighbours.max()

I admit that selecting the four values manually and recombining them is not very elegant (particularly if you would like to scale that to more than two dimensions). I don't see a general way to retrieve both neighbours at the same time using sel.

If you have a regularly spaced grid of coordinates, you can select all neighbours at the same time passing a slice to sel:

delta_lat = 10
delta_lon = 10
neighbours = data_xarray.sel(
    lat=slice(lat_search - delta_lat, lat_search + delta_lat),
    lon=slice(lon_search - delta_lon, lon_search + delta_lon),
)
Related