Cálculo de la reflectancia aparente en imágenes de satélite con R (Debian testing)

En el artículo anterior:

Analizando datos de Sensores remotos en R: paquete landsat (Debian testing)

refería el uso del lenguaje R para analizar datos de sensores remotos pero más bien se hacía énfasis en el modo de acceder a las aplicaciones (las que se probaron fueron muy sencillas y no tenían nada que ver con imágenes). Ahora, voy a exponer la forma de determinar la reflectancia aparente en imágenes de satélite con R en mi Debian testing. Para ello me basé en el artículo de Sarah Goslee y en la documentación de rgdal. Lo primero que hay que hacer es cargar la imagen en formato GeoTiff al espacio de trabajo en R pero en el formato que éste interpreta que es SpatialGridDataFrame. El comando ejecutado en R fue:

library(landsat)
b1<- readGDAL("b1_esc.tif")

La imagen en R es ahora el “objeto” b1 (una banda 1 escalada de 16 a 8 bits). De la metadata de la imagen para dicha banda obtuve:

Grescale = 0.6024314, Brescale = -1.5200000, ELEVACIÓN DEL SOL = 54.890000, FECHA DE ADQUISICIÓN DE LA IMAGEN = “1988-03-26”

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.

Obtenidos ya todos los parámetros necesarios, se ejecutó en el intérprete el siguiente comando:

b1.ar <- radiocorr(b1, Grescale = 0.6024314, Brescale = -1.5200000, sunelev = 54.89, edist = ESdist("1988-03-26"), Esun = 1983, method = "apparentreflectance")

La imagen en formato SpatialGridDataFrame (objeto b1.ar) necesita ahora ser convertida nuevamente en *.tif. Para ello usé el siguiente comando de rgdal:

writeGDAL(b1.ar,"/home/zeito/b1_esc_ar.tif",drivername="GTiff")

La imagen, cargada en QGIS con ensanchamiento de contraste a MinMax:

preservó su georreferenciación.

La imagen con valores de reflectancia aparente la transformé en ASCII con gdal_translate y un extracto del archivo es éste:

ncols        5055
nrows        5158
xllcorner    600915.000000000000
yllcorner    879945.000000000000
cellsize     30.000000000000
 0.0028731462080031633377 0.0028731462080031633377 0.0017131654312834143639 0.0017131654312834143639 0.0028731462080031633377 0.0028731462080031633377 0.0017131654312834143639 0.0017131654312834143639 0.0028731462080031633377 0.0028731462080031633377 0.0028731462080031633377 0.0017131654312834143639 0.0028731462080031633377 0.0017131654312834143639 0.0028731462080031633377 0.0028731462080031633377 0.0028731462080031633377 0.0051931077614426612854 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0051931077614426612854 0.0040331268683075904846 0.0051931077614426612854 0.0051931077614426612854 0.0051931077614426612854 0.0040331268683075904846 0.0051931077614426612854 0.0051931077614426612854 0.0051931077614426612854 0.0040331268683075904846 0.0040331268683075904846 0.0051931077614426612854 0.0063530886545777320862 0.0063530886545777320862 0.0063530886545777320862 0.0051931077614426612854 0.0063530886545777320862 0.0063530886545777320862 0.0063530886545777320862 0.0051931077614426612854 0.0051931077614426612854 0.0063530886545777320862 0.0051931077614426612854 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0028731462080031633377 0.0028731462080031633377 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0040331268683075904846 0.0051931077614426612854 0.0063530886545777320862 0.0063530886545777320862 0.0063530886545777320862 ...

obtenido simplemente para calibrar el procedimiento en QGIS.

Esta entrada fue publicada en Debian testing, GDAL, Lenguaje R, Teledetección. Guarda el enlace permanente.

12 respuestas a Cálculo de la reflectancia aparente en imágenes de satélite con R (Debian testing)

  1. Pingback: Cálculo de la reflectancia aparente en imágenes de satélite con QGIS |

  2. Pingback: Máscara de nubes usando lenguaje R y QGIS |

  3. Adhemar dijo:

    Hola José,

    Muchas gracias por responder a mi comentario. Por su puesto que antes de calcular la “reflectancia aparente” he revisado todos los artículos que has posteado respecto al tema, y si no hubiera sido así, tal vez no me hubiera enterado que este tipo de calculo es posible hacerlo en R.
    Bueno, te dejo el código que tengo para R:
    (Es una imagen LandsatTM 5, adquirida el 14-07-1986, esta es la banda 2)

    1. b2 <- readGDAL("LT50010721986195XXX02_B2_19s.tif")

    2. b2.ar <- radiocorr(b2, Grescale = 1.3220, Brescale = -4.16220, sunelev = 34.52449119, edist = ESdist("1986-07-14"), Esun = 1795, method = "apparentreflectance")

    3. writeGDAL(b2.ar, "b2_1986_001072_19s_apprefl.tif", drivername="GTiff")

    Es cierto, talvez no esté utilizando bien el metadato de la imagen (en este caso tengo un archivo MTL.txt), sobre todo tengo mis dudas en los valores de Grescale y Brescale.

    Los valores de la metadata que utilizo son:

    GROUP = IMAGE_ATTRIBUTES
        CLOUD_COVER = 0.00
        IMAGE_QUALITY = 9
        SUN_AZIMUTH = 46.81770736
        SUN_ELEVATION = 34.52449119
        GROUND_CONTROL_POINTS_MODEL = 161
        GEOMETRIC_RMSE_MODEL = 4.010
        GEOMETRIC_RMSE_MODEL_Y = 3.104
        GEOMETRIC_RMSE_MODEL_X = 2.539
        GROUND_CONTROL_POINTS_VERIFY = 3212
        GEOMETRIC_RMSE_VERIFY = 0.194
        GEOMETRIC_RMSE_VERIFY_QUAD_UL = 0.203
        GEOMETRIC_RMSE_VERIFY_QUAD_UR = 0.204
        GEOMETRIC_RMSE_VERIFY_QUAD_LL = 0.195
        GEOMETRIC_RMSE_VERIFY_QUAD_LR = 0.174
      END_GROUP = IMAGE_ATTRIBUTES
    

    Y para los valores de Grescale y Brescale.

    GROUP = RADIOMETRIC_RESCALING
        RADIANCE_MULT_BAND_1 = 0.671
        RADIANCE_MULT_BAND_2 = 1.322
        RADIANCE_MULT_BAND_3 = 1.044
        RADIANCE_MULT_BAND_4 = 0.876
        RADIANCE_MULT_BAND_5 = 0.120
        RADIANCE_MULT_BAND_6 = 0.055
        RADIANCE_MULT_BAND_7 = 0.066
        RADIANCE_ADD_BAND_1 = -2.19134
        RADIANCE_ADD_BAND_2 = -4.16220
        RADIANCE_ADD_BAND_3 = -2.21398
        RADIANCE_ADD_BAND_4 = -2.38602
        RADIANCE_ADD_BAND_5 = -0.49035
        RADIANCE_ADD_BAND_6 = 1.18243
        RADIANCE_ADD_BAND_7 = -0.21555
      END_GROUP = RADIOMETRIC_RESCALING
    

    Después de haber ejecutado el código obtengo una imagen de 216,6 MB (la imagen original ocupa un espacio de 54,2 MB).
    Para verificar los valores del rates utilizo GDAL:

    gdalinfo -stats b2_1986_001072_19s_apprefl.tif

    Band 1 Block=7801×1 Type=Float32, ColorInterp=Gray
    Minimum=-0.013, Maximum=1.062, Mean=0.088, StdDev=0.073
    Metadata:
    STATISTICS_MAXIMUM=1.0624703168869
    STATISTICS_MEAN=0.087802657331985
    STATISTICS_MINIMUM=-0.013282002881169
    STATISTICS_STDDEV=0.072701669810476

    Bueno, hasta qui es todo lo que utilizo.
    Gracias por revisar mi comentario.

    Saludos,
    Adhemar.

    • Búscate un archivo *.met que es donde están el GAIN y el BIAS. A eso se refiere el Grescale y el Brescale; respectivamente.

      • Adhemar dijo:

        Ok, ya he buscado este archivo *.met y no lo tengo, el archivo de la metadata de la imagen es este archivo MTL.txt (me parece que el metadato de las imágenes ya no vienen en *.met). En un principio estaba utilizando ENVI, y es en este programa donde he encontrado los valores de “ganancias y perdidas” GAIN y BIAS. He buscado estos valores en el archivo del matadato y son los único que coinciden con GROUP = RADIOMETRIC_RESCALING.

    • Adhemar dijo:

      Revisando mejor los metadatos _MTL.txt tengo esta columna de datos:

      GROUP = PRODUCT_PARAMETERS
      CORRECTION_GAIN_BAND_1 = “CPF”
      CORRECTION_GAIN_BAND_2 = “CPF”
      CORRECTION_GAIN_BAND_3 = “CPF”
      CORRECTION_GAIN_BAND_4 = “CPF”
      CORRECTION_GAIN_BAND_5 = “CPF”
      CORRECTION_GAIN_BAND_6 = “INTERNAL_CALIBRATION”
      CORRECTION_GAIN_BAND_7 = “CPF”
      CORRECTION_BIAS_BAND_1 = “CPF”
      CORRECTION_BIAS_BAND_2 = “CPF”
      CORRECTION_BIAS_BAND_3 = “CPF”
      CORRECTION_BIAS_BAND_4 = “CPF”
      CORRECTION_BIAS_BAND_5 = “CPF”
      CORRECTION_BIAS_BAND_6 = “CPF”
      CORRECTION_BIAS_BAND_7 = “CPF”

  4. Adhemar dijo:

    Y para GRASS utilizo lo siguiente:

    El código para GRASS es el siguiente:

    i.landsat.toar input_prefix=LT50010721986195XXX02_B output_prefix=LT50010721986195XXX02_toar_B metfile=D:\Imagenes_satelitales\ES_1986_L5\LT50010721986195XXX02\LT50010721986195XXX02_MTL.txt

    Según el manual de grass solo debo poner el nombre base de mis archivos de entrada y salida: (input_prefix=LT50010721986195XXX02_B output_prefix=LT50010721986195XXX02_toar_B),
    y después dar la dirección del archivo de “metadato”.

    Después que termina todo el cálculo ejecuto la función:

    r.info LT50010721986195XXX02_toar_b2

    y los valores que obtengo con grass son más o menos similares a los obtenidos con R

    Range of data: min = -0.00890854221307425 max = 1.04455794258934 |

    Debo decirte también que el calculo realizado con R lo hice con un sistema Windows 8 y OS X, los resultados son similares.

    Y para GRASS he utilizado Windows 8 (GRASS 7) y Ubuntu (GRASS 6.4)
    En todos obtengo los mismos valores.
    +—————————————————————————-+
    | Layer: LT50010721986195XXX02_toar_b2 Date: Tue Mar 25 17:51:24 2014 |
    | Mapset: L5_ES_1986_001072 Login of Creator: confade |
    | Location: Location_MDP_19N |
    | DataBase: D:\GrassData_MDP |
    | Title: ( LT50010721986195XXX02_toar_b2 ) |
    | Timestamp: none |
    |—————————————————————————-|
    | |
    | Type of Map: raster Number of Categories: 0 |
    | Data Type: DCELL |
    | Data Units: unitless Vertical datum: (none) |
    | Rows: 6941 |
    | Columns: 7801 |
    | Total Cells: 54146741 |
    | Projection: UTM (zone 19) |
    | N: -1813485 S: -2021715 Res: 30 |
    | E: 683415 W: 449385 Res: 30 |
    | Range of data: min = -0.00890854221307425 max = 1.04455794258934 |
    | |
    | Data Description: |
    | generated by i.landsat.toar |
    | |
    | Comments: |
    | Reflectance of Landsat-5 TM (method uncorrected) |
    | —————————————————————– |
    | Acquisition date (and time) ……….. 1986-07-14 (-0.0000 h) |
    | Production date ………………….. 2013-11-20 |
    | |
    | Earth-sun distance (d) ……………. -0.0000000 |
    | Sun elevation (and azimuth) ……….. 0.00000 (0.00000) |
    | Digital number (DN) range …………. 0 to 0 |
    | Calibration constants (Lmin to Lmax) .. +0.00000 to +0.00000 |
    | DN to Radiance (gain and bias) …….. +0.00000 and +0.00000 |
    | Mean solar irradiance (ESUN) ………. 0.000 |
    | Radiance to Reflectance (divide by) … +0.00000 |
    | —————————————————————— |
    | |
    | i.landsat.toar input_prefix=”LT50010721986195XXX02_B” output_prefix=\ |
    | “LT50010721986195XXX02_toar_b” metfile=”D:\Imagenes_satelitales\ES_1\ |
    | 986_L5\LT50010721986195XXX02\LT50010721986195XXX02_MTL.txt” method=”\ |
    | uncorrected” percent=0.01 pixel=1000 rayleigh=0.0 |
    | |
    +—————————————————————————-+

    Saludos,
    Adhemar.

  5. Ya encontré donde está el error. Los datos de Grescale están bien más no así los de Brescale. Por los datos de Grescale tu imagen fue adquirida entre Marzo de 1984 y Diciembre de 1991 pero procesada posterior al 2007. Por tanto, el Brescale tiene que ser de -2.84 y no de -4.16220 como tu haz colocado. Puedes corroborarlo en la Tabla IV de esta referencia:

    http://landsat.gsfc.nasa.gov/wp-content/uploads/2012/08/news_0077.pdf

    Saludos

  6. Adhemar dijo:

    Ok, voy a revisar el dato y te estaré comentando cuales son los resultados. Efectivamente la fecha de adquisición de la imagen es en 1984, pero la fecha de procesamiento de la imagen es del 2013, ha sido cuando hice el pedido de las imágenes.

    Saludos,
    Adhemar.

  7. Adhemar dijo:

    Bueno, ya modifiqué el valor del Brescale (-2.84) y el problema aún persiste.

    • Volví a recrear el procedimiento, ahora con unas imágenes de Perú, y las reflectividades aparentes oscilaron entre 0-1. Esta fue la metadata que tuve que ubicar para ejecutar radiocorr:

      Grescale = 0.766
      Brescale = -1.52
      SUN_ELEVATION = 42.53582961
      DATE_ACQUIRED = 2005-07-16
      Esun = 1983
      

      Yo he obtenido reflectividades aparentes con R y usando el raster calculator de QGIS y dan exactamente lo mismo; variando entre 0-1.

  8. Adhemar dijo:

    Muchas gracias José por tu ayuda. No sé que es lo que estoy haciendo mal, pero los resultados me siguen arrojando valores negativos y en algunas bandas mayores a 1.

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