Cálculo de la reflectancia aparente en imágenes de satélite con QGIS

En el artículo anterior señalé un método para determinar la reflectancia aparente en imágenes de satélite con R. A pesar de la potencia del lenguaje, que lo hace un procedimiento expedito, puede no parecer atractivo a muchos porque se maneja a través de líneas de comando; tanto en Linux como en Windows. Por esa razón, me dispuse a encontrar una vía alternativa usando la calculadora raster de QGIS; donde los valores de los números digitales (ND) obtenidos en R los empleé para calibrar el procedimiento en QGIS.

En un artículo precedente, referí que la determinación de la radiancia en el tope de la atmósfera (L) es el paso previo para la determinación de la reflectancia aparente en imágenes de satélite. La fórmula para el cálculo es la siguiente:

donde ρAS es la reflectancia aparente, d es la distancia Tierra-Sol (expresada en unidades astronómicas), Esun es la constante atmosférica extra-solar (Wm-2µm-1) y θz es el ángulo cenital (que se calcula a partir del ángulo de elevación del sol, θs, mediante θz = 90 – θs).

El valor de d, para la fecha de adquisición de la imagen (“1988-03-26″), lo determiné con R (comando landsat::ESdist(“1988-03-26”)). Se obtuvo un valor de 0.9971208. El Esun, constante atmosférica extra-solar (Wm-2µm-1), la obtuve de una tabla de Chander et al. (2009).

        Landsat 5     Landsat 7
Band 1    1983           1997
Band 2    1795           1812
Band 3    1539           1533
Band 4    1028           1039
Band 5     219.8          230.8
Band 7      83.49          84.90
Band 8      NA           1362

En mi caso, como es una banda 1 (Landsat 5), el valor de Esun = 1983. Finalmente, de la metadata se extrae que el ángulo de elevación solar es de 54.89 grados.

Todos estos valores combinados en la expresión de arriba, señalan que la constante:

K = π . d2 / Esun . Cos(90 – θs) = 0.00192549860081712

es el número que se empleará para multiplicar los valores de L para así obtener los de ρAS, es decir, se multiplica dicha banda (b1_esc_rad.tif) por 0.00192549860081712 con la calculadora raster de QGIS.

Como ya disponemos de los valores de L, basta multiplicarlos por 0.00192549860081712 para verificar si éstos son iguales a los que se obtuvieron con R. Efectivamente, la transformación pixel a pixel fue hecha de la manera esperada. A modo de ejemplo, para el primer pixel:

1.492156982421875 * 0.00192549860081712 = 0.00287314618185282

que concuerda, a partir de la décima cifra decimal (precisión que se usó para estimar la constante con LibreOffice), con el valor de 0.0028731462080031633377 obtenido con R para el primer pixel (e idéntico al que se encuentra en la imagen ASCII de dicha banda).

Por tanto, con QGIS, se carga la banda (b1_esc_rad.tif), se selecciona Raster -> Calculadora Raster, se introduce la fórmula, se coloca un nombre adecuado (b1_esc_ar.tif en este caso) y click para aceptar.

Esta entrada fue publicada en Software Libre, Teledetección. Guarda el enlace permanente.

28 respuestas a Cálculo de la reflectancia aparente en imágenes de satélite con QGIS

  1. Hola, tengo una inquietud, las unidades de los ángulos que las funciones del QGIS requiere son grados o radianes?

    Gracias.

  2. dario cantor dijo:

    estos valores de reflectivida aprente calculados son valores entre 0.01 y 0.04, como debo hacer para dejarlos en una imagen de 8 bits en valores de 0 a 255.

    saludos y gracias por el post

    • Usa gdal_translate en línea de comandos. Con man gdal_translate averiguas los parámetros. Fíjate lo que dice el extracto para el “scale” (lo que te interesa):

      -scale [src_min src_max [dst_min dst_max]]:
                 Rescale the input pixels values from the range src_min to src_max
                 to the range dst_min to dst_max. If omitted the output range is 0
                 to 255. If omitted the input range is automatically computed from
                 the source data.
      

      Si omites el intervalo de salida es de 0 a 255 y si omites el intervalo de entrada automáticamente calcula los datos a partir del ráster fuente, es decir, precisamente lo que necesitas.

      Aquí tienes una orientación:

      Escalamiento de imágenes georreferenciadas con gdal_translate (GDAL/OGR)

  3. daniel escargas dijo:

    Excelente aporte, me queda una duda respecto al angulo solar, ¿que valor tomar? porque en el encabezado de landsat tengo SUN_AZIMUTH = 60.0117949 SUN_ELEVATION = 47.4211734 y por otro lado me leí que use la calculadora http://solardat.uoregon.edu/cgi-bin/SolarPositionCalculator.cgi la que me arroja 76.99 y arriba refieres a restarle el valor a 90 con lo que tengo otro valor que tampoco coincide con ninguno de los anteriores. Agradecere me heches una mano

    • La calculadora arroja muchas cosas y tu no especificas. El coseno que hay que calcular es el del ángulo cenital. Por lo que señalas sería:

      cos(90 – 47.4211734) = cos(42,5788266) = 0,736347173

      • daniel escargas dijo:

        Muchas gracias, me estaba haciendo un embrollo pero esta respuesta viene a afirmar mi la idea que me venia formando. Gracias

  4. Mario Castañeda dijo:

    Hola, muy interesante la información. Tengo una duda, cuales son las unidades en las que se expresa la reflectancia??

    Saludos

  5. roman dijo:

    hola.. buena informacion..una pregunta..al convertir los ND, a reflectancia.. varios valores.. me salen negativos..y si los valores de reflectancia estan entre.. 0 y 1 q puedo hacer??.. por ejemplo.. el valor de reflectancia de un pixel me sale -0,002… debooo.. considerarlo como positivo o q?? ayuda por favor

    • Haz una estadística del ráster para que evalúes que fracción de los píxeles mantienen esa condición. Si es pequeña y de ese orden de magnitud (-0.002) puedes considerarlos cero y con álgebra de mapas convertirlos a ese valor (if mapa < 0 then mapa==0). Tampoco sería mala idea que revisaras si las fórmulas están bien escritas en el ráster calculator y probaras "manualmente" (en una hoja de cálculo o con una calculadora) con dos o tres ND seleccionados al azar si produce los valores esperados de reflectancia.

  6. Roman dijo:

    gracias.. y buenoo si he revisado mis formulas.. pero no hallo error,,.. peroo ahora veo q el signo negativo lo vengo jalando desde la convercion de ND a valores de radiancia.. siendo el “Lmin” = -1,52 (por ejemplo para la banda 1 del Landsat5), entonces.. todo ND<2 en el momento de la transfomacionn a valores de radiancia dara negativo… la formula q estoy utilizando para la radiancia es:

    L= ((Lmax-Lmin)/(Qmax-Qmin)) * ND + Lmin

    y para mi ejemplo.. el Lmax-Lmin/Qmax-Qmin = 0,762823529, entonces.el ND debera ser siempre mayo a 1 para q me de un valor positivoo… ya q si tengo un ND=1 el valor de radiancia me saldra L=0,757

    o es q ya en esta conversion de ND a radiancia(L) debo considerar q valor negativo de radiancia es igual a 0 (-0,757=0 ) .. por favor me espero tu comentariooo.. muchas gracias de antemano..

  7. Juan Carlos dijo:

    Hola que tal, gracias por tu información José. Me podrias pasar la bibliografia (articulos o libros) de donde sacaste toda esta información por favor. Buen dia

  8. paola tovar dijo:

    cual fue la formula que utilizo para calcular el parametro d para la fecha de adqusicion

  9. paola tovar dijo:

    estoy trabajando con imagenes spot 5 pero para calcular la radiancia no se que formula utilizar por que en los metadatos que me proporcionaron solo tengo el valor calibracion absoluta gains y estan en unidades w/m2/sr/um y dice que los offset son cero. No se que procedimiento seguir si tengo que convertir este valor a wm-2sr-1um-1 y cual formula utilizo

    • Esto confunde un poco porque existen varias fórmulas para los cálculos pero sólo una que se ajusta a la metadata. En la lámina 8 de esta referencia:

      http://www.rsgis.ait.ac.th/~honda/textbooks/advrs/DN2Reflectance_2p.pdf

      de acuerdo a lo que refieres, si tu offset es cero entonces multiplica el GAIN por los DN y tendrás los valores de radiancia. El offset es igual al bias y éste, a su vez, igual al Lmin. Después usa la fórmula que ya sabes para transformar en reflectancia absoluta.

      Toda esta “complicación” es por razones de almacenamiento. Es más barato hacerlo en Bytes (palabras de 8 bits). Dos elevado a la 8 es 256, por tanto, los DN varían entre 0 y 255.

  10. paola tovar dijo:

    muchas gracias por tu ayuda, la pregunta surgió porque estoy haciendo la corrección radiometrica en el model maker de erdas y me salio un error, entonces buscando los posibles errores busque formulas de radiancia para spot 5 y encontre este documento http://scielo.unam.mx/pdf/agro/v45n1/v45n1a10.pdf, donde dividen los ND/GAINS. otro error que pensé que podría ser son las unidades ya que la radiancia en unidades es Wm-2sr-1um-1 y yo tengo los de los GAINS en Wm2srum. El parametro de ESUN los descargue de esta pagina http://www.astrium-geo.com/es/918-resolucion-y-modos-espectrales y vienen en estas unidades (Wm2um). no se si el problema es de unidades y tengo que realizar la conversion la verdad no se como hacerla si tengo que pasar los GAINS y lo ESUN a Wm-2sr-1um-1 no se como hacer la conversión o si ese no es el problema. Gracias por tu ayuda no soy muy experta estoy aprendiendo y se consiguen muchas guias de lansadt y no muchas de spot.

  11. Mariana dijo:

    Podrías, José volver a cargar el link http://www.rsgis.ait.ac.th/~honda/textbooks/advrs/DN2Reflectance_2p.pdf nuevamente. Porque no puedo acceder a él y estoy necesitando calcular la reflectancia de una imagen SPOT 5 para luego, poder calcular indices de vegetación. Gracias!!! Muy buen post!

  12. Zenón Rizo Fernández dijo:

    Hola José Guerrero. Tengo una duda. Estoy calculando la reflectancia, con la misma formula, de unas imagenes pero me resulta con valores negativos en zonas de la imagen y, con valores positivos, en las zonas negras que contienen a dicha imagen. Es correcto o a que se debe?

    • No es correcto. Eso se debe a que no estás usando la metadata correspondiente a tus imágenes. Si son Landsat 5,7 hay un investigador de la NASA que publicó un artículo donde tiene allí resumidas todas las posibles combinaciones. Para otro tipo de imágenes habría que buscar en Internet.

  13. Hola buen día, necesito ayuda con urgencia. Tengo una imagen RapidEye a la que ya le calcule la radiancia y necesito calcular la reflectancia; estoy utilizando la misma formula, sin embargo los resultados me dan todos negativos, analice la formula parte por parte para ver en que parte me genera ese signo y encontré que el coseno de 26.73138 (angulo resultante de restar 90-angulo de elevación solar) es negativo…. que puedo hacer… es correcto que de este resultado?

  14. Pingback: Transformar día juliano en día calendario mediante Python | El Blog de José Guerrero

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

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