ELEVACION SARAVENA

INSTALAMOS LOS PAQUETES NECEASARIOS

install.packages("rgdal")

install.packages("raster")

install.packages("elevatr")

install.packages("rasterVis")

install.packages("rgl")

LLAMAMOS LIBRERIAS

library(rasterVis)

library(raster)
library(rgl)
library(rgdal)

library(elevatr)

REVISAMOS EL CONTENIO DE LA CARPETA

list.files("C:/Users/sebas/OneDrive/Documentos/R/R ARTURO/ARAUCA/81_ARAUCA/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"     "MGN_DPTO_POLITICO.dbf"    
 [3] "MGN_DPTO_POLITICO.prj"     "MGN_DPTO_POLITICO.sbn"    
 [5] "MGN_DPTO_POLITICO.sbx"     "MGN_DPTO_POLITICO.shp"    
 [7] "MGN_DPTO_POLITICO.shp.xml" "MGN_DPTO_POLITICO.shx"    
 [9] "MGN_MPIO_POLITICO.cpg"     "MGN_MPIO_POLITICO.dbf"    
[11] "MGN_MPIO_POLITICO.prj"     "MGN_MPIO_POLITICO.sbn"    
[13] "MGN_MPIO_POLITICO.sbx"     "MGN_MPIO_POLITICO.shp"    
[15] "MGN_MPIO_POLITICO.shp.xml" "MGN_MPIO_POLITICO.shx"    

LEEMOS EL ARCHIVO DE FORMA PARA EL RASTER

(munic <- shapefile("C:/Users/sebas/OneDrive/Documentos/R/R ARTURO/ARAUCA/81_ARAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 7 
extent      : -72.36662, -69.42756, 6.036228, 7.104381  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                           MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,      Shape_Area 
min values  :         81,      81001,     ARAUCA,                                 1959,  945.13540708,      2017,     ARAUCA, 1.29412603421, 0.0772141267598 
max values  :         81,      81794,       TAME, Decreto Nal 677 de Abril 13 de  1987, 5787.94213047,      2017,     ARAUCA, 4.66801606458,  0.471626116664 

LE PEDIMOS LOS ATRIBUTOS DEL MUNICIPIO

head(munic)

SELECCIONAMOS EL MUNICIPIO DE INTERES

Saravena <- munic[munic$MPIO_CNMBR=="SARAVENA",]
plot(Saravena, main="SARAVENA", axes=TRUE)
plot(munic, add=TRUE)

invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

elevation <- get_elev_raster(Saravena, z = 9)

VERIFIQUEMOS QUE HAY DENTRO DEL OBJETO

elevation

class      : RasterLayer 
dimensions : 1550, 1550, 2402500  (nrow, ncol, ncell)
resolution : 0.00137, 0.00136  (x, y)
extent     : -72.42872, -70.30522, 5.609792, 7.717792  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -369.7594, 5313.876  (min, max)

GRAFICAMOS LA ELEVACION DEL DEPARTAMENTO

plot(elevation, main="DEM DESCARGADO (METROS)")
plot(Saravena, add=TRUE)

RECORTAR AHORA LOS DATOS CORRESPONDIENTES A SARAVENA

elev_crop = crop(elevation, Saravena)
plot(elev_crop, main="DEM RECORTADO")
plot(Saravena, add=TRUE)

VERIFIQUEMOS QUE TIENE EL NUEVO OBJETO

elev_crop

class      : RasterLayer 
dimensions : 235, 314, 73790  (nrow, ncol, ncell)
resolution : 0.00137, 0.00136  (x, y)
extent     : -72.09992, -71.66974, 6.746752, 7.066352  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 159.8875, 2746.603  (min, max)

QUEREMOS QUE AHORA REPROYECTAR EL MUNICIPIO TENIENDO EN CUENTA OTRO SISTEMA DE CORDENADAS QUE TOMA COMO REFERENCIA BOGOTA

spatialref <- "+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_defs" 

pr3 <- projectExtent(elev_crop, spatialref)
res(pr3) <- 100
rep_elev <- projectRaster(elev_crop, pr3)

LO QUE HAY ALLI

rep_elev

class      : RasterLayer 
dimensions : 356, 477, 169812  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1218525, 1266225, 1238252, 1273852  (xmin, xmax, ymin, ymax)
crs        : +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_defs 
source     : memory
names      : layer 
values     : 165.2396, 2766.869  (min, max)

REPROYECCION DEL MUNICIPIO

(rep_Saravena = spTransform(Saravena,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 1218645, 1266131, 1238210, 1273788  (xmin, xmax, ymin, ymax)
crs         : +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_defs 
variables   : 9
names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                       MPIO_CRSLC,   MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,      Shape_Area 
value       :         81,      81736,   SARAVENA, Decreto Nal 204 de Feb 3 de 1976, 945.13540708,      2017,     ARAUCA, 1.29412603421, 0.0772141267598 

plot(rep_elev, main="DEM REPROYECTADO")
plot(rep_Saravena, add=TRUE)

GUARDAMOS EL DEM

writeRaster(rep_elev, filename="C:/Users/sebas/OneDrive/Documentos/R/R ARTURO/ARAUCA/rep_Saravena_elev.tif", datatype='INT4S', overwrite=TRUE)

SACAMOS LAS ESTADISTICAS DEL DEM

hist(rep_elev)

promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion  <- cellStats(rep_elev, 'sd')

CREAMOS UN VECTOR UNIDIMENCIONAL

metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)

(df_estadisticas <- data.frame(metricas, valores))

OBTENCION DE VARIABLES GEOMORFOMETRICAS: Pemdiente, Aspecto y sombreado.

slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)

GRAFICA ELEVACION

plot(rep_elev,main="DEM PARA SARAVENA (METROS)", col=terrain.colors(25,alpha=0.7))

GRAFICAMOS LA PENDIENTE

plot(slope,main="PENDIENTE PARA SARAVENA (GRADOS)", col=topo.colors(25,alpha=0.7))

GRAFICAMOS ASPECTO

plot(aspect,main="ASPECTO PARA SARAVENA (GRADOS)", col=rainbow(25,alpha=0.7))

PALETA COMBINADA

plot(hill,
        col=grey(1:100/100), 
        legend=FALSE,         
        main="DEM PARA SARAVENA",
        axes=FALSE)         
plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.35), add=TRUE)

MAPAS DE DATOS DE ELEVACION CON SOMBREADO DE RAYOS

install.packages("rayshader")

library(rayshader)

COMVERTIR EL DEM EN UNA MATRIZ

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 477x356."

elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

AGREGAR CAPA DE AGUAL

elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_map()

TRAZADO DE RAYOS DESDE LA DIRECCION DEL SOL

elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  plot_map()

OTRA FORMA DE VISUALIZACION

install.packages("jpeg")

library(jpeg)

Función de sombreado de mapeo de entorno esférico

getv=function(i,a,s){
  ct = dim(i)[1:2]/2
  sx = values(s)/90 * ct[1]
  sy = values(s)/90 * ct[2]
  a = values(a) * 0.01745
  px = floor(ct[1] + sx * -sin(a))
  py = floor(ct[2] + sy * cos(a))
  
  
  template = brick(s,s,s)
  values(template)=NA
   
  cellr = px + py * ct[1]*2
  cellg = px + py * ct[1]*2 + (ct[1]*2*ct[2]*2)
  cellb = px + py * ct[1]*2 + 2*(ct[1]*2*ct[2]*2)
  
  template[[1]] = i[cellr]
  template[[2]] = i[cellg]
  template[[3]] = i[cellb]
  
  template = template * 256
  
  template
}

map = readJPEG("C:/Users/sebas/OneDrive/Documentos/R/R ARTURO/ARAUCA/9pvbHjN.jpg")
out = getv(map, aspect, slope)

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

plotRGB(out, main = "ELEVACION EN SARAVENA")