[GIS] Indexing coordinates in order to get a value by coordinates from xarray

netcdfpythonspatio-temporal-dataxarray

Firstly, I am a super rookie and just recently started working with GIS. Therefore I am not sure how to ask/word the question, or even if it makes sense. So bear with me.

I have some NetCDF files, consisting some spatiotemporal data. When I import one of them, it looks like this:

<xarray.Dataset>
Dimensions:      (bnds: 2, time: 31, x: 720, y: 938)
Coordinates:
    lon          (y, x) float64 ...
    lat          (y, x) float64 ...
  * time         (time) datetime64[ns] 1999-01-01T11:30:00 ... 1999-01-31T11:30:00
Dimensions without coordinates: bnds, x, y
Data variables:
    time_bnds    (time, bnds) datetime64[ns] ...
    datum        (time) float64 ...
    temperature  (time, y, x) float32 ...
Attributes:
    CDI:            Climate Data Interface version 1.7.0 (http://mpimet.mpg.d...
    Conventions:    CF-1.4
    source:         surface and satellite observations, cosmo_090213_4.8_clm17
    institution:    Deutscher Wetterdienst
    title:          Temperature daily gridded dataset
    project_id:     TRY-advancement
    realization:    v1.0
    contact:        Stefan Kraehenmann, stefan.kraehenmann@dwd.de
    creation_date:  2016-01-22 21:35:27
    CDO:            Climate Data Operators version 1.7.0 (http://mpimet.mpg.d...
    history:        Wed Jun 22 09:48:12 2016: ncatted -a history,global,d,c, ...

What I want to do with this data is, I would like to call a function with parameters latitude and longitude, and get the temperature of that point. However as far as I understood, .sel() function can not help me since coordinates are only indexed(?) on time, not lat and long, from what I can see from the (*) sign near the coordinate time.

How would I somehow index it on lat-long as well, or attain the same functionality from some other function?

You can find example nc files here: https://opendata.dwd.de/climate_environment/CDC/grids_germany/daily/Project_TRY/air_temperature_mean/

I can just get the lat and lon from a point/cell, but can not do otherwise. I would like to use nearest neighbor or some similar method to retrieve a value (in this case, temp) of a point, given the coordinates.

Best Answer

You can use xarray.where. See: http://xarray.pydata.org/en/stable/indexing.html

xds = xarray.open_dataset("TT_199501_daymean.nc")
xds.where((xds.lon==5.8252) & (xds.lat==46.9359), drop=True)

enter image description here

UPDATE:

You can find the nearest neighbor using a KDTree. See: https://stackoverflow.com/questions/10818546/finding-index-of-nearest-point-in-numpy-arrays-of-x-and-y-coordinates

Then, you can select the point using .isel(...).

UPDATE: https://github.com/xarray-contrib/xoak

Related Question