💡 MATERIALES PARA EL EJERCICIO
► escena_landsat. Esta escena corresponde una escena del satélite Landsat 8.
Existen varias plataformas web a través de las que se puede descargar escenas de Landsat. Por ejemplo, se pueden citar las siguientes:
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.
Para descargar escenas de Landsat, debe registrarse (de forma gratuita) en el servidor https://earthexplorer.usgs.gov/.
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).
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.
Después de completar el registro, inicie con login. La aplicación pedirá su nombre de usuario (username) y contraseña (password)
EarthExplorer proporciona cuatro pestañas en el procedimiento de búsqueda para guiarlo a través de su solicitud de búsqueda:
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’.
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:
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:
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.
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:
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.
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.
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):
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:
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.
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)")
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):
Recibirá un archivo comprimido, que contiene todas las bandas espectrales como archivos geotiff georreferenciados.
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:
En la siguiente pantalla, presione Enviar pedido:
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.
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:
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:
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.
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
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
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"
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"
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
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)
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
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:
class
: el tipo de datos (en lenguaje R).
dimensions
: las dimensiones, según el número filas
(nrow), columnas (ncol) y de capas (nlyr).
resolution
: la resolución espacial, en este caso en
metros (30).
extent
: los límites (según el sistema de proyección;
dado que es UTM, serían metros) de la imagen.
`coord. ref.: sistema de coordenadas de referencia, en este caso WGS 84 / UTM zone 30N (EPSG:32630).
source
: ubicación del fichero.
names
: nombre(s) de la banda importada. Este nombre
puede ser cambiado por el investigador.
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.
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.
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,
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