Datos de elevación en R

Paula Virguez Gómez

28 de marzo de 2020

1.Introducción

Instalar paquetes

install.packages("rgdal")
install.packages("raster")
install.packages("elevatr")
install.packages("rasterVis")
install.packages("rgl")

library(rasterVis)

Loading required package: raster
Loading required package: lattice
Loading required package: latticeExtra

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

library(elevatr)

2. Obtención datos de elevación Raster Revisar el contenido de la carpeta del dpartamento deamazonas del DANE.

list.files("G:/Geomatica/Datos_de _elevacion_en_R/Amazonas_datos_DANE/91_AMAZONAS/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"    

Ahora, leer el shaoe usando la función del paquete raster:

(munic <- shapefile("G:/Geomatica/Datos_de _elevacion_en_R/Amazonas_datos_DANE/91_AMAZONAS/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 11 
extent      : -74.39634, -69.3955, -4.229406, 0.09725686  (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  :         91,      91001, EL ENCANTO, Decreto 106 de Enero 18 de 1984, 1432.79497455,      2017,   AMAZONAS, 1.96885902115, 0.116086252919 
max values  :         91,      91798,  TARAPACÁ,     ORD 24 DE AGOSTO 01 DE 1997, 16866.1197714,      2017,   AMAZONAS, 10.6360457915,  1.36686628423 

Atributos de munic:

head(munic)

Seleccionar solamente la capital ldel departamento

ltc <- munic[munic$MPIO_CNMBR=="LETICIA",]
plot(ltc, main="Leticia", axes=TRUE)
plot(munic, add=TRUE)

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

Descargar los datos de elevación:

elevacion <- get_elev_raster(ltc, z=9)

Verificar que hay dentro del objeto de elevación:

elevacion

class      : RasterLayer 
dimensions : 2057, 1550, 3188350  (nrow, ncol, ncell)
resolution : 0.00137, 0.00137  (x, y)
extent     : -71.02247, -68.89897, -5.622611, -2.804521  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -1267.138, 1273.577  (min, max)

plot(elevacion, main="DEM descargado (metros)")
plot(ltc, add=TRUE)

3.Recorte de datos para que coincidan con la extension del area de estudio

writeRaster(elevacion, filename = ".dem/elev_z9.tif", dataType('FLT4S', overwrite=TRUE))

elev_crop = crop(elevacion, ltc)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(ltc, add=TRUE)

Ver nuevo proyecto

elev_crop

class      : RasterLayer 
dimensions : 895, 475, 425125  (nrow, ncol, ncell)
resolution : 0.00137, 0.00137  (x, y)
extent     : -70.36487, -69.71412, -4.229321, -3.003171  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -607.5793, 435.8793  (min, max)

4. 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"

Reproyectar los datos de elevación de las coordenadas geograficas WGS84 en la zona MAGNA de Colombia Bogotá

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

rep_elev

class      : RasterLayer 
dimensions : 1362, 731, 995622  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1412458, 1485558, 22789.2, 158989.2  (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     : -470.7503, 401.2846  (min, max)

Ahora reproyectamos el Spatial PolygonsDataFrame que representa la capital de amazonas:

(rep_ltc = spTransform(ltc,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 1412656, 1485526, 22875.46, 158875.8  (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       :         91,      91001,    LETICIA, Decreto 352 de Feb 20 de 1964, 6183.12917992,      2017,   AMAZONAS, 3.88000064403, 0.500701815381 

plot(rep_elev,main="Modelo de elevación digital reproyectado")
plot(rep_ltc, add=TRUE)

writeRaster(rep_elev, filename = "./rep_ltc_elev.tif", datatype = 'INT4S', overwrite = TRUE)

5. Estadísticas básicas de los datos de elevación

Exploración de las estadísticas de los DEM

hist(rep_elev)

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

Crear vector unidimensional

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

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

6. 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 elevación

plot(rep_elev,main="DEM para Leticia (metros)", col=terrain.colors(25,alpha=0.7))

Plot pendiente

plot(slope,main="Pendiente para Leticia (grados)", col=rainbow(25,alpha = 0.7))

Plot Orientación

plot(aspect,main="Exposición para Leticia (grados)", col=rainbow(25,alpha = 0.7))

Plot combinado

plot(hill,
     col=grey(1:100/100),
     legend=FALSE,
     main="DEM para Leticia",
     axes=FALSE)

plot(rep_elev,
     axes=FALSE,
     col=rainbow(12,alpha = 0.35), add=TRUE)

**7. Mapeo de datos de elevación con rayshader

install.packages("rayshader")

library(rayshader)

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 731x1362."

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

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

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

8.Otra forma de vizualizació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[cellr]
  template[[3]] = i[cellb]
  
  template =  template * 256
  
  template
}

map = readJPEG("G:/Geomatica/Datos_de _elevacion_en_R/Amazonas_datos_DANE/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 Leticia")