1
votes

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.

1

1 Answers

1
votes

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.