Análisis espacial 3
Noé Osorio
2025-03-22
Orientación
Aspecto deriva la orientación de cada celda de una superficie ráster la orientación identifica curvatura: calcula la curvatura de una superficie ráster y opcionalmente incluye la curvatura del perfil y la del plano La herramienta Curvatura ajsuta un plano a las nueve celdas locales, pero un plano puede no ser una buena descripción del paisaje, y puede exagerar las variaciones naturales de interés valores de curvatura
Cuenvas de visión
determina las ubicaciones de la superficie ráster visibles para un conjunto0 de entidades del observador, a partir de donde está localizado se determina donde es visibley dónde no es visible
VIsibilidad
Determina las ubicaciones d elas superficie ráster visibles para un conjunto de características del observador o identifica que puntos son los buenos para ver
información solicitada
Puntos de observació: Puntos de observaión tiene como objetivo exposición visual de la superficie
- offseta indica la distancia en unidades de superficie que serán añadidas al valor z del punto de observación respectivo, Si Offseta está disponible
AZIMUTH: angulo horizontal de escaneo, el barrido se realiza en sentido de agujas del reloj
- vert1 limite inferior
- ver2 limite superior
Descargar un modelo de terreno
CEM INEGI
Al final descargamos tlaxcala, la razón: porque está más pequeño
Tiene una proyección lon/lat GCS_MEXICO_ITRF_2008
cem = terra::rast("../01_input/análisis espacial/CEM_V3_20170619_R15_E29_TIF/Tlaxcala_r15m.tif")
cem## class : SpatRaster
## dimensions : 4492, 7797, 1 (nrow, ncol, nlyr)
## resolution : 0.0001388889, 0.0001388889 (x, y)
## extent : -98.7084, -97.62548, 19.10507, 19.72896 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat GCS_MEXICO_ITRF_2008
## source : Tlaxcala_r15m.tif
## name : Tlaxcala_r15m
## min value : 2025
## max value : 4408se reproyecta el raster a 32614 y se visualiza con la función “plet”
utm_zona14n <- "+proj=utm +zone=14 +datum=WGS84 +units=m +no_defs"
cem<- terra::project(cem, utm_zona14n)
terra::plet(cem)Aspecto
El aspecto es lo que nos indica la dirección de la pendiente del raster, usando las funciones de terrain con el parámetro V es el aspecto, el rango es de 0 a 360 donde indica el ángulo en funciones norte, sur este oeste etc, como rosa de vientos
orientacion <- terra::terrain(x = cem, v = "aspect")
terra::plet(orientacion)Sombreado genera una visión
El hillshade es una técnica de visualización que simula la iluminación de una superficie considerando la posición del sol, lo que ayuda a resaltar el relieve del terreno.
Parámetros:
- slope
- aspect
- altitude: ángulo de elevación del sol (0-90 grados)
- azimuth: dirección de la fuente de luz (0-360 grados, 0=norte, 90=este)
- z: factor de exageración vertical (1 = sin exageración)
angle = 45, # ángulo de elevación del sol
direction = 315, # dirección de la fuente de luz (noroeste)
el 45 y 315 en priogramas como qgis y arcgis vienen definidos ya o sea se pueden cambiar
hillshade <- terra::shade(
slope = terra::terrain(cem, "slope"),
aspect = terra::terrain(cem, "aspect"),
angle = 45, # ángulo de elevación del sol
direction = 315, # dirección de la fuente de luz (noroeste)
normalize = TRUE
)
terra::plet(hillshade)Curvas de nivel
seleccionar modelo de terreno, y generar las curvas cada 200 metros
library(terra)## Warning: package 'terra' was built under R version 4.4.3## terra 1.8.29library(ggplot2)
library(tidyterra)## Warning: package 'tidyterra' was built under R version 4.4.3##
## Adjuntando el paquete: 'tidyterra'## The following object is masked from 'package:stats':
##
## filterintervalo = 200 #
min_elev <- minmax(cem)[1]
max_elev <- minmax(cem)[2]
niveles <- seq(
from = ceiling(min_elev / intervalo) * intervalo, # Redondear al siguiente intervalo
to = floor(max_elev / intervalo) * intervalo, # Redondear al anterior intervalo
by = intervalo
)
# Extraer curvas de nivel
curvas_nivel <- terra::as.contour(
cem,
levels = niveles
)
# Convertir a objeto sf para manipulación adicional
curvas_nivel_sf <- sf::st_as_sf(curvas_nivel)
# 2. Visualizar las curvas de nivel
plot(cem, main = "Modelo Digital de Elevación")
plot(curvas_nivel, add = TRUE)
library(tidyterra)
ggplot() +
geom_spatraster(data = cem) +
scale_fill_whitebox_c() +
geom_sf(data = curvas_nivel_sf, aes(color = level), linewidth = 0.5) +
scale_color_whitebox_c(
palette = "muted",
na.value = "white"
)+
labs(
title = "Curvas de Nivel",
subtitle = paste("Intervalo:", intervalo, "metros")
) +
theme_minimal() +
coord_sf()## <SpatRaster> resampled to 501018 cells.
Puntos altos o de visión
con la librería map edit y la geometría de puntos de nivel seleccioné los puntos donde según yo era lo más alto, los puntos son los siguientes:
library(leaflet)
observacion1 = sf::st_read("../01_input/análisis espacial/puntos_altos.gpkg")## Reading layer `puntos_altos' from data source
## `C:\Users\DELL\Desktop\Dip_UNAM_2025\01_input\análisis espacial\puntos_altos.gpkg'
## using driver `GPKG'
## Simple feature collection with 5 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 536513.7 ymin: 2126671 xmax: 611605.5 ymax: 2171545
## Projected CRS: WGS 84 / UTM zone 14Nleaflet() %>%
addProviderTiles(providers$Esri.WorldTopoMap) %>%
addMarkers(data=observacion1 %>% sf::st_transform(crs=4326))View shed
identificar campo de visión a través de puntos se observación, o sea los puntos de acá arriba cual ver mejor.
vs <- viewshed(loc=sf::st_coordinates(observacion1[1,]),
cem
# observer=1.7
)terra::plot(vs)
terra::plot(observacion1, add = TRUE)## Warning in plot.sf(observacion1, add = TRUE): ignoring all but the first
## attribute
EEso es lo que ves desde ahí?
factor_agregacion = 5
# Realizar el resample usando aggregate
# fact: factor de agregación en x y y
# fun: función para agregar valores (mean, max, min, median, sum, etc.)
cem_agregado = terra::aggregate(
vs,
fact = factor_agregacion
)## |---------|---------|---------|---------|========================================= leaflet() %>%
addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
addRasterImage(x = cem_agregado) %>%
addMarkers(data=observacion1[1,] %>% sf::st_transform(crs=4326))O sea queda pendiente hacer el análisis con los otros puntos
Area solar radiation
se usará un área más pequeña que será con la zona de la malinche
poligono_chico = sf::st_read("../01_input/análisis espacial/poligono_chico.gpkg")## Reading layer `poligono_chico' from data source
## `C:\Users\DELL\Desktop\Dip_UNAM_2025\01_input\análisis espacial\poligono_chico.gpkg'
## using driver `GPKG'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 587661 ymin: 2118953 xmax: 620230.5 ymax: 2144426
## Projected CRS: WGS 84 / UTM zone 14Nterra::plet(terra::vect(poligono_chico))## Warning: PROJ: proj_create_from_database: Cannot find proj.db (GDAL error 1)A ese polígono le cortaré el pedazo del raster, con la librería map edit hice el polígono y la función mask permite hacer el corte al polígono de referencia,
el primer objeto es cem que es el model de elevacióny en el mask usé la función vect para pasarlo a un objeto terra porque sino no sale o no sé si sí funciones con sf, probaré
cortado = terra::mask(x = cem,
mask = terra::vect(poligono_chico))
terra::plet(
cortado
)sí, sí funcionó, pensé que terra solo interactuaba con objetos terra pero igual permite sf
cortado = terra::mask(x = cem,
mask = poligono_chico)
terra::plet(
cortado
)Ahora con el objeto cortado se analizará el area solar radiation ahora yaque está cortado haremos lo de radiación y espero que no se tarde. del 1 al 2 de enero del 2025
library(solaR)## Warning: package 'solaR' was built under R version 4.4.3## Cargando paquete requerido: zoo## Warning: package 'zoo' was built under R version 4.4.3##
## Adjuntando el paquete: 'zoo'## The following object is masked from 'package:terra':
##
## time<-## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Cargando paquete requerido: lattice## Cargando paquete requerido: latticeExtra## Warning: package 'latticeExtra' was built under R version 4.4.3##
## Adjuntando el paquete: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layer## Time Zone set to UTC.##
## Adjuntando el paquete: 'solaR'## The following object is masked from 'package:terra':
##
## comparelibrary(solrad)## Warning: package 'solrad' was built under R version 4.4.3fecha_inicio = as.POSIXct("2025-01-01", tz = "UTC")
fecha_fin = as.POSIXct("2025-01-07 23:59:59", tz = "UTC")
fechas = seq(fecha_inicio, fecha_fin, by = "hour")
dem = cortado
dem_extent = ext(dem)
lon_central = (dem_extent[1] + dem_extent[2]) / 2
lat_central = (dem_extent[3] + dem_extent[4]) / 2pendiente <- terrain(dem, "slope")
orientacion <- terrain(dem, "aspect")
# Crear un objeto de parámetros solares con solaR para la ubicación central
params_solares <- calcSol(lat = 19.228438,
BTd = fechas,
keep.night = TRUE
)
params_solares## Object of class Sol
##
## Latitude: 19.2 degrees
##
## Daily values:
## Index decl eo
## Min. :2025-01-01 00:00:00 Min. :-0.4006 Min. :1.035
## 1st Qu.:2025-01-02 12:00:00 1st Qu.:-0.3982 1st Qu.:1.035
## Median :2025-01-04 00:00:00 Median :-0.3955 Median :1.035
## Mean :2025-01-04 00:00:00 Mean :-0.3952 Mean :1.035
## 3rd Qu.:2025-01-05 12:00:00 3rd Qu.:-0.3925 3rd Qu.:1.035
## Max. :2025-01-07 00:00:00 Max. :-0.3892 Max. :1.035
## ws Bo0d EoT
## Min. :-1.427 Min. :7320 Min. :-0.03057
## 1st Qu.:-1.426 1st Qu.:7341 1st Qu.:-0.02800
## Median :-1.425 Median :7364 Median :-0.02538
## Mean :-1.425 Mean :7366 Mean :-0.02532
## 3rd Qu.:-1.424 3rd Qu.:7390 3rd Qu.:-0.02269
## Max. :-1.423 Max. :7418 Max. :-0.01994
##
## Intradaily values:
## Index w aman
## Min. :2025-01-01 00:00:00 Min. :-2.9068 Min. :0.0000
## 1st Qu.:2025-01-02 17:45:00 1st Qu.:-1.3989 1st Qu.:0.0000
## Median :2025-01-04 11:30:00 Median : 0.1089 Median :0.0000
## Mean :2025-01-04 11:30:00 Mean : 0.1089 Mean :0.4583
## 3rd Qu.:2025-01-06 05:15:00 3rd Qu.: 1.6167 3rd Qu.:1.0000
## Max. :2025-01-07 23:00:00 Max. : 3.1265 Max. :1.0000
## cosThzS AlS AzS Bo0
## Min. :-0.9982 Min. :-1.5116 Min. :-1.36860 Min. : 0.0
## 1st Qu.:-0.7377 1st Qu.:-0.8298 1st Qu.:-1.14315 1st Qu.: 0.0
## Median :-0.1280 Median :-0.1284 Median : 0.09400 Median : 0.0
## Mean :-0.1265 Mean :-0.1876 Mean : 0.00506 Mean : 308.0
## 3rd Qu.: 0.4815 3rd Qu.: 0.5024 3rd Qu.: 1.12535 3rd Qu.: 681.3
## Max. : 0.7483 Max. : 0.8455 Max. : 1.36988 Max. :1058.7
## rd rg
## Min. :0.00000 Min. :0.00000
## 1st Qu.:0.00000 1st Qu.:0.00000
## Median :0.00000 Median :0.00000
## Mean :0.04181 Mean :0.04140
## 3rd Qu.:0.09307 3rd Qu.:0.08693
## Max. :0.14321 Max. :0.15452Lo que no entiendo es de dónde salen los datos