library(raster)
## Loading required package: sp
DEM <- raster(ncol=3, nrow=3, xmn=100, xmx=115, ymn=100, ymx=115)
ncell(DEM)
## [1] 9
#esta linea muestra la resolucion de las celdas en sus ejes X y Y
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)
# con esta instruccion, cambiamos el sistema de coordenadas de la matriz, para que use el sistema MAGNA para la zona de Bogota
crs(DEM) <- CRS('+init=epsg:3115')
(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)
plot(slope, main = "Pendiente")
text(slope)
# al igual que con la pendiente, realizamos el mismo procedimiento para el aspecto
(aspect=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)
plot(aspect, main = "Orientacion")
text(aspect)
# ahora realizaremos todo el procedimiento anterior, pero con una matriz 4x4 # empezamos primero por definir en un objeto las caracteristicas de la matriz
Ejemplo <- raster(ncol=4, nrow=4, xmn=100, xmx=120, ymn=100, ymx=120)
res(Ejemplo)
## [1] 5 5
Ejm_val<- c(50,45,50,48,30,29,30,29,10,9,9,10,25,23,19,21)
(values(Ejemplo)<-Ejm_val)
## [1] 50 45 50 48 30 29 30 29 10 9 9 10 25 23 19 21
plot(Ejemplo, main='ejemplo')
text(Ejemplo)
#cantidad de celdas en la matriz
ncell(Ejemplo)
## [1] 16
crs(Ejemplo) <- CRS('+init=epsg:3115')
(slope_Ejm = terrain(Ejemplo, 'slope', unit='degrees', neighbors=4))
## 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 : 31.31116, 76.29304 (min, max)
#imagen de la pendiente de las celdas mencionadas
plot(slope_Ejm, main='Pendiente de Ejemplo')
text(slope_Ejm)
# se repite el procedimiento anterior, pro en vez de la pendiente, ahora se calculara el aspecto de la misma
(aspect_Ejm=terrain(Ejemplo,'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 : 164.0546, 181.4688 (min, max)
plot(aspect_Ejm, main = "Orientacion de Ejemplo")
text(aspect_Ejm)