viernes, 17 de junio de 2016

Procesando EarthEnv-DEM90 en Python

EarthEnv-DEM90 es un dataset de elevaciones que cubre ~91% del planeta con una resolución aproximada de 90 metros.  Fue derivado de los productos CGIAR-CSI SRTM v4.1 y ASTER GDEM v2, y se encuentra descrito en Robinson et al. (2014).

Convenientemente, la descarga no llega a los 18GB (distribuida en 2016 bloques de 5°×5°).  Sin embargo, si pensamos descomprimirlo para procesarlo calculo que puede llegar a los 140GB… Asi que luego de buscar un poco encontre que GDAL permite cargar los tarball directamente y desde luego, la muy querida libreria rasterio hereda esta habilidad. Así que procesar estos archivos resulta bastante sencillo, tal y como se muestra en el siguiente script:


Haciendo algunas pruebas encontré que en los bordes de cada bloque rasterio retorna indices fuera de rango y tuve que corregirlos (lineas 11-14).  Esto no parece ser un problema de rasterio y se debe seguramente a que los bloques no inician exactamente en los valores indicados por el nombre del archivo, por ejemplo: el bloque "S05W075" no inicia realmente en (-5, -75) sino mas bien en (-5.0000000000000036, -74.99999999999997).

Por completitud, aquí esta la salida del script:

Everest: 8797 m a.s.l.
Aconcagua: 6869 m a.s.l.
Kilimanjaro: 5846 m a.s.l.
Edge case (off-by-one indices): 270 m a.s.l.

Alternativamente, si descomentamos la linea 12, veremos el error en los indices: aunque deberian estar en el rango (0, 5999), uno de los indices se desborda por encima y el otro por debajo:

Everest: 8797 m a.s.l.
Aconcagua: 6869 m a.s.l.
Kilimanjaro: 5846 m a.s.l.
==> [off-by-one indices: 6000, -1] <==
Edge case (off-by-one indices): 270 m a.s.l.

Finalmente, para no dejar el post sin una imagen, el circulo blanco en la siguiente imagen representa mi ubicación al momento de escribir :)  La imagen incluye un shapefile de la linea costera, que puede descargarse desde Natural Earth y utiliza el Basemap Toolkit de Matplotlib.

Usted esta en "EarthEnv-DEM90_N05W080" :)