diff --git a/book/espanol/00_Introduction_Setup/00_Getting_NASA_EarthData_Credentials.md b/book/espanol/00_Introduction_Setup/00_Getting_NASA_EarthData_Credentials.md new file mode 100644 index 0000000..18c7aa1 --- /dev/null +++ b/book/espanol/00_Introduction_Setup/00_Getting_NASA_EarthData_Credentials.md @@ -0,0 +1,19 @@ +# Cómo obtener las credenciales de NASA Earthdata + +En este cuaderno se describe el proceso para obtener las credenciales de NASA Earthdata. + +**Puedes completar este paso antes del día que se dicte el tutorial.** + +## Breve introducción + +El programa [Sistemas de datos científicos de la Tierra (ESDS, por sus siglas en inglés) de la NASA](https://www.earthdata.nasa.gov/) supervisa el ciclo de vida de los datos de ciencias de la Tierra de todas sus misiones de observación terrestre, desde su adquisición hasta su procesamiento y distribución. + +Para fines de esta guía, el sitio web NASA Earthdata es el punto de entrada que permite el acceso completo, gratuito y abierto a las colecciones de datos de ciencias de la Tierra de la NASA, con el fin de acelerar el progreso científico en beneficio de la sociedad. Para acceder a los datos mediante este portal, quienes deseen usarlos deben definir primero sus credenciales de acceso. Para crear una cuenta en EarthData, sigue estos pasos: + +- Ingresa al sitio web de EarthData de la NASA: [`https://www.earthdata.nasa.gov/`](https://www.earthdata.nasa.gov/). Luego, selecciona la opción "_Use Data_" en el menú, y luego "_Register_". Finalmente, ingresa a la siguiente dirección [`https://urs.earthdata.nasa.gov/`](https://urs.earthdata.nasa.gov/). + +![earthdata\_login](../assets/earthdata_login.png) (earthdata_inicio_de_sesión) + +- Selecciona la opción "_Register for a profile_" ("_Registrarse para obtener un perfil_") y elige un nombre de usuario y una contraseña. Necesitarás estos datos más adelante, así que elige aquellos que recuerdes bien. También tendrás que completar tu perfil para finalizar el registro, en el mismo se te pedirá información como correo electrónico, país y afiliación. Por último, elige "_Register for Earthdata Login_" ("_Registrarse para iniciar sesión en Earthdata_"). + +![earthdata\_perfil](../assets/earthdata_profile2.png) diff --git a/book/espanol/00_Introduction_Setup/01_Using_the_2i2c_Hub.md b/book/espanol/00_Introduction_Setup/01_Using_the_2i2c_Hub.md new file mode 100644 index 0000000..f035ab2 --- /dev/null +++ b/book/espanol/00_Introduction_Setup/01_Using_the_2i2c_Hub.md @@ -0,0 +1,82 @@ +# Uso del Hub de 2i2c + +Este cuaderno computacional contiene las instrucciones para iniciar sesión en la nube con ([JupyterHub](https://jupyter.org/hub)) plataforma proporcionada por [2i2c](https://2i2c.org) para este tutorial. + +**No podrás completar este paso hasta el día que inicie el tutorial (ese día recibirás la contraseña).** + +## 1. Iniciar sesión en el Hub de 2i2c + +Para iniciar sesión en el JupyterHub proporcionado por 2i2c, sigue estos pasos: + +1. **Navega hasta el Hub de 2i2c:** Tu navegador web debe apuntar a [este enlace] (https://climaterisk.opensci.2i2c.cloud). + +2. **Inicia sesión con tus credenciales:** + +- **Nombre de usuario:** Puedes elegir cualquier nombre de usuario que desees. Sugerimos que utilices tu nombre de usuario de GitHub para evitar conflictos. +- **Contraseña:** _Recibirás la contraseña el día que inicie el tutorial_. + +![2i2c\_login](../assets/2i2c_login.png) (2i2c_inicio_de_sesión) + +3. **Inicio de sesión:** + +El proceso de inicio de sesión puede tardar unos minutos, especialmente si es necesario crear un nuevo espacio de trabajo virtual solo para ti. + +![iniciar\_servidor2](../assets/start_server_2i2c.png) + +- **Qué esperar:** + +De manera predeterminada, al iniciar sesión en [`https://climaterisk.opensci.2i2c.cloud`](https://climaterisk.opensci.2i2c.cloud) se clonará automáticamente un repositorio para trabajar. Si el inicio de sesión es exitoso, verás la siguiente pantalla y estarás listo para empezar a trabajar. + +![entorno\_de\_trabajo\_jupyter\_lab](../assets/work_environment_jupyter_lab.png) + +**Notas:** Cualquier archivo en el que trabajes se mantendrá entre sesiones siempre y cuando uses el mismo nombre de usuario al iniciar sesión. + +## 2. Configuración del entorno en la nube para acceder a NASA EarthData desde Python + +Para acceder a los productos de NASA EarthData desde programas de Python o cuadernos computacionales de Jupyter, es necesario que guardes tus credenciales de NASA EarthData en un archivo especial llamado .netrc en tu directorio de inicio. + +- Puedes crear este archivo ejecutando en una terminal el script llamado `make_netrc.py`: + + ```bash + $ python make_netrc.py + ``` + + También puedes ejecutar este script dentro de este cuaderno de Jupyter ejecutando la celda Python siguiente (utilizando el comando mágico %run). + + Algunas advertencias: + + - El script no se ejecutará si ~/.netrc ya existe. Puedes borrar ese archivo o renombrarlo si quieres conservar las credenciales que contiene. + - El script solicitará tu nombre de usuario y contraseña de NASA EarthData, así que está atento si lo ejecutas desde un cuaderno de Jupyter. + +```bash +%run make_netrc.py +``` + +- De manera alternativa, puedes crear un archivo que se llame .netrc en tu carpeta de inicio (es decir, ~/.netrc) con el siguiente contenido: + ``` + machine urs.earthdata.nasa.gov login USERNAME password PASSWORD + ``` + Por supuesto, debes reemplazar el `USERNAME` y la `PASSWORD` en el archivo .netrc con los datos reales de tu cuenta NASA EarthData. Una vez que hayas guardado el archivo `.netrc` con tus credenciales correctas, es una buena práctica restringir el acceso al archivo: + ```bash + $ chmod 600 ~/.netrc + ``` + +## 3. Verificando el acceso a los productos de NASA EarthData + + + +Para asegurarte de que todo funciona correctamente, ejecuta el script `test_netrc.py`: + +```bash +$ python test_netrc.py +``` + +Una vez más, puedes ejecutarlo directamente desde un cuaderno computacional utilizando la celda de Python que se muestra a continuación: + + + +```bash +%run test_netrc.py +``` + +Si funcionó sin problemas, ¡ya está todo listo! ¡Ahora tienes todo lo que necesitas para explorar los datos de observación de la Tierra de la NASA mediante el portal EarthData! diff --git a/book/espanol/03_Using_NASA_EarthData/01_Using_OPERA_DIST_Products.md b/book/espanol/03_Using_NASA_EarthData/01_Using_OPERA_DIST_Products.md new file mode 100644 index 0000000..4506caa --- /dev/null +++ b/book/espanol/03_Using_NASA_EarthData/01_Using_OPERA_DIST_Products.md @@ -0,0 +1,338 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: "1.3" + jupytext_version: 1.16.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Utilización de los productos OPERA DIST + + + +## El proyecto OPERA + + + + + +
+ +
+ +Del proyecto (Observational Products for End-Users from Remote Sensing Analysis)](https://www.jpl.nasa.gov/go/opera) (Productos de observación para usuarios finales a partir del análisis por teledetección): + +> Iniciado en abril del 2021, el proyecto OPERA (Observational Products for End-Users from Remote Sensing Analysis) (Productos de observación para usuarios finales a partir del análisis por teledetección) del Jet Propulsion Laboratory recopila datos de instrumentos ópticos y de radar por satélite para generar seis conjuntos de productos: +> +> - un conjunto de productos sobre la extensión de las aguas superficiales a escala casi mundial +> - un conjunto de productos de alteración de la superficie casi mundial +> - un producto radiométrico corregido del terreno casi global +> - un conjunto de productos complejos Coregistered Single Look de Norteamérica +> - un paquete de productos de North America Displacement +> - un conjunto de productos de North America Vertical Land Motion + +Es decir, OPERA es una iniciativa de la NASA que toma, por ejemplo, datos de teledetección óptica o radar recopilados desde satélites, y genera una variedad de conjuntos de datos preprocesados para uso público. Los productos de OPERA no son imágenes de satélite sin procesar, sino el resultado de una clasificación algorítmica para determinar, por ejemplo, qué regiones terrestres contienen agua o dónde se ha desplazado la vegetación. Las imágenes de satélite sin procesar se recopilan a partir de mediciones realizadas por los instrumentos a bordo de las misiones de los satélites Sentinel-1 A/B, Sentinel-2 A/B y Landsat-8/9 (de ahí el término _HLS_ para "_Harmonized Landsat-Sentinel_" en numerosas descripciones de productos). + + + +--- + +## El producto OPERA Land Surface Disturbance (DIST) + + + +Uno de estos productos de datos OPERA es el producto _Land Surface Disturbance (DIST)_ (descrito con más detalle en la [OPERA DIST HLS product specification](https://d2pn8kiwq2w21t.cloudfront.net/documents/OPERA_DIST_HLS_Product_Specification_V1.pdf)) (especificación del producto OPERA DIST HLS). +El producto DIST cartografía la alteración de la vegetación (específicamente, la pérdida de cobertura vegetal por píxel HLS siempre que haya una disminución indicada) a partir de escenas del Harmonized Landsat-8 y el Sentinel-2 A/B (HLS). Una de las aplicaciones de estos datos es cuantificar los daños causados por los _incendios forestales_. El producto DIST_ALERT se publica a intervalos regulares (la misma cadencia de las imágenes HLS, aproximadamente cada 12 días en un determinado mosaico/región). El producto DIST_ANN resume las mediciones de las alteraciones a lo largo de un año. + +Los productos DIST cuantifican los datos de reflectancia de la superficie (SR) adquiridos a partir del Operational Land Imager (OLI) a bordo del satélite de teledetección Landsat-8 y del Multi-Spectral Instrument (MSI) a bordo del satélite de teledetección Sentinel-2 A/B. Los productos de datos HLS DIST son archivos de datos rasterizados, cada uno de ellos asociado a mosaicos de la superficie terrestre. Cada mosaico se representa mediante coordenadas cartográficas proyectadas alineadas con el [_Military Grid Reference System (MGRS)_](https://en.wikipedia.org/wiki/Military_Grid_Reference_System) (Sistema de Referencia de Cuadrículas Militares (MGRS)). Cada mosaico cubre 109.8 $km^2$ divididos en 3660 filas y 3660 columnas, con una separación de 30 metros entre píxeles, con mosaicos que se solapan con los vecinos en 4900 metros en cada dirección (los detalles se describen completamente en la [DIST product specification](https://d2pn8kiwq2w21t.cloudfront.net/documents/OPERA_DIST_HLS_Product_Specification_V1.pdf)) (especificación DIST del producto). + +Los productos OPERA DIST se distribuyen como [Cloud Optimized GeoTIFFs](https://www.cogeo.org/) (GeoTIFF optimizados para la nube). En la práctica, eso significa que las diferentes bandas se almacenan en archivos TIFF distintos. La especificación TIFF permite el almacenamiento de matrices multidimensionales en un único archivo. El almacenamiento de bandas distintas en archivos TIFF distintos permite que los archivos se descarguen de forma independiente. + + + +--- + + + +## Banda 1: Valor máximo de la anomalía de pérdida de vegetación (VEG_ANOM_MAX) + + + + + +Examinemos un archivo local con un ejemplo de datos DIST-ALERT. El archivo contiene la primera banda de datos de alteración: la _anomalía de pérdida máxima de vegetación_. Para cada píxel, se trata de un valor entre 0% y 100% que representa la diferencia porcentual entre la cobertura vegetal que se observa actualmente y un valor de referencia histórico. Es decir, un valor de 100 corresponde a una pérdida completa de vegetación en un píxel y un valor de 0 corresponde a ninguna pérdida de vegetación. Los valores de los píxeles se almacenan como enteros sin signo de 8 bits (UInt8) porque los valores de los píxeles solo deben oscilar entre 0 y 100. Un valor del píxel de 255 indica que faltan datos, es decir, que los datos HLS no pudieron destilar un valor máximo de anomalía en la vegetación para ese píxel. Por supuesto, el uso de datos enteros sin signo de 8 bits es mucho más eficiente para el almacenamiento y para la transmisión de datos a través de una red (en comparación con, por ejemplo, datos de punto flotante de 32 o 64 bits). + +Empecemos importando las bibliotecas necesarias. Observa que también estamos importando algunas clases de la biblioteca Bokeh para hacer que los gráficos interactivos sean un poco más atractivos. + + + +```python editable=true jupyter={"source_hidden": true} slideshow={"slide_type": ""} +# Notebook dependencies +import warnings +warnings.filterwarnings('ignore') +from pathlib import Path +import rioxarray as rio +import hvplot.xarray +from bokeh.models import FixedTicker +import geoviews as gv +gv.extension('bokeh') +``` + + + +Vamos a leer los datos de un archivo local `'OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-ANOM-MAX.tif'`. Antes de cargarlo, analicemos los metadatos incrustados en el nombre del archivo. + + + +```python jupyter={"source_hidden": true} +LOCAL_PATH = Path('..') / 'assets' / 'OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-ANOM-MAX.tif' +filename = LOCAL_PATH.name +print(filename) +``` + + + +Este nombre de archivo bastante largo incluye varios campos separados por caracteres de subrayado (`_`). Podemos utilizar el método `str.split` de Python para ver más fácilmente los distintos campos. + + + +```python jupyter={"source_hidden": true} +filename.split('_') # Use the Python str.split method to view the distinct fields more easily. +``` + + + +Los archivos del producto OPERA tienen un esquema de nombres particular (tal como se describe en la [DIST product specification](https://d2pn8kiwq2w21t.cloudfront.net/documents/OPERA_DIST_HLS_Product_Specification_V1.pdf) (especificación del producto DIST)). En la salida anterior, podemos extraer ciertos metadatos para este ejemplo: + +1. _Product_: `OPERA`; +2. _Level_: `L3` ; +3. _ProductType_: `DIST-ALERT-HLS` ; +4. _TileID_: `T10TEM` (cadena que hace referencia a un mosaico del [Military Grid Reference System (MGRS)](https://en.wikipedia.org/wiki/Military_Grid_Reference_System)) (Sistema de Referencia de Cuadrículas Militares (MGRS)); +5. _AcquisitionDateTime_: `20220815T185931Z` (cadena que representa una marca de tiempo GMT para la adquisición de los datos); +6. _ProductionDateTime_ : `20220817T153514Z` (cadena que representa una marca de tiempo GMT para cuando se generó el producto de los datos); +7. _Sensor_: `S2A` (identificador del satélite que adquirió los datos sin procesar: `L8` (Landsat-8), `S2A` (Sentinel-2 A) o `S2B` (Sentinel-2 B); +8. _Resolution_: `30` (por ejemplo, píxeles de longitud lateral $30\mathrm{m}$); +9. _ProductVersion_: `v0.1`; y +10. _LayerName_: `VEG-ANOM-MAX` + +Ten en cuenta que la NASA utiliza nomenclaturas convencionales como [Earthdata Search](https://search.earthdata.nasa.gov) para extraer datos significativos de los [SpatioTemporal Asset Catalogs (STAC)](https://stacspec.org/) ( Catálogos de Activos Espaciales y Temporales (STAC)). Más adelante utilizaremos estos campos, en particular los campos _TileID_ y _LayerName_, para filtrar los resultados de la búsqueda antes de recuperar los datos remotos. + + + + + +Vamos a subir los datos de este archivo local en un Xarray `DataArray` utilizando `rioxarray.open_rasterio`. Reetiquetaremos las coordenadas adecuadamente y extraeremos el CRS (sistema de referencia de coordenadas). + + + +```python jupyter={"source_hidden": true} +data = rio.open_rasterio(LOCAL_PATH) +crs = data.rio.crs +data = data.rename({'x':'easting', 'y':'northing', 'band':'band'}).squeeze() +``` + +```python jupyter={"source_hidden": true} +data +``` + +```python jupyter={"source_hidden": true} +crs +``` + + + +Antes de generar un gráfico, vamos a crear un mapa base utilizando mosaicos [ESRI](https://en.wikipedia.org/wiki/Esri). + + + +```python jupyter={"source_hidden": true} +# Creates basemap +base = gv.tile_sources.ESRI.opts(width=1000, height=1000, padding=0.1) +``` + + + +También utilizaremos diccionarios para capturar la mayor parte de las opciones de trazado que utilizaremos más adelante junto con `.hvplot.image`. + + + +```python jupyter={"source_hidden": true} +image_opts = dict( + x='easting', + y='northing', + rasterize=True, + dynamic=True, + frame_width=500, + frame_height=500, + aspect='equal', + cmap='hot_r', + clim=(0, 100), + alpha=0.8 + ) +layout_opts = dict( + xlabel='Longitude', + ylabel='Latitude' + ) +``` + + + +Por último, utilizaremos el método `DataArray.where` para filtrar los píxeles que faltan y los píxeles que no han tenido cambios en la vegetación, y modificaremos las opciones en `image_opts` y `layout_opts` a valores particulares para este conjunto de datos. + + + +```python jupyter={"source_hidden": true} +veg_anom_max = data.where((data>0) & (data!=255)) +image_opts.update(crs=data.rio.crs) +layout_opts.update(title=f"VEG_ANOM_MAX") +``` + + + +Esto nos permitirá generar un gráfico significativo. + + + +```python jupyter={"source_hidden": true} +veg_anom_max.hvplot.image(**image_opts).opts(**layout_opts) * base +``` + + + +En el gráfico resultante, los píxeles blancos y amarillos corresponden a regiones en las que se ha producido cierta deforestación, pero no mucha. Por el contrario, los píxeles oscuros y negros corresponden a regiones que han perdido casi toda la vegetación. + + + +--- + +## Banda 2: Fecha de alteración inicial de la vegetación (VEG_DIST_DATE) + + + +Los productos DIST-ALERT contienen varias bandas (tal como se resume en la [DIST product specification](https://d2pn8kiwq2w21t.cloudfront.net/documents/OPERA_DIST_HLS_Product_Specification_V1.pdf)) (especificación del producto DIST). La segunda banda que analizaremos es la _fecha de alteración inicial de la vegetación_ en el último año. Esta se almacena como un número entero de 16 bits (Int16). + +El archivo `OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-DIST-DATE.tif` se almacena localmente. La [DIST product specification](\(https://d2pn8kiwq2w21t.cloudfront.net/documents/OPERA_DIST_HLS_Product_Specification_V1.pdf\)) (especificación del producto DIST) describe cómo utilizar las convenciones para la denominación de archivos. Aquí destaca la _fecha y hora de adquisición_ `20220815T185931`, por ejemplo, casi las 7 p.m. (UTC) del 15 de agosto del 2022. + +Cargaremos y reetiquetaremos el `DataArray` como antes. + + + +```python jupyter={"source_hidden": true} +LOCAL_PATH = Path('..') / 'assets' / 'OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-DIST-DATE.tif' +data = rio.open_rasterio(LOCAL_PATH) +data = data.rename({'x':'easting', 'y':'northing', 'band':'band'}).squeeze() +``` + + + +En esta banda en particular, el valor 0 indica que no ha habido alteraciones en el último año y -1 es un valor centinela que indica que faltan datos. Cualquier valor positivo es el número de días desde el 31 de diciembre del 2020 en los que se midió la primera alteración en ese píxel. Filtraremos los valores no positivos y conservaremos estos valores significativos utilizando `DataArray.where`. + + + +```python jupyter={"source_hidden": true} +veg_dist_date = data.where(data>0) +``` + + + +Vamos a examinar el rango de valores numéricos en `veg_dist_date` utilizando DataArray.min`and`DataArray.max`. Ambos métodos ignorarán los píxeles que contengan `nan\` ("Not-a-Number") al calcular el mínimo y el máximo. + + + +```python jupyter={"source_hidden": true} +d_min, d_max = int(veg_dist_date.min().item()), int(veg_dist_date.max().item()) +print(f'{d_min=}\t{d_max=}') +``` + + + +En este caso, los datos significativos se encuentran entre 247 y 592. Recuerda que se trata del número de días transcurridos desde el 31 de diciembre del 2020, cuando se observó la primera alteración en el último año. Dado que estos datos se adquirieron el 15 de agosto del 2022, los únicos valores posibles estarían entre 227 y 592 días. Así que debemos recalibrar el mapa de colores en el gráfico + + + +```python jupyter={"source_hidden": true} +image_opts.update( + clim=(d_min,d_max), + crs=data.rio.crs + ) +layout_opts.update(title=f"VEG_DIST_DATE") +``` + +```python jupyter={"source_hidden": true} +veg_dist_date.hvplot.image(**image_opts).opts(**layout_opts) * base +``` + + + +Con este mapa de colores, los píxeles más claros mostraron algunos signos de deforestación hace cerca de un año. Por el contrario, los píxeles negros mostraron deforestación por primera vez cerca del momento de adquisición de los datos. Por tanto, esta banda es útil para seguir el avance de los incendios forestales a medida que arrasan los bosques. + + + +--- + + + +## Banda 3: Estado de alteración de la vegetación (VEG_DIST_STATUS) + + + + + +Por último, veamos una tercera banda de la familia de productos DIST-ALERT denominada _estado de alteración de la vegetación_. Estos valores de píxel se almacenan como enteros de 8 bits sin signo. Solo hay 6 valores distintos almacenados: + +- **0:** Sin alteración +- **1:** Alteración provisional (**primera detección**) con cambio en la cubierta vegetal < 50% +- **2:** Alteración confirmada (**detección recurrente**) con cambio en la cubierta vegetal < 50% +- **3:** Alteración provisional con cambio en la cubierta vegetal ≥ 50% +- **4:** Alteración confirmada con cambio en la cubierta vegetal ≥ 50% +- **255**: Datos no disponibles + +El valor de un píxel se marca como cambiado provisionalmente cuando la pérdida de la cobertura vegetal (alteración) es observada por primera vez por un satélite. Si el cambio se vuelve a notar en posteriores adquisiciones HLS sobre dicho píxel, entonces el píxel se marca como confirmado. + + + + + +Podemos usar un archivo local como ejemplo de esta capa/banda particular de los datos DIST-ALERT. El código es el mismo que el anterior, pero observa que: + +- los datos filtrados reflejan los valores de píxel significativos para esta capa (por ejemplo, `data<0` and `data<5`), y +- los límites del mapa de colores se reasignan en consecuencia (es decir, de 0 a 4). + + + +```python jupyter={"source_hidden": true} +LOCAL_PATH = Path('..') / 'assets' / 'OPERA_L3_DIST-ALERT-HLS_T10TEM_20220815T185931Z_20220817T153514Z_S2A_30_v0.1_VEG-DIST-STATUS.tif' +data = rio.open_rasterio(LOCAL_PATH) +data = data.rename({'x':'easting', 'y':'northing', 'band':'band'}).squeeze() +``` + +```python jupyter={"source_hidden": true} +veg_dist_status = data.where((data>0)&(data<5)) +image_opts.update(crs=data.rio.crs) +``` + +```python jupyter={"source_hidden": true} +layout_opts.update( + title=f"VEG_DIST_STATUS", + clim=(0,4), + colorbar_opts={'ticker': FixedTicker(ticks=[0, 1, 2, 3, 4])} + ) +``` + +```python jupyter={"source_hidden": true} +veg_dist_status.hvplot.image(**image_opts).opts(**layout_opts) * base +``` + + + +Este mapa de colores continuo no resalta especialmente bien las características de este gráfico. Una mejor opción sería un mapa de colores _categórico_. Veremos cómo conseguirlo en el próximo cuaderno (con los productos de datos OPERA DSWx). + + + +--- diff --git a/book/espanol/03_Using_NASA_EarthData/03_Using_PySTAC.md b/book/espanol/03_Using_NASA_EarthData/03_Using_PySTAC.md new file mode 100644 index 0000000..e2b5bb6 --- /dev/null +++ b/book/espanol/03_Using_NASA_EarthData/03_Using_PySTAC.md @@ -0,0 +1,575 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: "1.3" + jupytext_version: 1.16.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Uso de la API PySTAC + + + +En el [Earthdata Search website](https://search.earthdata.nasa.gov) (sitio web Earthdata Search) de la NASA se puede buscar gran cantidad de datos. El enlace anterior se conecta a una interfaz gráfica de usuario para buscar en [SpatioTemporal Asset Catalogs (STACs)](https://stacspec.org/) (catálogos de activos espaciotemporales (STAC)) al especificar una área de interés (AOI) y una _ventana temporal_ o un _intervalo de fechas_. + +Por motivos de reproducibilidad, queremos ser capaces de buscar en los catálogos de activos mediante programación. Aquí es donde entra en juego la biblioteca [PySTAC](https://pystac.readthedocs.io/en/stable/). + + + +--- + +## Esquema de las etapas del análisis + + + +- Identificación de los parámetros de búsqueda + - AOI, ventana temporal + - Endpoint, proveedor, identificador del catálogo ("nombre corto") +- Obtención de los resultados de la búsqueda + - Instrospección, análisis para identificar características, bandas de interés + - Envolver los resultados en un DataFrame para facilitar la exploración +- Explorar y refinar los resultados de la búsqueda + - Identificar los gránulos de mayor valor + - Filtrar los gránulos extraños con una contribución mínima + - Reunir los gránulos filtrados correspondientes en un DataFrame + - Identificar el tipo de resultado que se quiere generar +- Procesamiento de los datos para generar resultados relevantes + - Descargar los gránulos relevantes en Xarray DataArray, apilados adecuadamente + - Llevar a cabo los cálculos intermedios necesarios + - Ensamblar los fragmentos de datos relevantes en la visualización + + + +--- + +## Identificar los parámetros de búsqueda + +### Definir el AOI y el rango de fechas + + + +Empezaremos tomando en cuenta un ejemplo concreto. [Heavy rains severely impacted Southeast Texas in May 2024](https://www.texastribune.org/2024/05/03/texas-floods-weather-harris-county/) (Fuertes lluvias afectaron gravemente al sureste de Texas en mayo del 2024), provocando [flooding and causing significant damage to property and human life](https://www.texastribune.org/series/east-texas-floods-2024/) (inundaciones y causando importantes daños materiales y humanos). + +Como es usual, se requieren ciertas importaciones relevantes. Las dos primeras celdas son familiares (relacionadas con las herramientas de análisis y la visualización de los datos que ya se examinaron). La tercera celda incluye importaciones de la biblioteca `pystac_client` y de la biblioteca `gdal`, seguidas de algunos ajustes necesarios para utilizar la [GDAL (the Geospatial Data Abstraction Library)](https://gdal.org) (Biblioteca de Abstracción de Datos Geoespaciales). Estos detalles en la configuración permiten que las sesiones de tu cuaderno interactúen sin problemas con las fuentes remotas de datos geoespaciales. + + + +```python jupyter={"source_hidden": true} +from warnings import filterwarnings +filterwarnings('ignore') +# data wrangling imports +import numpy as np +import pandas as pd +import xarray as xr +import rioxarray as rio +import rasterio +``` + +```python jupyter={"source_hidden": true} +# Imports for plotting +import hvplot.pandas +import hvplot.xarray +import geoviews as gv +from geoviews import opts +gv.extension('bokeh') +``` + +```python jupyter={"source_hidden": true} +# STAC imports to retrieve cloud data +from pystac_client import Client +from osgeo import gdal +# GDAL setup for accessing cloud data +gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/.cookies.txt') +gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/.cookies.txt') +gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR') +gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF, TIFF') +``` + + + +A continuación, definiremos los parámetros de búsqueda geográfica para poder recuperar los datos correspondientes a ese evento de inundación. Esto consiste en especificar un _área de interés (AOI)_ y un _intervalo de fechas_. + +- El AOI se especifica como un rectángulo de coordenadas de longitud-latitud en una única 4-tupla con la forma + $$({\mathtt{longitude}}_{\mathrm{min}},{\mathtt{latitude}}_{\mathrm{min}},{\mathtt{longitude}}_{\mathrm{max}},{\mathtt{latitude}}_{\mathrm{max}}),$$ + por ejemplo, las coordenadas de la esquina inferior izquierda seguidas de las coordenadas de la esquina superior derecha. +- El intervalo de fechas se especifica como una cadena de la forma + $${\mathtt{date}_{\mathrm{start}}}/{\mathtt{date}_{\mathrm{end}}},$$ + donde las fechas se especifican en el formato estándar `YYYY-MM-DD`. + + + +```python jupyter={"source_hidden": true} +# Center of the AOI +livingston_tx_lonlat = (-95.09,30.69) # (lon, lat) form +``` + + + +Escribiremos algunas funciones cortas para encapsular la lógica de nuestros flujos de trabajo genéricos. Para el código de investigación, estas se colocarían en archivos de módulos Python. Por conveniencia, incrustaremos las funciones en este cuaderno y otras para que puedan ejecutarse correctamente con dependencias mínimas. + + + +```python jupyter={"source_hidden": true} +# simple utility to make a rectangle with given center of width dx & height dy +def make_bbox(pt,dx,dy): + '''Returns bounding-box represented as tuple (x_lo, y_lo, x_hi, y_hi) + given inputs pt=(x, y), width & height dx & dy respectively, + where x_lo = x-dx/2, x_hi=x+dx/2, y_lo = y-dy/2, y_hi = y+dy/2. + ''' + return tuple(coord+sgn*delta for sgn in (-1,+1) for coord,delta in zip(pt, (dx/2,dy/2))) +``` + +```python jupyter={"source_hidden": true} +# simple utility to plot an AOI or bounding-box +def plot_bbox(bbox): + '''Given bounding-box, returns GeoViews plot of Rectangle & Point at center + + bbox: bounding-box specified as (lon_min, lat_min, lon_max, lat_max) + Assume longitude-latitude coordinates. + ''' + # These plot options are fixed but can be over-ridden + point_opts = opts.Points(size=12, alpha=0.25, color='blue') + rect_opts = opts.Rectangles(line_width=0, alpha=0.1, color='red') + lon_lat = (0.5*sum(bbox[::2]), 0.5*sum(bbox[1::2])) + return (gv.Points([lon_lat]) * gv.Rectangles([bbox])).opts(point_opts, rect_opts) +``` + +```python jupyter={"source_hidden": true} +AOI = make_bbox(livingston_tx_lonlat, 0.5, 0.25) +basemap = gv.tile_sources.OSM.opts(width=500, height=500) +plot_bbox(AOI) * basemap +``` + + + +Agreguemos un intervalo de fechas. Las inundaciones ocurrieron principalmente entre el 30 de abril y el 2 de mayo. Estableceremos una ventana temporal más larga que cubra los meses de abril y mayo. + + + +```python jupyter={"source_hidden": true} +start_date, stop_date = '2024-04-01', '2024-05-31' +DATE_RANGE = f'{start_date}/{stop_date}' +``` + + + +Por último, crearemos un un diccionario `search_params` que almacene el AOI y el intervalo de fechas. Este diccionario se utilizará para buscar datos en los STAC. + + + +```python jupyter={"source_hidden": true} +search_params = dict(bbox=AOI, datetime=DATE_RANGE) +print(search_params) +``` + +--- + +## Obtención de los resultados de búsqueda + +### Ejecución de una búsqueda con la API PySTAC + + + +Para iniciar una búsqueda de datos se necesitan tres datos más: el _Endpoint_ (una URL), el _Proveedor_ (una cadena que representa una ruta que extiende el Endpoint) y los _Identificadores de la colección_ (una lista de cadenas que hacen referencia a catálogos específicos). Generalmente, debemos experimentar con [Earthdata Search website](https://search.earthdata.nasa.gov) (sitio web Earthdata Search) de la NASA para determinar correctamente estos valores para los productos de los datos específicos que queremos recuperar. El [NASA CMR STAC GitHub repository also monitors issues](https://github.com/nasa/cmr-stac/issues) (repositorio de GitHub de la NASA CMR STAC también supervisa los problemas) relacionados con la API para las consultas de búsqueda de EarthData Cloud. + +Para la búsqueda de productos de datos DSWx que queremos ejecutar, estos parámetros son los que se definen en la siguiente celda de código. + + + +```python jupyter={"source_hidden": true} +ENDPOINT = 'https://cmr.earthdata.nasa.gov/stac' # base URL for the STAC to search +PROVIDER = 'POCLOUD' +COLLECTIONS = ["OPERA_L3_DSWX-HLS_V1_1.0"] +# Update the dictionary opts with list of collections to search +search_params.update(collections=COLLECTIONS) +print(search_params) +``` + + + +Una vez que se definieron los parámetros de búsqueda en el diccionario de Python `search_params`, podemos instanciar un `Cliente` y buscar en el catálogo espacio-temporal de activos utilizando el método `Client.search`. + + + +```python jupyter={"source_hidden": true} +catalog = Client.open(f'{ENDPOINT}/{PROVIDER}/') +search_results = catalog.search(**search_params) +print(f'{type(search_results)=}\n',search_results) +``` + + + +El objeto `search_results` que se obtuvo al llamar al método `search` es del tipo `ItemSearch`. Para recuperar los resultados, invocamos el método `items` y lanzamos el resultado como una `list` de Python que enlazaremos con los gránulos identificadores. + + + +```python jupyter={"source_hidden": true} +%%time +granules = list(search_results.items()) +print(f"Number of granules found with tiles overlapping given AOI: {len(granules)}") +``` + + + +Analicemos el contenido de la lista `granules`. + + + +```python jupyter={"source_hidden": true} +granule = granules[0] +print(f'{type(granule)=}') +``` + +```python jupyter={"source_hidden": true} +granule +``` + + + +El objeto `granule` tiene una representación de salida enriquecida en este cuaderno de Jupyter. Podemos ampliar los atributos en la celda de salida haciendo clic en los triángulos. + +![](../assets/granule_output_repr.png) + +El término _gránulo_ se refiere a una colección de archivos de datos (datos ráster en este caso), todos ellos asociados a datos sin procesar adquiridos por un satélite concreto en una fecha y hora fijas sobre un mosaico geográfico concreto. Hay una gran variedad de atributos interesantes asociados con este gránulo. + +- properties['datetime']: una cadena que representa la hora de adquisición de los datos de los archivos de datos ráster de este gránulo, +- properties['eo:cloud_cover']: el porcentaje de píxeles oscurecidos por nubes y sombras de nubes en los archivos de datos ráster de este gránulo, y +- `assets`: un `dict` de Python cuyos valores resumen las bandas o niveles de los datos ráster asociados con este gránulo. + + + +```python jupyter={"source_hidden": true} +print(f"{type(granule.properties)=}\n") +print(f"{granule.properties['datetime']=}\n") +print(f"{granule.properties['eo:cloud_cover']=}\n") +print(f"{type(granule.assets)=}\n") +print(f"{granule.assets.keys()=}\n") +``` + + + +Cada objeto en `granule.assets` es una instancia de la clase `Asset` que tiene un atributo `href`. Es el atributo `href` el que nos indica dónde localizar un archivo GeoTiff asociado con este activo de este gránulo. + + + +```python jupyter={"source_hidden": true} +for a in granule.assets: + print(f"{a=}\t{type(granule.assets[a])=}") + print(f"{granule.assets[a].href=}\n\n") +``` + + + +Además, el `Item` tiene un atributo `.id` que almacena una cadena. Al igual que ocurre con los nombres de archivo asociados a los productos OPERA, esta cadena `.id` contiene el identificador de un mosaico geográfico MGRS. Podemos extraer ese identificador aplicando manipulaciones de cadenas de Python al atributo `.id` del gránulo. Hagámoslo y almacenemos el resultado en `tile_id`. + + + +```python jupyter={"source_hidden": true} +print(granule.id) +tile_id = granule.id.split('_')[3] +print(f"{tile_id=}") +``` + +--- + +### Resumiendo los resultados de la búsqueda en un DataFrame + + + +Los detalles de los resultados de la búsqueda son complicados de analizar de esta manera. Vamos a extraer algunos campos concretos de los gránulos obtenidos en un Pandas `DataFrame` utilizando una cómoda función de Python. Definiremos la función aquí y la reutilizaremos en cuadernos posteriores. + + + +```python jupyter={"source_hidden": true} +# utility to extract search results into a Pandas DataFrame +def search_to_dataframe(search): + '''Constructs Pandas DataFrame from PySTAC Earthdata search results. + DataFrame columns are determined from search item properties and assets. + 'asset': string identifying an Asset type associated with a granule + 'href': data URL for file associated with the Asset in a given row.''' + granules = list(search.items()) + assert granules, "Error: empty list of search results" + props = list({prop for g in granules for prop in g.properties.keys()}) + tile_ids = map(lambda granule: granule.id.split('_')[3], granules) + rows = (([g.properties.get(k, None) for k in props] + [a, g.assets[a].href, t]) + for g, t in zip(granules,tile_ids) for a in g.assets ) + df = pd.concat(map(lambda x: pd.DataFrame(x, index=props+['asset','href', 'tile_id']).T, rows), + axis=0, ignore_index=True) + assert len(df), "Empty DataFrame" + return df +``` + + + +Invocar `search_to_dataframe` en `search_results` codifica la mayor parte de la información importante de la búsqueda como un Pandas `DataFrame`, tal como se muestra a continuación. + + + +```python jupyter={"source_hidden": true} +df = search_to_dataframe(search_results) +df.head() +``` + + + +El método `DataFrame.info` nos permite examinar el esquema de este DataFrame. + + + +```python jupyter={"source_hidden": true} +df.info() +``` + + + +Vamos a limpiar el DataFrame de resultados de búsqueda. Esto podría estar incrustado en una función, pero, vale la pena saber cómo se hace esto interactivamente con Pandas. + +En primer lugar, para estos resultados, solo es necesaria una columna `Datetime`, podemos eliminar las demás. + + + +```python jupyter={"source_hidden": true} +df = df.drop(['start_datetime', 'end_datetime'], axis=1) +df.info() +``` + + + +A continuación, arreglaremos el esquema del `DataFrame` `df` convirtiendo las columnas en tipos de datos sensibles. También será conveniente utilizar la marca de tiempo de la adquisición como índice del DataFrame. Hagámoslo utilizando el método `DataFrame.set_index`. + + + +```python jupyter={"source_hidden": true} +df['datetime'] = pd.DatetimeIndex(df['datetime']) +df['eo:cloud_cover'] = df['eo:cloud_cover'].astype(np.float16) +str_cols = ['asset', 'href', 'tile_id'] +for col in str_cols: + df[col] = df[col].astype(pd.StringDtype()) +df = df.set_index('datetime').sort_index() +``` + +```python jupyter={"source_hidden": true} +df.info() +``` + + + +Al final esto da un DataFrame con un esquema conciso que puede utilizarse para manipulaciones posteriores. Agrupar los resultados de la búsqueda STAC en un Pandas `DataFrame` de forma sensata es un poco complicado. Varias de las manipulaciones anteriores podrían haberse incluido en la función `search_to_dataframe`. Pero, dado que los resultados de búsqueda de la API STAC aún están evolucionando, actualmente es mejor ser flexible y utilizar Pandas de forma interactiva para trabajar con los resultados de búsqueda. Veremos esto con más detalle en ejemplos posteriores. + + + +--- + +## Explorar y refinar los resultados de búsqueda + + + +Si examinamos la columna numérica `eo:cloud_cover` del DataFrame `df`, podemos recopilar estadísticas utilizando agregaciones estándar y el método `DataFrame.agg`. + + + +```python jupyter={"source_hidden": true} +df['eo:cloud_cover'].agg(['min','mean','median','max']) +``` + + + +Observa que hay varias entradas `nan` en esta columna. Las funciones de agregación estadística de Pandads suelen ser "`nan`-aware", esto significa que ignoran implícitamente las entradas `nan` ("Not-a-Number") al calcular las estadísticas. + + + +### Filtrado del DataFrame de búsqueda con Pandas + + + +Como primera operación de filtrado, vamos a mantener solo las filas para las que la cobertura en la nube es inferior al 50%. + + + +```python jupyter={"source_hidden": true} +df_clear = df.loc[df['eo:cloud_cover']<50] +df_clear +``` + + + +Para esta consulta de búsqueda, cada gránulo DSWX comprende datos ráster para diez bandas o niveles. Podemos ver esto aplicando el método Pandas `Series.value_counts` a la columna `asset`. + + + +```python jupyter={"source_hidden": true} +df_clear.asset.value_counts() +``` + + + +Vamos a filtrar las filas que corresponden a la banda `B01_WTR` del producto de datos DSWx. El accesorio de Pandas `DataFrame.str` hace que esta operación sea sencilla. Llamaremos al `DataFrame` filtrado `b01_wtr`. + + + +```python jupyter={"source_hidden": true} +b01_wtr = df_clear.loc[df_clear.asset.str.contains('B01_WTR')] +b01_wtr.info() +b01_wtr.asset.value_counts() +``` + + + +También podemos ver que hay varios mosaicos geográficos asociados a los gránulos encontrados que intersecan el AOI proporcionado. + + + +```python jupyter={"source_hidden": true} +b01_wtr.tile_id.value_counts() +``` + + + +Recuerda que estos códigos se refieren a mosaicos geográficos MGRS especificados en un sistema de coordenadas concreto. Como identificamos estos códigos en la columna `tile_id`, podemos filtrar las filas que corresponden, por ejemplo, a los archivos recopilados sobre el mosaico T15RUQ del MGRS: + + + +```python jupyter={"source_hidden": true} +b01_wtr_t15ruq = b01_wtr.loc[b01_wtr.tile_id=='T15RUQ'] +b01_wtr_t15ruq +``` + + + +Ahora tenemos un `DataFrame` `b01_wtr_t15ruq` mucho más corto que resume las ubicaciones remotas de los archivos (por ejemplo, GeoTiffs) que almacenan datos ráster para la banda de aguas superficiales `B01_WTR` en el mosaico MGRS `T15RUQ` recopilados en varias marcas de tiempo que se encuentran dentro de la ventana de tiempo que especificamos. Podemos utilizar este DataFrame para descargar esos archivos para su análisis o visualización. + + + +--- + +## Procesamiento de datos para obtener resultados relevantes + +### Apilamiento de los datos + + + +Tenemos un `DataFrame` que identifica archivos remotos específicos de datos ráster. El siguiente paso es combinar estos datos ráster en una estructura de datos adecuada para el análisis. El Xarray `DataArray` es adecuado en este caso. La combinación puede generarse utilizando la función `concat` de Xarray. La función `urls_to_stack` en la siguiente celda es larga pero no complicada. Toma un `DataFrame` con marcas de tiempo en el índice y una columna etiquetada `href` de las URL, lee los archivos asociados a esas URL uno a uno, y apila las matrices bidimensionales relevantes de datos ráster en una matriz tridimensional. + + + +```python jupyter={"source_hidden": true} +def urls_to_stack(granule_dataframe): + '''Processes DataFrame of PySTAC search results (with OPERA tile URLs) & + returns stacked Xarray DataArray (dimensions time, latitude, & longitude)''' + + stack = [] + for i, row in granule_dataframe.iterrows(): + with rasterio.open(row.href) as ds: + # extract CRS string + crs = str(ds.crs).split(':')[-1] + # extract the image spatial extent (xmin, ymin, xmax, ymax) + xmin, ymin, xmax, ymax = ds.bounds + # the x and y resolution of the image is available in image metadata + x_res = np.abs(ds.transform[0]) + y_res = np.abs(ds.transform[4]) + # read the data + img = ds.read() + # Ensure img has three dimensions (bands, y, x) + if img.ndim == 2: + img = np.expand_dims(img, axis=0) + lon = np.arange(xmin, xmax, x_res) + lat = np.arange(ymax, ymin, -y_res) + bands = np.arange(img.shape[0]) + da = xr.DataArray( + data=img, + dims=["band", "lat", "lon"], + coords=dict( + lon=(["lon"], lon), + lat=(["lat"], lat), + time=i, + band=bands + ), + attrs=dict( + description="OPERA DSWx B01", + units=None, + ), + ) + da.rio.write_crs(crs, inplace=True) + stack.append(da) + return xr.concat(stack, dim='time').squeeze() +``` + +```python jupyter={"source_hidden": true} +%%time +stack = urls_to_stack(b01_wtr_t15ruq) +``` + +```python jupyter={"source_hidden": true} +stack +``` + +### Producir una visualización de los datos + +```python jupyter={"source_hidden": true} +# Define a colormap with RGBA tuples +COLORS = [(150, 150, 150, 0.1)]*256 # Setting all values to gray with low opacity +COLORS[0] = (0, 255, 0, 0.1) # Not-water class to green +COLORS[1] = (0, 0, 255, 1) # Open surface water +COLORS[2] = (0, 0, 255, 1) # Partial surface water +``` + +```python jupyter={"source_hidden": true} +image_opts = dict( + x='lon', + y='lat', + project=True, + rasterize=True, + cmap=COLORS, + colorbar=False, + tiles = gv.tile_sources.OSM, + widget_location='bottom', + frame_width=500, + frame_height=500, + xlabel='Longitude (degrees)', + ylabel='Latitude (degrees)', + title = 'DSWx data for May 2024 Texas floods', + fontscale=1.25 + ) +``` + + + +Trazar las imágenes totalmente puede consumir mucha memoria. Vamos a utilizar el método Xarray `DataArray.isel` para extraer un trozo del arreglo `stack` con menos píxeles. Esto permitirá un renderizado y desplazamiento rápidos. + + + +```python jupyter={"source_hidden": true} +view = stack.isel(lon=slice(3000,None), lat=slice(3000,None)) +view.hvplot.image(**image_opts) +``` + +```python jupyter={"source_hidden": true} +stack.hvplot.image(**image_opts) # Construct view from all slices. +``` + + + +Antes de continuar, recuerda apagar el kernel de este cuaderno para liberar memoria para otros cálculos. + + + +--- + + + +Este cuaderno proporciona principalmente un ejemplo para ilustrar el uso de la API PySTAC. + +En cuadernos posteriores, utilizaremos este flujo de trabajo general: + +1. Establecer una consulta de búsqueda mediante la identificación de un _AOI_ particular y un _intervalo de fechas_. +2. Identificar un _endpoint_, un _proveedor_ y un _catálogo de activos_ adecuados, y ejecutar la búsqueda utilizando `pystac.Client`. +3. Convertir los resultados de la búsqueda en un Pandas DataFrame que contenga los principales campos de interés. +4. Utilizar el DataFrame resultante para filtrar los archivos de datos remotos más relevantes necesarios para el análisis y/o la visualización. +5. Ejecutar el análisis y/o la visualización utilizando el DataFrame para recuperar los datos requeridos. + + diff --git a/book/espanol/04_Case_Studies/00_Template.md b/book/espanol/04_Case_Studies/00_Template.md new file mode 100644 index 0000000..7534e0c --- /dev/null +++ b/book/espanol/04_Case_Studies/00_Template.md @@ -0,0 +1,232 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: "1.3" + jupytext_version: 1.16.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# Plantilla para el uso de la nube EarthData + +## Esquema de los pasos para el análisis + + + +- Identificación de los parámetros de búsqueda + - AOI, ventana temporal + - Endpoint, proveedor, identificador del catálogo ("nombre corto") +- Obtención de los resultados de la búsqueda + - Instrospección, análisis para identificar características, bandas de interés + - Envolver los resultados en un DataFrame para facilitar la exploración +- Explorar y refinar los resultados de la búsqueda + - Identificar los gránulos de mayor valor + - Filtrar los gránulos extraños con una contribución mínima + - Reunir los gránulos filtrados relevantes en un DataFrame + - Identificar el tipo de resultado que se quiere generar +- Procesamiento de los datos para generar resultados relevantes + - Descargar los gránulos relevantes en Xarray DataArray, apilados adecuadamente + - Llevar a cabo los cálculos intermedios necesarios + - Ensamblar los fragmentos de datos relevantes en la visualización + + + +--- + +### Importaciones preliminares + +```python jupyter={"source_hidden": true} +from warnings import filterwarnings +filterwarnings('ignore') +# data wrangling imports +import numpy as np +import pandas as pd +import xarray as xr +import rioxarray as rio +import rasterio +``` + +```python jupyter={"source_hidden": true} +# Imports for plotting +import hvplot.pandas +import hvplot.xarray +import geoviews as gv +from geoviews import opts +gv.extension('bokeh') +``` + +```python jupyter={"source_hidden": true} +# STAC imports to retrieve cloud data +from pystac_client import Client +from osgeo import gdal +# GDAL setup for accessing cloud data +gdal.SetConfigOption('GDAL_HTTP_COOKIEFILE','~/.cookies.txt') +gdal.SetConfigOption('GDAL_HTTP_COOKIEJAR', '~/.cookies.txt') +gdal.SetConfigOption('GDAL_DISABLE_READDIR_ON_OPEN','EMPTY_DIR') +gdal.SetConfigOption('CPL_VSIL_CURL_ALLOWED_EXTENSIONS','TIF, TIFF') +``` + +### Utilidades prácticas + +Estas funciones podrían colocarse en archivos de módulos para proyectos de investigación más desarrollados. Para fines didácticos, se incluyen en este cuaderno. + +```python jupyter={"source_hidden": true} +# simple utility to make a rectangle with given center of width dx & height dy +def make_bbox(pt,dx,dy): + '''Returns bounding-box represented as tuple (x_lo, y_lo, x_hi, y_hi) + given inputs pt=(x, y), width & height dx & dy respectively, + where x_lo = x-dx/2, x_hi=x+dx/2, y_lo = y-dy/2, y_hi = y+dy/2. + ''' + return tuple(coord+sgn*delta for sgn in (-1,+1) for coord,delta in zip(pt, (dx/2,dy/2))) +``` + +```python jupyter={"source_hidden": true} +# simple utility to plot an AOI or bounding-box +def plot_bbox(bbox): + '''Given bounding-box, returns GeoViews plot of Rectangle & Point at center + + bbox: bounding-box specified as (lon_min, lat_min, lon_max, lat_max) + Assume longitude-latitude coordinates. + ''' + # These plot options are fixed but can be over-ridden + point_opts = opts.Points(size=12, alpha=0.25, color='blue') + rect_opts = opts.Rectangles(line_width=0, alpha=0.1, color='red') + lon_lat = (0.5*sum(bbox[::2]), 0.5*sum(bbox[1::2])) + return (gv.Points([lon_lat]) * gv.Rectangles([bbox])).opts(point_opts, rect_opts) +``` + +```python jupyter={"source_hidden": true} +# utility to extract search results into a Pandas DataFrame +def search_to_dataframe(search): + '''Constructs Pandas DataFrame from PySTAC Earthdata search results. + DataFrame columns are determined from search item properties and assets. + 'asset': string identifying an Asset type associated with a granule + 'href': data URL for file associated with the Asset in a given row.''' + granules = list(search.items()) + assert granules, "Error: empty list of search results" + props = list({prop for g in granules for prop in g.properties.keys()}) + tile_ids = map(lambda granule: granule.id.split('_')[3], granules) + rows = (([g.properties.get(k, None) for k in props] + [a, g.assets[a].href, t]) + for g, t in zip(granules,tile_ids) for a in g.assets ) + df = pd.concat(map(lambda x: pd.DataFrame(x, index=props+['asset','href', 'tile_id']).T, rows), + axis=0, ignore_index=True) + assert len(df), "Empty DataFrame" + return df +``` + +```python jupyter={"source_hidden": true} +# utility to process DataFrame of search results & return DataArray of stacked raster images +def urls_to_stack(granule_dataframe): + '''Processes DataFrame of PySTAC search results (with OPERA tile URLs) & + returns stacked Xarray DataArray (dimensions time, latitude, & longitude)''' + + stack = [] + for i, row in granule_dataframe.iterrows(): + with rasterio.open(row.href) as ds: + # extract CRS string + crs = str(ds.crs).split(':')[-1] + # extract the image spatial extent (xmin, ymin, xmax, ymax) + xmin, ymin, xmax, ymax = ds.bounds + # the x and y resolution of the image is available in image metadata + x_res = np.abs(ds.transform[0]) + y_res = np.abs(ds.transform[4]) + # read the data + img = ds.read() + # Ensure img has three dimensions (bands, y, x) + if img.ndim == 2: + img = np.expand_dims(img, axis=0) + lon = np.arange(xmin, xmax, x_res) + lat = np.arange(ymax, ymin, -y_res) + bands = np.arange(img.shape[0]) + da = xr.DataArray( + data=img, + dims=["band", "lat", "lon"], + coords=dict( + lon=(["lon"], lon), + lat=(["lat"], lat), + time=i, + band=bands + ), + attrs=dict( + description="OPERA DSWx B01", + units=None, + ), + ) + da.rio.write_crs(crs, inplace=True) + stack.append(da) + return xr.concat(stack, dim='time').squeeze() +``` + +--- + +## Identificación de los parámetros de búsqueda + +```python jupyter={"source_hidden": true} +AOI = ... +DATE_RANGE = ... +``` + +```python jupyter={"source_hidden": true} +# Optionally plot the AOI +``` + +```python jupyter={"source_hidden": true} +search_params = dict(bbox=AOI, datetime=DATE_RANGE) +print(search_params) +``` + +--- + +## Obtención de los resultados de la búsqueda + +```python jupyter={"source_hidden": true} +ENDPOINT = ... +PROVIDER = ... +COLLECTIONS = ... +# Update the dictionary opts with list of collections to search +search_params.update(collections=COLLECTIONS) +print(search_params) +``` + +```python jupyter={"source_hidden": true} +catalog = Client.open(f'{ENDPOINT}/{PROVIDER}/') +search_results = catalog.search(**search_params) +``` + +```python +df = search_to_dataframe(search_results) +df.head() +``` + +Limpiar el DataFrame `df` de forma que tenga sentido (por ejemplo, eliminando columnas/filas innecesarias, convirtiendo columnas en tipos de datos fijos, estableciendo el índice, etc.). + +```python +``` + +--- + +## Explorar y refinar los resultados de la búsqueda + + + +Consiste en filtrar filas o columnas adecuadamente para limitar los resultados de la búsqueda a los archivos de datos ráster más adecuados para el análisis y/o la visualización. Esto puede significar enfocarse en determinados mosaicos geográficos, bandas específicas del producto de datos, determinadas fechas/marcas de tiempo, etc. + + + +```python +``` + +--- + +## Procesamiento de los datos para generar resultados relevantes + +Esto puede incluir apilar matrices bidimensionales en una matriz tridimensional, mosaico de imágenes ráster de mosaicos adyacentes en un solo mosaico, etc. + +```python +``` + +--- diff --git a/book/espanol/intro.md b/book/espanol/intro.md new file mode 100644 index 0000000..493dde7 --- /dev/null +++ b/book/espanol/intro.md @@ -0,0 +1,17 @@ +# Determinación de riesgos climáticos con NASA Earthdata Cloud + +![](assets/banner.jpg) + +Repositorio de Github: https://github.com/ScienceCore/climaterisk + +Desarrollado de forma conjunta por + +![](assets/MD_logo.png) + +y + +![](assets/2i2c_logo.png) + +como un módulo del ScienceCore para + +![](assets/TOPS.png) diff --git a/espanol/CONDUCT.md b/espanol/CONDUCT.md new file mode 100644 index 0000000..37ef03e --- /dev/null +++ b/espanol/CONDUCT.md @@ -0,0 +1,8 @@ +# Código de conducta + +Como proyecto entre 2i2c y MetaDocencia, el proyecto ScienceCore:climaterisk está sujeto a las políticas siguientes: + +- [Código de conducta de 2i2c](https://compass.2i2c.org/code-of-conduct) +- [Pautas de Convivencia de MetaDocencia](https://www.metadocencia.org/pdc) + +Puede seguirse el procedimiento de notificación de cualquiera de las dos políticas. diff --git a/espanol/CONTRIBUTING.md b/espanol/CONTRIBUTING.md new file mode 100644 index 0000000..a627abf --- /dev/null +++ b/espanol/CONTRIBUTING.md @@ -0,0 +1,13 @@ +# Contribuciones + +Este repositorio ha sido creado mediante una colaboración entre 2i2c y MetaDocencia. Para hacer un cambio en este repositorio, por favor genera un _pull request_. + +Si el cambio es menor(como un arreglo de formato o una corrección gramatical), se invita a quienes proporcionan mantenimiento a auto-mergear el cambio sin revisión (omitiendo las protecciones de la rama). + +Si el cambio no es trivial, solicita una revisión a otro miembro del equipo de ScienceCore:climaterisk. + +Una vez que el cambio haya sido aprobado, el PR podrá fusionarse. Si el PR proviene de un miembro del equipo de ScienceCore:climaterisk, normalmente se fusionará por sí mismo. Si el RP proviene de alguien externo al equipo ScienceCore:climaterisk, cualquier miembro de ese equipo puede hacer la fusión o merge. + +## Código de conducta + +Ten en cuenta que este tutorial ScienceCore: Determinación de riesgos con NASA Earthdata Cloud se publica con un [Código de conducta de contribuidores](CONDUCT.md). Al contribuir con este proyecto, te comprometes a cumplir sus términos.