Corrección atmosférica absoluta: sustracción de objetos oscuros (DOS) en imágenes de satélite con lenguaje R

El método de sustracción de objetos oscuros (DOS) asume que si hay áreas en una imagen con valores muy bajos de reflectancia, cualquier reflectancia aparente debería ser producto de la dispersión (scatering) atmosférica y esta información puede ser usada para calibrar el resto de la imagen (Chavez, 1988; 1989). Los píxeles más oscuros pueden ser seleccionados examinando el histograma de los ND de una imagen o estableciendo un umbral, como por ejemplo, “los valores de ND más bajo encontrados por lo menos en n píxeles” o algún otro criterio apropiado al tamaño de la imagen que está siendo analizada. El valor de ND elegido, el “Starting Haze Value” (SHV), se convierte en radiancia (L) mediante cualquiera (la que se ajuste a la metadata de la imagen) de las siguientes ecuaciones:

Como es improbable que la mayoría de las imágenes contengan píxeles que sean verdaderamente negros, entonces la corrección se aplica asumiendo un porcentaje de reflectancia del 1 % en esas áreas; tal como lo señalan las siguientes expresiones:

La forma más simple de DOS convierte los valores de Lhaze calculados en reflectancias en el sensor:

y los sustrae de la imagen completa; también convertida en reflectancias en el sensor.

Como ejemplo para probar el método DOS, usé una serie de 7 bandas Landsat 5, las cuales fueron “recortadas” con un script de bash para ajustarlas todas por igual a un área que incluyera nubes, agua (lago de Valencia), ciudad y zonas montañosas. Para determinar el SHV de inicio se usó la primera banda (b1.tif). El procedimiento completo en R:

> b1<-rgdal::readGDAL("b1.tif")
> SHV<- table(b1@data[,1])
> SHV <- min(as.numeric(names(SHV)[SHV > 1000]))
> SHV

[1] 69

El valor mínimo obtenido, para cúmulos superiores a 1000 píxeles, fue de 69. Ahora hay que calcular nuevos valores de SHV para cada banda, según el procedimiento de Chavez (1988) y con la función DOS de landsat, los cuales están correlacionados según las condiciones atmosféricas. Para ello se requiere la metadata del archivo *.met que acompaña las imágenes (satélite Landsat 5, Gain = 0.6024314, Bias = -1.5200000, elevación del sol = 52.300000, fecha de adquisición de la imagen = 1986-03-13). El proceso completo se reseña a continuación:

> b1.DOS <- landsat::DOS(sat = 5, SHV = SHV, SHV.band = 1, Grescale =0.6024314, Brescale = -1.5200000, sunelev = 52.300000, edist = landsat::ESdist("1986-03-13"))
> b1.DOS <- b1.DOS[["DNfinal.mean"]]
> b1.DOS
         coef-4    coef-2   coef-1 coef-0.7 coef-0.5
band1 68.160162 68.160162 68.16016 68.16016 68.16016
band2 39.451800 51.756131 59.36949 61.87529 63.60689
band3 21.663036 37.967268 50.75640 55.42699 58.78936
band4 10.175615 24.934890 40.87729 47.58544 52.69736
band5  3.013091  8.194174 21.81642 30.37980 38.10901
band7  2.673985  5.670024 16.89510 25.19077 33.23688

La salida de arriba muestra los posible valores de SHV calculados para 5 coeficientes de dispersión; con los más pequeños (-4) describiendo las condiciones más claras y -0.5 representando bruma extrema. Como la banda 1 fue la empleada para determinar los SHV no se corrige con los coeficientes de dispersión. En términos generales, Chavez (1989) sugiere que valores de SHV menores o iguales a 55 corresponden a atmósfera muy clara (-4); de 56-75 clara (-2); de 76-95 moderada (-1); de 95-115 brumosa (-0.7); y > 115 muy brumosa (-0.5). El valor de SHV=69 para la banda 1 del día 13 de Marzo señala que la atmósfera era clara por lo que se selecciona toda la columna de coef -2 para ser usada en la función DOS() para las bandas restantes. En R ésto se hace de la siguiente manera:

b1.DOS<-b1.DOS[,2]
b2.DOSrefl <- landsat::radiocorr(b2, Grescale =0.6024314, Brescale = -1.5200000, sunelev = 52.300000, edist = landsat::ESdist("1986-03-13"), Esun =1795, Lhaze=b1.DOS[2], method="DOS")
rgdal::writeGDAL(b2.DOSrefl,"b2_dosrefl.tif", drivername="GTiff")

donde b2_dosrefl.tif es la banda 2 cuyos píxeles contienen información de reflectividad aparente para la cual se han sustraído los objetos oscuros. Un procedimiento similar se realiza para las bandas restantes.

About these ads
Esta entrada fue publicada en Lenguaje R, SIG, Teledetección. Guarda el enlace permanente.

6 respuestas a Corrección atmosférica absoluta: sustracción de objetos oscuros (DOS) en imágenes de satélite con lenguaje R

  1. Felicidades nuevamente, por introducir técnicamente en R las correcciones atmosféricas como DOS. También hay otras opciones como el algoritmo “6S code” de Vermote(?) la cual podemos acceder a este modulo desde GRASS-GIS, supongo que es posible ejecutarlo en “R” con spgrass6. Saludos

  2. daniel escargas dijo:

    Amigo, vengo aprendiendo mucho de tus posts, me surge una duda, veo que la funcion de Reflectancia en sensor que propones es la que propone Chander donde dice que “d” es la distancia “sol-tierra” en unidades astronomicas (que sale de una tabla en base al dia juliano), la duda es que Chuvieco propone para el calculo de Reflectancia Aparente (que entiendo que es lo mismo) la misma funcion donde en lugar de “d^2″ pone “k” y “k = (1+0,0167 (sen(2π (D – 93,5)/365)))^2″ donde D = Dia Juliano.
    Los valores que obtengo con el segundo metodo son practicamente 10 veces menores. Yo quiero los valores de reflectancia para calcular IVs (NDVI, SAVI y ARVI).
    Desde ya mil Gracias
    Daniel

  3. daniel escargas dijo:

    El tema de los valores de Chander y Chuvieco esta solucionado, resulta que la tercera edicion de Chuvieco tiene un error de parentesis de en la formula de K (abre 4 parentesis y cierra 3) y yo puse mal ese parentesis faltante, ahora consulte la segunda edicion y la funcion es “k = 1+0,0167 (sen(2π (D – 93,5)/365))^2″.
    Lo que hace mas similares los valores obtenidos por Chander y Chuvieco.
    Nuevamente gracias y disculpas por la molestia

  4. Andrés dijo:

    Buenos días. Me estoy iniciando en el tema del análisis de imágenes satelitales y he utilizado tu otro post en el cual calculas la temperatura al tope atmosférico mediante R. La cuestión es que como tal lo dije anteriormente la temperatura que obtengo dista mucho de ser la Temp. Superficial. Por esa razón supongo que debo hacer la corrección atmosférica. Al ingresar a la metadata de mis imagenes landsat TM5, descargadas de glovis, no existe el parametro GAINS para empezar a calcular L.

    Otra cuestión que me parece importante aclarar es que estoy trabajando con la banda infrarrojo térmica, por esa razón no se muy bien que arreglos debo hacer para obtener la temperatura superficial. Desde ya agradezco con anticipación la ayuda.

    Saludos

Deja un comentario

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s