Collection 2 images

Landsat bands

Landsat bands

Primera pregunta:

¿Que productos mantener en Nor?

Tareas

• Selección Escena de referencia

Se han buscado imágenes Landsat 8 y 9 posteriores a la fecha de la primeras imágenes de Landsat 9 (2021-10-31) que estuvieran completamente libres de nubes.

Una vez seleccionadas las escenas se procedió a ver el histograma de las bandas 1 y 2 (Coastal blue y blue).

Las siguientes escenas cumplían los requisitos:

LC08_L2SP_202034_20220802_20220806_02_T1
LC08_L2SP_202034_20220818_20220824_02_T1
LC08_L2SP_202034_20230805_20230812_02_T1
LC09_L2SP_202034_20230712_20230714_02_T1
LC09_L2SP_202034_20230728_20230802_02_T1

• Selección Escena de referencia

Finalmente la decisión ha sido la escena Landsat 8 del 2 de agosto de 2022.

Se ha preferido una Landsat 8 porque parecen radiométricamente más coherentes que Landsat 9 y porque están disponibles desde 2013, habiendo por tanto muchas más escenas en ese periodo.

     20220802l8oli_202_34

• Selección Escena de referencia

Imagen de referencia:

1. ¿Se aplica la nueva imagen a toda la serie histórica?
2. ¿Se normaliza la nueva con la antigua escena de referencia (o viceversa)?
3. ¿Se usa una escena de referencia por sensor y/o periodo?

• Selección Escena de referencia

Nomenclatura:

1. ¿Ori por SR?
2. ¿Geo y Rad existen pero no se guardan?
3. ¿Nor pasa a Nor2?
4. ¿Se mantiene el antiguo formato?

• Selección Escena de referencia

Formato:

1. Valores: Float/Integer (ahora Float 0-1)
2. COG/GTIFF

Productos

Productos

El código está hecho y funcionando, a falta de normalizar las escenas para probarlo con las imágenes normalizadas

El mayor esfuerzo de tiempo se ha empleado en la mejora de la máscara de máxima inundación. Se trata de un enfoque dinámico basado tanto en donde puede haber agua como en donde no puede haber agua:

• Zonas por debajo de 5 m s.n.m.
• Zonas con Σ Fmask >= 5
• Zonas con NDVI alto todo el año
• Zonas con altura de la vegetación elevada (no podemos ver debajo de los árboles)
• Zonas con mucha pendiente
• Zonas con sombra orográfica
• Zonas de nubes y sombra de nubes (*)

Productos

Productos

Productos

P10 >= 0.4

Mean >= 0.4

Productos Cob Veg

Productos Slope

Productos (zona No Inundable)

Productos (DTM <= 5 m s.n.m.)

Productos (Fmask)

Productos (Fmask)

Productos (zona inundable)

Productos (Water Mask)

• ¿ Que sumamos a Fmask la máscara de Fmask o la máscara de agua resultante?
• ¿Hacemos otro sumatorio con nuestras máscaras de agua?
• ¿Merece la pena calcular esos parámetros por ciclos?

Otras cuestiones

• ¿Actualización de máquina virtual?
• Landsat 7 gapfill
• Nuevos productos (Fenología, GPP, etc...)
• Metadatos

No Data:               -9999

Tierra limpia:                0

Agua:                              1

Nubes:                           2

Sombra Nubes:            3

Sombra Orográfica:    4

Valores Propuestos

Protocolo V2. Estructura

__init__.py

download.py

protocolov2.py

productos.py

hidroperiodo.py

utils.py

run_download.sh

run_hydroperiod.sh

https://github.com/Digdgeo/ProtocoloV2

Protocolo V2. Estructura

protocolov2.py

Métodos

__init__(ruta_escena)

get_hillshade()

get_cloud_pn()

remove_masks()

projwin()

coef_sr_st()

normalize()

nor1(banda, mascara, std=1)

nor2l8(banda, slope, intercept)

run()

El proceso de normalización comienza pasando la ruta a la escena original de Landsat (Colection 2 Level 2).

El método run() se encarga de ejecutar todo el proceso. En primer lugar se realiza la descarga del quicklook de la escena y se añade la entrada a la base de datos (MongoDB). 

A continuación se genera el Hillshade con los datos solares obtenidos del MTL, posteriormente se reclasifica y recorta Fmask con el shape del Parque Nacional. Pasamos a recortar la escena con el shape del WRS 202-34 y a aplicar los coeficientes para tener las reflectividades en 0-1.

Finalmente pasamos a la normalización, controlado el flujo por normalize, nor1 para obetner los parámetros y nor2 para aplicarlos y guardar las bandas normalizadas

Protocolo V2. Estructura

productos.py

Métodos

__init__(ruta_nor)

ndvi()

ndwi()

mndwi()

flood()

turbidity()

depth()

get_flood_surface()

run()

La obtención de productos comienza instanciando la clase pasando la ruta a la escena normalizada. 

Los productos a fecha de hoy son NDVI, NDWI, MNDWI, Flood, Turbidity y Depth, otros*. Existe otro método get_flood_surface() que es el que obtiene la superficie inundada de los recintos principales de la marisma**.

El método run() va ejecutándolos todos y guardando la información en la base de datos

**¿Polígonos de censo aéreo?

*¿Otros productos?

***LST***

Protocolo V2. Estructura

Base de datos (MongoDB)

Protocolo V2. Estructura

hidroperiodo.py

funciones

get_escenas_values()

get_hydroperiod()

get_products()

Gini() & IRT()**

El hidroperiodo no es una clase sino un script con 3 funciones.

La primera get_escenas_values() se encarga de calcular los días "que vale" cada escena del ciclo hidrológico (1 de septiembre - 31 de agosto). 

La segunda get_hydroperiod() calcula el número de días que un pixel ha estado inundado o seco o ha sido NoData.

La tercera get_products() genera los rasters con los hidroperiodos y el número de días válidos del ciclo para cada pixel. 

Protocolo V2. Estructura

Protocolo V2. Estructura

Base de datos (MongoDB)

Protocolo V2. Estructura

utils.py

funciones

enviar_correo(destinatarios, asunto, cuerpo, archivo_adjunto)

proceso_finalizado(info_escena, archivo_adjunto)

prepare_hydroperiod(productos_dir, output_dir, ciclo_hidrologico, umbral_nubes)

En utils.py se engloban las dos funciones para mandar el correo con la información de la escena y una función para preparar los datos para ejecutar el hidroperiodo. 

Los mails se mandan desde la cuenta de gmail del LAST lastebd.protocolo.vlab@gmail.com

al terminar de procesar una escena. Informa sobre la cobertura de nubes, la superficie inundada  y adjunta el quicklook de la misma. 

La función del hidroperiodo se encarga de copiar las máscaras de agua agrupadas por ciclo

Protocolo V2. Estructura

download.py

funciones

download_landsat_scenes(username,
password,
latitude,
longitude,
days_back=15,
end_date=None, 
process=True,
max_cloud_cover=100, 
output_dir='path/to/save/rars')

Es donde sucede "la magia". Se trata del script que se encarga de buscar escenas nuevas para el rango de fechas seleccionado, cotejarlas con la base de datos y en caso de no estar, realizar la descarga y llamar a las clases Landsat() y Product() para procesar la escena. Si todo va bien, termina enviando un mail llamando a las funciones de utils.py

Protocolo V2. Estructura

run_download.sh --> download_landsat_scenes()

run_hydroperiod.sh --> prepare_hydroperiod(), get_escenas_values(), get_hydroperiod() y get_products()

Son dos archivos bash que son los que se añaden al crontab para que se ejecuten una vez a la semana (run_download.sh) y una vez al año (run_hydroperiod, que se ejecutará cada 15 de octubre)

Protocolo V2. Estructura

Archivos en /nor

Archivos en /pro

• bandas (blue-swir2).tif
• gráfico regresión bandas.png
• coeficientes.txt
• fmask.tif
• hillshade.tif

• NDVI.tif
• NDWI.tif
• MNDWI.tif
• Flood.tif
• Depth.tif
• Turbidity.tif
• LST.tif
• superficie_inundada.csv

Archivos en /hyd

• máscaras de agua.tif
• hydroperiod.itf
• valid_days.tif
• IRT.tif

/hidroperiodoYYYY-YYYY/

Protocolo V2. Estructura

Protocolo V2. Flood

Protocolo V2. Flood

La máscara basada en Fmask, aunque mejora bastante a la anterior, sigue sin ser perfecta. 

Buscar una metodología que no depende de ningún tipo de máscara de zonas inundables, sino de definir en que zonas no podríamos detectar agua aunque la hubiera. 

Criterios:

Swir1 <= 0.12

slope > 8 (**embalses NDWI_scene & MNDWI_scene**)

hillshade < p30

ndvip10 > 0.3 & ndvimean > 0.5

cobveg > 75

ndvi_scene > 0.6 & dtm > 2.5

fmask_scene clouds & shadows

sum(reclass((fmask_scene, ndwi_scene, mndwi_scene)) >= 2

Protocolo V2. Flood

Protocolo V2. Flood

Protocolo V2. Flood

Protocolo V2. Flood

Protocolo V2. Flood

• En la marisma y Santa Olalla las máscaras de agua son prácticamente iguales, aunque hay veces que la máscara de nubes de la Colección 2 Level 2 juega malas pasadas (algunas zonas de aguas muy someras las marca como nube, nieve, etc...). Y por lo general marca más nubes.
• La v2 puede tener más errores de falsos positivos (no tanto más, como en zonas distintas (plásticos negros y algunas parcelas de cutlivos en lugar de en los pinares y en las sombras orográficas de la v1)
• Las sombras de nubes no enmascaradas correctamente siguen dando errores (al igual que en la v1)
• Se ha mejorado mucho en las sombras orográficas, en zonas de vegetación y en la costa. Al igual que en algunos embalses y muchas balsas de riego que no estaban en la antigua máscara de máxima inundación

Conclusiones comparativa v1-v2

Protocolo V2. Issues

• Las escenas de la colección 2 nivel 2 no son recomendables para trabajar con agua. Desde la propia USGS lo indican y han sacado un nuevo producto con otra corrección más enfocada al agua que sirven bajo demanda. 
• Hay valores muy extraños que hacen que tanto el NDVI como la turbidez den valores muy raros en esos lugares. 
• En tierra y en la marisma funcionan bien, aunque sí hay a veces valores anómalos en los embalses 
• No parece afectar a la mascara de agua ni por tanto al hidroperiodo
• Se podría probar la colección 2 nivel 1 para ver si se evitan esos errores

Issues Colection 2 Level 2

Protocolo V2. Issues

Actualización 25/10/2024

• Procesada por completo la versión 2 (colección 2 Level 2)
• 1093 escenas 1984-presente (TM, ETM+ y OLI)
• 748 productos  máscaras de agua (744 máscaras de agua)

• Se ejecuta en la máquina CentOS y guarda directamente en el "nuevo Saturno"

• 413 TM (Landsat 5)
• 370 ETM+ (Landsat 7)
• 310 OLI (Landsat 8, 9)

• Desviaciones del Hidroperiodo con respecto a la media (1984/85-2023/24)

• Se ha añadido un script al final del proceso para exportar la base de datos a json para que pueda ser abierta fuera de un entorno de MongoDB ("/datos_last/mongo_data")

• FMASK falla, especialmente en las TM y en situaciones de inundación pecularies, como exceso de sedimentos/turbidez o aguas muy someras. 
• Hay otra máscara que se encuentra en las escenas ETM+ y TM (no en las OLI) que parece que funciona mejor
• Gapfill (se están procesando las escenas nuevamente rellenando los gapfills con GDAL (funciona muy muy aceptablemente bien)
• La idea sería procesar de nuevo las TM con la máscara alternativa a FMASK y volver a correr productos e hidroperiodos

• Conclusiones Versión 2

Los productos de inundación no son perfectos, especialmente en zonas de sombras que no quedan enmascaradas por las sombras en el Hillshade (ej. zonas urbanas) y en algunas zonas de invernaderos (sobre todo plásticos negros)

La idea es ir mejorando el producto, aún así, no funciona peor que el anterior y, a riesgo de ser crucificado, creo que tanto la inundación como el hidroperiodo han mejorado con respecto a la versión 1.

Se han normalizado más escenas que en la versión 1. Al ser la descarga automática, se procesan (y hasta normalizan!) escenas que de haberlas comprobado manualmente probablemente se hubieran descargado

¡Necesidad de validación (productos y normalización)!

¡Necesidad de nuevos datos de campo!

• Próximos pasos

1. Guardar los datos del censo aéreo en postgreSQL (Eduardo)
2. Método Inundación con DTM
3. Geoserver (COG?)
4. Sentinel 2 (¿Misma metodología?)
5. Actualización del Servidor (cloudlast02? Rocky Linux?)

Actualización Protocolo v2 14/11/2024

Se va a procesar nuevamente la colección 2 (~1 semana)

para incluirles Máscara de nubes alternativa a TM y ETM+ y relleno de gapfill a ETM+

- Se ha incluido el envío de mail en caso de que no se normalice la escena (añadiendo porcentaje de superficies por recintos)

- Se ha añadido el método apply_gapfill con GDAL

- Se ha añadido el procesado con la máscara de nubes alternativa a Fmask en las escenas TM y ETM+

- Pruebas (exitosas) con Geoserver (Ya se podrían mandar los rasters al terminar cada escena)

Actualización Protocolo v2 14/11/2024

A decidir:

1) De verdad gapfill?

2) Estructura y nomenclatura en /pro

3) Umbral de nubes para el hidroperiodo (automático?)

4) Data type. Propuesta seguir en float y eliminar /ori (748 gb) y los rasters auxiliares (~250 gb) de hidroperiodo

5) COG!?

Actualización Protocolo v2 14/11/2024

A realizar:

* Cruce con lagunas Carola y Censo aéreo

* Exportar a PostgreSQL

* Flood DTM

* Línea de costa

* Read the Docs

* Videotutorial