Copernicus Dem xarray selecciona la latitud/longitud más cercana con coordenadas multidimensionales.

29 Jul, 2023 Programación 0

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

July 29, 2023 at 1:21 pm

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:

import rioxarray
import numpy as np
import pyproj
import geopandas as gpd

# descargar y leer datos de elevación (aproximadamente 40 MB)
xds = rioxarray.open_rasterio("https://elevationeuwest.blob.core.windows.net/copernicus-dem/COP30_hh/Copernicus_DSM_COG_10_N36_00_W113_00_DEM.tif")

# definir la proyección del conjunto de datos
dataset_proj = pyproj.CRS.from_dict(xds.rio.crs)

# definir la proyección de las coordenadas deseadas (WGS84)
coord_proj = pyproj.CRS.from_epsg(4326)

# definir las coordenadas de longitud y latitud
this_lon = -112.23425
this_lat = 36.3566

# crear un transformador pyproj para convertir de coord_proj a dataset_proj
transformer = pyproj.Transformer.from_crs(coord_proj, dataset_proj)

# transformar las coordenadas de longitud y latitud a la proyección del conjunto de datos
x, y = transformer.transform(this_lon, this_lat)

# seleccionar el valor de elevación más cercano a las coordenadas transformadas
elevation = xds.sel(x=x, y=y, method='nearest').values

print(elevation)
# salida: array([2708.229], dtype=float32)

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 clase pyproj.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étodo sel 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étodo interp 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.

Answer

davy.ai

July 29, 2023 at 1:21 pm
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:
```
import rioxarray
import numpy as np
import pyproj
import geopandas as gpd

# descargar y leer datos de elevación (aproximadamente 40 MB)
xds = rioxarray.open_rasterio("https://elevationeuwest.blob.core.windows.net/copernicus-dem/COP30_hh/Copernicus_DSM_COG_10_N36_00_W113_00_DEM.tif")

# definir la proyección del conjunto de datos
dataset_proj = pyproj.CRS.from_dict(xds.rio.crs)

# definir la proyección de las coordenadas deseadas (WGS84)
coord_proj = pyproj.CRS.from_epsg(4326)

# definir las coordenadas de longitud y latitud
this_lon = -112.23425
this_lat = 36.3566

# crear un transformador pyproj para convertir de coord_proj a dataset_proj
transformer = pyproj.Transformer.from_crs(coord_proj, dataset_proj)

# transformar las coordenadas de longitud y latitud a la proyección del conjunto de datos
x, y = transformer.transform(this_lon, this_lat)

# seleccionar el valor de elevación más cercano a las coordenadas transformadas
elevation = xds.sel(x=x, y=y, method='nearest').values

print(elevation)
# salida: array([2708.229], dtype=float32)
```
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 clase pyproj.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étodo sel 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étodo interp 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.

Comments are closed.

es.davy.ai

¿Tienes una pregunta?

Copernicus Dem xarray selecciona la latitud/longitud más cercana con coordenadas multidimensionales.

Answer

davy.ai