Lat/long in 2d array - nearest pixel indexing

Question

My dataset xyz_data = xr.open_dataset("xyz.nc") has following properties:

xarray.Dataset
Dimensions:
time: 2967x: 61y: 91
Coordinates:
lat
(x, y)
float32
...
lon
(x, y)
float32
...
time
(time)
datetime64[ns]
2021-12-01 ... 2021-12-31T23:45:00
Data variables:
xyz_probability
(time, x, y)
int8
...
Indexes:
time
PandasIndex

I am using this code to find the index for the pixel nearest to my location of interest.

lat = xyz_data['lat']
lon = xyz_data['lon']

target_lon = 8.117500
target_lat = 48.59111

lon_diff = np.abs(lon - target_lon)
lat_diff = np.abs(lat - target_lat)

nearest_lon_idx = np.unravel_index(lon_diff.argmin(), lon.shape)
nearest_lat_idx = np.unravel_index(lat_diff.argmin(), lat.shape)

print(nearest_lon_idx)
print(nearest_lat_idx)

# The  output of print command is (9, 21) and (25, 25) 

xyz_filtered = xyz_data.isel(time=slice(None), x= [nearest_lon_idx[0]], y=[nearest_lat_idx[1]])

Problem: xyz_filtered doesn't give me the co-ordinates nearest to my target location. I suppose it has to do with indexing but I can't solve it. Suggestions would be much appreciated!

I also followed https://gis.stackexchange.com/questions/357026/indexing-coordinates-in-order-to-get-a-value-by-coordinates-from-xarray thread but I get an error which says:

ValueError: Cannot apply_along_axis when any iteration dimensions are 0

Suggestions would be much appreciated!

Answer 1

I went through your code an i realized that you are finding the minimum absolute difference for each latitude and longitude independently. This won't guarantee the minimum combined distance between the target and the grid point.

What you can do is:

1. You can calculate the euclidean distance between the target location an each grid point and find the minimum index of the combined distance.

   import numpy as np
   
   lon_diff = lon - target_lon
   lat_diff = lat - target_lat
   dist = np.sqrt(lon_diff**2 + lat_diff**2)   
   
   nearest_idx = np.unravel_index(dist.argmin(), dist.shape)
   

2. You can also use both nearest_idx[0] and nearest_idx[1] to index the filtered dataset:

   xyz_filtered = xyz_data.isel( time=slice(None), x=[nearest_idx[0]], y=[nearest_idx1] )

3. But what i would recommand you is that instead of manually calculating distances, you sshould use apply_ufunc of xarray

your code should contain this snipped:

def euclidean_distance(lon1, lat1, lon2, lat2):
    return np.sqrt((lon1 - lon2)**2 + (lat1 - lat2)**2)

xyz_filtered = xyz_data.apply_ufunc(
    euclidean_distance,
    input_vars=["lon", "lat"],
    kwargs={"lon2": target_lon, "lat2": target_lat},
    vectorize=True,
).argmin(dim=("x", "y"))

xyz_filtered = xyz_data.isel(time=slice(None), x=xyz_filtered.x, y=xyz_filtered.y)

This approach uses vectorized operations, directly finding the index of the minimum distance.