PREPROCESAMIENTO DE IMÁGENES DE SATÉLITE: CORRECCIONES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla y Francisco Conde
2023-04-17
- 1 INTRODUCCION
- 2 CORRECCIÓN RADIOMÉTRICA.
- 2.1 Lectura del fichero de metadatos.
- 2.2 Importación de ficheros y creación
de un objeto
RasterStacka partir de los metadatos. - 2.3 Cálculo de la radiancia espectral de la parte superior de la atmósfera (TOA) manualmente.
- 2.4 Cálculo de la radiancia espectral
de la parte superior de la atmósfera (TOA) de manera automática con
RStoolbox. - 2.5 Conversión en reflectancia en lo alto de la atmósfera (“Top-of-the-atmosphere reflectance”, también “apparent reflectance”)
- 2.6 Exportar los resultados a un raster multibanda.
1 INTRODUCCION
Cuando se descargan las imágenes de los servidores como Earth Explorer, las imágenes sólo han sido sometidas a un preprocesamiento limitado. Este preprocesamiento suele incluir correcciones para la geometría del sensor y distorsiones geométricas inducidas por el terreno (comúnmente conocidas como ortorrectificación). Sin embargo, antes de utilizar esos datos, deben considerarse aplicar otro conjunto de correcciones adicionales.
Los sensores a borde de los satélites registran la intensidad de la radiación electromagnética emitida o reflejada en forma de números digitales (DN). La resolución radiométrica de Landsat5 es 8 bits, mientras que la de Landsat 8 es 16 bits. Esto significa que los valores de píxeles proporcionados como números enteros van de 1 a 255 en el caso de Landsat 5 (y entre 1 y 65535, respectivamente, para Landsat 8). Los DN de cada banda son una trasformación lineal de la radiación real registrada por el sensor para reducir el tamaño de los archivos.
Pero sus imágenes incorporan distorsiones causadas por los propios sensores, por el sol, la atmosféra o la topografía, entre otros elementos. Por lo tanto, debe realizarse un preprocesamiento incial para minimizar estos efectos, que en la mayoría de los casos, y antes de la transformación y/o clasificación de imágenes, se centrará en la corrección radiométrica y geométrica.
El propósito de la corrección radiométrica es corregir aquellos errores procedentes del propio sensor y de la atmósfera.
La corrección geométrica es un procedimiento para minimizar distorsiones debidas a variaciones en la geometría sensor-Tierra. Esta corrección ajusta las imágenes a las coordenadas del mundo real (p. ej., latitud y longitud) sobre la supeficie terrestre.
A continuación se presentarán diversas opciones para la calibración y corrección de imágenes de satélite. En algunos casos, esta fase puede omitirse; por ejemplo, para la clasificación de una sola escena es aceptable el uso directo de los DN. Pero si el objetivo es comparar los valores espectrales de diferentes escenas o la misma escena en varios momentos del tiempo, se debe aplicar la corrección de reflectancia de la parte superior de la atmósfera y, cuando sea posible, corregir atmosféricamente los datos de reflectancia en la superficie.
:::
Acceso a datos corregidos
El USGS también proporciona imágenes Landsat corregidas, procesadas a nivel de reflectancia en superficie. Puede ser conveniente usar directamente estos datos sin pasar por la fase citada anteriormente, ya que probablemente el producto final sea superior al que se puede obtener a partir de los métodos señalados en líneas precedentes, y además ahorran tiempo de procesamiento y espacio en el disco duro. Hay que tener en cuenta, sin embargo, que la fiabilidad de esa corrección depende de la presencia de la llamada “vegetación densa oscura”, por lo que no funcionan tan bien en áreas con poca o ninguna cubierta vegetal.
:::
2 CORRECCIÓN RADIOMÉTRICA.
La corrección radiométrica se realiza para eliminar o minimizar la contaminació producida por la atmósfera, ya que esta afecta a la radiancia registrada por el satélite a través de los procesos de dispersión, absorción y refracción. La decisión de realizar la corrección radiométrica depende de:
La naturaleza del problema en estudio;
El tipo de imágenes de sensores remotos disponibles;
La cantidad de datos atmosféricos in situ disponibles;
La precisión de la información de naturaleza biofísica que se obtendrá de las imágenes.
Existen cierto número de modelos de corrección radiométrica disponibles, unos de carácter absoluto y otros relativos (Du et al. 2002).
Los modelos absolutos son los más complejos, dado que se utilizan datos atmosféricos in situ y modelos de transferencia radiativa. Si bien estos modelos son los más precisos, dependen de la disponibilidad de los datos atmosféricos in situ, que no siempre disponibles. Ejemplos de estos modelos son el 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) y el MODTRAN (Moderate Resolution Atmospheric Radiance and Transmittance).
Los modelos relativos o basado en imágenes dependen de la información de la imagen digital, por lo que no requiere datos atmosféricos in situ. Se basan en la suposición de que las diferentes bandas de una misma imagen podrían compensar las distorsiones generadas por la atmósfera. Ejemplos de ello son el ajuste de histograma de una única imagen, o la normalización de imágenes de fechas múltiples a partir de histogramas o de regresión.
En esencia, el procedimiento comprende transformar los números digitales escalados en unidades de intensidad como la radiancia o la reflectancia. Este proceso de conversión de los números digitales en radiancia/reflectancia requiere información específica procedente de cada sensor, aplicando unos coeficientes de escalado radiométrico. En el caso de las imágenes de Landsat, los coeficientes de escalado radiométrico están incluidos en un archivo de metadatos (archivo .MTL), que también contiene las constantes necesarias para convertir los datos TIRS a la temperatura de brillo del satélite. No obstante,
De los productos Landsat ofreciedos por el USGS, el nivel 1 es el que contiene números digitales (DN) escalados (por lo general, enteros sin signo de 8 o 16 bits), que son los que habrán de ser convertidos en reflectancia en la parte superior de la atmósfera o en el sensor utilizando los metadatos proporcionados junto con el producto. El producto de nivel 2 son datos ya corregidos atmosféricamente (reflectancia de superficie, temperatura de la superficie terrestre) que están “listos para usar”, ya que no tiene que implementar ninguna corrección.
En este capítulo, realizamos la corrección radiométrica para Landsat utilizando los modelos disponibles en el paquete RStoolbox (Leutner y Horning 2016). Los detalles acerca de este paquete de R pueden ser consultados en la página web https://bleutner.github.io/RStoolbox/.
Proceso de calibración y corrección radiométrica
Para ello, el primer paso es cargar tanto el paquete raster como el paquete RStoolbox para iniciar el proceso.
library(raster) ## Loading required package: splibrary(RStoolbox) 2.1 Lectura del fichero de metadatos.
Iniciamos el procedimiento estableciendo la carpeta donde se encuentran todos los materiales necesarios para realizar la actividad
setwd("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/")
getwd()## [1] "D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite"Las imágenes de Landsat se acompañan de un conjunto de metadatos que proporcionan información esencial acerca del tiempo de adquisición, el nivel de procesamiento o la proyección de datos. Más importante aún, los metadatos también incluyen coeficientes de calibración y/o reescalado. Los metadatos están almacenados en un archivo de texto que termina en “_MTL.txt”, que se puede abrir con cualquier editor de texto.
En algunos casos, se pueden leer estos parámetros manualmente
abriendo el archivo de texto, aunque también puede ser importado
directamente a R con la función greadLines() (o
XML::xmlParse() para archivos XML). Sin embargo, la función
RStoolbox::readMeta() está diseñada para leer el fichero de
metadatos de Landsat:
fichero_mtl <-"D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/correccion_radiometrica/LC08_L1TP_201032_20210723_20210729_01_T1_MTL.txt"
fichero_metaData <- readMeta(fichero_mtl)El resumen del objeto metaData es ofrecer una visión
general de la información disponible para la escena correspondiente
(igualmente, la función str() se utiliza para ver qué
información está almacenada en dicho objeto).
summary(fichero_metaData)## Scene: LC82010322021204LGN00
## Satellite: LANDSAT8
## Sensor: OLI_TIRS
## Date: 2021-07-23
## Path/Row: 201/32
## Projection: +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
##
## Data:
## FILES QUANTITY CATEGORY
## B1_dn LC08_L1TP_201032_20210723_20210729_01_T1_B1.TIF dn image
## B2_dn LC08_L1TP_201032_20210723_20210729_01_T1_B2.TIF dn image
## B3_dn LC08_L1TP_201032_20210723_20210729_01_T1_B3.TIF dn image
## B4_dn LC08_L1TP_201032_20210723_20210729_01_T1_B4.TIF dn image
## B5_dn LC08_L1TP_201032_20210723_20210729_01_T1_B5.TIF dn image
## B6_dn LC08_L1TP_201032_20210723_20210729_01_T1_B6.TIF dn image
## B7_dn LC08_L1TP_201032_20210723_20210729_01_T1_B7.TIF dn image
## B9_dn LC08_L1TP_201032_20210723_20210729_01_T1_B9.TIF dn image
## B10_dn LC08_L1TP_201032_20210723_20210729_01_T1_B10.TIF dn image
## B11_dn LC08_L1TP_201032_20210723_20210729_01_T1_B11.TIF dn image
## B8_dn LC08_L1TP_201032_20210723_20210729_01_T1_B8.TIF dn pan
## QA_BQA LC08_L1TP_201032_20210723_20210729_01_T1_BQA.TIF dn qa
##
## Available calibration parameters (gain and offset):
## dn -> radiance (toa)
## dn -> reflectance (toa)
## dn -> brightness temperature (toa)El volumen de información que contiene es elevado. Ya se vió su contenido en el primer apartado. Además, también incluye la información necesaria para llevar a cabo el proceso de corrección, como es:
Ángulo del azimuth del sol.
Elevación del sol.
Distancia de la Tierra al Sol.
Valores “offset” y “gain” de cada banda para su conversión en radiancia, reflectancia y temperatura de brillo (sólo en las bandas térmicas)
2.2 Importación de
ficheros y creación de un objeto RasterStack a partir de
los metadatos.
A continuación, debe importase el fichero de metadatos con la función
RStoolbox::stackMeta(). Esta función facilita el
pre-procesado porque establece un enlace entre las bandas de la escena y
sus correspondientes parámetros de metadatos:
Esta función lee los ficheros de metadatos MTL o XML y, a
continuación, carga todos los ficheros creando un
RasterStack. Los argumentos pueden consultarse en la
siguiente página
web:
file Path to Landsat MTL metadata (_MTL.txt) file or an Landsat CDR xml metadata file (*.xml). quantity
quantity: qué se devuelve. Opciones: números digitales (‘dn’), reflectancia en lo alto de la atmósfera (‘tre’), reflectancia en superficie (‘sre’), temperatura de brillo (‘bt’), índice espectral (‘index’), todas las categorías (‘all’).
category categoría de datos a devolver (‘image’: una imagen, ‘pan’: imagen pancromática, ‘index’: índice multibandas, ‘qa’ bandas de calidad, ‘all’: todas las categorías).
allResolutions Si TRUE apilarían todas las capas, independientemente de su resolución.
lsat <- stackMeta(fichero_mtl,
quantity = "dn", # Parámetro números digitales
category = "image",
allResolutions = FALSE) Por si acaso, se define el valor no DATA en el fichero (DN de los valores ausentes es 0)
NAvalue(lsat) ## [1] -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -Inf -InfAtributos del RasterStack
lsat ## class : RasterStack
## dimensions : 7911, 7791, 61634601, 10 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_dn, B2_dn, B3_dn, B4_dn, B5_dn, B6_dn, B7_dn, B9_dn, B10_dn, B11_dn
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535Resumen de los valores
summary(lsat) ## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (0.16% of all cells)## B1_dn B2_dn B3_dn B4_dn B5_dn B6_dn B7_dn B9_dn B10_dn B11_dn
## Min. 0 0 0.00 0 0.00 0 0 0 0 0.0
## 1st Qu. 0 0 0.00 0 0.00 0 0 0 0 0.0
## Median 10883 10318 10068.00 10327 16049.00 16147 12224 5078 30912 27692.5
## 3rd Qu. 12016 11802 12250.75 13739 19024.75 20728 16029 5103 32902 29229.0
## Max. 48846 49505 47859.00 50522 54234.00 36620 31795 13825 37899 33112.0
## NA's 0 0 0.00 0 0.00 0 0 0 0 0.02.3 Cálculo de la radiancia espectral de la parte superior de la atmósfera (TOA) manualmente.
La radiancia espectral en el sensor (también llamada radiancia TOA -Top of the atmosphere-) es la cantidad de energía recibida por el sensor por segundo por estereorradián (sr) por metro cuadrado. Para calcular este valor necesitamos aplicar parámetros de calibración específicos del sensor y de la banda a los DN. Normalmente, estos parámetros consisten en:
Un término multiplicativo (a menudo denominado ganancia -“gain”-).
Un término aditivo (“offset”).
Ambos términos son diferentes para cada banda y deben aplicarse a cada pixel siguiendo la fórmula:
\[\lambda= M_L * Q_{Cal} + A_L\]
En esta ecuación,
\[M_L\] es el factor de reescalado multiplicativo específico de cada banda (en el fichero
.MTLaparece como MULTIPLICATIVE_ADD_BAND_x, donde x es el número de banda).\[A_L\] es el factor de reescalado aditivo específico de cada banda (en el fichero
.MTLaparece como RADIANCE_ADD_BAND_x, donde x es el número de banda).\[Q_{Cal}\] son los números digitales (DN), tal y cómo se captura la imagen sin procesar de Landsat.
Para conocer los valores específicos de estos factores de conversión podemos solicitarlos mediante la siguiente sintaxis:
radParameters <- fichero_metaData$CALRAD
radParameters ## offset gain
## B1_dn -60.83240 0.01216600
## B2_dn -62.29312 0.01245900
## B3_dn -57.40257 0.01148100
## B4_dn -48.40508 0.00968100
## B5_dn -29.62150 0.00592430
## B6_dn -7.36660 0.00147330
## B7_dn -2.48294 0.00049659
## B8_dn -54.78123 0.01095600
## B9_dn -11.57675 0.00231540
## B10_dn 0.10000 0.00033420
## B11_dn 0.10000 0.00033420Una vez disponemos de ellos, se puede aplicar a todas las bandas
lsat_rad <- lsat * radParameters$gain + radParameters$offset ## Warning in a[] <- e2: número de items para para sustituir no es un múltiplo de
## la longitud del reemplazo
## Warning in a[] <- e2: número de items para para sustituir no es un múltiplo de
## la longitud del reemplazo2.4 Cálculo de la
radiancia espectral de la parte superior de la atmósfera (TOA) de manera
automática con RStoolbox.
La función RStoolbox::radCor() es la herramienta
utilizada para implementar los diferentes métodos de corrección de las
imágenes de satélite.
Para ello, se seguirán varios pasos, que se explicitan en el
argumento method:
Conversión de valores de DN a radiación en el sensor (argumento
method = rad).Conversión a reflectancia TOA (es decir, reflectancia aparente; también conocida como TOA, top-of-the-atmosphere reflectance) utilizando el modelo de reflectancia aparente (argumento
method = apref).Conversión a reflectancia superficial comparando los resultados de varios métodos :
• El método de sustracción simple de objetos oscuros (simple dark object substraction, SDOS), que requiere previamente la corrección de la neblina (argumento
method = sdos).• El método de sustracción de objetos oscuros (dark object substraction, DOS), que asume que la dispersión es más alta en la banda azul y disminuye hacia el NIR (argumento
method = dos).• Por último, la corrección radiométrica utilizando el modelo COST (estimación del coseno de transmitancia atmosférica, argumento
method = costz).
Además, la función RStoolbox::radCor() incluye otros
argumentos de interés:
metaData: es el fichero de metadatos del que se extraen los factores de conversión.bandSet: permite seleccionar un conjunto específico de bandas.atmosphere: suministra información a la función sobre las condiciones atmosféricas o la cantidad de nubes durante el proceso de adquisión de la imagen.
lsat_rad <- radCor(lsat, # Objeto ráster
method = "rad", # Conversión de DN a radiancia
metaData = fichero_metaData, # Fichero de metadatos (*.MTL)
bandSet = c("B1_dn", "B2_dn", "B3_dn", "B4_dn", "B5_dn", "B6_dn", "B7_dn"), # Bandas
atmosphere = "veryClear")Comprobamos los atributos del nuevo fichero, por ejemplo el nuevo rango de valores para cada banda:
lsat_rad ## class : RasterStack
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_tra, B2_tra, B3_tra, B4_tra, B5_tra, B6_tra, B7_tra
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 736.46641, 746.64483, 695.00477, 586.03926, 358.62750, 89.18612, 30.06109La radiancia (la radiación que emite la Tierra, es decir, la que registra el sensor) se mide en vatios por metro cuadrado por estereorradián (\(Wm^2 sr^{(-1)}\)).Como consecuencia, los datos han pasado de números enteros a números flotantes:
dataType(lsat[[1]])## [1] "INT2U"dataType(lsat_rad[[1]])## [1] "FLT4S"A continuación, representaremos gráficamente una imagen en falso color de esta nueva imagen
plotRGB(lsat_rad,
r=5, g=4, b=3,
axis = TRUE,
stretch="lin")## Warning in plot.window(...): "axis" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "axis" is not a graphical parameter## Warning in title(...): "axis" is not a graphical parameter## Warning in graphics::rasterImage(z, bb[1], bb[3], bb[2], bb[4], interpolate =
## interpolate, : "axis" is not a graphical parameter
2.5 Conversión en reflectancia en lo alto de la atmósfera (“Top-of-the-atmosphere reflectance”, también “apparent reflectance”)
A continuación, se convertirá la radiación del sensor en una reflectancia espectral aparente, que tiene en cuenta los cambios temporales en la iluminación solar debido a la geometría de la Tierra y el Sol.
La reflectancia es una proporción adimensional (sin unidades) entre la emitancia de la radiación de un objeto y la irradiancia (Mather y Koch 2011). La emitancia de la radiación se refiere a la energía que se aleja de la superficie (por ejemplo, la energía térmica emitida por la Tierra), mientras que la irradiancia es la energía radiante que alcanza una superficie (Mather y Koch 2011). Si normalizamos la radiancia respecto a la irradiación solar aproximada en un día dado del año, obtenemos la reflectancia en la parte superior de la atmósfera (Top-of-atmosphere-reflectance), también conocida como reflectancia aparente, medida en unidades (radiancia * radiance-1) que van de 0 a 1.
Para su cálculo también se usa la función
RStoolbox::radCor() si bien el argumento
method se rellena con apref.
lsat_apref<- radCor(lsat, # Objeto raster
method = "apref", # Conversión de DN to reflectancia en lo alto de la atmósfera
metaData = fichero_metaData, # Fichero de metadatos
bandSet = c("B1_dn", "B2_dn", "B3_dn", "B4_dn", "B5_dn", "B6_dn", "B7_dn"), # Bandas
atmosphere = "veryClear")Comprobación de los atributos del objeto corregido
lsat_apref## class : RasterStack
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_tre, B2_tre, B3_tre, B4_tre, B5_tre, B6_tre, B7_tre
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1Además de calcular la reflectancia, la función radCor()
también puede convertir las bandas térmicas en temperatura de
brillo (en grados Kelvin), que es la temperatura que tendría un
cuerpo negro en función de la radiación medida. La temperatura de brillo
es una aproximación a la temperatura de la superficie terrestre. Para la
temperatura real de la superficie terrestre, la temperatura de brillo
debe ajustarse para la emisividad de la cubierta terrestre, que es
diferente para diferentes para agua, bosque etc.
2.5.1 Conversión a reflectancia de superficie.
La reflectancia TOA no tiene en cuenta ningún efecto atmosférico y es simplemente la relación entre la radiación en el sensor con la radiación entrante del sol, donde se consideran únicamente la distancia y el ángulo del sensor al objetivo. La radiación emitida por el sol llega a la atmósfera de la Tierra relativamente inalterada pero, a medida que atraviesa la atmósfera, se dispersa y es absorbida por los gases atmosféricos y los aerosoles, modificando sustancialmente el espectro de radiación. Adicionalmente, una vez que la radiación solar se refleja desde la superficie terrestre y vuelve al espacio, tiene que atravesar otra vez la atmósfera antes de alcanzar al sensor. A todo ello debe añadirse que también la radiación que no procede directamente del Sol o de la Tierra, sino que se origina en el entorno o en la atmósfera misma, llamada radiación de trayectoria (Path radiation).
El proceso de eliminación de estos efectos atmosféricos se denomina “corrección atmosférica”. Los métodos para la corrección atmosférica van desde simples correcciones basadas en imágenes hasta simulaciones altamente complejas de transferencia radiativa en la atmósfera.
El enfoque más simple se llama sustracción de objetos oscuros (DOS). Se basa en la idea de que siempre hay unos pocos píxeles que casi no contienen radiación (con una reflectividad de alrededor del 1%) y, por lo tanto, deberían ser oscuros. Sin embargo, debido a que hay una fuerte dispersión en las longitudes de onda cortas (azul, verde), cualquier radiación captada por esos píxeles oscuros se debe a la dispersión atmosférica (también llamada “haze”. En consecuencia, con una simple corrección atmosférica podemos restar ese efecto provocado por el “haze”, lo que mejoraría la nitidez de la imagen.
Se supone que cualquier radiación que se origine en los píxeles oscuros se debe a la neblina atmosférica y no a la reflectancia de la superficie. Los valores de neblina son estimaciones de la radiación de la trayectoria, que se pueden sustraer en un procedimiento de sustracción de objetos oscuros (DOS). La neblina atmosférica afecta casi exclusivamente al rango de longitud de onda visible. Por lo tanto, por lo general, solo querrá estimar la neblina en las bandas azul, verde y roja, pero ocasionalmente también en la banda NIR.
El enfoque de sustracción de objetos oscuros se basa en la estimación de la neblina atmosférica basada en píxeles oscuros. Se supone que los píxeles oscuros tienen una reflectancia cero, de ahí el nombre. Luego se supone además que cualquier radiación que se origine de tales píxeles oscuros se debe a la neblina atmosférica y no a la reflectancia de la superficie misma. Para la sustracción de objetos oscuros (DOS), usamos la función estimar Haze () para estimar el valor de píxel del número digital (DN) para objetos oscuros en el rango de longitud de onda visible y luego usamos los valores para convertir la reflectancia del sensor en reflectancia de la superficie.
Los métodos de corrección atmosférica (sdos, dos y costz) se aplican a la región óptica (solar) del espectro y no afectan a la banda térmica.
Existen diferentes métodos para realizar la eliminación del “haze”. Los siguientes métodos están disponibles en la función radCor():
dos - Sustracción de objetos oscuros siguiendo a Chávez (1989) costz - Sustracción de objetos oscuros siguiendo a Chávez (1996) sdos - Sustracción simple de objetos oscuros.
Para la sustracción de objetos oscuros (DOS), primero tenemos que
estimar el valor de píxel del número digital (DN) de los objetos oscuros
para el rango de longitud de onda visible usando la función estimar
Haze(), luego usar estos valores para convertir la
reflectancia del sensor en reflectancia de la superficie.
2.5.1.1 Método sdos
En general, se supone que cualquier radiación que se origine en píxeles oscuros (es decir, con reflectancia cercana a cero) se debe a la neblina atmosférica (“atmospheric haze”) y no a la reflectancia de la superficie. Por lo tanto, los valores correspondientes a esta neblina son estimaciones de la radiación de la trayectoria. SDOS requiere la estimación de los valores de turbidez en las bandas visible y NIR. Esto se debe a que la neblina atmosférica afecta principalmente al rango de longitud de onda visible (bandas azul, verde y roja). Los siguientes son pasos clave para el modelo SDOS.
El más simple (argumento method = "sdos") requiere
obtener un valor de “haze” por banda; este valor se sustrae de todos los
píxeles en la imagen. Para estimar el valor del “haze” se usa la función
RStoolbox::estimateHaze(), tomando los DN y especificando
la proporción de píxeles oscuros que esperamos en el image (argumento
darkProp). La proporción de “haze” sólo se estima en las
bandas 1-4 (longitud de onda visible y NTR), ya que su efecto en
longitudes de onda > 800 nm es insignificante.
En primer lugar, se efectuará una estimación automática de la neblina utilizando la función RStoolbox::estimarHaze(), que estima el valor (en términos de DN) de los píxeles correspondientes a los objetos oscuros para las bandas visible y NIR. La sintaxis es la siguiente:
x Un objeto raster
hazeBands incluye el número de bandas o el nombre de las bandas para las que se quiere estimar la neblina atmosférica.
darkProp es un número con la proporción de píxeles que se estima son oscuros.
hazeDN <- estimateHaze(lsat, # Objeto ráster
hazeBands = 2:5, # Bandas a calcular
plot = TRUE, # Muestra un histograma y los valores de "haze"
darkProp=0.01) # Proporción de píxeles para estimar el "haze"
A continuación, se efectua la corrección radiométrica utilizando el método SDOS, en el que el “valor de neblina inicial” (Start Haze Valu) en el que se basará la corrección atmosférica se ha calculado anteriormente.
lsat_sref_sdos <- radCor(lsat, # Objeto ráster
metaData = fichero_metaData, # Fichero de metadatos
method = "sdos", # Método de corrección
hazeValues = hazeDN, # Objeto con la estimación del haze
hazeBands = 2:5, # Bandas a calcular
atmosphere = "veryClear") # Características atmosféricasComprobación de los atributos del objeto corregido con el método SDOS
lsat_sref_sdos## class : RasterStack
## dimensions : 7911, 7791, 61634601, 10 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre, B9_sre, B10_bt, B11_bt
## min values : 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 147.5171, 141.6706
## max values : 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 0.231737, 322.676098, 318.302567Visualización
plotRGB(lsat_sref_sdos, r=5, g=4, b=3, stretch="lin")
2.5.1.2 Método dos
Otro método (method = "dos") se basa en un modelo de
dispersión atmosférica que asume que la dispersión es más alta en las
longitudes de onda azules y disminuye exponencialmente hacia el NIR;
así, estimando el “haze” en la banda azul se puede calcular para el
resto de bandas. A diferencia del enfoque simple “sdos” este método
requiere parámetros específicos del sensor, es decir, sólo es válido
sólo para el sensor que desea corregir.
El procedimiento DOS utiliza un modelo de reducción de la neblina atmosférica, que asume que la dispersión es más alta en la longitud de onda azul y disminuye hacia el NIR. Los pasos para realizarlo son los son los siguientes.
De igual manera, los números digitales de todas las capas raster en
el raster stack pueden convertirse en valores de reflectancia en
superficie mediante el argumento dos. Este método calcula
la reflectancia en superficie usando la ténica de substracción de
objetos oscuros (“Darck Object substraction) para corrección
del”scattering” de la atmósfera, diseñada por Chavez (1989).
lsat_srf_dos <- radCor(lsat, # Objeto ráster
method = "dos", # Método de corrección
metaData = fichero_metaData, # Fichero de metadatos
darkProp = 0.01,
bandSet = c("B1_dn", "B2_dn", "B3_dn", "B4_dn", "B5_dn", "B6_dn", "B7_dn"), # Bandas
atmosphere = "veryClear")Comprobación de los atributos del objeto corregido con el método DOS. Obsérvese que, en comparación con los resultados del modelos anterior (el SDOS) existe un ligero cambio en los valores de reflectancia.
lsat_srf_dos## class : RasterStack
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1Visualización
plotRGB(lsat_srf_dos, r=5, g=4, b=3, stretch="lin")
2.5.1.3 Método costz (Cosine estimation of atmospheric transmittance)
A continuación, realicemos la corrección radiométrica especificando method = “costz”, acompañado del objeto en el que se ha calculado el valor de la neblina atmosférica (hazeValues = hazeDN) y de las bandas que se van a corregir.
lsat_srf_costz <- radCor(lsat, # Objeto ráster
method = "costz", # Método de corrección
metaData = fichero_metaData, # Fichero de metadatos
darkProp = 0.01,
bandSet = c("B1_dn", "B2_dn", "B3_dn", "B4_dn", "B5_dn", "B6_dn", "B7_dn"), # Bandas
atmosphere = "veryClear")Comprobación de los atributos de las imágenes Landsat corregidas con el método COST
lsat_srf_costz## class : RasterStack
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1El resultado muestra que hay un cambio significativo en los valores máximos de reflectancia de las bandas 1–3.
Visualización
plotRGB(lsat_srf_costz, r=5, g=4, b=3, stretch="lin")
2.6 Exportar los resultados a un raster multibanda.
A continuación se muestran las bandas del rasterStack
con su denominación, que se mantendrá para su posterior grabación.
lsat_srf_costz## class : RasterStack
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## names : B1_sre, B2_sre, B3_sre, B4_sre, B5_sre, B6_sre, B7_sre
## min values : 0, 0, 0, 0, 0, 0, 0
## max values : 1, 1, 1, 1, 1, 1, 1Para grabar por separado las diferentes bandas de un objeto
rasterStack podemos utilizar la siguiente sintaxis:
writeRaster(lsat_srf_costz,
filename= "LC08_L1TP_201032_20210723_20210729_SR.tif",
bylayer=TRUE, # Crea un fichero separado para cada capa
format="GTiff",
overwrite=TRUE)