Manejo de imágenes en R
1 Manejo de rasters en R
1.1 Descargar imagen
Para descargarla desde mi googleDrive vamos a utilizar los siguientes comandos
## Loading required package: sp
library(curl)
#Ubicación para guardar la imagen descargada
location <- "D:/Descargas/imagen1.tif"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=17PyQnEIICpjNtPP3v59gLEaoY3zzzEcR"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl_download(dwnld_URL,
destfile = location)
#Leer imagen
imagen1 <- raster(location)
1.2 Lectura de imágenes
Recordar que:
Cargar paquete raster
Escribir completa la ubicación del archivo
Incluir extension
Para ubicaciones de archivos “/” en lugar de “’'”
A continuación, se puede graficar como cualquier gráfica en R
Sin embargo, aquí nada más se cargó una banda. Para cargar una imagen con más de una banda hay que utilizar otra función.
También se puede utilizar la función brick
Regresamos a la imagen importada con stack
1.3 Características del objeto
Ver las características del objeto ¿Qué es cada cosa?
## class : RasterStack
## dimensions : 81, 103, 8343, 10 (nrow, ncol, ncell, nlayers)
## resolution : 0.0001796631, 0.0001796631 (x, y)
## extent : -101.2353, -101.2168, 19.64364, 19.65819 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : B2, B3, B4, B5, B6, B7, B8, B8A, B11, B12
## min values : 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535, 65535
## [1] 10
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
Más características
## [1] 0.0001796631 0.0001796631
## [1] 8343
## [1] 81 103 10
Más características
## [1] "B2" "B3" "B4" "B5" "B6" "B7" "B8" "B8A" "B11" "B12"
## class : Extent
## xmin : -101.2353
## xmax : -101.2168
## ymin : 19.64364
## ymax : 19.65819
¿Qué es esto?
## class : RasterLayer
## band : 1 (of 10 bands)
## dimensions : 81, 103, 8343 (nrow, ncol, ncell)
## resolution : 0.0001796631, 0.0001796631 (x, y)
## extent : -101.2353, -101.2168, 19.64364, 19.65819 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : D:/Descargas/imagen1.tif
## names : B2
## values : 0, 65535 (min, max)
¿y esto?
## [1] "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U" "INT2U"
## [10] "INT2U"
2 Sentinel-2
Imágenes mutiespectrales con 13 bandas de distinta resolución espacial. Por eso, en este stack sólo se tienen las imágenes
2.1 Visualización
También se puede ver la imagen con otros paquetes
## Loading required package: lattice
## Loading required package: latticeExtra
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:latticeExtra':
##
## layer
Otra forma de mostrar imagenes
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
## Warning: Raster values found that are outside the range [0, 2100]
2.2 Selección de ciertas bandas
Se puede realizar de dos maneras
3 Manejo de DEM
3.1 DEM
Primero cargar la imagen y calcular el NDVI
library(raster)
#Ubicación para guardar la imagen descargada
location_dem <- "D:/Descargas/img_DEM.tif"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=1qOJxEaAFgYMr8PQNuIKJgcTYq_SyvhO0"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl_download(dwnld_URL,
destfile = location_dem)
dem<-raster(location_dem)
plot(dem)
## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## [1] 0.0002694946 0.0002694946
## [1] 55 69 1
3.1.1 Calculo de pendiente
3.1.2 Calculo de aspecto
4 Calculo de índices
4.1 Indices de vegetación
4.1.1 Simple ratio
Otras formas de realizar lo mismo
4.1.3 EVI
EVI <- 2.5 * (imagen1$NIR - imagen1$R) / ((imagen1$NIR + 6 * imagen1$R - 7.5 * imagen1$B + 1))
plot(EVI)
No se alcanza a distinguir nada. Se ven valores de -40 a 20. Probablemente hicimos algo mal. A veces es mejor pasar de la escala 0 - 10000 a una de 0 - 1 antes de hacer los cálculos de algunos de estos índices. Entonces dividimos cada banda entre 10000.
4.2 Índices de humedad
4.2.1 NDWI
Para calcular este índice requerimos otras bandas: B8 y B11, así que volvemos a cargar la imagen.
4.3 Índices de cicatrices de incendios
4.4 Índices n-dimensionales
4.4.1 Tasseled Cap
#Tasselled Cap - Brightness
#Ahorita no se va a calcular porque no tenemos la B10
# Brightness <- 0.3037*imagen1[[1]]+0.2793*imagen1[[2]] +0.4743*imagen1[[3]]+0.5585*imagen1[[4]]+ 0.5082*B10+0.1863*imagen1[[10]]
#Tasselled Cap - Greeness
Greeness <- (-0.2848*imagen1[[1]])+(-0.2435*imagen1[[2]])+(-0.5436*imagen1[[3]]) +0.7243*imagen1[[8]]+0.0840*imagen1[[9]]+(-0.1800*imagen1[[10]])
#Tasselled Cap - Wetness
Wetness <- 0.1509*imagen1[[1]]+0.1973*imagen1[[2]]+0.3279*imagen1[[3]] +0.3406*imagen1[[8]]+(-0.7112*imagen1[[9]])+(-0.4572*imagen1[[10]])
plot(Greeness)
5 Calculo de texturas
5.1 GLCM
Este método se basa en la matriz de co-ocurrencias de tonos de gris (GLCM). Calcula diferentes métricas en ciertas direcciones. Es un método de “moving window”
Primero cargar la imagen y calcular el NDVI
library(raster)
library(glcm)
imagen1<-stack(location)
imagen1<-subset(imagen1, c(1:3,8))
names(imagen1)<-c("B","G","R","NIR")
NDVI <- (imagen1$NIR - imagen1$R) / (imagen1$NIR + imagen1$R)
5.1.1 Calculo de GLCM con dirección
Aquí podemos indicar el tamaño de la ventana que vamos a utilizar para calcular la textura. Para que sólo un pixel obtenga el valor de la textura de la ventana hay que utilizar un tamaño non. El parámetro shift va decir cuántos pixeles se va a mover en X y en Y. En este caso el shift indica que se va a comparar cada pixel dentro de la ventana con el pixel que esté a una distancia de un pixel a la derecha y uno arriba.
5.1.2 GLCM sin dirección
Debido a que se puede calcular la textura en 4 direcciones, para obtener una medida sin efecto de dirección se ponen las cuatro posibilidades: (0,1); (1,1); (1,0); (1,-1)
5.2 FOTO
Este método primero realiza una transformada de Fourier de la imagen y luego realiza una ordenación de estos datos. Por eso se llama Ordenación de la Transformada de Fourier (FOTO). La idea de este método es que caracteriza la textura a partir de la frecuencia de las ondas dominantes (r-spectrum). De tal manera, si una ventana presenta una textura dominante a una distancia de varios pixeles, estará caracterizado por ondas de frecuencia corta (pocos ciclos por km). En cambio, si una ventana presenta una textura dominante a una distancia de pocos pixeles, estará caracterizado por ondas de frecuencia largas (muchos ciclos por km). Este método se calcula por áreas de la imagen que correponden al tamaño de la ventana.
Este método sólo permite calcularse para ventanas cuadradas. Por ello, sólo se indica el lado de la ventana cuadrada.
</div>
</div>