Experimento de tiempos de desplazamiento
Contenido:- Preparación de datos
- Crear un raster de costes
- Crear raster base y puntos de velocidad
- Crear polígonos de Voronoi a partir de los puntos de velocidades
- Calcular tiempos de desplazamiento para todos los puntos (parallel)
Esta entrada de blog tiene intenciones oscuras, ya que solo es un experimento para resolver un problema de desplazamiento con las siguientes condiciones:
- Hay una serie de puntos que representan puntos de salida.
- Hay una capa de vías que contienen un atributo de velocidad del tramo.
- Hay una capa de ríos, curro cruzamiento representa una modificación en la velocidad.
- Hay una capa de polígonos dentro de los cuales la velocidad se ve modificada.
Estos ingredientes conforman una de las preguntas sin respuesta del GIS Stack Exchange. . Esta es concretamente la cuestión a la que este experimento intenta dar respuesta: Constructing a Voronoi diagram using a complicated travel time metric.
Antes de empezar con proceso se destacará que actualmente hay herramientas potentisimas, libres y de pago para realizar este tipo de análisis, sin necesidad de inferir datos ni cacharrear con los rasters y capas de polígonos Voronoi. Aquí os dejo algunas fuentes de interés sobre este problema:
- 7 herramientas de análisis de redes y rutas en QGIS
- Crear áreas de tiempo de recorrido
- Open Source Routing Machine
Preparación de datos
En primer lugar cargaremos las librerías a utilizar. Leeremos los datos con las capas indicadas anteriormente. Además se les asignaran los valores de velocidad de forma aleatoria. Para este ejercicio nos vale cualquier valor.
library(sf)
library(dplyr)
library(mapview)
library(htmlwidgets)
library(leaflet)
library(raster)
library(ggplot2)
# Read raw data
locs <- st_read("./data/20230126/data-raw/LEIP/LEIP.shp")
roads <- st_read("./data/20230126/data-raw/Carreteras_IMD/Carreteras_IMD.shp")
terrain <- st_read("./data/20230126/data-raw/AEIP/AEIP.shp")
rios <- st_read("./data/20230126/data-raw/Rios/Rios.shp")
# Create random speeds for features
r <- r %>% mutate(SPEED = sample(50, size = nrow(r))) %>% st_cast("LINESTRING")
w <- w %>% mutate(SPEED = sample(500, size = nrow(w))/10) %>% st_cast("LINESTRING")
t <- t %>% mutate(SPEED = sample(50, size = nrow(t)))
l1 <- mapview(l %>% dplyr::select(-ID), alpha=0, leyend = FALSE)
l2 <- mapview(r, zcol = "SPEED")
l3 <- mapview(w, lwd=0.5)
l4 <- mapview(t, alpha.regions = 0.2, zcol = "SPEED")
Crear un raster de costes
En este caso, el proceso es el siguiente:
- Las carreteras a puntos.
- Los rios a puntos.
- Las áreas a puntos (a partir de un sampleado regular dentro del contorno de los polígonos).
wpoints <- w %>% st_cast("POINT")
rpoints <- r %>% st_cast("POINT")
tpoints <- st_sample(t, size = 500, type = "regular") %>% st_as_sf() %>% st_join(t, join = st_nearest_feature)
# merge all points with speeds
speedpoints <- bind_rows(wpoints, rpoints, tpoints) %>% dplyr::select(SPEED)
head(as.data.frame(speedpoints %>% dplyr::select(-geometry)))
## SPEED geometry
## 343 28.8 POINT (602942 4789750)
## 343.1 28.8 POINT (602888 4789734)
## 343.2 28.8 POINT (602851 4789719)
## 343.3 28.8 POINT (602827 4789711)
## 343.4 28.8 POINT (602805 4789698)
## 343.5 28.8 POINT (602743 4789644)
Crear raster base y puntos de velocidad
Se creará un raster con una resolución de 50 m para realizar los cálulos más rápido.
#función para empty raster (raster base)
library(raster)
library(fasterize)
rasterbase <- function(layer, res){
r <- raster(extent(layer), res)
res(r) <- c(res, res)
crs(r) <- crs(layer)
values(r) <- NA
return(r)
}
raster <- rasterbase(l, 50)
Crear polígonos de Voronoi a partir de los puntos de velocidades
Se crearán los polígonos a partir de los puntos con los valores de velocidad, y se rasterizará el conjunto de polígonos voronoi para crear un mapa raster de velocidad.
g = st_combine(st_geometry(speedpoints))
x <- st_voronoi(g, bOnlyEdges = FALSE)
v = st_as_sf(st_collection_extract(x)) %>%
st_join(speedpoints, join = st_nearest_feature) %>%
st_make_valid() %>%
st_cast("POLYGON") %>%
st_crop(l)
speeds <- fasterize(v, raster, field = "SPEED")
Ahora, siguiendo el proceso realizar ya en un anterior posts (Shortest path with Tobler function), se realizarán los raster de coste de desplazamiento y tiempo acumulado de desplazamiento.
library(gdistance)
tr <- transition(speeds, transitionFunction=mean, directions=8)
trc <- geoCorrection(tr, type="c", multpl=FALSE, scl=TRUE)
p <- l %>% dplyr::filter(ID == "1905")
ac <- accCost(trc, st_coordinates(p))
par(mfrow=c(1,2))
plot(raster(trc))
points(as_Spatial(p))
text(st_coordinates(p), labels = "start", pos = 4)
plot(ac)
points(as_Spatial(p))
text(st_coordinates(p), labels = "start", pos = 4)
Isocronas
library(rasterVis)
levels <- seq(min(values(ac)), max(values(ac)), 2)
rasterVis::levelplot(ac, par.settings = BuRdTheme, margin=FALSE,
contour=TRUE, at = levels, lty=3,
labels = list(cex = 0.5), label.style = 'align',
pretty = TRUE)
Polígonos de Voronoi con tiempos de desplazamiento (respuesta a GIS Stack Exchange)
ttras <- round(ac, digits = 1)
polygonizator <- function(raster){
sf::st_as_sf(stars::st_as_stars(raster),
as_points = FALSE,
merge = TRUE)
}
ttpol <- polygonizator(ttras) %>% rename(TIME = layer)
Recordatorio: Esta capa de tiempo de viaje es una aproximación al problema. No es válida por varias razones:
- Como la velocidad de desplazamiento procede de un raster burdo, se supone que el acceso se hace igual por todo el territorio, lo cual es incorrecto pues realmente no se puede caminar por todas partes.
- Las velocidades de desplazamiento se han asignado a las entidades vectoriales originales de forma aleatoria, por lo que no se puede extraer conclusiones de la aproximación.
- Se han mantenido los puntos extraidos de velocidad en carretera y no se han reescrito totalmente con las velocidades asignadas a las áeras. De esta forma se mantienen ambos valores, lo cual genera zonas de desplazamiento raro, pues puede haber puntos próximos con velocidades muy dispares.
Se podría volver a generar el diagrama de voronoi a partir de estos resultados, aplicando las herramientas vistas anteriormente. Habría que samplear puntos en la capa y volver a construir el diagrama. Realmente no aportaría nueva información y reduciría la compresión del problema, pero ahí va:
ttpoints <- rasterToPoints(ttras, spatial=TRUE)
ttpoints <- ttpoints %>%
st_as_sf() %>%
st_transform(st_crs(l)) %>%
st_make_valid() %>%
unique() %>%
dplyr::rename(TIME=layer)
# Reducir el conjunto de datos para facilitar el cálculo
if(nrow(ttpoints) > 2000){ttpoints <- sample_frac(ttpoints, 1/(nrow(ttpoints)/2000))}
g = st_combine(st_geometry(ttpoints)) %>% st_drop_geometry()
x <- st_voronoi(g, bOnlyEdges = FALSE)
v = st_as_sf(st_collection_extract(x)) %>%
st_join(ttpoints, join = st_nearest_feature) %>%
st_make_valid() %>%
st_cast("POLYGON") %>%
st_crop(l)
library("ggspatial")
ggplot() +
geom_sf(data = v, aes(fill = TIME)) +
geom_sf(data = p, size = 5, shape = 3)+
geom_label(data = p, aes(p$geometry[[1]][1], p$geometry[[1]][2], label = "Start"), size = 5)+
scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
xlab("x") + ylab("y")+
theme_minimal()
Calcular tiempos de desplazamiento para todos los puntos (parallel)
Llegados a este punto, vamos a formalizar una función con la que calcular los tiempos de desplazamiento para cada uno de los puntos del territorio. Haremos que de la función se derive 1 única capa de información geográfica (ojo, que si el proceso es enorme puedes meterte en problemas). Por otra parte, se dejará fuera del proceso la reconversión a voronoi.
# Función individual
traveltime <- function(p, trc){
# p <- p <- l[1,]
ac <- accCost(trc, st_coordinates(p))
ttras <- round(ac)
ttpol <- polygonizator(ttras) %>% rename(TIME = layer)
ttpol$ID <- p$ID
return(ttpol)
}
# lsel <- l[1:20,]
lsel <- l
library(parallel)
cl <- parallel::makeCluster(detectCores(), type="FORK")
doParallel::registerDoParallel(cl)
TTIMES <- parLapply(cl, 1:nrow(lsel), function(x){traveltime(lsel[x,], trc)})
stopCluster(cl)
Unir todos los resultados en una única capa, guardar y mostrar.
results <- bind_rows(TTIMES)
length(unique(results$ID))
## [1] 184
Seleccionar punto para mostrar resultados sobre uno de los lugares calculados.
p1 <- l %>% dplyr::filter(ID == "1905")
t1 <- results %>% dplyr::filter(ID == "1905")
