- Reproyectar datos Raster
3. Reproyeccion de datos raster en R
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
Introduccion
Algunas veces se encuentran datos raster que no se alinean cuando se trazan o analizan. Los datos raster que no se alinean es por que se encuentran en diferentes sistemas de coordenadas (CRS), por tanto algunas veces es necesario reproyectar los rasters utilizando la funcion projectRaster() del paquete raster
- Cargar librerias
library(raster)## Warning: package 'raster' was built under R version 3.4.4## Loading required package: splibrary(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-7library(sp)
library(ggplot2)- Importar datos Raster DEM y DEM hillshade
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")zona <- raster("Landsat_recorte_B2.TIF")
zona## class : RasterLayer
## dimensions : 2618, 3571, 9348878 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 616785, 723915, 397515, 476055 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Landsat_recorte_B2.TIF
## names : Landsat_recorte_B2
## values : 0, 35668 (min, max)- Convertir los datasets en dataframe
DEM_zona_df <- raster::as.data.frame(DEM, xy = TRUE)
str(DEM_zona_df)## 'data.frame': 1068028 obs. of 3 variables:
## $ x : num 1023028 1023119 1023211 1023302 1023393 ...
## $ y : num 979178 979178 979178 979178 979178 ...
## $ DEM_zona: int NA NA NA NA NA NA NA NA NA NA ...Hill_df <- raster::as.data.frame(Hillshade, xy = TRUE)
str(Hill_df)## 'data.frame': 1068028 obs. of 3 variables:
## $ x : num 1023028 1023119 1023211 1023302 1023393 ...
## $ y : num 979178 979178 979178 979178 979178 ...
## $ Hillshade_DEM1: int NA NA NA NA NA NA NA NA NA NA ...zona_df <- raster::as.data.frame(zona, xy = TRUE)
str(zona_df)## 'data.frame': 9348878 obs. of 3 variables:
## $ x : num 616800 616830 616860 616890 616920 ...
## $ y : num 476040 476040 476040 476040 476040 ...
## $ Landsat_recorte_B2: num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...Crear una capa del mapa DEM sobre el hillshade
ggplot() +
geom_raster(data = Hill_df,
aes(x = x, y = y,
alpha = Hillshade_DEM1)
) +
geom_raster(data = DEM_zona_df ,
aes(x = x, y = y,
fill = DEM_zona,
alpha=0.8)
) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10)) +
guides(fill = guide_colorbar()) +
scale_alpha(range = c(0.25, 0.65), guide = "none") +
coord_equal()
Ploteo DEM de la zona
ggplot() +
geom_raster(data = DEM_zona_df,
aes(x = x, y = y,
fill = DEM_zona)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10))
Ploteo hillshade de la zona
ggplot() +
geom_raster(data = Hill_df,
aes(x = x, y = y,
alpha = Hillshade_DEM1))
Tanto el Hillshade como el DEM sen encuentran en las mismas proyecciones
Ploteo de la zona de estudio
ggplot() +
geom_raster(data = zona_df,
aes(x = x, y = y,
alpha = Landsat_recorte_B2))
A pesar de que el Hillshade y el DEM se encuentran en la misma proyeccion, la imagen de la zona de estudio no esta en el mismo sistema de coordenadas (CRS) y por tanto es importante revisar el sistema de coordenadas para reprojectar alguno de los raster para que queden todas las capas alineadas.
- Visualizacion sistema de coordenadas
crs(DEM)## CRS arguments:
## +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666
## +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defscrs(Hillshade)## CRS arguments:
## +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666
## +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defscrs(zona)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0Reproyeccion de un raster
Se puede tuilizar la funcion projectRaster() para poner un raster en un nuevo CRS. La reproyeccion funciona uniamente cuando se tiene ya definido un sistema de coordenadas para el objeto raster y este se desea reproyectar. La funcion no puede ser usada si no hay un sistema de coordenadas definido.
Cuando se reproyecta un raster, es moverlo de una grilla a otra
Para utilizar la funcion projectRaster, se deben definir dos cosas: 1. el objeto que se quiere reproyectar 2. el sistema de coordenadas al cual quiero reproyectar.
La sintaxis de la funcion es projectRaster(RasterObject, crs = CRSToReprojectTo)
Para el caso particular se va amanterner el CRS de la imagen, es decir se va a reproyectar el DEM y el hillshade. La funcion projectRaster() se utiliza sobre el objeto raster y no sobre el data.frame()
DEM_zona_proj <- projectRaster(DEM,
crs = crs(zona))
Hill_proj <- projectRaster(Hillshade,
crs = crs(zona))crs(DEM_zona_proj)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0crs(Hill_proj)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0crs(zona)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0Si se compara, ahora todos los datos tienen el mismo sistema de coordenadas.
- Comparacion del contenido de uno de los rasters con su reproyeccion, donde se observa que los valores cambian al reproyectar.
extent(DEM_zona_proj)## class : Extent
## xmin : 624874.7
## xmax : 724355.9
## ymin : 396921.8
## ymax : 487930.7extent(Hill_proj)## class : Extent
## xmin : 624874.7
## xmax : 724355.9
## ymin : 396921.8
## ymax : 487930.7Resolucion del Raster
Se puede comparar la resolucion del DEM y hillshade reproyectados versus el dato original
res(DEM)## [1] 91.12718 91.12718res(Hillshade)## [1] 91.12718 91.12718res(DEM_zona_proj)## [1] 91.1 91.1res(Hill_proj)## [1] 91.1 91.1En el hillshade la resolucion es diferente, pero eso representa el mismo dato. Se le puede decir a R que obligue al raster recien reproyectado que tenga una resolucion de 1m * 1m agregando una linea de codigo res=, dentro de la funcion projectRaster()
Una vez reproyectados, es posible plotear los nuevos datos en conjunto. Para plotear con ggplot() es necesario crear un nuevo dataframe con los nuevos rasters reproyectados
DEM_proj_df <- raster::as.data.frame(DEM_zona_proj, xy = TRUE)
Hill_proj_df <- raster::as.data.frame(Hill_proj, xy = TRUE)ggplot() +
geom_raster(data = DEM_proj_df,
aes(x = x, y = y,
alpha = DEM_zona)
) +
geom_raster(data = Hill_proj_df ,
aes(x = x, y = y,
fill = Hillshade_DEM1,
alpha=0.8)
) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10))