EXPLORACION DE IMAGENES SATELITALES
Zuly Olaya
10/3/2021
La imagen satelital “Landsat 8 "LC08_L1TP_022031_20200831_20200906_01_T1.TIF”, es la zona sur del lago Michigan, cubre parcialmente los estados de Illionis, Indiana y Michigan.
A continuación se describe la nomenclatura del nombre de la imagen satelital.
- La composición del nombre: LC08_L1TP_022031_20200831_20200906_01_T1.TIF
- L: Landsat
- C: Sensor de captura de imagen: C(OLI/TIRS), entonces es OLI / TIRS combinado con Landsat 8
- 08: Generación del satélite: Landsat 8
- L1TP: Nivel de procesado de la imagen según su corrección geométrica
- 022031: Path: 022 y Row:031
- 20200831: Fecha de adquisición de la imagen: 31 de agosto de 2020
- 20200906: Fecha de procesado de la imagen: 6 de septiembre de 2020
- 01: Numero de colección
- T1 Categoría de colección T1(Tier 1), presentan corrección en la precisión y radiometría
- Tiff: Formato de la imagen
Se cargan las librerias necesarias para reaizar los diferentes procesos.
Sys.setenv("PATH"=sprintf("%s;C:\\Users\\me\\AppData\\Local\\Programs\\MiKTeX 2.9\\miktex\\bin\\x64\\",Sys.getenv("PATH")))
library(raster)## Loading required package: splibrary(rgdal)## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/Antonio/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Antonio/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Users/Antonio/Documents/R/win-library/4.0/rgdal/projlibrary(ggplot2)
library(terra)## terra version 1.0.10##
## Attaching package: 'terra'## The following object is masked from 'package:rgdal':
##
## projectlibrary(rgeos)## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-5
## Polygon checking: TRUElibrary(sp)
library(corrplot)## corrplot 0.84 loadedSys.which("pdflatex")## pdflatex
## "C:\\PROGRA~2\\MIKTEX~1.9\\miktex\\bin\\pdflatex.exe"Se define la ruta donde se encuentra la imagen satelital
Se visualiza la ruta anteriormente escrita
getwd()## [1] "C:/Users/Antonio/Desktop/PERCEPCION_REMOTA/TALLER/LABORATORIO_2"PROPIEDADES DE LA IMAGEN
Se descarga cada banda y se guarda en una variable SpatRaster.
# Ultra azul: Estudios costeros y de aerosoles
b1 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF')
# Azul:Cartografía batimétrica, que distingue el suelo de la vegetación y la vegetación caducifolia de la vegetación de coníferas.
b2 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF')
# Verde: Destaca los picos de máxima vegetación, que son útiles para evaluar el vigor de las plantas.
b3 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF')
# Rojo: Distingue las laderas de vegetación
b4 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B4.TIF')
# Infrarrojo Cercano (NIR): Destaca el contenido de biomasa y las costas
b5 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B5.TIF')
# Infrarrojo de Onda Corta 1 (SWIR 1):Distingue la humedad del suelo y de la vegetación; penetra a través de nubes finas
b6 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B6.TIF')
# Infrarrojo de Onda Corta 2 (SWIR 2):Mejora de la lectura de la humedad del suelo y la vegetación y la penetración a través de nubes finas
b7 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B7.TIF')
# Pancromatico: Resolución de 15 metros, definición de imagen más nítida
b8 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B8.TIF')
# Cirrus: Mejor detección de la contaminación en cirros
b9 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B9.TIF')
# Sensor Térmico Infrarrojo 1 (TIRS 1):Resolución de 100 metros, mapeo térmico y humedad estimada del suelo
b10 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B10.TIF')
# Sensor Térmico Infrarrojo 2 (TIRS 2):Resolución de 100 metros, mapeo térmico y humedad estimada del suelo
b11 <- raster('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_BQA.TIF')Ya almacenadas las bandas individualmente se imprimen sus variables. Ejemplo banda B3 (verde):
+Clase: Raster +Dimensiones: Conformada por filas, columnas y numero de celdas +Resolución: Indica que el valor de pixel de la imagen tiene un cubrimiento de 30*30 metros en terreno +Extensión de la imagen: Los valores de pixel extremos en X (máximo y mínimo) y (máximo y mínimo) +Sistema de referencia: Proyección UTM zona 16 con Datum WGS 84 +Ruta de almacenamiento: Ruta de guardado +Nombre de la banda: Nombre original de la banda +Valores máximos y mínimos contenidos en la imagen: La imagen tiene un valor mínimo de 0 y un valor máximo de 65535.
b3## class : RasterLayer
## dimensions : 7871, 7751, 61008121 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## source : C:/Users/Antonio/Desktop/PERCEPCION_REMOTA/TALLER/LABORATORIO_2/LC08_L1TP_022031_20200831_20200906_01_T1/LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B3
## values : 0, 65535 (min, max)INFORMACIÓN ESTADISTICA DE LA IMAGEN
Se accede a varias propiedades de la banda. Ejemplo: Banda 3. Verde
# Numero de celdas
ncell(b3)## [1] 61008121# Dimensiones
dim(b3)## [1] 7871 7751 1# Resolución
res(b3)## [1] 30 30# Número de Layers en la imagen
nlayers(b3)## [1] 1# Comparación de bandas
compareRaster(b3,b2)## [1] TRUELa comparación de bandas permite identificar si tienen la misma extensión, número de filas y columnas, proyección, resolución y origen.
## histograma Banda 3,
hist(b3, maxcell=5000000, plot=TRUE)## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.## Warning in plot.window(xlim, ylim, "", ...): "maxcell" is not a graphical
## parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "maxcell" is not a graphical parameter## Warning in axis(1, ...): "maxcell" is not a graphical parameter## Warning in axis(2, ...): "maxcell" is not a graphical parameter
Se crea un RasterStack, un objeto con varias capas las cuales pueden ser de archivos separados a partir de los objetos RasterLayer (banda única) existentes.
s <- stack(b4, b3, b2)
s## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B4, LC08_L1TP_022031_20200831_20200906_01_T1_B3, LC08_L1TP_022031_20200831_20200906_01_T1_B2
## min values : 0, 0, 0
## max values : 65535, 65535, 65535Se crea un RasterStack usando los nombres de archivo con la ayuda de paste0.
Bandas <- paste0('LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B', 1:7, ".tif")
Bandas## [1] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B1.tif"
## [2] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B2.tif"
## [3] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B3.tif"
## [4] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B4.tif"
## [5] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B5.tif"
## [6] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B6.tif"
## [7] "LC08_L1TP_022031_20200831_20200906_01_T1//LC08_L1TP_022031_20200831_20200906_01_T1_B7.tif"Se crea un RasterStack con los elementos de la lista y se comprueban sus propiedades.
landsat <- rast(Bandas)
landsat## class : SpatRaster
## dimensions : 7871, 7751, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## sources : LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF
## ... and 4 more source(s)
## names : LC08_~T1_B1, LC08_~T1_B2, LC08_~T1_B3, LC08_~T1_B4, LC08_~T1_B5, LC08_~T1_B6, ...MAPAS Y COMPUESTOS DE ÚNICA BANDA
Se muestran las capas de manera individual de un RasterStack de la imagen multiespectral. Ejemplo de las bandas Azul, Verde, Rojo y NIR
par(mfrow = c(2,2))
plot(b2, main = "Azul", col = gray(0:100 / 100))
plot(b3, main = "Verde", col = gray(0:100 / 100))
plot(b4, main = "Rojo", col = gray(0:100 / 100))
plot(b5, main = "NIR", col = gray(0:100 / 100))
Las leyendas pueden variar entre 0 y 1. Observe la diferencia en el sombreado y el rango de leyendas entre las diferentes bandas. Esto se debe a que las diferentes características de la superficie reflejan la radiación solar incidente de manera diferente. Cada capa representa la cantidad de radiación solar incidente que se refleja para un rango de longitud de onda particular. Por ejemplo, la vegetación refleja más energía en NIR que otras longitudes de onda y, por lo tanto, parece más brillante. Por el contrario, el agua absorbe la mayor parte de la energía en la longitud de onda NIR y parece oscura.
No obtenemos tanta información de estas parcelas en escala de grises; a menudo se combinan en un “compuesto” para crear tramas más interesantes. Puede obtener más información sobre los compuestos de color en la teledetección aquí y también en la sección a continuación.
Tomado de: https://rpubs.com/ials2un/image_exploration_terra
LANDSAT EN COLOR VERDADERO.
Para componer una imagen Landsat 8 en color verdadero el orden de almacenamiento de sus bandas: +Rojo =Banda 4 +Verde= Banda 3 +Azul= Banda 2 Para esto se agrupan por medio de función stack y se guarda en la variable ImagenCV. Se verifica la estructura de los datos alamcenados en ImagenCV.
ImagenCV <- stack(b4, b3, b2)
ImagenCV## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B4, LC08_L1TP_022031_20200831_20200906_01_T1_B3, LC08_L1TP_022031_20200831_20200906_01_T1_B2
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(ImagenCV, r=1, g=2, b=3, scale=65535, axes=TRUE, stretch="lin", main="Landsat 8 en Color Verdadero")
LANDSAT EN FALSO COLOR
La banda 5 NIR, al ingresar por el canal 1( Rojo ) delata la vegetación facilitando su visualización con tonalidades rojizas.
ImagenFC <- stack(b5, b4, b3)
ImagenFC## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B5, LC08_L1TP_022031_20200831_20200906_01_T1_B4, LC08_L1TP_022031_20200831_20200906_01_T1_B3
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(ImagenFC, axes=TRUE, stretch="lin", main="Landsat Falso Color")
COMBINACIONES DE BANDA
Se grafican las bandas en las que existe mayor interacción de energía con el elemento así:
# Matriz de 2*2
par(mfrow = c(2,2))
ImagenZU <- stack(b7, b6, b4)
ImagenZU## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B7, LC08_L1TP_022031_20200831_20200906_01_T1_B6, LC08_L1TP_022031_20200831_20200906_01_T1_B4
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(ImagenZU, r=1, g=2, b=3, scale=65535, axes=TRUE, stretch="lin", main="Zonas Urbanas")
Imagenveg <- stack(b6, b5, b4)
Imagenveg## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B6, LC08_L1TP_022031_20200831_20200906_01_T1_B5, LC08_L1TP_022031_20200831_20200906_01_T1_B4
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(Imagenveg, r=1, g=2, b=3, scale=65535, axes=TRUE, stretch="lin", main="Vegetación")
Imagenagua <- stack(b5, b6, b4)
Imagenagua## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B5, LC08_L1TP_022031_20200831_20200906_01_T1_B6, LC08_L1TP_022031_20200831_20200906_01_T1_B4
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(Imagenagua, r=1, g=2, b=3, scale=65535, axes=TRUE, stretch="lin", main="Agua")
Imagenagric <- stack(b6, b5, b2)
Imagenagric## class : RasterStack
## dimensions : 7871, 7751, 61008121, 3 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## names : LC08_L1TP_022031_20200831_20200906_01_T1_B6, LC08_L1TP_022031_20200831_20200906_01_T1_B5, LC08_L1TP_022031_20200831_20200906_01_T1_B2
## min values : 0, 0, 0
## max values : 65535, 65535, 65535plotRGB(Imagenagric, r=1, g=2, b=3, scale=65535, axes=TRUE, stretch="lin", main="Agricultura")
Falso color para Zonas urbanas: Landsat 8 (7,6,4) Las áreas urbanas se reflejan en tonos “magentas”, como las que se observan al lado izquierdo cerca al lago Michigan, se observan muchas zonas urbanas distribuidas por toda la imagen, los cuerpos se ven reflejados con un azul oscuro y las zonas verdes se ven reflejadas con un tono verde claro.
Vegetación: Landsat 8 (6,5,4) La vegetación se resalta con un color verde claro, el cual se encuentra distribuido en la imagen, indicando que contiene una gran cantidad de parques, zonas verdes o bosques, las zonas urbanas se ven reflejadas de color morado claro y el agua con azul oscuro.
Agua: Landsat 8 (5,6,4) Los cuerpos de agua se reflejan con color azul oscuro como se evidencia en el lago Michigan, también se ven cuerpos de agua en la zona oriente que a diferencia de las anteriores composiciones no se evidenciaba, las zonas urbanas se reflejan de color celeste y las zonas verdes de color naranja.
Agricultura: Landsat 8 (6,5,2) Las zonas agrícolas se ven representadas por una tonalidad ver brillante, en la imagen se detectan por franjas, las zonas urbanas ya no se detectan y se resaltan los cuerpos de agua con una tonalidad azul oscura.
##SUBCONJUNTOS Y RENOMBRAR BANDAS Se seleccionan capas específicas (bandas) mediante la funcion subset o mediante la indexación. En este ejemplo se muestran los dos casos con las tres primeras bandas
landsatsub1 <- subset(landsat, 1:3)
# igual
landsatsub2 <- landsat[[1:3]]
landsatsub2## class : SpatRaster
## dimensions : 7871, 7751, 3 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## sources : LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF
## names : LC08_L1TP_~6_01_T1_B1, LC08_L1TP_~6_01_T1_B2, LC08_L1TP_~6_01_T1_B3Se identifica los números de bandas almacenados en las variables, landsat, landsatsub1 y landsatsub2
nlyr(landsat)## [1] 7nlyr(landsatsub1)## [1] 3nlyr(landsatsub2)## [1] 3Se extrae las bandas que se van a usar, para esto usamos subset de la variable “landsat”:
landsat <- subset(landsat, 1:7)
landsat## class : SpatRaster
## dimensions : 7871, 7751, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## sources : LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF
## LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF
## ... and 4 more source(s)
## names : LC08_~T1_B1, LC08_~T1_B2, LC08_~T1_B3, LC08_~T1_B4, LC08_~T1_B5, LC08_~T1_B6, ...Se almacenan los nombres de las 7 bandas
names(landsat)## [1] "LC08_L1TP_022031_20200831_20200906_01_T1_B1"
## [2] "LC08_L1TP_022031_20200831_20200906_01_T1_B2"
## [3] "LC08_L1TP_022031_20200831_20200906_01_T1_B3"
## [4] "LC08_L1TP_022031_20200831_20200906_01_T1_B4"
## [5] "LC08_L1TP_022031_20200831_20200906_01_T1_B5"
## [6] "LC08_L1TP_022031_20200831_20200906_01_T1_B6"
## [7] "LC08_L1TP_022031_20200831_20200906_01_T1_B7"names(landsat) <- c('Ultra Azul', 'Azul', 'Verde', 'Rojo', 'NIR', 'SWIR1', 'SWIR2')
names(landsat)## [1] "Ultra Azul" "Azul" "Verde" "Rojo" "NIR"
## [6] "SWIR1" "SWIR2"SUB CONJUNTO ESPACIAL O RECORTE
Para realizar el recorte se utiliza el archivo vector “Illionis_utm_zona16.shp” que pertenece al estado de Illionis en Estados Unidos", con el fin de tener un área delimitada de trabajo. Para esto se carga el shapefile por medio de la ruta de almacenamiento.
knitr::opts_chunk$set(echo = TRUE)
setwd("C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2")Se visualiza la ruta del shapefile
getwd()## [1] "C:/Users/Antonio/Desktop/PERCEPCION_REMOTA/TALLER/LABORATORIO_2"shape<-("ILLLIONIS_utm_zone16.shp")
shape## [1] "ILLLIONIS_utm_zone16.shp"Se guarda el shapefile del estado de Illionis en la variable “Zona_estudio”
zona_estudio<-readOGR(".","ILLLIONIS_utm_zone16")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Antonio\Desktop\PERCEPCION_REMOTA\TALLER\LABORATORIO_2", layer: "ILLLIONIS_utm_zone16"
## with 1 features
## It has 3 fieldsSe identifica que el shapefile “Zona de estudio” y la “imagen landsat” tengan el mismo sistema de coordenadas. Se verifica la clase del shapefile en este caso es polígono
class(zona_estudio)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"La proyección
proj4string(zona_estudio)## Warning in proj4string(zona_estudio): CRS object has comment, which is lost in
## output## [1] "+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs""Si el shapefile tuviera un sistema de coordenadas diferente, se modificarian las propiedades de esta forma: En utm16 se define como zona=16 y en WGS84 = se definen las coordenadas geograficas WGS84.
###utm16 <- "+proj=utm +zone=16 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs" #utm zona 16 surcon datum WGS84
wgs84 <- "+proj=longlat +ellps=WGS84" #coordenadas geograficas WGS84
proj4string(zona_estudio) <- CRS(utm16)###Propiedades de “zona de estudio”
zona_estudio## class : SpatialPolygonsDataFrame
## features : 1
## extent : 115181.3, 455925, 4095660, 4712751 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## variables : 3
## names : STATE_NAME, STATE_FIPS, STATE_ABBR
## value : Illinois, 17, ILCon la función ext(extensión) se identifica las dimensiones de la imagen satelital.
ext(landsat)## SpatExtent : 403485, 636015, 4504485, 4740615 (xmin, xmax, ymin, ymax)Se realiza lo mismo para la “zona de estudio”
ext(zona_estudio)## SpatExtent : 115181.346239442, 455924.962844811, 4095659.50650682, 4712750.858519 (xmin, xmax, ymin, ymax)Con la función crop se realiza el corte entre la imagen “landsat” y “zona de estudio”
(landsatcrop <- crop(landsat, zona_estudio))## class : SpatRaster
## dimensions : 6942, 1748, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 403485, 455925, 4504485, 4712745 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## source : spat_mVy1JQtlnFXglQC.tif
## names : Ultra Azul, Azul, Verde, Rojo, NIR, SWIR1, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 29467, 31077, 31398, 34441, 39142, 36252, ...Se visualiza el corte
plotRGB(landsatcrop, r=4, g=3, b=2, axes=TRUE, stretch="lin", main="Recorte Landsat RGB432 color verdadero")
###Guardar los resultados en Disco
Para guardar la imagen recortada (raster) en el disco se usa la función writeRaster. Admiten varios tipos de archivos. Usaremos el formato GeoTiff de uso común. Si bien se conserva el orden de las capas, los nombres de las capas lamentablemente se pierden en el formato GeoTiff.
writeRaster(landsatcrop, filename="C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\ILLLIONIS_recorte.tiff",
datatype="INT1U", overwrite=TRUE)###Relación entre bandas
Un scatterplot matrix o diagrama de dispersión son utiles para explorar las relaciones que existen entre capas ráster. Esto se puede hacer con la función pairs(Pares) de paquete Terra. EL siguiente gráfico representa los valores digitales en la longitud de onda azul frente a los valores digitales en la longitud de onda verde.
pairs(landsatcrop[[2:3]], main = "Azul Vs Verde")
En el gráfico se evidencia altas correlaciones entre las bandas azul y verde. Si se usan las dos bandas no se perdería mucha información por su correlación.
pairs(landsatcrop[[4:5]], main = "Rojo vr NIR")
En este ejemplo se evidencia la correlación entre la banda NIR y la banda rojo, los puntos no se ven de forma alargada sino triangular, con una correlación alta-media.
Se grafican todas las bandas para poder visualizar su correlación.
pairs(landsatcrop[[1:7]], main = "Comparación de bandas")
Se observa en general que todas las bandas tienen una buena correlación, las más altas se reflejan entre las bandas ultra-azul, azul, verde y rojo y la menor correlación entre bandas se comparan con la banda NIR.
###EXTRAER VALORES DE PIXEL
Para extraer valores de pixel de las celdas raster para ubicaciones geográficas específicas.
La función de extracción se utiliza para obtener valores ráster en las ubicaciones de otros datos espaciales. Se utilizara como insumo las coberturas de la tierra de “Global Land Cover 2000”que se pueden descargar de la siguiente página: https://forobs.jrc.ec.europa.eu/products/glc2000/products.php.
(landcover <- shapefile("C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_2\\illionis_coberturas_.shp"))## class : SpatialPolygonsDataFrame
## features : 107
## extent : 406408.1, 455802.1, 4547467, 4667565 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs
## variables : 6
## names : OBJECTID_1, OBJECTID, gridcode, USER_CLASS, Shape_Leng, Shape_Area
## min values : 1, 1, 11, Consolidated Rock Sparse Vegetation, 3289.23631651, 504451.735068
## max values : 99, 99, 9, Water bodies, 558670.755792, 2069140874.67Se crean automaticamente 300 puntos aleatorios dentro del shape
ptsamp <- spsample(landcover, 300, type = "random")## Warning in proj4string(obj): CRS object has comment, which is lost in outputSe agrega la calse de coberturas a los puntos creados
ptsamp$Clase <- over(ptsamp, landcover)$USER_CLASSptsamp <- vect(ptsamp)
# Usamos las coordenadas x-y para extraer los valores espectrales de las ubicaciones.
xy <- as.matrix(geom(ptsamp)[,c('x','y')])
df <- extract(landsat, xy)
# Para ver algunos de los valores digitales
head(df)## Ultra Azul Azul Verde Rojo NIR SWIR1 SWIR2
## 1 9660 8736 8199 6889 24952 12193 7928
## 2 11303 10665 9837 9949 11897 10201 9131
## 3 9959 9064 9257 7915 22521 12798 8631
## 4 9773 8837 8230 7311 17537 12004 8562
## 5 12319 11827 11237 11285 13544 12738 10917
## 6 9678 8697 8153 6768 27540 12831 8043Perfiles espectrales
Un gráfico del espectro (todas las bandas) para píxeles que representan ciertas características de la superficie terrestre (por ejemplo, agua) se conoce como perfil espectral. Dichos perfiles demuestran las diferencias en las propiedades espectrales de varias características de la superficie terrestre y constituyen la base para el análisis de imágenes. Los valores espectrales se pueden extraer de cualquier conjunto de datos multiespectrales utilizando la función de extracción. En el ejemplo anterior, extrajimos valores de datos Landsat para las muestras. Estas muestras incluyen: bosques, arbustos, pastizales, tierras de cultivo, agua, barbecho, construidos y abiertos. Primero calculamos los valores medios de reflectancia para cada clase y cada banda.
ms <- aggregate(df, list(ptsamp$Clase), mean)
# en lugar de la primera columna, usamos nombres de fila
rownames(ms) <- ms[,1]
ms <- ms[,-1]
ms <- ms * 2E-05
head(ms) ## Ultra Azul
## Consolidated Rock Sparse Vegetation 0.2190600
## Cropland 0.2024378
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.2144800
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.1966900
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.2244000
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2529200
## Azul
## Consolidated Rock Sparse Vegetation 0.2036200
## Cropland 0.1852471
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.1993733
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.1776625
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.2108858
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2443800
## Verde
## Consolidated Rock Sparse Vegetation 0.1928100
## Cropland 0.1747399
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.1901667
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.1665100
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.1995033
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2382000
## Rojo
## Consolidated Rock Sparse Vegetation 0.1855900
## Cropland 0.1586798
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.1856400
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.1448875
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.1976775
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2450000
## NIR
## Consolidated Rock Sparse Vegetation 0.3044900
## Cropland 0.3902626
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.3485200
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.4160675
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.3036233
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2740600
## SWIR1
## Consolidated Rock Sparse Vegetation 0.2505900
## Cropland 0.2568075
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.2929867
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.2375600
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.2660967
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2808400
## SWIR2
## Consolidated Rock Sparse Vegetation 0.2002200
## Cropland 0.1800748
## Temperate or Sub-polar Broadleaved Deciduous Forest - Closed Canopy 0.2340667
## Temperate or Sub-polar Mixed Broadleaved or Needleleaved Forest - Closed Canopy 0.1594075
## Temperate or Sub-polar Needleleaved Evergreen Forest - Closed Canopy 0.2073033
## Temperate or Subpolar Needleleaved Evergreen Shrubland - Open Canopy 0.2418400# Se crea un vector de color para las clases de cobertura terrestre para usarlo en el trazado
mycolor <- c("#B2182B", "#D6604D", "#F4A582", "#FDDBC7", "#D1E5F0", "#92C5DE", "#4393C3" ,"#2166AC", "#9ea832", "#32a873", "#3292a8", "#5532a8", "#9e32a8",
"#a8325d", "#0ce2ed", "#32a8a0", "#3267a8", "#324ea8")
# Se transforma ms de un data.frame a una matriz
ms <- as.matrix(ms)
# Primero se crea una celda vacía
plot(0, ylim=c(0,0.8), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Reflectancia")
# add the different classes
for (i in 1:nrow(ms)){
lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Titulo
title(main="Perfiles espectrales para Landsat", font.main = 2)
# Legenda
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
El perfil espectral la muestra similitud en la reflectancia de las diferentes características que hay en la superficie de la tierra (o por encima de ella). Agua en el grafico como (Water bodies o cuerpos de agua) muestra una reflexión relativamente alta en todas las bandas.
“Bosque templado o subpolar mixto de hojas anchas o de hojas acuosas: dosel cerrado”, tienen la reflectancia mas alta entre las bandas 4 y 6.
La “Vegetación escasa de rocas consolidadas” es la que menos reflectancia en todas las bandas
##Este archivo esta basado en el la publicación del profesor Ivan lizarazo que se puede cnsultar en: https://rpubs.com/ials2un/image_exploration_terra
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252 LC_CTYPE=Spanish_Spain.1252
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] corrplot_0.84 rgeos_0.5-5 terra_1.0-10 ggplot2_3.3.3 rgdal_1.5-23
## [6] raster_3.4-5 sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 highr_0.8 compiler_4.0.4 pillar_1.5.1
## [5] tools_4.0.4 digest_0.6.27 evaluate_0.14 lifecycle_1.0.0
## [9] tibble_3.1.0 gtable_0.3.0 lattice_0.20-41 pkgconfig_2.0.3
## [13] rlang_0.4.10 DBI_1.1.1 yaml_2.2.1 xfun_0.21
## [17] withr_2.4.1 stringr_1.4.0 dplyr_1.0.5 knitr_1.31
## [21] generics_0.1.0 vctrs_0.3.6 grid_4.0.4 tidyselect_1.1.0
## [25] glue_1.4.2 R6_2.5.0 fansi_0.4.2 rmarkdown_2.7
## [29] purrr_0.3.4 magrittr_2.0.1 scales_1.1.1 codetools_0.2-18
## [33] htmltools_0.5.1.1 ellipsis_0.3.1 assertthat_0.2.1 colorspace_2.0-0
## [37] utf8_1.1.4 stringi_1.5.3 munsell_0.5.0 crayon_1.4.1