Copernicus Dem xarray selecciona la latitud/longitud más cercana con coordenadas multidimensionales.
Esta pregunta está relacionada, pero de alguna manera todavía necesito ayuda para que funcione.
https://stackoverflow.com/questions/58758480/xarray-select-nearest-lat-lon-with-multi-dimension-coordinates
import rioxarray
import numpy as np
import geopandas as gpd
import cartopy.crs as ccrs
# descargar y leer los datos de elevación (aproximadamente 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")
# ahora deseo encontrar la elevación en las siguientes coordenadas:
this_lon = -112.23425
this_lat = 36.3566
# puedo obtener una elevación cercana redondeando las coordenadas:
xds.loc[dict(x=-112.2, y=36.4)].values
# array([2708.229], dtype=float32)
# pero dado que los datos tienen una cuadrícula de 30 metros, debería
# poder ser más preciso que redondear las coordenadas
# Si uso las coordenadas exactas, obtengo un error porque no están
# en los índices:
xds.loc[dict(x=-112.23425, y=36.3566)].values
# KeyError: -112.23425
He intentado usar cartopy, pero esto falla:
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
La documentación menciona que “Las instancias de Copernicus DEM están disponibles en Coordenadas Geográficas; el datum horizontal de referencia es el World Geodetic System 1984 (WGS84-G1150; EPSG 4326)”, pero no sé cómo utilizar esta información.
davy.ai
El problema aquí es que las coordenadas x e y en el conjunto de datos xarray no están etiquetadas como latitud y longitud, sino como coordenadas cartesianas en metros. Por lo tanto, al intentar seleccionar una coordenada específica usando
xds.loc[dict(x=-112.23425, y=36.3566)]
, obtenemos un error de clave porque estas coordenadas exactas no existen en los índices x e y.Para seleccionar el valor de elevación más cercano a una latitud y longitud específicas, primero debemos convertir las coordenadas de latitud y longitud a las coordenadas cartesianas en metros que están presentes en el conjunto de datos. Podemos usar la biblioteca
pyproj
para hacer esto, como se muestra en el siguiente ejemplo:En este código, primero definimos la proyección del conjunto de datos utilizando el atributo
rio
del conjunto de datos xarray. También definimos la proyección de las coordenadas deseadas como EPSG:4326 (WGS84). Luego, usamos la clasepyproj.Transformer
para crear un transformador que pueda convertir las coordenadas de longitud y latitud a la proyección del conjunto de datos. Finalmente, usamos el métodosel
del conjunto de datos xarray para seleccionar el valor de elevación más cercano a las coordenadas transformadas.Se debe tener en cuenta que usamos el argumento
method='nearest'
para seleccionar el valor que está más cerca de las coordenadas transformadas. Esto es necesario porque los índices x e y no son matrices continuas, por lo que no podemos usar el métodointerp
para interpolar los valores de elevación entre los puntos de la cuadrícula. En cambio, seleccionamos el punto más cercano y asumimos que el valor de elevación en ese punto es representativo del área circundante. Si se necesitan valores de elevación más precisos, se debe utilizar un conjunto de datos diferente con una mayor resolución espacial.