I have a netcdf file that has the following (example) variables:
latitude
longitude
temperature
The dimensions are in [x, y] (pixel coordinates) primarily because both the latitude and longitude are both 2-D arrays of an irregular grid.
I want to extract pixel values of temperature in eg. 53.55, 3.5 (lat/lon degree decimal coordinates). Typically, for a 1-D array of latitude/longitude, I would be able to use numpy.argmin() to find the index of both lat/lon and thus, the pixel coordinate for the temperature variable.
Alternatively, in xarray, I could use eg.
import xarray as xr
ds = open_dataset(some_file)
ds.close()
ds.temperature.sel(lat=53.55, lon=3.5, method='nearest')
Except my dimensions are now in (x, y). Perhaps due to an insufficient knowledge of the netcdf data format, but I have trouble coming up with ways to extract the data I need. How can I extract the pixel value I need? Let me know if I can clarify the question better.
You can still use argmin() if you first calculate the distance of each (2D) grid point to your requested location. Small example:
import xarray as xr
import numpy as np
f = xr.open_dataset('tauu.his.NETHERLANDS.NL_W.DOWA_40h12tg2_fERA5_WF2019_fix.201601.nc')
# 2D latitude and longitude
lons = f['lon'].values
lats = f['lat'].values
# Goal latitude and longitude
lon = 4.
lat = 52.
# Calculate distance (in degrees..) for all grid points
distance = np.sqrt( (lons-lon)**2 + (lats-lat)**2 )
# `argmin` on a 2D (or 3D, 4D, ..) field gives the index in the flattened array:
ji = distance.argmin()
# "Unravel" the index:
j,i = np.unravel_index(ji, lons.shape)
print(lons[j,i], lats[j,i]) # Prints: 3.989169 52.00158
In this case, I simply used the Euclidean distance in degrees. You can always replace that by something more fancy or accurate, like e.g. the distance in spherical coordinates.