Cuando descargamos datos raster o vectoriales es importante conocer el sistema de referencia en el que se encuentra, en r la funcion que utilizamos para conocer esto es crs() Cuando tenemos datos que se encuentran en un sist de referencia diferente al que queremos trabajar, podemos reproyectar colocando el sist. indicado

Mediante un dem vamos a crear un hillshade, y vamos a utilizar la banda 5 de una imagen Landsat8

#librerias
library(raster)
## 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: D:/Lore/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: D:/Lore/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-5
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.4.4
#ubicacion del directorio
setwd("D:\\maestria\\percepcion remota avanzada\\ejerciciosR")
dem <- raster("dem.tif")
b5 <- raster("b5.tif")
summary(b5)
## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.38% of all cells)
##            b5
## Min.     5843
## 1st Qu. 13208
## Median  15133
## 3rd Qu. 17591
## Max.    51722
## NA's        0
#para crear un hillshade necesitamos crear primero un slope y aspect

slope <- terrain(dem, opt='slope', units='radians', filename = '')
slope
## class       : RasterLayer 
## dimensions  : 210, 249, 52290  (nrow, ncol, ncell)
## resolution  : 345.5232, 345.5232  (x, y)
## extent      : 652842.1, 738877.3, 9649416, 9721976  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : slope 
## values      : 0, 0.6802791  (min, max)
plot(slope)

aspect <- terrain(dem, opt='aspect', units='radians', filename = '')
plot(aspect)

hillshade <- hillShade(slope= slope, aspect = aspect)
plot(hillshade)

Vamos a utilizar las libreria de ggplot2 para mejorar la visualización del mapa. Vamos a colocar los dos raster en un mismo mapa. Al tener sistema de coordenadas diferente vemos que no podemos visualizar nuestro mapa porque se encuentran en diferentes pocisiones

# Lo primero es tranformar el RasterLayer a un data.frame
df_hillshade <- as.data.frame(hillshade, xy= TRUE)
df_hillshade <- na.omit(df_hillshade)
df_hillshade
str(df_hillshade)
## 'data.frame':    25441 obs. of  3 variables:
##  $ x    : num  678929 678584 678929 679275 678238 ...
##  $ y    : num  9721458 9721112 9721112 9721112 9720767 ...
##  $ layer: num  0.717 0.716 0.726 0.733 0.718 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:26849] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:26849] "1" "2" "3" "4" ...
df_b5 <- as.data.frame(b5, xy= TRUE)
df_b5 <- na.omit(df_b5)
str(df_b5)
## 'data.frame':    3544468 obs. of  3 variables:
##  $ x : num  677790 677820 677850 677880 677910 ...
##  $ y : num  -276060 -276060 -276060 -276060 -276060 ...
##  $ b5: int  33786 35432 35774 35978 36287 35801 31765 31431 33083 34475 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:3678199] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:3678199] "1" "2" "3" "4" ...
ggplot() +
  geom_raster(data = df_hillshade, 
              aes(x = x, y = y, 
                  alpha = layer)
  ) +
  geom_raster(data = df_b5 , 
              aes(x = x, y = y, 
                  fill = b5,
                  alpha=0.8)
  ) + 
  scale_fill_gradientn(name = "combinacion", colors = terrain.colors(10)) +
  guides(fill = guide_colorbar()) +
  scale_alpha(range = c(0.25, 0.65), guide = "none") +
  coord_equal()

Visulizamos los mapas individualmente

par(mfrow = c(2, 1))
ggplot() +
  geom_raster(data = df_b5,
              aes(x = x, y = y,
                  fill = b5)) +
  scale_fill_gradientn(name = "niveles digitales", colors = terrain.colors(10))

ggplot() +
  geom_raster(data = df_hillshade,
              aes(x = x, y = y,
                  fill = layer)) +
  scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10))

Vamos a reproyectar el hillshade con el sistema de referncia que tiene el raster b5. Para eso vamos primero vamos a conocer el sistema de referencia de cada uno. y utilizar el de uno de ellos para reproyectar el otro

#hillshade =  17s
#b5= 17n
crs(hillshade)
## CRS arguments:
##  +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0
crs(b5)
## CRS arguments:
##  +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
hillshade_proj <- projectRaster(hillshade, crs = crs(b5))
hillshade_proj
## class       : RasterLayer 
## dimensions  : 220, 259, 56980  (nrow, ncol, ncell)
## resolution  : 346, 346  (x, y)
## extent      : 651112.1, 740726.1, -352414, -276294  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.1357706, 0.9789219  (min, max)

como vemos los dos rasters se encuentran en el mismo sistemas de coordenadas

df_hillshade_proj <- as.data.frame(hillshade_proj, xy= TRUE)
df_hillshade <- na.omit(df_hillshade)



ggplot() +
  geom_raster(data = df_hillshade_proj, 
              aes(x = x, y = y, 
                  alpha = layer)
  ) +
  geom_raster(data = df_b5 , 
              aes(x = x, y = y, 
                  fill = b5,
                  alpha=0.8)
  ) + 
  scale_fill_gradientn(name = "ND", colors = terrain.colors(10)) +
  guides(fill = guide_colorbar()) +
  scale_alpha(range = c(0.25, 0.65), guide = "none") +
  coord_equal()