Elevación Caquetá

Objetivo

El objetivo de este trabajo es obtener, procesar y visuaizar el modelo digital de la elevación del departamento de Caquetá - Colombia.

Introducción a elevatr

Los datos de elevación juegan un papel importante aplicaciones como hidrografía, modelado geográfico, visualización, entre otros. En R, el paquete elevatr es útil ya que permite el acceso a los datos de elevación desde las API web.

library(rasterVis)

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

library(elevatr)

Obtener datos de elevación raster.

Usando get_elev_raster () para acceder a los Terrain Tiles en AWS.

list.files("/cloud/project/Caquetá/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"    
 [2] "MGN_DPTO_POLITICO.dbf"    
 [3] "MGN_DPTO_POLITICO.prj"    
 [4] "MGN_DPTO_POLITICO.sbn"    
 [5] "MGN_DPTO_POLITICO.sbx"    
 [6] "MGN_DPTO_POLITICO.shp"    
 [7] "MGN_DPTO_POLITICO.shp.xml"
 [8] "MGN_DPTO_POLITICO.shx"    
 [9] "MGN_MPIO_POLITICO.cpg"    
[10] "MGN_MPIO_POLITICO.dbf"    
[11] "MGN_MPIO_POLITICO.prj"    
[12] "MGN_MPIO_POLITICO.sbn"    
[13] "MGN_MPIO_POLITICO.sbx"    
[14] "MGN_MPIO_POLITICO.shp"    
[15] "MGN_MPIO_POLITICO.shp.xml"
[16] "MGN_MPIO_POLITICO.shx"    

(munic <-  shapefile("/cloud/project/Caquetá/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 16 
extent      : -76.30622, -71.25385, -0.70584, 2.964148  (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  :         18,      18001,    ALBANIA, Decreto 1678 de Septiembre 7 de 1967,  403.35827433,      2017,    CAQUETÁ,  1.1128288978, 0.0327393194813 
max values  :         18,      18860, VALPARAÍSO,  Ordenanza 3 de Noviembre 12 de 1985, 42312.7997376,      2017,    CAQUETÁ, 20.0458589874,   3.43579026269 

¿Cuáles son los atributos del objeto munich?

head(munic)

Seleccionemos solo la ciudad capital de este departamento.

medallo <- munic[munic$MPIO_CNMBR=="FLORENCIA",]
plot(medallo, main="Florencia", axes=TRUE)
plot(munic, add=TRUE)

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

elevation <- get_elev_raster(medallo, z = 8)

Downloading DEMs [========>------------------]  33% eta:  1s
Downloading DEMs [=============>-------------]  50% eta:  1s
Downloading DEMs [=================>---------]  67% eta:  0s
Downloading DEMs [=====================>-----]  83% eta:  0s
Downloading DEMs [===========================] 100% eta:  0s
Merging DEMs
Reprojecting DEM to original projection
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfNote: Elevation units are in meters.
Note: The coordinate reference system is:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

¿Qué hay dentro del objeto de elevación?

elevation

class      : RasterLayer 
dimensions : 1544, 1033, 1594952  (nrow, ncol, ncell)
resolution : 0.00275, 0.00275  (x, y)
extent     : -75.95125, -73.1105, -1.420879, 2.825121  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -611.0625, 3796.604  (min, max)

####Recorte los datos de elevación para que coincidan con la extensión del área de estudio.

plot(elevation, main="This the downloaded DEM [meters]")
plot(medallo, add=TRUE)

writeRaster(elevation, filename="./dem/elev_z8.tif", datatype='INT4S', overwrite=TRUE)

elev_crop = crop(elevation, medallo)
plot(elev_crop, main="Cropped Digital elevation model")
plot(medallo, add=TRUE)

elev_crop

class      : RasterLayer 
dimensions : 318, 169, 53742  (nrow, ncol, ncell)
resolution : 0.00275, 0.00275  (x, y)
extent     : -75.8055, -75.34075, 1.320871, 2.195371  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 216.2132, 3375.029  (min, max)

Reproyectar los datos de elevación

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)

rep_elev

class      : RasterLayer 
dimensions : 968, 518, 501424  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 807662.2, 859462.2, 637831.4, 734631.4  (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     : 215.1738, 3368.633  (min, max)

(rep_medallo = spTransform(medallo,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 807669.7, 859391.2, 637973.3, 734552.1  (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       :         18,      18001,  FLORENCIA, Decreto 642 de Junio 17 de 1912, 2547.63842143,      2017,    CAQUETÁ, 2.94250780503, 0.206927769189 

plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_medallo, add=TRUE)

Estadísticas básicas de datos de elevación

hist(rep_elev)

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

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

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

Obtención de variables geomorfométricas.

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

plot(rep_elev,main="DEM for Caquetá [meters]", col=terrain.colors(25,alpha=0.7))

plot(slope,main="Slope for Caquetá [degrees]", col=topo.colors(25,alpha=0.7))

plot(aspect,main="Aspect for Caquetá [degrees]", col=rainbow(25,alpha=0.7))

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

Mapeo de datos de elevación con sombreador de rayos

install.packages("rayshader")

library(rayshader)

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 518x968."

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

Calculating Surface Normal [---------------] ETA:  0s
Calculating Surface Normal [==============-] ETA:  0s
Calculating Surface Normal [===============] ETA:  0s
                                                     

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

Calculating Surface Normal [---------------] ETA:  0s
Calculating Surface Normal [===============] ETA:  0s
                                                     

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

Calculating Surface Normal [---------------] ETA:  0s
Calculating Surface Normal [===============] ETA:  0s
                                                     

Raytracing [-------------------------------] ETA:  0s
Raytracing [====---------------------------] ETA:  1s
Raytracing [=====--------------------------] ETA:  1s
Raytracing [======-------------------------] ETA:  1s
Raytracing [=======------------------------] ETA:  1s
Raytracing [========-----------------------] ETA:  1s
Raytracing [=========----------------------] ETA:  1s
Raytracing [==========---------------------] ETA:  1s
Raytracing [===========--------------------] ETA:  1s
Raytracing [============-------------------] ETA:  1s
Raytracing [=============------------------] ETA:  1s
Raytracing [==============-----------------] ETA:  1s
Raytracing [===============----------------] ETA:  1s
Raytracing [================---------------] ETA:  1s
Raytracing [=================--------------] ETA:  1s
Raytracing [==================-------------] ETA:  1s
Raytracing [===================------------] ETA:  1s
Raytracing [====================-----------] ETA:  1s
Raytracing [=====================----------] ETA:  1s
Raytracing [=====================----------] ETA:  0s
Raytracing [======================---------] ETA:  0s
Raytracing [=======================--------] ETA:  0s
Raytracing [========================-------] ETA:  0s
Raytracing [=========================------] ETA:  0s
Raytracing [==========================-----] ETA:  0s
Raytracing [===========================----] ETA:  0s
Raytracing [============================---] ETA:  0s
Raytracing [=============================--] ETA:  0s
Raytracing [==============================-] ETA:  0s
Raytracing [===============================] ETA:  0s
                                                     

Otra forma de visualización

install.packages("jpeg")

library(jpeg)

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("./9pvbHjN.jpg")

out = getv(map, aspect, slope)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

plotRGB(out, main = "Supposedly pretty mountains in Caquetá")