Este codigo fue tomado de los tutoriales publicados por Data Carpentry, realizado por Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul, Joseph Stachelek como parate del curso Introduccion a los Datos Geoespaciales con R.
A menudo se necesitan combinar valores y realizar calculos entre rasters para crear nuevos rasters de salida. Esto incluye como restar un raster de otro usando matematica basica entre rasters y la funcion overlay(). Tambien incluye como extraerlos valores de pixeles de un conjunto de ubicaciones, por ejemplo el buffer de una region alrededor alrededor de la ubicacion delas parcelas de un sitio en campo.
library(raster)
## Warning: package 'raster' was built under R version 3.4.4
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.2-16, (SVN revision 701)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
## Path to GDAL shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/Natalia/Documents/R/win-library/3.4/rgdal/proj
## Linking to sp version: 1.2-7
library(sp)
library(ggplot2)
DEM <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\DEM_zona.tif")
# import DTM hillshade
Hillshade <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Hillshade_DEM1.tif")
# Imagen B2
Imagen <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Landsat_recorte_B2.tif")
DEM_zona_df <- raster::as.data.frame(DEM, xy = TRUE)
DEM_zona_df <- na.omit(DEM_zona_df)
str(DEM_zona_df)
## 'data.frame': 553959 obs. of 3 variables:
## $ x : num 1051460 1051460 1051369 1051460 1051369 ...
## $ y : num 979087 978996 978905 978905 978814 ...
## $ DEM_zona: int 2956 2839 2706 2690 2638 2633 2604 2560 2560 2540 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:514069] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..- attr(*, "names")= chr [1:514069] "1" "2" "3" "4" ...
Hill_df <- raster::as.data.frame(Hillshade, xy = TRUE)
Hill_df <- na.omit(Hill_df)
str(Hill_df)
## 'data.frame': 549211 obs. of 3 variables:
## $ x : num 1051460 1051369 1051460 1051278 1051369 ...
## $ y : num 978723 978632 978632 978541 978541 ...
## $ Hillshade_DEM1: int 7 98 48 111 112 56 67 63 12 6 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:518817] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..- attr(*, "names")= chr [1:518817] "1" "2" "3" "4" ...
Imagen_df <- raster::as.data.frame(Imagen, xy = TRUE)
Imagen_df <- na.omit(Imagen_df)
str(Imagen_df)
## 'data.frame': 5012512 obs. of 3 variables:
## $ x : num 635820 635850 635880 635910 635940 ...
## $ y : num 476040 476040 476040 476040 476040 ...
## $ Landsat_recorte_B2: num 7796 7920 7951 8034 8240 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:4336366] 1 2 3 4 5 6 7 8 9 10 ...
## .. ..- attr(*, "names")= chr [1:4336366] "1" "2" "3" "4" ...
Los calculos entre dos o mas rasters se realiza para crear una nueva salida. Por ejemplo si se esta interesado en mapear la altura de los arboles en toda una zona, se puede calcular la diferencia entre el modelo digital de superficie (DSM) y el modelo digital del terreno (DTM). como resultado el dataset se refiere al modelo de altura del canopi (CHM) y representa la altura actual de los arboles, edificios, etc, sin la influencia de la elevacion del suelo.
ggplot() +
geom_raster(data = DEM_zona_df,
aes(x = x, y = y,
fill = DEM_zona)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10))
ggplot() + geom_raster(data = DEM_hill_df, aes(x = x, y = y, alpha = Hillshade_DEM1)) + scale_alpha(range = c(0.15, 0.65), guide = “none”)
ggplot() +
geom_raster(data = Hill_df, aes(x = x, y = y, alpha = Hillshade_DEM1)) +
scale_alpha(range = c(0.15, 0.65), guide = "none")
Se puede calcular la diferencia entre dos rasters en dos formas diferentes:
Substrayendo los dos rasters en R utilizando raster math,
Ppara un proceso mas eficiente, particularmente si los rasters son grandes y el calculo se vuelve mas complejo, se utiliza la funcion overlay()
Es posible hacer calculos raster por sustracción, adicion, multiplicacion, etc entre dos rasters. En el mundo geoespacial, se le da el nombre de “raster math”
Resta del Modelo Digital del Terreno de el modelo digital de elevacion para crear el modelo de altura del canopy. Despues de la sustraccion, se crea un dataframe para plotear con ggplot.
CHM <- DEM - Hillshade
CHM_df <- raster::as.data.frame(CHM, xy = TRUE)
CHM_df <- na.omit(CHM_df)
Plot del modelo de altura del canopy (CHM)
ggplot() +
geom_raster(data = CHM_df ,
aes(x = x, y = y,
fill = layer)
) +
scale_fill_gradientn(name = "Altura Canopy", colors = terrain.colors(10))
Distribucion de los valores de la nueva capa creada de altura del canopy
ggplot(CHM_df) +
geom_histogram(aes(layer))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
El rango esta entre 0 y 2600
Raster math, es una aproximacion apropiada para el calculo raster si:
Sin embargo, raster math es poco eficiente en calculos mas complejos o archivos muy grandes. La funcion overlay() es mas eficiente cuando:
La funcion overlay() toma dos o mas rasters y aplica una funcion utilizando metodos de procesamiento eficientes. La sintaxis es:
outputRaster <- overlay(raster1, raster2, fun=functionName)
Si los rasters estan almacenados como RasterStack o RasterBrick en R se debe usar calc(). Ya que overlay() no funciona con raster stack
Hora se va a realizar el mismo calculo de resta en que se utilizó raster math, pero ahora utilizando la funcion overlay().
CHM_ov <- overlay(DEM, Hillshade,
fun = function(r1, r2) { return( r1 - r2) })
plot(CHM_ov)
Ahora se necesita convertir el nuevo objeto en un data frame para plotear con ggplot
CHM_ov_df <- raster::as.data.frame(CHM_ov, xy = TRUE)
CHM_ov_df <- na.omit(CHM_ov_df)
ggplot() +
geom_raster(data = CHM_ov_df,
aes(x = x, y = y,
fill = layer)
) +
scale_fill_gradientn(name = "Altura canopy", colors = terrain.colors(10))
Para exportar el raster creado como un GeoTIFF se utiliza la funcion writeRaster().
Cuando se escribe un objeto raster a un archivo GeoTIFF, el nombre queda nombre-tiff. La funcion writeRaster() escribe por defecto el archivo de salida en su directorio de trabajo a menos que se especifique una ruta de archivo diferente.
Se especifica el formato de salida (“GTiff”), los Na (NAflag =-9999). Tambien se debe especificar a R sobre reescribir sobre un archivo con el mismo nombre (overwrite =TRUE)
writeRaster(CHM_ov, "CHM.tiff",
format="GTiff",
overwrite=TRUE,
NAflag=-9999)