Sublime's custom image

Caso de estudio

Monitoreo de Humedales Altoandinos: hacia la protección integral de nuestro ecosistema

La Superintendencia Del Medio Ambiente realizó un estudio, que terminó con la clausura de pozos en uno de estos humedales en la región de Atacama. Dicho estudio es replicado aquí con las herramientas del DataCube del DO. Si quieres mayor información al respecto, puedes revisar el siguiente artículo.

Definir el Área de estudio

En ese caso, se definirá con un rectángulo la zona que abarca el humedal cuyo pozo fue clausurado y algunos humedales circundantes para poder compararlos. Se le puede dar un vistazo a la zona en el siguiente mapa (el cual es interactivo: zoom in/out y mover).

Se quieren las bandas azul, verde, roja e infrarroja, para armar los índices más conocidos. También se solicitan las bandas de calidad, para filtrar las nubes y otros artefactos.

La resolución espacial que se pide es de 30x30 metros y queremos que la imagen esté en WGS84 UTM 19 SUR (EPSG: 32719).

Los productos a utilizar son Landsat 5 y 7, que son fue la selección realizada por la SMA debido a la compatibilidad radiométrica entre estos sensores (ligeramente diferente de Landsat 8). Hay otroas productos disponibles (y se irán agregando otros), los cuáles pueden ser listados con dc.list_products().

Por cada uno de los productos definidos, se hace una búsqueda y luego se unen en un único dataset

Podemos ver el rango de fechas disponibles por sensor:

Enmascaramiento de valores no válidos

En este primer paso, se aplicará una máscara para la nubosidad, utilizando la banda de calidad de Landsat. Para Landsat 8, un resumen de los valores apropiados puede ser encontrado aquí. Si bien estas imágenes corresponden a Landsat 5 y 7, se aplican las mismas considereacinoes.

Siguiendo esa guía como paso, se fijan los valores válidos (libres de nubes) los que pertenezcan a la siguiente lista. Si la banda de calidad tiene alguno de estos valores, se considera que dicho pixel no tiene nubes.

Se puede apreciar del objeto, que se tienen 871 mediciones combinadas entre landsat 5 y 7.

Cabe destacar que el Cubo es "lazy", es decir, el código no se ejecuta inmediatamente y sólo se harán los cálculos cuando sea requerido por alguna operación o cuando el usuario lo explicite.

Calculando índices

Los índices de vegetación son muy utilizados en el monitreo del estado de la vegetación, entre otros fenómenos de interés. Uno de los más usados para ver el estado de vigor (y uno de las más antiguos), es el NDVI que utiliza las bandas rojo e infra-rojo cercano.

Se usará este índice a modo de mostración, pero se podría utilizar cualquier otro.

Nota importante: cada banda del producto está en valores enteros (int16), para indicar la reflectancia. Se guardan de esta manera para disminuir el tamaño del archivo, pero antes de realizar cálculos, deben ser escalados de vuelta a valores decimales contenidos entre [0, 1]. Este factor de escalamiento es de 0.0001 y está especificado en la documentación.

Y se agrega el NDVI como una variable extra al Cubo. De paso, es un buen momento para pedirle al Cubo que ejecute el código "lazy", principalmente porque este es un punto en el procesamiento que es la base de lo que se viene más adelante y por si desea hacer algunas visualizaciones.

Revisando sólo las zonas de interés (humedales)

Es necesario descargar los archivos vectoriales que contienen los humedales. En este caso, se usarán los del ministerio de medio ambiente, para el año 2015, los que pueden ser visualizados y descargados en el visor oficial. Para ir directo a la descarga, se puede usar este otro sitio del ministerio.

Una vez descargado el archivo, se puede cargar directamente sin descromprimir:

Como el archivo contiene los humedales para todo Chile, es necesario cortar el vector y conservar sólo la zona de estudio.

Y finalmente, se realiza en buffer a los humedales de la zona en estudio. En este caso, se hace un buffer de 100 metros, pero se puede optar por cualquier otro radio (incluso uno "negativo", que haría el buffer hacia el interior, que es utilizado a veces para minimizar el efecto borde). Este nuevo vector, será utilizado para extraer las estadísticas zonales por cada entidad.

Agrupar datos por tiempo

En primera instancia, calcular la mediana entre los meses de diciembre, enero, febrero y marzo de cada año. Primero, se arma un vector con las fechas de interés que están en el cubo, y se genera un subset

Es necesario generar una variable por temporada. Diciembre (el mes 12), pertenece al año anterior, por lo que no se puede utilizar el año directamente. Se toma la coordenada de base time y se le adicionan 31 días (aquí nos interesa el año, no el mes, por lo tanto no hay problema al hacerlo de esta forma).

En la zona están presentes algunas escenas que están "movidas" y que hay que retirar del estudio (ya fueron informadas al USGS y se espera que un futuro sean corregidas).

La nueva coordenada de tiempo (año) ha sido añadida satisfactoriamente, y se utiliza como variable de agrupamiento para calcular la mediana de la temporada (ignorando valores ausentes).

Extracción de datos por zona (zonal statitics)

Se itera por cada año y se extraen las estadísticas. Hay que procurar que tanto el archivo vectorial como el xarray tengan el mismo CRS.

Y finalmente, se visualizan los resultados. Primero se genera un gráfico dinámico y luego se genera un mapa con algunos de los NDVI de resultado (mediana por temporada) para apreciar su evolución en el tiempo.