💡 MATERIALES PARA EL EJERCICIO

escena_landsat. Esta escena corresponde una escena del satélite Landsat 8.

script_laboratorio_2_landsat

DESCARGA DE IMÁGENES LANDSAT

Plataformas para la descarga de imágenes

Existen varias plataformas web a través de las que se puede descargar escenas de Landsat. Por ejemplo, se pueden citar las siguientes:

Descarga de imágenes mediante la plataforma del USGS

El USGS pone gratuitamente a disposición de la comunidad científica tanto imágenes históricas como otras más recientes a través del sitio web https://earthexplorer.usgs.gov/. Se mostrará a continuación el procedimiento para la descarga en esta plataforma.

Registro

Para descargar escenas de Landsat, debe registrarse (de forma gratuita) en el servidor https://earthexplorer.usgs.gov/.

Pantalla con login USGS EarthExplorer
Pantalla con login USGS EarthExplorer

Una vez que haya completado el registro, debe recibir un correo electrónico en su dirección de correo electrónico proporcionada para confirmar su cuenta.

Después de la confirmación, será redirigido a la página de inicio de sesión donde deberá completar su nombre de usuario y contraseña (o hacer clic en “iniciar sesión” en la barra de menú del encabezado).

La interfaz de usuario

EarthExplorer admite la búsqueda en línea en bases de datos integrales, visualizaciones de vista rápida, exportación de metadatos y descarga de datos de ciencias de la tierra de los archivos del Servicio Geológico de EE. UU. (USGS).

Puede acceder al USGS EarthExplorer a través de la siguiente URL utilizando cualquier navegador web: https://earthexplorer.usgs.gov/

Una vez que se hace clic, se debe cargar la interfaz gráfica de usuario (GUI) principal de EarthExplorer, que se compone de tres elementos clave:

  • Barra de menú de encabezado: Botones para los servicios de inicio de sesión y registro, así como funcionalidades de ayuda, RSS y comentarios. Después de iniciar sesión, puede guardar y cargar consultas aquí

  • Barra lateral de búsqueda de datos: los componentes de búsqueda se dividen en cuatro pestañas y le permiten ingresar criterios de búsqueda, seleccionar conjuntos de datos para consultar, ingresar criterios adicionales y revisar los resultados en una ventana tabular

  • Vista de imagen con elementos de navegación: componentes integrados de Google Maps para visualizar resultados de búsqueda, con herramientas de navegación estándar de Google Maps, es decir, acercar/alejar, vista de la calle (esquina inferior derecha) e información de coordenadas de la posición actual del cursor (esquina superior derecha) ). Puede alternar entre la vista de imágenes satelitales y la vista de datos GIS seleccionando el botón adecuado en la esquina superior izquierda.

Interface de usuario del USGS EarthExplorer
Interface de usuario del USGS EarthExplorer

Selección de imágenes

Después de completar el registro, inicie con login. La aplicación pedirá su nombre de usuario (username) y contraseña (password)

Pantalla inicial tras entrar con nombre usuario y contraseña
Pantalla inicial tras entrar con nombre usuario y contraseña

EarthExplorer proporciona cuatro pestañas en el procedimiento de búsqueda para guiarlo a través de su solicitud de búsqueda:

EarthExplorer search tabs
EarthExplorer search tabs

El siguiente paso es definir los criterios de búsqueda, como la región de interés (ROI), las fechas preferidas y el número de escenas preferidas. Estas preferencias se pueden establecer en la pestaña ‘Criterios de búsqueda’.

Localización

Hay varias otras formas de definir el área temática con mayor precisión, por ejemplo, mediante archivos de forma, características, áreas predefinidas o kml. El ROI puede ser definido por

  • Escribir dirección/nombre del lugar.

  • Dibujando interactivamente el ROI haciendo clic en el mapa. Para obtener las coordenadas precisas, haga clic en Agregar coordenadas y agregue las coordenadas.

  • Proporcionando un archivo .kml o .shp del ROI.

La forma más sencilla de definir el AOI es utilizar Google Maps integrado escribiendo una dirección o el nombre de un lugar, por ejemplo, “Madrid, Casa de Campo”, haga clic y elija una de las sugerencias solicitadas. También son posibles las designaciones para la longitud geográfica y la latitud de la posición deseada, por ejemplo, “52.5194, 13.4067”.

Además, existe una notación global utilizada para catalogar los datos de Landsat, llamada Sistema de referencia mundial (WRS), en la que Landsat 8 sigue al WRS-2. Este sistema divide la superficie de la Tierra en las geometrías de registro de las adquisiciones de Landsat. Hay un mapa general WRS-2 y un convertidor de ruta/fila a latitud/longitud WRS-2 proporcionado por el USGS. Especialmente el convertidor ayuda a encontrar todas las combinaciones posibles de Path-Row para su AOI, que es 192/024 y 193/23 en nuestro caso para Berlín. Ingrese uno de esos pares de Ruta/Fila y haga clic para agregar el centro adecuado de la adquisición de Landsat como una coordenada a su interfaz de Google Map:

EarthExplorer search tabs
EarthExplorer search tabs

Para eliminar cualquier coordenada dada, presione la cruz roja.

Otra manera fácil de definir su AOI es simplemente hacer clic con el botón izquierdo dentro del mapa, lo que automáticamente agrega una coordenada para una búsqueda de una sola coordenada. Al definir dos puntos en el mapa, realizará una búsqueda de línea, lo que da como resultado todos los productos de datos que intersecan la línea. Al definir tres o más coordenadas, se muestra automáticamente un polígono, formando su AOI:

Ejemplo de uso de Path and Row para la búsqueda de una localización
Ejemplo de uso de Path and Row para la búsqueda de una localización
Polígono definido por cuatro coordenadas haciendo click en el botón izquerdo (interfaz Google Map)
Polígono definido por cuatro coordenadas haciendo click en el botón izquerdo (interfaz Google Map)

Finalmente, también proporcionando un archivo .kml o .shp del ROI. Haga clic en Cargar KML/Shapefile, seleccione Shapefile y haga clic en Seleccionar archivo. Localice el archivo de formas y haga clic en Abrir. El archivo de forma debe estar en formato archivado y no tener más de 30 nodos. Haga clic en Cerrar y aparecerá el ROI con nodos en el mapa.

Periodo temporal.

El criterio de fecha se puede dar haciendo clic en la pestaña “Rango de fechas”. Defina las fechas de búsqueda Desde y Hasta. Los meses preferidos también pueden ser definido

Para que su consulta sea más concreta, puede definir el lapso de tiempo dentro del cual desea obtener datos en la parte inferior de la pestaña Criterios de búsqueda. Simplemente configure la fecha de inicio y la fecha de finalización, así como todos los meses deseados, en nuestro caso, junio, julio y agosto:

Selección del periodo temporal de interés
Selección del periodo temporal de interés

La pestaña Conjunto de datos clasifica los conjuntos de datos en colecciones de datos similares. Hay una estructura de árbol dinámica, que le permite expandir/contraer productos presionando los signos más y menos al lado. Como puede ver, hay una gran cantidad de datos para elegir. Después de definir los criterios de búsqueda, haga clic en la pestaña “Conjuntos de datos” para seleccionar los datos. Aquí se incluye una lista de conjuntos de datos disponibles de forma gratuita, incluidos Landsat, Sentinel, SRTM DEM, MODIS, AVHRR, etc. Estaremos descargando los últimos datos Landsat-8 OLI/TIRS.

Selección del producto Landsat
Selección del producto Landsat

Criterios adicionales (opcional): esta pestaña ayuda a reducir aún más los resultados de su consulta de búsqueda al definir criterios de búsqueda adicionales, por ejemplo, la proporción permitida de cobertura de nubes sobre la tierra, cobertura absoluta de nubes, día o noche, así como restricciones de ruta/fila. Además, puede usar la identificación única del producto (por ejemplo, LC08_L1TP_192023_20170830_20170914_01_T1) para encontrar escenas individuales específicas. Restrinja la cobertura de nubes terrestres a “Menos del 10 %”. Debido a que Landsat solo pasa sobre una punto de planeta cada 16 días, es bastante probable que coincida con algún día nuboso. Por lo tanto, encontrar una escena sin nubes es bastante complicado, especialmente durante el invierno y en climas como el de la propia Cantabria.

Pestaña para seleccionar criterios adicionales: cubierta de nubes
Pestaña para seleccionar criterios adicionales: cubierta de nubes

Para visualizar las escenas disponibles con criterios de búsqueda definidos, haga clic en la pestaña Resultados. Cuando esté configurado, haga clic en en la parte inferior de la barra lateral de búsqueda de datos o en la pestaña Resultados en la parte superior para ejecutar su investigación. Lo más probable es que obtenga cuatro productos de datos como resultado (dependiendo de la forma de su AOI):

Resultados de la búsqueda
Resultados de la búsqueda

Cada producto recibe una identificación única, así como un tiempo de adquisición y la ruta y fila WGS-2.

Además, hay una serie de controles de superposición y descarga que puede elegir para cada escena:

Controles Overlay y download
Controles Overlay y download
  • Show footprint: muestra el contorno de una escena en Google Map en color verdadero.

  • Show Browse Overlay: muestra una imagen de vista previa de la escena.

  • Compare Browse: active este botón en varias escenas y luego superpóngalas aquí

  • Show Metadata and Browse: muestra la imagen de exploración y los metadatos completos de la escena seleccionada

  • Download Options: permite a los usuarios registrados descargar los datos seleccionados.

  • Add to Bulk Download: permite a los usuarios registrados descargar masivamente los datos seleccionados

  • Order Scene: permite a los usuarios registrados ordenar o solicitar un procesamiento especializado de productos

  • Exclude Scene from Results: elimina la escena en particular de la ventana de resultados actual.

Visualización de las escenas disponibles

Antes de la descarga definitiva de las escenas, podría ser conveniente una primera selección. A pesar de su relativamente baja resolución temporal, la cantidad de imágenes comienza a ser importante, por lo que es conveniente seleccionar las de mayor interés. Esto se puede conseguir fácilmente buscan y descargando un fichero *.csv que proporciona la misma página web. Este fichero contiene la lista de escenas de Landsat que se han seleccionado. Dicho fichero debe descargarse en el disco duro y descomprimirse. Con la primera línea de código se identifican todos los ficheros que cumplen las condiciones (argumento pattern) en el directorio de trabajo.

La función RStoolbox::readEE() proporciona una tabla con los resultados de la búsqueda y agrega columns útiles para su posterior representación gráfica.

library(RStoolbox)
## Warning: package 'RStoolbox' was built under R version 4.3.3
## This is version 1.0.0 of RStoolbox
ee <- readEE(system.file("external/EarthExplorer_LS8.txt", package = "RStoolbox"))

Si quisiéramos identificar qué escenas seleccionadas tienen una cobertura nubosa inferior al 20% el código sería el siguiente:

ee[ee$Cloud.Cover < 20,]
##    Landsat.Scene.Identifier WRS.Path WRS.Row Data.Category Cloud.Cover
## 6     LC82240632015157LGN00      224      63       NOMINAL       13.20
## 27    LC82240632014186LGN00      224      63       NOMINAL       15.94
## 40    LC82240632013327LGN00      224      63       NOMINAL       13.19
## 46    LC82240632013215LGN00      224      63       NOMINAL       14.72
## 49    LC82240632013167LGN00      224      63       NOMINAL        6.35
##    Station.Identifier Day.Night Data.Type.Level.1 Date.Acquired Sun.Elevation
## 6                 LGN       DAY               L1T    2015/06/06      51.98871
## 27                LGN       DAY               L1T    2014/07/05      51.01591
## 40                LGN       DAY              L1GT    2013/11/23      61.83434
## 46                LGN       DAY               L1T    2013/08/03      54.43226
## 49                LGN       DAY               L1T    2013/06/16      51.64735
##    Sun.Azimuth Geometric.RMSE.Model.X Geometric.RMSE.Model.Y
## 6     43.59274                  6.244                  5.663
## 27    44.79093                  5.452                  5.420
## 40   126.95937                  0.000                  0.000
## 46    51.68815                  5.631                  5.840
## 49    42.59918                  6.377                  5.862
##               Display.ID           Ordering.ID
## 6  LC82240632015157LGN00 LC82240632015157LGN00
## 27 LC82240632014186LGN00 LC82240632014186LGN00
## 40 LC82240632013327LGN00 LC82240632013327LGN00
## 46 LC82240632013215LGN00 LC82240632013215LGN00
## 49 LC82240632013167LGN00 LC82240632013167LGN00
##                                                                Download.Link
## 6  http://earthexplorer.usgs.gov/download/options/4923/LC82240632015157LGN00
## 27 http://earthexplorer.usgs.gov/download/options/4923/LC82240632014186LGN00
## 40 http://earthexplorer.usgs.gov/download/options/4923/LC82240632013327LGN00
## 46 http://earthexplorer.usgs.gov/download/options/4923/LC82240632013215LGN00
## 49 http://earthexplorer.usgs.gov/download/options/4923/LC82240632013167LGN00
##    Browse.Link       Date Doy Year Satellite Num
## 6           NA 2015-06-06 157 2015       LS8   8
## 27          NA 2014-07-05 186 2014       LS8   8
## 40          NA 2013-11-23 327 2013       LS8   8
## 46          NA 2013-08-03 215 2013       LS8   8
## 49          NA 2013-06-16 167 2013       LS8   8

Una vez creado el objeto ee, se pueden dibujar las escenas de Landsat disponibles con el paquete ggplot2. Si juega con diferentes niveles de cobertura de nubes máxima, se puede conocer qué imágenes están realmente disponibles y cuáles son las variaciones estacionales de las mismas. Es evidente que en nuestras latitudes las posibilidades de obtener imágenes libres de nubes son menores durante el invierno.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
ggplot(ee) + 
     geom_segment(aes(x = Date, xend = Date, y = 0, yend = 100 - Cloud.Cover, 
     col = as.factor(Year))) +
        scale_y_continuous(name = "Scene quality (% clear sky)")

Descarga del fichero con las escenas seleccionadas

Haga clic en el icono Descargar de la escena seleccionada y seleccione el producto Geo Tiff de nivel 1 para descargar el archivo . Archivos tif de todas las bandas. Para descargar una sola escena de Nivel 1 (consulte el capítulo Landsat 8), haga clic en Opciones de descarga y elija el último elemento, que debe ser el archivo más grande (aproximadamente 700-900 MB para Landsat 8):

Opciones de descarga
Opciones de descarga

Recibirá un archivo comprimido, que contiene todas las bandas espectrales como archivos geotiff georreferenciados.

Añadir escenas al carrito
Añadir escenas al carrito

Cuando haya terminado, haga clic en Cesta de artículos en la barra de menú del encabezado para enviar su pedido. Verá una lista de todas las escenas seleccionadas. Confirme su selección presionando Proceder al pago:

Revisión de la orden
Revisión de la orden

En la siguiente pantalla, presione Enviar pedido:

Enviar la petición
Enviar la petición

Se le dará una identificación de pedido única y se le enviará un correo electrónico de confirmación a su dirección de correo electrónico.

Confirmación de la orden
Confirmación de la orden

Todos los pedidos de nivel 2 enviados a través de ESPA se procesan en un plazo de 2 a 5 días, según el tamaño del pedido y la acumulación que ya se encuentra en el sistema. Se enviará un segundo correo electrónico de confirmación cuando los productos estén listos para descargar. A partir de este momento, todas las escenas permanecerán disponibles durante 10 días. Haga clic en la URL de estado del pedido en su segundo correo electrónico de confirmación, que lo redirige al siguiente sitio web:

Order overview given by the link of the second email
Order overview given by the link of the second email

Haga clic en su ID de pedido, que lo lleva al sitio de descarga. Simplemente haga clic en el enlace de descarga para obtener sus datos:

Ventana final
Ventana final
Nomenclatura de una imagen
Nomenclatura de una imagen

El ID de una escena se interpreta así:

► En el primer conjunto la letra L indica Landsat, seguido del sensor, que puede ser: C = OLI/TIRS combinados, 0 = OLI solamente, T = TIRS solamente. Los dos números siguientes indican el satélite, donde 09 = Landsat 9.

► El segundo conjunto hace referencia al nivel de corrección del procesamiento, que puede ser L1TP- Nivel 1 de precisión del terreno (corregido) o L1GT-Nivel 1 Terreno sistemático (corregido).

► El conjunto sucesor tiene 6 números, los primeros 3 indican el Path, y los 3 últimos el Row; para este ejemplo el Path/Row 008055.

PROCESAMIENTO BÁSICO DE ESCENAS LANDSAT CON R

La escena que servirá para adiestrar a los alumnos en el manejo de las escenas de Landsat es una escena perteneciente a la colección 2 Nivel 1. Los datos son número digitales (DN) con un formato de 16-bit. Estos datos pueden ser convertidos a reflectancia en el techo de la atmósfera (TOA, Top of Atmosphere) o radiancia usando los factores de escalado radiométrico que aparecen en el fichero de metadatos.

El producto Landsat descargado de EarthExplorer es un fichero comprimido .tar. Es recomendable crear un nuevo directorio de trabajo, donde se descomprimirá el citado fichero. Este fichero comprimido tiene un formato *.tar, cuyo contenido puede extraerse utilizando un programa como 7-zip en un entorno Windows, pero también es posible ejecutar todo el proceso en R, siguiendo las siguientes instrucciones.


En la actividad anterior se ha creado la estructura de carpetas y se han descomprimido los ficheros que contienen las imágenes. Se puede comprobar que se ha realizado correctamente listando todos los ficheros del directorio de trabajo.

list.files("D:/G174_2025/lab_2/ejercicio_landsat/")

Una vez decomprimidos, cada una de las bandas landsat individuales de una escena es un fichero en formato GeoTIFF.

La escena descargada como ejemplo corresponde al satélite Landsat 8 (LC08), colección 2 nivel 1 (L1TP), path 201 y fila 032 (201032) tomada el 23 de julio (07) de 2021 (20210723) Tier 1. Contiene tres tipos de archivos:

.TIF: Archivos de imágenes.

.txt: Archivos de Metadatos.

.xml: Archivos de Metadatos en formato xml.

json: Metadatos en formato json

jpeg: Archivo para la previsualización de la escena

Archivos con imágenes

Dentro del conjunto de ficheros en formato .tif se encuentran los ficheros correspondientes a las 11 bandas espectrales disponibles en esa plataforma. Sus características son las siguientes:

  • Banda 1: Permite visualizar imágenes de aguas poco profundas y partículas finas como polvo y humo. Denominación de Intervalo: Aerosol / Costero. Ancho (µm): 0.43 - 0.45

  • Banda 2 : Permite detectar cuerpos de agua (mapeo batimétrico), delimitar costas, diferenciación entre suelo y vegetación, diferenciación entre la vegetación conífera y decidua, detección de rasgos urbanos, vías y construcciones. Denominación de Intervalo: Azul. Ancho (µm): 0.45 - 0.51

  • Banda 3: Permite evaluar el vigor de la vegetación sana, diferenciar tipos de rocas, delinear aguas poco profundas, medir la calidad de agua (discrimina sedimentos en suspensión), rasgos urbanos y de infraestructura. Denominación de Intervalo: Verde. Ancho (µm): 0.53 - 0.59

  • Banda 4: Permite determinar la absorción de clorofila, por ello es muy útil para la clasificación de la cubierta vegetal, contrasta áreas con y sin vegetación, delimita áreas agrícolas y urbanas. Denominación de Intervalo: Rojo. Ancho (µm): 0.64 - 0.67

  • Banda 5: Las plantas saludables lo reflejan. Permite el cálculo de biomasa. Al compararlo con otras bandas, se obtienen índices como NDVI, que permite medir la salud de la planta con mayor precisión. Además, permite delimitar costas, para diferenciación suelos–cultivos y suelos–agua, para geomorfología, suelos y geología. Denominación de Intervalo: Infrarrojo Cercano (NIR). Ancho (µm): 0.85 - 0.88

  • Banda 6: Permite diferenciar la tierra húmeda de la seca, diferenciar entre nubes, nieve y hielo. Denominación de Intervalo: Infrarrojo de onda corta 1 (SWIR - 1). Ancho (µm): 1.57 - 1.65

  • Banda 7: Es útil para geología, puesto que rocas y suelos que parecen similares en otras franjas a menudo tienen fuertes contrastes en esta banda. Además mejora la determinación de contenidos de humedad en suelos y vegetación. Denominación de Intervalo: Infrarrojo de onda corta 2 (SWIR - 2). Ancho (µm): 2.11 - 2.29

  • Banda 8: Actúa como una película en blanco y negro, ya que en lugar de coleccionar colores visibles por separado, los combina en un solo canal. Es la más nítida de todas las bandas. Se utiliza para crear un mayor contraste entre áreas con y sin cubierta vegetal. Denominación de Intervalo: Pancromático (Pan). Ancho (µm): 0.50 - 0.68

  • Banda 9: Es la banda que muestra menos. Cubre una porción muy fina de longitudes de onda. Permite detectar nubes, está diseñada especialmente para cirrus. Los cirrus debido a sus bordes suaves son difíciles de detectar, y pueden esconder información importante. Denominación de Intervalo: Cirrus. Ancho (µm): 1.36 - 1.38

  • Banda 10: es la banda térmica. Denominación de Intervalo: Infrarrojo térmico 1 (TIR - 1). Ancho (µm): 10.6 - 11.19

  • Banda 11: Empleado para el mapeo termal mejorado y estimación de humedad del suelo. Denominación de Intervalo: Infrarrojo térmico 2 (TIR - 2). Ancho (µm): 11.5 - 12.51

Bandas de Landsat 8
Bandas de Landsat 8

Archivos adicionales

Además, aparecen otra serie de ficheros GeoTIFF adicionales con información más detallada:

  • Ficheros con la iluminación solar y el coeficiente angular del sensor (_SAA, _SZA, _VAA y _VZA).

  • Fichero de control de calidad (_QA)

Por último, también se incluye otros ficheros como

  • Ficheros con metadatos (_MTL)

  • Fichero con los coeficientes angulares (_ANG)

Se puede comprobar esta relación de listando todos los ficheros del directorio de trabajo.

list.files("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/")
##  [1] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"         
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"          
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"         
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"         
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"          
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"          
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"          
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"          
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"          
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"          
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"          
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"          
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.json"        
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"         
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"         
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF"    
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"   
## [18] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"         
## [19] "LC08_L1TP_201032_20210723_20210729_02_T1_stac.json"       
## [20] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"         
## [21] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_large.jpeg"
## [22] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_small.jpeg"
## [23] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"         
## [24] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"

Como alternativa, también se puede utilizar el siguiente código.

dir(path = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/")
##  [1] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"         
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"          
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"         
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"         
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"          
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"          
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"          
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"          
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"          
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"          
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"          
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"          
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.json"        
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"         
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"         
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF"    
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"   
## [18] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"         
## [19] "LC08_L1TP_201032_20210723_20210729_02_T1_stac.json"       
## [20] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"         
## [21] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_large.jpeg"
## [22] "LC08_L1TP_201032_20210723_20210729_02_T1_thumb_small.jpeg"
## [23] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"         
## [24] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"

También podemos crear listas por tipo de archivo:

tifs <- dir(path = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/", 
            pattern ="TIF")
tifs
##  [1] "LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF"       
##  [2] "LC08_L1TP_201032_20210723_20210729_02_T1_B10.TIF"      
##  [3] "LC08_L1TP_201032_20210723_20210729_02_T1_B11.TIF"      
##  [4] "LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF"       
##  [5] "LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF"       
##  [6] "LC08_L1TP_201032_20210723_20210729_02_T1_B4.TIF"       
##  [7] "LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF"       
##  [8] "LC08_L1TP_201032_20210723_20210729_02_T1_B6.TIF"       
##  [9] "LC08_L1TP_201032_20210723_20210729_02_T1_B7.TIF"       
## [10] "LC08_L1TP_201032_20210723_20210729_02_T1_B8.TIF"       
## [11] "LC08_L1TP_201032_20210723_20210729_02_T1_B9.TIF"       
## [12] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.TIF" 
## [13] "LC08_L1TP_201032_20210723_20210729_02_T1_QA_RADSAT.TIF"
## [14] "LC08_L1TP_201032_20210723_20210729_02_T1_SAA.TIF"      
## [15] "LC08_L1TP_201032_20210723_20210729_02_T1_SZA.TIF"      
## [16] "LC08_L1TP_201032_20210723_20210729_02_T1_VAA.TIF"      
## [17] "LC08_L1TP_201032_20210723_20210729_02_T1_VZA.TIF"
txts <- dir(path = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/", 
            pattern ="txt")
txts
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1_ANG.txt"
## [2] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt"
xmls <- dir(path = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/", 
            pattern ="xml")
xmls
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1_MTL.xml"
Raster monobanda vs raster multibanda
Raster monobanda vs raster multibanda
https://www.neonscience.org/resources/learning-hub/tutorials/dc-multiband-rasters-r

Archivos con metadatos

Son fundamentalmente 2, uno en formato .MTL y otro con formato .XML. Estos ficheros contienen la siguiente información:

  • Nombre del satélite y del sensor.

  • Identificador (ID) de la escena.

  • Fecha de adquisición y procesado de la escena.

  • Número de pasada y de fila.

  • Proyección.

Además, también incluye la información necesaria para llevar a cabo el proceso de corrección, como es:

  • Ángulo del azimuth del sol.

  • Elevación del sol.

  • Distancia de la Tierra al Sol.

  • Bandas disponibles.

  • Resolución espacial y radiométrica de cada banda.

Es importante tener acceso a los metadatos de la imagen para extraer información de la escena. R cuenta con una familia de funciones para la carga de archivos ASCII. La siguiente configuración nos permite cargar los datos existente en el fichero *_MTL.txt`

parametros <- read.csv(file = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_MTL.txt",
                       header = FALSE,
                       sep = "=", 
                       stringsAsFactors = FALSE, 
                       strip.white = TRUE)

names(parametros) <- c('id','valor')

Para visualizar el contenido de dicho fichero se utiliza la función View().

View(parametros)

Del conjunto de datos disponibles, podemos revisar algunos interesantes para el análisis posterior:

• Porcentaje de cobertura Nubosa de la imagen:

cloud <- parametros[which(parametros$id == "CLOUD_COVER"), 2]
cloud
## [1] "2.97"

• Porcentaje de cobertura nubosa sobre superficie terrestre:

cloud_land <- parametros[which(parametros$id == "CLOUD_COVER_LAND"), 2]
cloud_land
## [1] "2.97"

• Fecha de toma de la imagen:

fecha <- parametros[which(parametros$id == "DATE_ACQUIRED"), 2]
fecha
## [1] "2021-07-23"

• Nombre oficial de la imagen; se usa el índice [1] ya que existe varias veces en el archivo:

id <- parametros[which(parametros$id == "LANDSAT_PRODUCT_ID"), 2][1]
id
## [1] "LC08_L1TP_201032_20210723_20210729_02_T1"

Archivos de control de calidad

Uno de los ficheros más importantes del conjunto es la imagen de control de calidad de los píxeles. Esta banda, en formato *.tif, presenta la etiqueta QA_PIXEL. Para abrir y leer este archivos necesitamos cargar el paquete terra()

library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.21
qa_px <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL.tif")

En relación al contenido de la banda podemos revisar sus valores y revisar su significado

unique(qa_px)
##    LC08_L1TP_201032_20210723_20210729_02_T1_QA_PIXEL
## 1                                                  1
## 2                                              21762
## 3                                              21824
## 4                                              21890
## 5                                              21952
## 6                                              22018
## 7                                              22080
## 8                                              22146
## 9                                              22280
## 10                                             23826
## 11                                             23888
## 12                                             24082
## 13                                             24144
## 14                                             54534
## 15                                             54596
## 16                                             54790
## 17                                             54852
## 18                                             55052
## 19                                             56598
## 20                                             56660
## 21                                             56854
## 22                                             56916
Valores de la banda de calidad
Valores de la banda de calidad

Esos datos los podríamos utilizar como máscara para futuros procesos:

mascara <- qa_px
mascara[mascara != 21824] <- NA                 # Si el valor del pixel no es igual a  21824 se transforma en valor NA

Esta máscara se puede

plot(mascara)

Archivo de imagen color natural

Los archivos de imagen de color natural están disponibles directamente como archivos GeoTIFF (.tif) que no requieren procesamiento adicional. Puede ser representada mediante el paquete magick que tiene capacidad para leer y modificar (rotar, escalar, recortar, aplicar efectos, etc) diferentes formatos de imagen (png, jpeg…). A continuación, con la función magick::image_read() se representará gráficamente esa imagen en color verdadero.

library(magick)
## Warning: package 'magick' was built under R version 4.3.3
## Linking to ImageMagick 6.9.12.98
## Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fontconfig, x11
landsat_color_natural <- image_read("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_thumb_large.jpeg")

landsat_color_natural

IMPORTAR IMÁGENES DE LANDSAT EN R

A partir de este momento, nos centraremos en el manejo de las imágenes de satélite que corresponden a las bandas de Landsat 8. Para importar cada banda de manera individual (en este caso, la banda NIR) la sintaxis es la siguiente

imagen <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF")

Al citar el nombre del objeto, R devuelve un resumen con información sobre dicho objeto.

imagen
## class       : SpatRaster 
## dimensions  : 7911, 7791, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source      : LC08_L1TP_201032_20210723_20210729_02_T1_B5.TIF 
## name        : LC08_L1TP_201032_20210723_20210729_02_T1_B5

La información proporcionada (visto en un apartado anterior) es la siguiente:

La importanción individual de diferentes bandas de la escena se puede realizar varias maneras. Una de ellas es utilizar la función rast() de manera consecutiva.

temp1 <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF") 
temp2 <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF") 
temp3 <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF") 

No obstante, lo habitual es convertir esos ficheros en un archivo multibanda, mediante la función concatenar c(). Esta última función también puede ser aplicada a imágenes ya cargadas tanto como las almacenadas en el disco duro.

multibanda <- rast(c(temp1, temp2, temp3))
multibanda
## class       : SpatRaster 
## dimensions  : 7911, 7791, 3  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630)

Una vez ilustrado este ejemplo, y dado que no necesitaremos más los objetos anteriores, se eliminan del Global Environment.

rm(temp1, temp2, temp3, multibanda)

El procedimiento de importanción anterior es farrogoso y propenso a errores, dado que es fácil equivocarse por la complejidad de los nombres de los ficheros originales. Existen varias alternativas para solventar este problema.

• Una primera es crear una lista con las bandas, aprovechando que comparten la misma denominación (salvo el número de banda) y que estan alojadas en la misma carpeta. Una vez creada la lista, será cargados con la función rast().

lista.escena.completa <- list.files("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/",               # Directorio de los datos
                              pattern = ".*B[1234567]\\.tif$",                  # El signo $ lista todos los ficheros que acaban igual
                              ignore.case=TRUE,                                 # Sensible a mayúsculas y minúsculas
                              full.names = TRUE)
escena.completa <- rast(lista.escena.completa)
escena.completa
## class       : SpatRaster 
## dimensions  : 7911, 7791, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## sources     : LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF  
##               ... and 4 more sources
## names       : LC08_~T1_B1, LC08_~T1_B2, LC08_~T1_B3, LC08_~T1_B4, LC08_~T1_B5, LC08_~T1_B6, ...

• Otro procedimiento alternativo para crear esa lista es el siguiente

lista.escena.completa2 <- paste0("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/datos/ejercicio_landsat/LC08_L1TP_201032_20210723_20210729_02_T1_B", 
                  1:7,
                  pattern = ".TIF")
escena.completa2 <- rast(lista.escena.completa2)
escena.completa2
## class       : SpatRaster 
## dimensions  : 7911, 7791, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## sources     : LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF  
##               ... and 4 more sources
## names       : LC08_~T1_B1, LC08_~T1_B2, LC08_~T1_B3, LC08_~T1_B4, LC08_~T1_B5, LC08_~T1_B6, ...

En ambos casos, sólo se han cargado las primeras 7 bandas, dado que las bandas 8, 9, 10 y 11 tienen una resolución espacial diferente al resto.

Además, y dado que los nombres de banda son bastante largos y redundantes, es conveniente acortarlos, pero sin que pierdan significado. El uso de nombres simples facilita el acceso a las capas individuales. Existen varias posibilidades para realizar el cambio.

• Por un lado, se pueden cambiar utilizando la función names() de acuerdo con la siguiente sintaxis:

names(escena.completa) <- gsub(pattern = "LC08_L1TP_201032_20210723_20210729_02_T1",      # conjunto de caracteres que se desea sustituir
                                     replace ="",                                          # caracteres que sustituirán a los anteriores (ninguno)
                                     x <- names(escena.completa))                          # vector o dataframe que sustituirán a los anteriores
escena.completa
## class       : SpatRaster 
## dimensions  : 7911, 7791, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## sources     : LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF  
##               ... and 4 more sources
## names       : _B1, _B2, _B3, _B4, _B5, _B6, ...

• Con una sintaxis más fácil, es posible también reducir el nombre de las bandas así.

names(escena.completa) <- paste("b",
                                1:7,
                                sep="")
escena.completa
## class       : SpatRaster 
## dimensions  : 7911, 7791, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## sources     : LC08_L1TP_201032_20210723_20210729_02_T1_B1.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B2.TIF  
##               LC08_L1TP_201032_20210723_20210729_02_T1_B3.TIF  
##               ... and 4 more sources
## names       : b1, b2, b3, b4, b5, b6, ...

DIFERENCIA ENTRE EL NIVEL 1 Y EL NIVEL 2:

Como ya se ha señalado, la misma escena de Landsat puede ser descargada con un nivel de procesamiento 1 ó 2. En las fases iniciales del curso se trabajará con datos de Nivel 1, pero más adelante se utilizarán los datos de Nivel 2. En caso de elegir Landsat para realizar uno de los trabajos individuales los alumnos deberán usar siempre las escenas pertenecientes a:

Collection 2

Level 2

Tier 1

La estructura de ficheros dentro de este Nivel 2 presenta algunas diferencias con respecto a la estructura del Nivel 1, pero estos cambios no son demasiado relevantes en este curso. Sí lo es el tipo de información ofrecida.

Estructura ficheros nivel 2
Estructura ficheros nivel 2

Los productos de Nivel 1 son imágenes almacenadas en forma de Números digitales (DN) escalados (generalmente números enteros sin signo de 8 o 16 bits), que deben ser covertidos en reflectancia en el techo de la atmósfera (TOA) o radiancia usando un conjunto de metadatos proporcionados en el mismo fichero. Los productos de Nivel 2 (reflectancia superficial y temperatura de superficie) son obtenidos a partir de los productos de Nivel 1, corrigiéndose las variaciones en la radiancia medida por el sensor debido a ciertos procesos atmosféricos, con lo que más fácil la comparación de las imágenes tomadas sobre la misma área geográfica en diferentes momentos del tiempo.

Una circunstancia relevante al utilizar los datos de la Colección 2 es la necesidad de transformar los valores originales de los píxeles en valores de reflectancia reales aplicando una fórmula matemática que incluye un factor de escala y un número real.

La fórmula sería \[sr = (DN * 0.0000275) - 0.2\]

Supongamos que su píxel de la Colección 2, Nivel 2, Nivel 1 tiene un DN de 15432. Para calcular la reflectancia real de la superficie:

data_sr <- (15432 * 0.0000275) - 0.2
data_sr

Esto significa que el valor original (DN) equivale a una reflectancia de 0.22438 (o sea, el 22.438%). Hay que recordar que la reflectancia de superficie es la fracción de la radiación solar incidente que proviene de la superficie terrestre.

En el caso de la temperatura de superficie se debe recurrir al mismo procedimiento. Este producto se genera a partir de las bandas infrarrojas térmicas de Landsat Collection 2 Level 1 utilizando la reflectancia de la parte superior de la atmósfera (TOA), la temperatura de brillo (BT), el radiómetro de emisión y reflexión térmica avanzada del espacio (ASTER), la base de datos de emisividad global (GED), datos del índice de vegetación de diferencia normalizada (NDVI) de ASTER y perfiles atmosféricos de altura geopotencial, humedad específica y temperatura del aire extraídos de los datos de reanálisis.

Para transformar los ND en ºC se recurre a la siguiente transformación: \[tem ºC = (DN * 0.00341802 + 149) - 273.15\]

Por ejemplo un ND de 46432 unidades equivale a una temperatura en ºC de

lst <- (46432 * 0.00341802 + 149) - 273.15

Recuérdese también que se pueden encontrar las nomenclaturas de RT (Real Time), T1 (Tier 1) o T2 (Tier 2). Los archivos Tier 1 presentan corrección en la precisión y radiometría. Los archivos Tier 2 no presentan correcciones en la geometría debido a la imprecisión de la órbita, por lo que su nivel de precisión es menor.

Exportación e importación de imágenes en formato ráster.

Las funciones habituales de R para guardar datos (como save() o saveRDS ()) no son apropiadas para trabajar con objetos en formato ráster. Esto se debe a que, por lo general, los objetos ráster creados en R realmente no almacenan datos, sólo algunos metadatos y la ruta a los archivos originales. Esto es un “atajo” para el manejo de ficheros que “pesan” mucho y pueden causar problemas si la memoria disponible en un ordenador personal no es suficiente. Una vez que se realiza un cálculo sobre un objeto ráster, R carga los trozos o mosaicos más pequeños en la memoria, realiza los cálculos y vuelve a escribir los resultados en un archivo. Solo los rásteres pequeños se leerán completamente en la memoria. Pero este diseño tiene un inconveniente: si los archivos se cambian de directorio o se eliminan del directorio temporal, el enlace se romperá y el objeto raster quedará inutilizable.

Para guardar objetos raster en el disco duro es necesario utilizar la función terra::writeRaster(), con la que es posible elegir diferentes opciones para almacenar los datos (más información en ?writeFormats). Si se intercambian frecuentemente los ficheros con otros software SIG es recomendable el uso del formato GeoTIFF. Sin embargo, téngase en cuenta que los nombres de las bandas no se guardan bajo este formato,

  • Si se intercambian frecuentemente los ficheros con otros software SIG es recomendable el uso del formato GeoTIFF. Sin embargo, téngase en cuenta que los nombres de las bandas no se guardan bajo este formato, por lo que, en caso de archivos multicapa tipo RasterBrick se deberá recordar el orden de las bandas. Por otro lado, leer y escribir archivos GeoTiff en R es más lento que en el formato ráster .grd.

por lo que, en caso de archivos multicapa tipo RasterBrick se deberá recordar el orden de las bandas. Por otro lado, leer y escribir archivos GeoTiff en R es más lento que en el formato ráster .grd.

writeRaster(escena.completa,
            filename = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/resultados/multibanda_landsat.tif",
            filetype = "GTiff",
            overwrite = TRUE)

⚠️ Argumento = overwrite

Otra precaución es la inclusión del argumento overwrite, que sobreescribe la última versión del objeto ráster. Si deseamos sobreescribir, pero olvidamos incluir este argumento en la sintaxis, proporcionará un mensaje de error. No obstante, hay que tener cuidado, precisamente por su capacidad para reemplazar versiones anteriores del mismo raster.

Podemos volver a importar este fichero para comprobar que todo se ha realizado correctamente.

multibanda_tif <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/resultados/multibanda_landsat.tif")
multibanda_tif
## class       : SpatRaster 
## dimensions  : 7911, 7791, 7  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 316485, 550215, 4345785, 4583115  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 30N (EPSG:32630) 
## source      : multibanda_landsat.tif 
## names       :    b1,    b2,    b3,    b4,    b5,    b6, ... 
## min values  :  8486,  7594,  7003,  5998,  5207,  4755, ... 
## max values  : 65535, 57121, 65535, 65535, 65535, 65535, ...

En el caso de preferir el formato .grd la sintaxis sería la siguiente (¡no realizar inmediatamente, pues la grabación ocupa mucho tiempo!)

writeRaster(multibanda, 
            filename = "D:/G174_2025/LABORATORIO_2_Imagenes_satelite/resultados/multibanda_landsat.grd",
            overwrite = TRUE)

En el caso de importar el fichero *.grd se mantiene el nombre original de las capas, lo cual es recomendable.

multibanda_grd <- rast("D:/G174_2025/LABORATORIO_2_Imagenes_satelite/resultados/multibanda_multibanda.grd")
multibanda_grd