Protocolo Automático para la normalización de imágenes Landsat

Protocolo Automático

  • Importación
  • Reproyección
  • Corrección Radiométrica
  • Normalización

Python

MongoDB

Prerrequisitos

Miramón

Importación?

Instanciación

Al trabajar con clases debemos "instanciar" o crear de la clase. En este paso hacemos muchas cosas que son necesarias para la correcta ejecución del Protocolo:

def __init__(self, ruta):
        
        
        '''Instanciamos la clase con la escena que queramos, hay que introducir la ruta a la carpeta en ori
        y de esa ruta el constructor obtiene el resto de rutas que necesita para ejecutarse'''
                
        self.ruta_escena = ruta
        self.ori = os.path.split(ruta)[0]
        self.escena = os.path.split(ruta)[1]
        self.raiz = os.path.split(self.ori)[0]
        self.geo = os.path.join(self.raiz, 'geo')
        self.rad = os.path.join(self.raiz, 'rad')
        self.nor = os.path.join(self.raiz, 'nor')
        self.data = os.path.join(self.raiz, 'data')
        self.mimport = os.path.join(self.ruta_escena, 'miramon_import')
        if not os.path.exists(self.mimport):
            os.makedirs(self.mimport)
        self.bat = os.path.join(self.ruta_escena, 'import.bat')
        self.bat2 = os.path.join(self.rad, 'importRad.bat')
        self.equilibrado = os.path.join(self.data, 'equilibrada.img')
        self.noequilibrado = os.path.join(self.data, 'MASK_1.img')
        self.parametrosnor = {}
        self.iter = 1
        self.cloud_mask = None
        for i in os.listdir(self.ruta_escena):
            if i.endswith('MTL.txt'):
                mtl = os.path.join(self.ruta_escena,i)
                arc = open(mtl,'r')
                for i in arc:
                    if 'LANDSAT_SCENE_ID' in i:
                        usgs_id = i[-23:-2]
                    elif 'CLOUD_COVER' in i:
                        cloud_scene = float(i[-6:-1])
        arc.close()
        #print "El porcentaje de nubes en la escena es de " + str(cloud_scene)
        
        self.newesc = {'_id': self.escena, 'usgs_id': usgs_id, 'Clouds': {'cloud_scene': cloud_scene},\
                       'Info': {'Tecnico': 'LAST-EBD Auto', 'Iniciada': time.ctime(),'Pasos': {'geo': '', 'rad': '', 'nor': ''}}}
        
        # Conectamos con MongoDB
        connection = pymongo.MongoClient("mongodb://localhost")

        # DB teledeteccion, collection landsat
        
        db=connection.teledeteccion
        landsat = db.landsat
        
        try:
        
            landsat.insert_one(self.newesc)
        
        except Exception as e:
            
            landsat.update_one({'_id':self.escena}, {'$set':{'Info.Iniciada': time.ctime()}})
  • Definimos las rutas
  • Creamos variables que se utilizarán en posteriores métodos
  • Insertamos el registro en la base de datos

 

Protocolo Automático

Importación

Métodos

  • fmask()
  • fmask_legend()
  • mascara_cloud_pn()
  • get_cloud_pn()
  • createG_bat()
  • callG_bat()
  • get_kl_csw()
  • remove_masks()
  • Máscara de Nubes
  • Añadirle leyenda
  • Obtener cobertura de nubes PN
  • Importar la escena con Miramón
  • Obtener Kls

Máscara de  Nubes

Idealmente usamos Fmask. 

  • Se trata de un algoritmo diseñado por Zhe Zhu (EROS, United States Geological Services) y Curtis E. Woodcock, Boston University).
  • Aplicable a las Landsat 4-8 y Sentinel 2 .
  • Muy contrastado y utilizado, que en teoría sustituirá a la Banda de calidad de Landsat (algún día...).

 

 

 

 

Pero... A veces hay fallos

Alrededor de un 5% de las imágenes dan error al ejecutar Fmask. En esas ocasiones se utiliza la banda de calidad que viene con las Landsat 8 (la que será sustituida por Fmask y se aplicará a toda la serie histórica de Landsat).

En esas ocasiones usaremos la librería de Python pymasker. Diseñada para extraer información de la banda de calidad de Landsat 8 (también funciona con MODIS!)

a = os.system('C:/Cloud_Mask/Fmask 1 1 0 50')
a
if a == 0:
    self.cloud_mask = 'Fmask'
    print 'Mascara de nubes (Fmask) generada en ' + str(t-time.time()) + ' segundos'
else:
    print 'comenzando BQA'
    for i in os.listdir(self.ruta_escena):
        if i.endswith('BQA.TIF'):
            masker = landsatmasker(os.path.join(self.ruta_escena, i))
            mask = masker.getcloudmask(confidence.high, cirrus = True, cumulative = True)
            masker.savetif(mask, os.path.join(self.ruta_escena, self.escena + '_Fmask.TIF'))
    self.cloud_mask = 'BQA'
    print 'Mascara de nubes (BQA) generada en ' + str(t-time.time()) + ' segundos'
                           
except Exception as e:

print "Unexpected error:", type(e), e
  • La máscara de nubes se recorta con un shape del Parque Nacional de Doñana y se calcula el porcentaje de nubes y/o sombra de nubes sobre el mismo. 
  • Este dato se pasará a la base de datos. En la base de datos se almacenan:
    1. Nubes en Escena
    2. Nubes en PN 
    3. Umbral usado (Fmask). 
  • Se añade al hdr de la máscara de nubes la información de clases {1: Agua, 2: Sombra de Nubes...}. En este paso solo al hdr (Fmask), en /nor tanto al hdr como al doc.

 

 

 

 

Importación a Miramón

  • Miramón se emplea en 2 pasos (Importación y Corrección Radiométrica).
  • Se llama mediante bat a los módulos del programa tifimg y corrad.
  • Tanto el rel como los docs se van modificando, tratados como archivos de texto plano, adecuándose a los valores que deben llevar en cada paso.
def createI_bat(self):
        
        '''-----\n
        Este metodo crea un archivo bat con los parametros necesarios para realizar la importacion'''

        #estas son las variables que necesarias para crear el bat de Miramon
        tifimg = 'C:\\MiraMon\\TIFIMG'
        num1 = '9'
        num2 = '1'
        num3 = '0'
        salidapath = self.mimport #aquí va la ruta de salida de la escena
        dt = '/DT=c:\\MiraMon'

        for i in os.listdir(self.ruta_escena):
            if i.endswith('B1.TIF'):
                banda1 = os.path.join(self.ruta_escena, i)
            elif i.endswith('MTL.txt'):
                mtl = "/MD="+self.ruta_escena+"\\"+i
            else: continue

        lista = [tifimg, num1, banda1,  salidapath, num2, num3, mtl, dt]
        print lista

        batline = (" ").join(lista)

        pr = open(self.bat, 'w')
        pr.write(batline)
        pr.close()


def callI_bat(self):
    
    '''-----\n
    Este metodo llama ejecuta el bat de la importacion. Tarda entre 7 y 21 segundos en importar la escena'''

    #import os, time
    ti = time.time()
    a = os.system(self.bat)
    a
    if a == 0:
        print "Escena importada con éxito en " + str(time.time()-ti) + " segundos"
    else:
        print "No se pudo importar la escena"
    #borramos el archivo bat creado para la importación de la escena, una vez se ha importado ésta
    os.remove(self.bat)

Obtención de Kls

Obtención automática del objeto oscuro de cada banda (P.S. Chaves, 1988)

  • Obtener un DTM de cada escena
  • Generar el hillshade a partir del DTM de la escena, teniendo en cuenta la altura y el azimuth solar obtenidos del archivo de metadatos original _MTL.txt
  • Obtener la máscara de cuerpos de agua (Fmask o BQA)
  • En caso de Landsat 7 se aplicará también un Gapmask

 

 

 

 

Finzalizada la Importación:

  • Máscara de nubes (Fmask o BQA)
  • GeoTIF originales y los doc, img y rel de Miramón
  • Kl.rad modificado con los valores de la escena
  • Histogramas de cada banda generados y guardados en /rad

Protocolo Automático

Reproyección

Reproyectamos la escena a ED 1950 UTM huso 30 mediante la misma reproyección de 3 parámetros usada en el Protocolo Manual

Métodos

  • reproject()
  • copyDocG()
  • modifyDocG()
  • modifyRelG()
  • modify_hdr_rad_pro()
  • Reproyectamos las bandas desde las GeoTIF originales a img + hdr con Gdalwarp
  • Copiamos y modificamos los doc y el rel desde MiramonImport
['gdalwarp', '-s_srs', '"+proj=utm +zone=29 +datum=wgs84 +units=m"', '-t_srs', 
'"+proj=utm +zone=30 +ellps=intl +towgs84=-84,-107,-120,0,0,0,0 +units=m +no_defs"', 
'-r', 'cubic', '-te', '78000 4036980 340020 4269000', '-tr', '30 30', '-of', 'ENVI']

Finalizada la Reproyección:

  • Todas las bandas convertidas a img + hdr
  • docs y rel para cada banda
  • Máscara de nubes reproyectada y guardada en la carpeta de la escena en nor

Protocolo Automático

Corrección Radiométrica

Métodos

  • copy_files_GR()
  • createR_bat()
  • callR_bat()
  • cleanR()
  • modify_hdr_rad()
  • correct_sup_inf()
  • clean_correct_R()
  • modify_hdr_rad_255()
  • translate_bytes_gdal()
  • modifyRelRad()
  • re_clean_R()
  • modifyDocR()
  • modify_hdr_rad_pro()
  • Copiamos desde /geo
  • Creamos y ejecutamos el bat
  • Borramos los archivos copiados
  • Corregimos ref 0/100
  • Convertimos de float32 a byte
  • Volvemos a eliminar los archivos sobrantes
  • modificamos docs y rel
  • Necesitamos tener en la carpeta raíz /rad:
    1. sindato.img (DTM plano)
    2. kl.rad (que será el que llamé Miramón, luego, se copiará con el nombre de la escena + '_kl.rad' en el directorio de la escena en /rad
  • El método no cambia entre Landsat 7-8, ya que en ambos casos la salida es en float32 (reales en porcentaje)
  • Todas las reconversiones y la lectura y creación de rasters a partir de arrays se hace con GDAL y/o Rasterio

Finalizada la Corrección Radiométrica:

  • Bandas (1-9) en reflectancia en superficie sin normalizar 
  • Formato byte convertido desde el float32 del Corrad
  • Archivo de kl copiado en la carpeta de la escena en /rad y añadido también a la base de datos

Protocolo Automático

Normalización

Métodos

  • normalize()
  • nor1()
  • nor2()
  • copyDocR()
  • modifyRelNor()
  • fmask_legend()
  • fmask_binary()
  • modify_hdr_rad_pro()
  • Normalizar la escena con la escena de referencia 
  • Criterios R > 0.85 y nº pixeles por PIA >= 10
  • OLS ref/actual (SciPy)
  • Añadir la leyenda a la máscara de nubes (hdr y doc)

Normalización. 

  • Se va leyendo las bandas como arrays y se van aplicando máscaras con las Áreas Pseudo Invariante que corresponda.
  • Una vez aplicadas todas las máscaras se realiza una primera regresión entre los valores válidos de que han quedado en las PIAs y los valores de la escena de referencia
  • Se eliminan los píxeles cuyo residuo es superior a 11 y se vuelve a realizar una regresión con los restantes. Si cumple los objetivos de R > 0.85 y nº píxeles por PIA >= 10 se produce la normalización (Rasterio)
  • Sí... También está contemplado nor1bis para ir contemplando criterios más laxos, de cara a poder normalizar más imágenes (máscara equilibrada y std de 22 y 33) 

 

 

 

 

Finalizada la Normalización:

  • Máscara de nubes (Fmask o BQA)
  • Máscara de nubes CM (valores 0 y 1)
  • img + hdr + doc de las bandas de la 2 a la 7 que se hayan podido normalizar (junto con el rel)
  • Archivo txt de los coeficientes de normalización
  • Registro en la base de datos de los coeficientes de normalización y del número de píxeles por PIA
{ 
    "_id" : "20140812l8oli202_34", 
    "Info" : {
        "Iniciada" : "Tue Nov 24 01:08:18 2015", 
        "Pasos" : {
            "rad" : {
                "Fecha" : "Tue Nov 24 01:17:22 2015", 
                "Kl-Values" : {
                    "b4" : NumberInt(5762), 
                    "b5" : NumberInt(4674), 
                    "b6" : NumberInt(4203), 
                    "b7" : NumberInt(4453), 
                    "b1" : NumberInt(9125), 
                    "b2" : NumberInt(8083), 
                    "b3" : NumberInt(6664), 
                    "b9" : NumberInt(4968)
                }, 
                "Corrad" : "True"
            }, 
            "geo" : {
                "Georef" : "True", 
                "Fecha" : "Tue Nov 24 01:14:02 2015"
            }, 
            "nor" : {
                "Normalize" : "True", 
                "Fecha" : "Tue Nov 24 01:18:01 2015", 
                "Nor-Values" : {
                    "b4" : {
                        "Parametros" : {
                            "slope" : 1.119179679624905, 
                            "r" : 0.9946001583686911, 
                            "intercept" : -0.818075379766352, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(56020)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(711), 
                            "Embalses" : NumberInt(2222), 
                            "Mineria" : NumberInt(163), 
                            "Mar" : NumberInt(49501), 
                            "Urbano-1" : NumberInt(628), 
                            "Urbano-2" : NumberInt(817), 
                            "Aeropuertos" : NumberInt(354), 
                            "Pastizales" : NumberInt(113), 
                            "Pinar" : NumberInt(1511)
                        }
                    }, 
                    "b5" : {
                        "Parametros" : {
                            "slope" : 1.4735184393927765, 
                            "r" : 0.9944821925427072, 
                            "intercept" : -8.002146493486533, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(54093)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(206), 
                            "Embalses" : NumberInt(1487), 
                            "Mineria" : NumberInt(180), 
                            "Mar" : NumberInt(49500), 
                            "Urbano-1" : NumberInt(424), 
                            "Urbano-2" : NumberInt(671), 
                            "Aeropuertos" : NumberInt(287), 
                            "Pastizales" : NumberInt(121), 
                            "Pinar" : NumberInt(1217)
                        }
                    }, 
                    "b6" : {
                        "Parametros" : {
                            "slope" : 1.3603867421884617, 
                            "r" : 0.9971326638449448, 
                            "intercept" : -5.338541949997632, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(55311)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(698), 
                            "Embalses" : NumberInt(1953), 
                            "Mineria" : NumberInt(194), 
                            "Mar" : NumberInt(49501), 
                            "Urbano-1" : NumberInt(599), 
                            "Urbano-2" : NumberInt(741), 
                            "Aeropuertos" : NumberInt(332), 
                            "Pastizales" : NumberInt(117), 
                            "Pinar" : NumberInt(1176)
                        }
                    }, 
                    "b7" : {
                        "Parametros" : {
                            "slope" : 1.5470858289888685, 
                            "r" : 0.9961333630704381, 
                            "intercept" : -1.4178915745814127, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(55864)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(683), 
                            "Embalses" : NumberInt(2284), 
                            "Mineria" : NumberInt(200), 
                            "Mar" : NumberInt(49501), 
                            "Urbano-1" : NumberInt(555), 
                            "Urbano-2" : NumberInt(713), 
                            "Aeropuertos" : NumberInt(398), 
                            "Pastizales" : NumberInt(108), 
                            "Pinar" : NumberInt(1422)
                        }
                    }, 
                    "b2" : {
                        "Parametros" : {
                            "slope" : 1.897435537990908, 
                            "r" : 0.9796848721026005, 
                            "intercept" : 0.931654933456219, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(53180)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(662), 
                            "Embalses" : NumberInt(1590), 
                            "Mineria" : NumberInt(80), 
                            "Mar" : NumberInt(47729), 
                            "Urbano-1" : NumberInt(489), 
                            "Urbano-2" : NumberInt(726), 
                            "Aeropuertos" : NumberInt(238), 
                            "Pastizales" : NumberInt(85), 
                            "Pinar" : NumberInt(1581)
                        }
                    }, 
                    "b3" : {
                        "Parametros" : {
                            "slope" : 1.4573014972232725, 
                            "r" : 0.9910393887505663, 
                            "intercept" : -2.521427924507977, 
                            "iter" : NumberInt(1), 
                            "N" : NumberInt(55643)
                        }, 
                        "Tipo_Area" : {
                            "Arena" : NumberInt(571), 
                            "Embalses" : NumberInt(2111), 
                            "Mineria" : NumberInt(125), 
                            "Mar" : NumberInt(49499), 
                            "Urbano-1" : NumberInt(532), 
                            "Urbano-2" : NumberInt(810), 
                            "Aeropuertos" : NumberInt(287), 
                            "Pastizales" : NumberInt(115), 
                            "Pinar" : NumberInt(1593)
                        }
                    }
                }
            }
        }, 
        "Tecnico" : "LAST-EBD Auto", 
        "Finalizada" : "Tue Nov 24 01:18:04 2015"
    }, 
    "Clouds" : {
        "cloud_scene" : 0.1, 
        "umbral" : NumberInt(50), 
        "cloud_PN" : 0.0
    }, 
    "usgs_id" : "LC82020342014224LGN00"
}

GitHub

Repositorio del código

ToDo List

  • Días:
    1. Añadir Landsat 7 dentro de la clase
    2. Quitar del rel de /nor las bandas que no se normalizan
    3. Fichero js para servidor
  • Quien sabe:
    1. Artículo
    2. Normalizar Landsat 8 con Landsat 8
    3. Versión sin Miramón

Va a ser difícil porqué...

Queda mucho por hacer

Así que el aplauso... que sea pequeño

Copy of Protocolo Automático

By Diego García Díaz

Copy of Protocolo Automático

Presentación del Protocolo Automático para el Tratamiento de imágenes Landsat del Laboratorio de SIG y Teledetección de la Estación Biológica de Doñana.

  • 410

More from Diego García Díaz