R Notebook
28.03.2020
#librerias a utilizar
library(rgdal)## Loading required package: sp## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: C:/Users/juane/Documents/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: C:/Users/juane/Documents/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1library(raster)
library(elevatr)
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtralibrary(rgl)revisar el contenido del folder
list.files("D:/unal/geomatica/santander/marco geoestadisitico/68_SANTANDER/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"linea para leer el archivo .shp con una funcion de la libreria raster
(municipios<-shapefile("D:/unal/geomatica/santander/marco geoestadisitico/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))## class : SpatialPolygonsDataFrame
## features : 87
## extent : -74.52895, -72.47706, 5.707536, 8.14501 (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 : 68, 68001, AGUADA, 1539, 19.69444164, 2017, SANTANDER, 0.270018178002, 0.00160984367921
## max values : 68, 68895, ZAPATOCA, Ordenanza 33 de 1968, 3174.27867176, 2017, SANTANDER, 4.36038638488, 0.259459176725atributos del objeto municipios
head(municipios)## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC
## 0 68 68001 BUCARAMANGA 1623
## 1 68 68013 AGUADA 1775
## 2 68 68020 ALBANIA Ordenanza 33 de 1919
## 3 68 68051 ARATOCA 1750
## 4 68 68077 BARBOSA Ordenanza 30 del 25 de Abril de 1936
## 5 68 68079 BARICHARA 1799
## MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
## 0 152.90497 2017 SANTANDER 0.6922752 0.012513526
## 1 75.23125 2017 SANTANDER 0.4758098 0.006146093
## 2 166.21697 2017 SANTANDER 0.8761299 0.013570466
## 3 169.79155 2017 SANTANDER 0.6746922 0.013882031
## 4 46.66489 2017 SANTANDER 0.2703415 0.003810882
## 5 137.27581 2017 SANTANDER 0.5610888 0.011223345ploteo de la capital del departamento
bucara<-municipios[municipios$MPIO_CNMBR=="BUCARAMANGA",]
plot(bucara, main="Bucaramanga", axes=TRUE)
plot(municipios, add=TRUE)
invisible(text(coordinates(municipios), labels=as.character(municipios$MPIO_CNMBR), cex=0.8))
# descarga de los datos de elevacion ubicado en la pagina Amazon Web Services
elevation <- get_elev_raster(bucara, z = 10)## Merging DEMs## Reprojecting DEM to original projection## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0atributos del objeto recien creado
elevation## class : RasterLayer
## dimensions : 1033, 1545, 1595985 (nrow, ncol, ncell)
## resolution : 0.000687, 0.000682 (x, y)
## extent : -73.48, -72.41858, 6.661371, 7.365877 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 108.0003, 4506.756 (min, max)ploteo del objeto elevation y encima de este el ploteo de la capital del departamento
plot(elevation, main="This the downloaded DEM [meters]")
plot(bucara, add=TRUE)
# writeRaster(elevation, filename=“./elev_z8.tif”, datatype=‘INT4S’, overwrite=TRUE) #ploteo del modelo de elevacion pero zolo de lazona delimitada por el objeto bucara
elev_crop = crop(elevation, bucara)
plot(elev_crop, main="Cropped Digital elevation model")
plot(bucara, add=TRUE)
#atributos del objeto elev_crop
elev_crop## class : RasterLayer
## dimensions : 292, 183, 53436 (nrow, ncol, ncell)
## resolution : 0.000687, 0.000682 (x, y)
## extent : -73.17222, -73.0465, 7.067161, 7.266305 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 507.3208, 2260.499 (min, max)#reproyectacion de los datos de elevacion en un nuevo sistema de coordenadas (MAGNA-Colombia)
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 "ajuste de las caracterisicas del raster
pr3 <- projectExtent(elev_crop, spatialref)
# Adjust the cell size
res(pr3) <- 100
# now project
rep_elev <- projectRaster(elev_crop, pr3)#atributos del raster creado recien
rep_elev## class : RasterLayer
## dimensions : 221, 139, 30719 (nrow, ncol, ncell)
## resolution : 100, 100 (x, y)
## extent : 1099976, 1113876, 1273307, 1295407 (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 : 507.6196, 2258.253 (min, max)reproyeccion de los datos de la capital en un SpatialPolygonsDataFrame
(rep_bucara = spTransform(bucara,spatialref))## class : SpatialPolygonsDataFrame
## features : 1
## extent : 1100015, 1113881, 1273395, 1295432 (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 : 68, 68001, BUCARAMANGA, 1623, 152.90496962, 2017, SANTANDER, 0.692275167631, 0.0125135257452ploteo de los datos reproyectados
plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_bucara, add=TRUE)
writeRaster(rep_elev, filename="./rep_bucara_elev.tif", datatype='INT4S', overwrite=TRUE)estadisticas de los datos de reproyeccion
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))## metricas valores
## 1 mean 1136.9049
## 2 min 507.6196
## 3 max 2258.2531
## 4 std 360.6739estos objetos seran necesarios para las variables geomorfometricas
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)este ploteo muestra los datos de elevacion de la capital bajo el sistema de coordenadas MAGNA
plot(rep_elev,main="DEM for Bucaramanga [meters]", col=terrain.colors(25,alpha=0.7))
# este ploteo muestra el grado de las pendientes existentes en el tereno de la capital
plot(slope,main="Slope for Bucaramanga [degrees]", col=topo.colors(25,alpha=0.7))
#este ploteo muestra el aspecto de las pendientes existentes en el tereno de la capital
plot(aspect,main="Aspect for Bucaramanga [degrees]", col=topo.colors(25,alpha=0.7))
# este ploteo muestra la elevacion de las colinas en la capital, y encima de esto esta el ploteo de la elevacion propia de dicho terreno
plot(hill,
col=grey(1:100/100), # create a color ramp of grey colors for hillshade
legend=FALSE, # no legend, we don't care about the grey of the hillshade
main="DEM for Bucaramanga",
axes=FALSE) # makes for a cleaner plot, if the coordinates aren't necessary
plot(rep_elev,
axes=FALSE,
col=terrain.colors(12, alpha=0.35), add=TRUE) # color method
# esta libreria permitira tener diferentes visualizaciones basadas en datos de elevacion
library(rayshader)para realizar las visalizaciones, hay que convertir la DEM en una matriz, la primera imagen es una visualizacion general del terreno, la segunda muestra ademas el agua existente y la tercera muestra estas dos ultimas mas la sombra del sol sobre el terreno
elmat = raster_to_matrix(rep_elev)
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()
#install.packages("jpeg")
library(jpeg)# funcion hillshade para mapear un ambiente espacial
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")
# Do the mapping
out = getv(map, aspect, slope)
# Plot pretty mountains
plotRGB(out, main = "Supposedly pretty mountains in Bucaramanga")