selecting nearest lat/lon in dataset · general

This is probably really simple, but I am reading in a new dataset, and am unfamiliar with how to extract locations with lat/lon in the way I usually do with model output. I am using xarray.

xarray.Dataset
Dimensions: x: 761y: 761z: 30nbounds: 2
Coordinates:
x (x) float64 -3.04e+06 -3.032e+06 ... 3.04e+06
y (y) float64 -3.04e+06 -3.032e+06 ... 3.04e+06
lon (y, x) float64
lat (y, x) float64
z (z) float64 -30.0 -90.0 ... -1.71e+03 -1.77e+03
z_bnds (z, nbounds) float64
Data variables:
temperature (z, y, x) float64
Attributes: (0)

I want to find, say, temp at all depths at a the location nearest to a specified lat/lon.

The typical things I use don't work (either give error, or nonsense). For example:

temps.temperature.sel(lat=-65, lon=-110, method = 'nearest').isel(z=0)
or
xylar_temps.temperature.where(xylar_temps.lat==-65, xylar_temps.lon==-110).isel(z=0)

Philip Chmielowiec (Dec 20 2023 at 11:08):

You mentioned you are working on a Cartesian grid, so attempting to select using spherical (lat, lon) coordinates might be causing this unexpected behavior.

Philip Chmielowiec (Dec 20 2023 at 11:10):

I would start off by converting your lat and lon into their respective x y and z Cartesian representation and see if that gives your expected results when selecting.

Philip Chmielowiec (Dec 20 2023 at 11:11):

Katelyn FitzGerald (Dec 20 2023 at 16:56):

Mira Berdahl (Dec 20 2023 at 18:31):

@Philip Chmielowiec
The min/max range for lat is : -90 to -51.84956701
for lon, the range is -179.9 to 180.
(The grid is centered over Antarctica roughly).

I'm also not sure what you mean by converting the lat/lon onto the cartesian grid. They are already given in terms of the (y,x) grid.

Philip Chmielowiec (Dec 20 2023 at 18:32):

Got it, I've been working with unstructured grids mostly, so sometimes instead of lat lon our coordinates are in terms of x y z, so just wanted to make sure :)

Daniel Adriaansen (Dec 20 2023 at 18:57):

Hi Mira. What does your code return when this is used:
temps.temperature.sel(lat=-65, lon=-110, method = 'nearest').isel(z=0)

Mira Berdahl (Dec 20 2023 at 18:58):

Hi @Daniel Adriaansen
I get the following error:
KeyError: "no index found for coordinate 'lat'"

Daniel Adriaansen (Dec 20 2023 at 19:08):

I believe because lat/lon are not dimension coordinates, using sel() will not work in this case. Do your data have a specific projection on the Earth?

Mira Berdahl (Dec 20 2023 at 19:16):

I don't believe they are projected, its a regular grid.
Basically I want to be able to do the opposite of this, where isntead of selecting an x/y to get a lat/lon, I want to be able to do the opposite:

temps.temperature.isel(x=100,y=100, z=0)

xarray.DataArray'temperature'
array(0.379751)
Coordinates:
x () float64 -2.24e+06
y () float64 -2.24e+06
lon () float64 -135.0
lat () float64 -61.43
z () float64 -30.0
Attributes:
units :
degree_celsius
long_name :
corrected pot. temp

Daniel Adriaansen (Dec 20 2023 at 19:27):

OK thanks! If there are two-dimensional latitude/longitude arrays I think it should have some projection unless it is a latitude/longitude projection and instead of 1D lat/lon you have 2D lat/lon for some reason but I've never seen this and I think your dimensions would be lat/lon and not x/y. I'd double check the file header to see if there is grid information in there. Regardless, take a look here: https://stackoverflow.com/questions/58758480/xarray-select-nearest-lat-lon-with-multi-dimension-coordinates. There are several approaches mentioned, but the only thing I will caution is that simply taking a minimum difference between the latitude and longitude arrays may not give you the correct result if you do have projected data. I could take a look at a file if you're able to share, or provide header information here.

Mira Berdahl (Dec 20 2023 at 19:54):

Thanks very much, I'll check out the link you provided. If it helps, the header for the file is :

netcdf obs_temperature_1995-2017_8km_x_60m {
dimensions:
        x = 761 ;
        y = 761 ;
        z = 30 ;
        nbounds = 2 ;
variables:
        double x(x) ;
                x:_FillValue = nan ;
                x:units = "meters" ;
        double y(y) ;
                y:_FillValue = nan ;
                y:units = "meters" ;
        double temperature(z, y, x) ;
                temperature:_FillValue = nan ;
                temperature:units = "degree_celsius" ;
                temperature:long_name = "corrected pot. temp" ;
                temperature:coordinates = "lon lat" ;
        double lon(y, x) ;
                lon:_FillValue = nan ;
                lon:units = "degrees" ;
        double lat(y, x) ;
                lat:_FillValue = nan ;
                lat:units = "degrees" ;
        double z(z) ;
                z:_FillValue = nan ;
                z:units = "meters" ;
                z:bounds = "z_bnds" ;
                z:standard_name = "depth" ;
                z:positive = "up" ;
                z:axis = "Z" ;
        double z_bnds(z, nbounds) ;
                z_bnds:_FillValue = nan ;
                z_bnds:comment = "depth bounds" ;

// global attributes:
                :coordinates = "z_bnds" ;

Daniel Adriaansen (Dec 21 2023 at 16:53):

You're right, I see no grid information in there. Do you have any information about the model or observation dataset that this file came from?

Mira Berdahl (Dec 21 2023 at 23:08):

I believe I figured out a way, with help from the stackoverflow site in the thread above:

#  plot the verticies of each point
temps.temperature.isel(z=0).plot(x='lon', y='lat')

# I want to find the temp at a certain lat/lon point.
mylat = -65
mylon = -110

# First, find the index of the grid point nearest a specific lat/lon.
abslat = np.abs(temps.lat-mylat)
abslon = np.abs(temps.lon-mylon)
c = np.maximum(abslon, abslat)

([xloc], [yloc]) = np.where(c == np.min(c))

# Now I can use that index location to get the values at the x/y diminsion
temp_at_myLoc = temps.isel(x=yloc,y=xloc)
temp_at_myLoc

Stream: general

Topic: selecting nearest lat/lon in dataset

Mira Berdahl (Dec 20 2023 at 01:18):