This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

Install packages before load them

#Para Mac se necesita Xquartz. Instalar…

install.packages(c(‘elevatr’, ‘sf’, ‘leaflet’, ‘terra’, ‘MultiscaleDTM’, ‘exactextractr’))

library(elevatr) library(sf) library(leaflet) library(terra) library(MultiscaleDTM) #to use this one, it´s needed to install rlg package library(exactextractr)

Cargar capa

munic <- st_read(“~/Desktop/cordoba/Cordoba.gpkg”)

Conectar con Marco Geoestadistico Nacional de Referencia

(cordoba_crs = st_crs(st_geometry(munic)))

Area municipios de Córdoba

munic$area_m2 = sf::st_area(munic)

el área expresada en kilómetros cuadrados

munic\(area = munic\)area_m2/1000000

#redondear el area a dos digitos munic\(area = signif(munic\)area, digits = 6)

#llamar municipios munic

#nuevo objeto con centroides de municipios (centers <- st_centroid(munic))

#dividir cada uno en sus coordenadas x, y: centers\(x = st_coordinates(centers)[,1] centers\)y = st_coordinates(centers)[,2]

#shora se comprueba centers

#ahora se obtiene la DEM elevation <- get_elev_raster(munic, z = 10) Mosaicing & Projecting Note: Elevation units are in meters.

#objeto de elevación elevation

#escribir el DEM en su directorio local usando la función writeRaster writeRaster(elevation, “~/Desktop/cordoba/elev_cordoba_z10.tif”, overwrite=TRUE)

#convertir la elevación en un objeto ráster de Terra (elevation2 <- terra::rast(elevation))

#Visualizar los datos de elevación.Primero, necesitaremos una paleta de colores pal <- colorNumeric(c(“cyan”, “forestgreen”,“yellow”,“tan”,“orange”, “brown”), values(elevation), na.color = “transparent”)

#Ahora, reduciremos la resolución espacial de los datos (elevation3 <- terra::aggregate(elevation2, fact = 2))

#recortaremos el objeto de elevación3 para que se ajuste a los límites de cordoba elevation4 <- terra::crop(elevation3, munic, mask=TRUE)

#Ahora es el momento de trazar leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% + addPolygons(color = “white”, weight = 1.0, smoothFactor = 0.5, + opacity = 0.25, fillOpacity = 0.15, + popup = paste(“Municipio:”, munic\(mpio_cnmbr, "<br>", + "Km2: ", munic\)area, “
”)) %>% + addLabelOnlyMarkers(data = centers, + lng = ~x, lat = ~y, label = ~mpio_cnmbr, + labelOptions = labelOptions(noHide = TRUE, direction = ‘top’, textOnly = TRUE, textsize = “10px”)) %>% + addRasterImage(elevation4, colors = pal, opacity = 0.9) %>% + addLegend(pal = pal, values = values(elevation), + title = “Elevation data for Cordoba (mts)”)

#… sessionInfo()

#DEM usando terra (dem = terra::rast(“~/Desktop/cordoba/elev_cordoba_z10.tif”))

#ahora vamos a reducir la resolución del DEM para evitar problemas de memoria: dem2 = terra::aggregate(dem,2, “mean”)

#Ahora leeremos nuestros municipios usando la biblioteca (munic <- sf::st_read(“~/Desktop/cordoba/cordoba.gpkg”))

#Ahora, recortemos el objeto dem2 a los límites de nuestro departamento. dem3 = terra::crop(dem2,munic, mask=TRUE)

#se Modifican los siguientes dos fragmentos de código para transformar los datos de entrada al https://rodolfofrancoweb.com/sig/proyecciones_y_sistemas_de_coordenadas/magna-sirgas-origen-nacional/ (dem_plane = project(dem3, “EPSG:9377”))

#necesitamos realizar una transformación similar para el objeto vectorial (munic_plane = sf::st_transform(munic, “EPSG:9377”))

#Ahora, llamemos a la función SlpAsp (slp_asp = MultiscaleDTM::SlpAsp( + dem_plane, + w = c(3, 3), + unit = “degrees”, + method = “queen”, + metrics = c(“slope”, “aspect”), + na.rm = TRUE, + include_scale = FALSE, + mask_aspect = TRUE))

Se divide el objeto slp_asp en dos partes. Capa 1

(slope = subset(slp_asp, 1))

#¿Cuál es la distribución de los valores de pendiente? terra::hist(slope, main = “Cordoba’s slope”, xlab = “Slope (in degrees)”)

#Capa 2 (aspect =subset(slp_asp, 2))

#terra terra::hist(aspect, main = “Cordoba’s aspect”, xlab = “Aspect (in degrees)”)

#convertimos los grados de pendiente en porcentaje de pendiente (slope_perc = tan(slope(pi/180))100)

#terra terra::hist(slope_perc, main = “Cordoba’s slope”, xlab = “Slope (in percentage)”)

Reclassify the slope raster

#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)

m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3, 12, 25, 4, 25, 50, 5, 50, 75, 6, 75, 160, 7) m <- matrix(m, ncol=3, byrow = TRUE) rc <- classify(slope_perc, m, right=TRUE)

#Ahora, calculemos estadísticas zonales (munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, ‘mean’)) #### Warning in .local(x, y, …): Polygons transformed to raster CRS (EPSG:9377)

#pendiente media de córdoba hist(munic$mean_slope, main = “Cordoba’s mean slope”, xlab = “Slope (in percentage)”)

#hist(munic$class, main = “Cordoba’s reclassified slope”, xlab = “Slope (as a category)”) #terra::hist(slope_perc)

#Transformemos la pendiente nuevamente a coordenadas geográficas (rc.geo = project(rc, “EPSG:4326”))

#Esta es la pendiente porcentual Esta es la pendiente porcentual

#Primero, necesitaremos una paleta de colores para la pendiente

#Ahora es el momento de trazar leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% addPolygons(color = “gray”, weight = 1.0, smoothFactor = 0.5, opacity = 0.4, fillOpacity = 0.10, popup = paste(“Municipio:”, munic\(mpio_cnmbr, "<br>", "Slope class: ", munic\)class, “
”)) %>% addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>% addLegend(pal = palredgreen, values = values(slope.geo), title = “Terrain slope in Cordoba (%)”)

#pendiente del terreno en grados y agregue etiquetas que muestren los nombres de los municipios junto con sus correspondientes valores de pendiente media. ## Si la pendiente está en porcentaje, conviértela a grados slope.geo <- atan(slope.geo / 100) * (180 / pi)

Extraer valores de pendiente media para cada municipio

valores_pendiente <- extract(slope.geo, munic, fun = mean, na.rm = TRUE) ##Mensaje de advertencia: [extraer] transformando datos vectoriales al CRS del ráster munic$pendiente_media <- valores_pendiente

#hora es el momento de trazar

leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, zoom = 9) %>% addPolygons(color = “gray”, weight = 1.0, smoothFactor = 0.5, opacity = 0.4, fillOpacity = 0.10, popup = paste(“Municipio:”, munic\(mpio_cnmbr, "<br>", "Pendiente media: ", round(munic\)pendiente_media, 2), “°
”, “Slope class:”, munic$class, “
”)) %>% addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>% addLegend(pal = palredgreen, values = values(slope.geo), title = “Pendiente del terreno (°)”)

##Mensaje de advertencia:La capa sf tiene datos inconsistentes (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).Necesita ’+proj=longlat +datum=WGS84

#exolicación codigo Explicación del código: setView(-75.9, 8.7, zoom = 9): Ajusta la vista inicial del mapa (con las coordenadas lon y lat que desees y el nivel de zoom). addPolygons(): Dibuja los límites de los municipios y agrega un popup que muestra el nombre del municipio, la pendiente media en grados y la clase de pendiente. addRasterImage(): Añade el raster de pendiente al mapa, donde palredgreen es una paleta de colores (puedes elegir la que prefieras). addLegend(): Añade una leyenda para la pendiente del terreno en grados