Cálculo de temperaturas a partir de la banda térmica de Landsat usando lenguaje R

En el artículo “Analyzing Remote Sensing Data in R: The landsat Package” existe un procedimiento para determinar las temperaturas, en ºK, en el tope de la atmósfera empleando la función thermalband del paquete landsat para R. Emplea la banda 6 (bandas 61 y 62 en Landsat 7) que corresponde al infrarrojo térmico. Como es sencillo, me avoqué a la tarea de aplicarlo en mi sistema Debian, sin embargo, al principio, la aplicación me arrojó un error diciendo que existía la imposibilidad de acceder a una pequeña porción de la imagen. El error persistía a pesar de probar con varias bandas 6 procedentes de diferentes sectores de Venezuela.

Tratando de averiguar de que se trataba el error, abrí con gimp la imagen de una de las bandas 6 y pude observar que ésta incluía un pequeña imagen sobrepuesta (no sé conque fin) sobre ella. Lo mismo sucedía para las otras bandas térmicas. Sospechando que aquí estaba el origen del error hice finalmente un subset de una las bandas térmicas y corrió sin problemas la aplicación.

El subset de la imagen fue hecho con QGIS; es decir, Raster -> Extraction -> Clipper para luego seguir las instrucciones para el corte incorporando el nombre del archivo de salida (b6-nuev.tif), la proyección y la zona de corte (fácil de establecer con dragado del ratón). Posteriormente, se cargó la imagen *.tif al intérprete de R con:

library(landsat)
b6<- readGDAL("b6-nuev.tif")

El objeto b6 es la imagen en formato SpatialGridDataFrame. Para obtener la banda térmica (b6.thermal; formato SpatialGridDataFrame):

b6.thermal <- thermalband(b6, band = 6)

y para exportarla a formato GeoTiff:

writeGDAL(b6.thermal,"/home/zeito/b6_thermal.tif",drivername="GTiff")

La imagen desplegada con QGIS, equalizada por ensanchamiento de contraste, mostrando las estadísticas de la misma:

El raster con temperaturas en ºC se puede obtener usando la calculadora raster para restar a la banda térmica en ºK la cantidad de 273.15. Las estadísticas de la banda en ºC, obtenidas con gdalinfo, son las siguientes:

Driver: GTiff/GeoTIFF
Files: b6_thermal-C.tif
       b6_thermal-C.tif.aux.xml
Size is 4457, 4003
Coordinate System is:
PROJCS["WGS 84 / UTM zone 19N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-69],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32619"]]
Origin = (627960.000000000000000,1029600.000000000000000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  627960.000, 1029600.000) ( 67d50'5.56"W,  9d18'44.60"N)
Lower Left  (  627960.000,  909510.000) ( 67d50'17.72"W,  8d13'34.96"N)
Upper Right (  761670.000, 1029600.000) ( 66d37'4.57"W,  9d18'22.81"N)
Lower Right (  761670.000,  909510.000) ( 66d37'29.43"W,  8d13'15.74"N)
Center      (  694815.000,  969555.000) ( 67d13'44.25"W,  8d46'1.31"N)
Band 1 Block=4457x1 Type=Float32, ColorInterp=Gray
  Min=9.451 Max=66.599 
  Minimum=9.451, Maximum=66.599, Mean=28.833, StdDev=2.848
  NoData Value=-3.4028234663852886e+38
  Metadata:
    STATISTICS_MINIMUM=9.4505615234375
    STATISTICS_MAXIMUM=66.598785400391
    STATISTICS_MEAN=28.833118577415
    STATISTICS_STDDEV=2.8480360652522

Al principio, me desconcertaron los valores extremos. Cuando se despliega, por ejemplo, una imagen con mapa de color con 5 clases:

queda delimitada la muy pequeña zona que presenta los mayores valores de temperatura (incluye dos clases). Corresponde a un incendio de vegetación muy localizado que puede visualizarse fácilmente en las bandas del visible. La zona de menor temperatura en la imagen es muy probable que correspondan a las nubes que se observan en la parte superior derecha. Por otra parte, en las otras dos clases, la de menor temperatura delimita la vegetación, cuerpos de agua y nubes y la más cálida las zonas desprovistas de vegetación y la infraestructura civil.

Otra imagen pero con pseudocolor:

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

3 respuestas a Cálculo de temperaturas a partir de la banda térmica de Landsat usando lenguaje R

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

  2. Pingback: Cálculo de temperaturas a partir de la banda térmica de Landsat usando lenguaje R | Geoprocessing - Geoprocessamento | Scoop.it

  3. Pingback: Cómo solventar la falta del operador logaritmo en el raster calculator de QGIS con gdal_calc.py | 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