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.
El Blog de José Guerrero 
Pingback: Cálculo de la reflectancia aparente en imágenes de satélite con QGIS |
Pingback: Máscara de nubes usando lenguaje R y QGIS |
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_ATTRIBUTESY 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_RESCALINGDespué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.
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.
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”
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.
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:
Click to access news_0077.pdf
Saludos
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.
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:
Yo he obtenido reflectividades aparentes con R y usando el raster calculator de QGIS y dan exactamente lo mismo; variando entre 0-1.
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.