Copernicus Dem xarray select nearest lat/lon with multi-dimension coordinates

Question

This question is related, but somehow I still need some help to get this to work. xarray select nearest lat/lon with multi-dimension coordinates

import rioxarray
import numpy as np
import geopandas as gpd
import cartopy.crs as ccrs

# download and read elevation data (about 40MB)
xds = rioxarray.open_rasterio("https://elevationeuwest.blob.core.windows.net/copernicus-dem/COP30_hh/Copernicus_DSM_COG_10_N36_00_W113_00_DEM.tif")

# now I wish to find the elevation at the following coordinates:
this_lon = -112.23425
this_lat = 36.3566

# I can get elevation nearby by rounding the coordinates:
xds.loc[dict(x=-112.2, y=36.4)].values
# array([2708.229], dtype=float32)

# but since the data has a 30 meters grid, I should be able 
# to be more precise than rounding the coordinates
# If I use the exact coordinates I get an error since they are 
# not in the indexes:
xds.loc[dict(x=-112.23425, y=36.3566)].values
# KeyError: -112.23425

I have tried using cartopy, but this fails:

data_crs = ccrs.LambertConformal(central_longitude=-100)
x, y = data_crs.transform_point(-112.23425, 36.3566, src_crs=ccrs.PlateCarree())
xds.sel(x=x, y=y)
# KeyError: -1090022.066606806

The documentation mentions that "The Copernicus DEM instances are available in Geographic Coordinates; the horizontal reference datum is the World Geodetic System 1984 (WGS84-G1150; EPSG 4326)", but I do not know how to use this information.

You can think of [WGS84](https://en.m.wikipedia.org/wiki/World_Geodetic_System) as simply lat/lon - it’s a particular definition of a reference ellipsoid representing the earth but for most uses that’s what you’re looking for when using lat long data. Also, take a look at the `nearest` and `tolerance` arguments to [xarray.Dataset.sel](https://xarray.pydata.org/en/stable/user-guide/indexing.html#nearest-neighbor-lookups) — Michael Delgado, Jan 03 '22 at 15:26

Timothée HENRY · Accepted Answer · 2021-12-28T17:45:22.740

Actually the longitude and latitude are evenly space at 1/3600 degrees. We can see that with:

xds.x.values
array([-113.        , -112.99972222, -112.99944444, ..., -112.00083333,
       -112.00055556, -112.00027778])

(xds.x.values + 113)*3600
array([0.000e+00, 1.000e+00, 2.000e+00, ..., 3.597e+03, 3.598e+03,
   3.599e+03])

Hence we can find quickly the nearest points in the Copernicus DEM indexes as follows:

i = round((-112.23425 + 113) * 3600)
j = round((36.3566 - 36) * 3600)

lon_nearest = xds.x.values[i]
lat_nearest = xds.y.values[3600-j]

xds.loc[dict(x=lon_nearest, y=lat_nearest)].values

Copernicus Dem xarray select nearest lat/lon with multi-dimension coordinates

1 Answers1