Slope and aspect calculation

Cuando se trabaja con datos raster del tipo DEM es util poder calcular su pendiente y su aspecto para ello relaizmaos el siguiente proceso. Pirmero cargamos la libreria

library(raster)
dem <- raster(ncol=3, nrow=3, xmn=100, xmx=115, ymn=100, ymx=115)
ncell(dem)

[1] 9

#resolución de la celda
res(dem)

[1] 5 5

valores <- c(50, 45, 50, 30, 30, 30, 8, 10, 10)
(values(dem) <- valores)

[1] 50 45 50 30 30 30  8 10 10

plot(dem, main = "DEM")
text(dem)

#el codigo epsg 3115 es de MAGMA-SIRGAS/colombia west zone
crs(dem) <- CRS('+init=epsg:3115')

Slope and aspect calculation

(slope = terrain(dem, 'slope', unit='degrees', neighbors=8))

class      : RasterLayer 
dimensions : 3, 3, 9  (nrow, ncol, ncell)
resolution : 5, 5  (x, y)
extent     : 100, 115, 100, 115  (xmin, xmax, ymin, ymax)
crs        : +init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : slope 
values     : 75.25766, 75.25766  (min, max)

Slope calculation

plot(slope, main = "Pendiente")
text(slope)

(aspecto = terrain(dem, 'aspect', unit='degrees', neighbors=8))

class      : RasterLayer 
dimensions : 3, 3, 9  (nrow, ncol, ncell)
resolution : 5, 5  (x, y)
extent     : 100, 115, 100, 115  (xmin, xmax, ymin, ymax)
crs        : +init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : aspect 
values     : 180.7538, 180.7538  (min, max)

Aspect calculation

plot(aspecto, main = "Aspecto")
text(aspecto)

##Assigment Ahora utilicemos más valores para DEM y veremos como la pendiente y el aspecto cambian, esta vez utilizaremos un raster con 4 filas y 4 columnas

dem1 <- raster(ncol=4, nrow=4, xmn=100, xmx=120, ymn=100, ymx=120)
ncell(dem1)

[1] 16

res(dem1)

[1] 5 5

valores <- c(48, 50, 45, 50, 35, 30, 30, 30, 8, 8, 10, 10, 25, 15, 20, 20)
(values(dem1) <- valores)

 [1] 48 50 45 50 35 30 30 30  8  8 10 10 25 15 20 20

plot(dem1, main = "DEM")
text(dem1)

#el codigo epsg 3115 es de MAGMA-SIRGAS/colombia west zone
crs(dem1) <- CRS('+init=epsg:3115')

(slope = terrain(dem1, 'slope', unit='degrees', neighbors=8))

class      : RasterLayer 
dimensions : 4, 4, 16  (nrow, ncol, ncell)
resolution : 5, 5  (x, y)
extent     : 100, 120, 100, 120  (xmin, xmax, ymin, ymax)
crs        : +init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : slope 
values     : 48.92371, 75.91133  (min, max)

plot(slope, main = "Pendiente")
text(slope)

(aspecto = terrain(dem1, 'aspect', unit='degrees', neighbors=8))

class      : RasterLayer 
dimensions : 4, 4, 16  (nrow, ncol, ncell)
resolution : 5, 5  (x, y)
extent     : 100, 120, 100, 120  (xmin, xmax, ymin, ymax)
crs        : +init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666 +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
source     : memory
names      : aspect 
values     : 173.1572, 191.3099  (min, max)

plot(aspecto, main = "Aspecto")
text(aspecto)