Diego García Díaz
GIS & Remote Sensing. Python & Statistics. Surf & MTB :)
Protocolo automático para el tratamiento de imágenes Landsat y la generación de productos derivados
I meeting proyecto GOYAS (2/12/2024)
Protocolo automático para el tratamiento de imágenes Landsat y la generación de productos derivados
Javier Bustamante Díaz
Responsable científico
(Investigador Científico)
Ricardo Díaz-Delgado Hernández
Dr. Biología
Isabel Afán Asencio
Responsable Técnico
Dra. Biología
David Aragonés Borrego
Responsable Técnico
Ingeniero Forestal
Diego García Díaz
Geógrafo
Pedro Jesús Gómez Giráldez
Ingeniero de Montes
Gabriela Patricia Romero Olivos
Geógrafa
1- Productos
(y metadatos)
2 - Scripts
3 - Geoportal
Flood Mask
Turbidez
Profundidad
Hidroperiodo.
¿Qué podemos aportar al proyecto?
2 - Scripts de procesado
2 - Scripts
3 - Geoportal
class Landsat:
"""Clase para trabajar con las Landsat de la nueva Collection 2 Level 2 del USGS"""
def __init__(self, ruta_escena):
"""Inicializa un objeto Landsat basado en el path de la escena.
Args:
ruta_escena (str): Ruta al directorio de la escena Landsat.
"""
# Definimos los paths de entrada
self.ruta_escena = ruta_escena
self.escena = os.path.split(self.ruta_escena)[1]
self.ori = os.path.split(self.ruta_escena)[0]
self.base = os.path.split(self.ori)[0]
self.data = os.path.join(self.base, 'data')
self.geo = os.path.join(self.base, 'geo')
self.rad = os.path.join(self.base, 'rad')
self.nor = os.path.join(self.base, 'nor')
self.pro = os.path.join(self.base, 'pro')
# Definimos un grupo de variables de la escena
data_escena = self.escena.split('_')
self.sat = "L" + data_escena[0][-1]
sensores = {'L4': 'TM', 'L5': 'TM', 'L7': 'ETM+', 'L8': 'OLI', 'L9': 'OLI'}
self.sensor = sensores[self.sat]
self.nprocesado = data_escena[1]
self.path = data_escena[2][:3]
self.row = data_escena[2][-3:]
self.escena_date = data_escena[3]
self.escena_procesado_date = data_escena[4]
self.collection = data_escena[5]
self.tier = data_escena[6]
# Mascara de nuebes. Hay que hacerlo así porque sabe dios por qué ^!·/&"! los valores no son los mismos en OLI que en ETM+ y TM
if self.sensor == 'OLI':
self.cloud_mask_values = [21824, 21952]
else:
self.cloud_mask_values = [5440, 5504]
#Definimos nombre last
if self.sensor == 'ETM+':
self.last_name = (self.escena_date + self.sat + self.sensor[:-1] + self.path + '_' + self.row[1:]).lower()
else:
self.last_name = (self.escena_date + self.sat + self.sensor + self.path + '_' + self.row[1:]).lower()
# Definimos paths de salida
self.pro_escena = os.path.join(self.pro, self.last_name)
os.makedirs(self.pro_escena, exist_ok=True)
self.geo_escena = os.path.join(self.geo, self.last_name)
os.makedirs(self.geo_escena, exist_ok=True)
self.rad_escena = os.path.join(self.rad, self.last_name)
os.makedirs(self.rad_escena, exist_ok=True)
self.nor_escena = os.path.join(self.nor, self.last_name)
os.makedirs(self.nor_escena, exist_ok=True)
self.pro_escena = os.path.join(self.pro, self.last_name)
os.makedirs(self.pro_escena, exist_ok=True)
# Definimos máscaras a utilizar
self.equilibrado = os.path.join(self.data, 'Equilibrada.tif')
self.noequilibrado = os.path.join(self.data, 'NoEquilibrada.tif')
self.parametrosnor = {}
self.iter = 1
#CReamos un diccionario a partir del MTL
self.mtl = {}
for i in os.listdir(self.ruta_escena):
if i.endswith('MTL.txt'):
mtl = os.path.join(self.ruta_escena, i)
f = open(mtl, 'r')
for line in f.readlines():
if "=" in line:
l = line.split("=")
self.mtl[l[0].strip()] = l[1].strip()
#Creamos las variables con la franja del espectro de cada banda
if self.sat in ['L8', 'L9']:
for i in os.listdir(self.ruta_escena):
if i.endswith('.TIF'):
banda = os.path.splitext(i)[0].split('_')[-1]
if banda == 'B2':
self.blue = os.path.join(self.ruta_escena, i)
elif banda == 'B3':
self.green = os.path.join(self.ruta_escena, i)
elif banda == 'B4':
self.red = os.path.join(self.ruta_escena, i)
elif banda == 'B5':
self.nir = os.path.join(self.ruta_escena, i)
elif banda == 'B6':
self.swir1 = os.path.join(self.ruta_escena, i)
elif banda == 'B7':
self.swir2 = os.path.join(self.ruta_escena, i)
elif i.endswith('QA_PIXEL.TIF'):
self.qa = os.path.join(self.ruta_escena, i)
print(self.qa)
else:
continue
elif self.sat in ['L7', 'L5']:
for i in os.listdir(self.ruta_escena):
if i.endswith('.TIF'):
banda = os.path.splitext(i)[0].split('_')[-1]
if banda == 'B1':
self.blue = os.path.join(self.ruta_escena, i)
elif banda == 'B2':
self.green = os.path.join(self.ruta_escena, i)
elif banda == 'B3':
self.red = os.path.join(self.ruta_escena, i)
elif banda == 'B4':
self.nir = os.path.join(self.ruta_escena, i)
elif banda == 'B5':
self.swir1 = os.path.join(self.ruta_escena, i)
elif banda == 'B7':
self.swir2 = os.path.join(self.ruta_escena, i)
elif i.endswith('QA_PIXEL.TIF'):
self.qa = os.path.join(self.ruta_escena, i)
print(self.qa)
else:
continue
else:
print('No encuentro ninguna escena Landsat')
# Descargamos el quicklook de la escena
url_base = 'https://landsatlook.usgs.gov/gen-browse?size=rrb&type=refl&product_id={}'.format(self.mtl['LANDSAT_PRODUCT_ID'].strip('""'))
self.qk_name = os.path.join(self.ruta_escena, self.escena + '_Quicklook.jpeg')
if not os.path.exists(self.qk_name):
qk = open(self.qk_name, 'wb')
qk_open = urlopen(url_base)
urlimg = qk_open.read()
qk.write(urlimg)
qk.close()
print('QuicKlook descargado')
else:
print('El Quicklook ya estaba previamente descargado')
#print(url_base)
print('Landsat iniciada con éxito')
self.pn_cover = None
#Creamos el json para instarlo en la base de datos MongoDB
self.newesc = {'_id': self.last_name,
'usgs_id': self.mtl['LANDSAT_SCENE_ID'][1:-1],
'tier_id': self.mtl['LANDSAT_PRODUCT_ID'][1:-1],
'lpgs': self.mtl['PROCESSING_SOFTWARE_VERSION'][1:-1],
'category': self.mtl['COLLECTION_CATEGORY'][1:-1],
'Clouds': {'cloud_scene': float(self.mtl['CLOUD_COVER']),
'land cloud cover': float(self.mtl['CLOUD_COVER_LAND'])},
'Info': {'Tecnico': 'LAST-EBD Auto',
'Iniciada': datetime.now(), 'Pasos': {'rad': '', 'nor': ''}}}
try:
db.insert_one(self.newesc)
except Exception as e:
db.update_one({'_id': self.last_name}, {'$set':{'Info.Iniciada': datetime.now()}}, upsert=True)
#print("Unexpected error:", type(e)) #se Podria dar un error por clave unica, por eso en
#ese caso, lo que hacemos es actualizar la fecha en la que tratamos la imagen
print('Landsat instanciada y subida a la base de datos')¿Qué podemos aportar al proyecto?
3 - Geoportal:
presente y
futuro
3 - Geoportal
http://venus.ebd.csic.es/imgs/
¿Qué podemos aportar al proyecto?
Protocolo Automático v2
- Versión 2 del Protocolo automático para el tratamiento de imágenes Landsat.
- Desarrollado en el LAST-EBD como respuesta a un problema complejo, como es el estudio de la inundación en la marisma de Doñana.
- Las condiciones cambiantes de la marisma hacen muy difícil definir una metodología para estimar la cubierta de agua (en una serie temporal) por medio de índices de agua u otro tipo de algoritmos.
- Por ello, se decidió llevar a cabo una normalización (basada en Áreas Pseudo Invariantes (PIAs)) de las imágenes en base a la reflectividad de unas áreas en las que a priori debería ser homogenéa y constante a lo largo del tiempo (arenas, mar y embalses, cascos históricos, etc...)
- Basado en las campañas de campo realizadas por personal del LAST y de la Estación Biológica de Doñana.
- Reajustado dentro del Proyecto Ecopotential con unos 1440 puntos de campo
Máquina Virtual
- CentOS 7
- 64 bits
- 64 Gigas RAM
- 32 núcleos
- Conectado a un NAS (32 TB)
Características
Software
- MongoDB
- GDAL
- Git
- Anaconda (python 3.8)
- Numpy
- Rasterio
- Fiona
- Geopandas
- Scipy
download.py
protocolov2.py
productos.py
hidroperiodo.py
utils.py
run_download.sh
run_hydroperiod.sh
GitHub Repository
Código
Carpetas
/ori
/geo
/rad
/nor
/pro
/hyd
/data
/mongo
- Descarga Automática (SR-Colección 2 L2)
- Geo (recorte 202/34)
- Rad (aplicar coeficientes escalado)
- Normalización
- Productos
- Actualización de MongoDB y PostgreSQL
Busqueda y descarga
Hillshade
MDT
Geo & Rad
Productos
- NDVI
- NDWI
- MNDWI
- FLOOD
- DEPTH
- Turbidity
- Hydroperiod
Normalización
Escena
Landsat
Pias
PIAs Ref
FMASK
Data
No
GeoServer
ICTS
LAST
{"_id": "20240815l9oli202_34",
"usgs_id": "LC92020342024228LGN01",
"tier_id": "LC09_L1TP_202034_20240815_20240815_02_T1",
"lpgs": "LPGS_16.4.0",
"category": "T1",
"Clouds": {"cloud_scene": 0.06,
"land cloud cover": 0.06,
"cloud_PN": 0 },
"Info": {
"Tecnico": "LAST-EBD Auto",
"Iniciada": {"$date": {"$numberLong": "1724690285816"}},
"Pasos": {"rad": "",
"nor": {"Normalize": "True",
"Nor-Values": {
"blue":
{"Parametros":
{"slope": 1.020172035444719,
"intercept": -0.00711237939135306,
"std": 0.016095327213406563,
"r": 0.98906989518987,
"N": 54886,
"iter": 1},
"Tipo_Area":
{"Mar": 49347,
"Embalses": 793,
"Pinar": 1564,
"Urbano-1": 575,
"Urbano-2": 944,
"Aeropuertos": 547,
"Arena": 664,
"Pastizales": 266,
"Mineria": 186}},
"Productos":
["NDVI", "NDWI", "MNDWI",
{"Flood": {
"El Rincon del Pescador": 72.72,
"Marismillas": 0.99,
"Caracoles": 0,
"FAO": 14.94,
"Marisma Occidental": 1.26,
"Marisma Oriental": 6.57,
"Entremuros": 115.29} },
"Turbidity",
"Depth" ]}{"_id": "hidroperiodo_2023-2024",
"escenas": [
{"escena_id": "20240409l9oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 19896.300000000003},
{"escena_id": "20240503l8oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 11161.53},
{"escena_id": "20240511l9oli202_34",
"nubes_marismas": 8.31,
"ha_inundacion": 8354.789999999999},
{"escena_id": "20240527l9oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 3754.53},
{"escena_id": "20240604l8oli202_34",
"nubes_marismas": 0.06,
"ha_inundacion": 2661.84},
{"escena_id": "20240714l9oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 249.93},
{"escena_id": "20240722l8oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 215.1},
{"escena_id": "20240807l8oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 209.43},
{"escena_id": "20240815l9oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 211.76999999999998},
{"escena_id": "20231008l8oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 305.91},
{"escena_id": "20231024l8oli202_34",
"nubes_marismas": 7.69,
"ha_inundacion": 1131.8400000000001},
{"escena_id": "20231117l9oli202_34",
"nubes_marismas": 0,
"ha_inundacion": 408.51},
{"escena_id": "20231125l8oli202_34",
"nubes_marismas": 0.03,
"ha_inundacion": 421.83},
{"escena_id": "20231219l9oli202_34",
"nubes_marismas": 4.76,
"ha_inundacion": 589.5899999999999},
{"escena_id": "20231227l8oli202_34",
"nubes_marismas": 0.63,
"ha_inundacion": 512.8199999999999},
{"escena_id": "20240112l8oli202_34",
"nubes_marismas": 0.09,
"ha_inundacion": 696.69},
{"escena_id": "20240120l9oli202_34",
"nubes_marismas": 0.01,
"ha_inundacion": 1138.59},
{"escena_id": "20240221l9oli202_34",
"nubes_marismas": 7.73,
"ha_inundacion": 4288.32},
{"escena_id": "20240229l8oli202_34",
"nubes_marismas": 0.01,
"ha_inundacion": 3901.3200000000006},
{"escena_id": "20240316l8oli202_34",
"nubes_marismas": 6.05,
"ha_inundacion": 4150.98}]}Landsat
Hidroperiodo
Productos
-
dtm
-
fmask_escena
-
hillshade_escena
-
ndvi_escena
-
ndwi_escena
-
mndwi_escena
-
cob_veg
-
ndvi_p10
-
ndvi_mean
-
fmask
-
mndwi
-
ndwi
-
swir1
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
Flood
Depth
Hydroperiod
Valid days
Hydroperiod
nº días inundado / nº días válidos * 365
Anomalías
ciclo / media(84-24)
Conclusiones
- Metadatos: MongoDB pero hay que implementar un estandar (¿ayuda?)
- Productos: Máscara de inundación y productos derivados: Profundidad, Turbidez, Hidroperiodo (nº de días válidos representados por cada escena)
- El propio Protocolo puede ser un producto en sí mismo, al igual que el hidroperiodo
Gracias por su atención
deck
By Diego García Díaz
