In Robinlovelace/pctSeville:

---

title: "Exploring cycling potential in Seville with official data" author: "Robin Lovelace" date: "r Sys.Date()" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Vignette Title} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8}

---

The Seville data can be loaded as follows, using a function in the pctSeville package:

library(pctSeville)
library(stplanr)
if(grepl(pattern = "vign", x = getwd())) {
  indir = "../pctSeville-data/"
} else {
  indir = "pctSeville-data/"
}
sev_dat = load_city_sev(data_dir = indir)
library(tidyverse)
sev_dat$su09_grid_poblacion_amsevilla =
  rename(sev_dat$su09_grid_poblacion_amsevilla, POBLACION = Pob_Tot)

The results can be viewed as follow:

sapply(sev_dat, nrow)

library(sf)
plot(sev_dat$corredores_verdes[1])

Plot the results

The cycle paths, green corridors and 3000 m buffer are plotted below.

library(tmap)
tmap_mode("view")
sev_dat$corredores_verdes = st_zm(sev_dat$corredores_verdes)
qtm(sev_dat$buffer_3000_definitivo_sin_rodrigo) +
  qtm(sev_dat$corredores_verdes, lines.col = "green") +
  qtm(sev_dat$InventarioSevilla) +
  qtm(sev_dat$MetropolitanasSevilla)

The population and rail data are illustrated below.

qtm(sev_dat$su08_grid_poblacion_amsevilla, fill = "POBLACION") +
  # qtm(sev_dat$vc01_1_carretera_amsevilla) + # roads
  # qtm(sev_dat$vc03_ffcc_amsevilla) + # rail links
  qtm(sev_dat$vc03_ffcc_amsevilla) 

The zones area represented in grey:

sev_dat$todas_zonas$NOM_ARTICU = as.numeric(sev_dat$todas_zonas$NOM_ARTICU)
sev_dat$todas_zonas = dplyr::arrange(sev_dat$todas_zonas, NOM_ARTICU)
qtm(sev_dat$todas_zonas) +
  qtm(sev_dat$su08_grid_poblacion_amsevilla, fill = "POBLACION") +
  qtm(sev_dat$todas_zonas[c(5, 7, 14),], fill = "red") +
  tm_scale_bar()

Propensity to cycle analysis

To explore the propensity to cycle, we will use the example of the Dos Hermanas catchement.

From population to cyling propensity to stations. The stations in the study area can be found and plotted as follows:

sel_ex = sev_dat$todas_zonas$NOM_ARTICU == 14
expo = sev_dat$todas_zonas[sel_ex,]
expo_buff = st_buffer(expo, dist = 500)
library(osmdata)
q = add_feature(opq = opq(bbox = c(-6.08, 37.32, -6.03, 37.38)), key = "railway", value = "station")
res_osm = osmdata_sf(q)
stations = res_osm$osm_points
grep(pattern = "San Juan Alto", stations$name)

# failed attempt to find san juan alto
# stations_polypoints = st_centroid(res_osm$osm_polygons)
# num_nas = sapply(stations, function(x) sum(is.na(x)))
# names_keep = names(head(sort(num_nas), n = 4))
# stations = stations[names_keep]
# stations_polypoints = stations_polypoints[names_keep]
# stations = rbind(stations, stations_polypoints)

st_crs(stations)
stations = st_transform(stations, st_crs(expo))
plot(stations$geometry)
plot(expo, add = TRUE)
stations = stations[expo_buff,]
stations_buff = st_buffer(stations, dist = 500) %>% st_union()
# saveRDS(stations, "data/stations.Rds")
# saveRDS(stations_buff, "data/stations_buff.Rds")

The likely trips to the station can be visualised as follows, for each of the population cells:

pop = sev_dat$su09_grid_poblacion_amsevilla[expo,]
pop_cents = st_centroid(pop)
pop_cents = pop_cents[expo,]
# write_sf(st_transform(pop_cents, 4326), "data/pop_cents.geojson") # only run in interactive mode
# write_sf(st_transform(pop, 4326), "data/pop.geojson") # only run in interactive mode
pop_cents = pop_cents %>% filter(!st_within(., stations_buff, sparse = F))
station = st_centroid(expo)
grid_num = 3
l1 = st_linestring(rbind(st_coordinates(stations[1,]), st_coordinates(pop_cents[grid_num,])))
l1 = st_sf(st_sfc(l1), crs = st_crs(pop))
# routing
l1w = as(st_transform(l1, crs = 4326), "Spatial")
lp = line2df(l1w)
r = viaroute(lp$fy, lp$fx, lp$ty, lp$tx, profile = "bicyce")
r = viaroute2sldf(r)
qtm(expo) +
  qtm(pop, fill = "POBLACION") +
  qtm(pop_cents) +
  qtm(stations, symbols.size = 5) +
  qtm(l1) +
  qtm(r, lines.col = "green") +
  tm_scale_bar()

# coefficients of trips
ck = c(0.556, 0.247, 0.137)
pop_cents$POBLACION[grid_num] *
  ck[1] # how many people travel from A to B

# d = rep(seq(from = 100, to = 3000, by = 100), 3)
# dd = rep(c(0.0002, 0.0005, 0.001), each = 3000 / 100)
# pcycle = exp(d * -dd)
# dd = as.character(dd)
# # plot(d, pcycle)
# dd_df = data_frame(`Distance (m)` = d, `Proportion cycling` = pcycle, beta = dd)
# ggplot(dd_df, aes(`Distance (m)`, `Proportion cycling`, colour = beta)) +
#   geom_line() +
#   ylim(c(0, 1))
# filter(dd_df, `Distance (m)` == 500 | `Distance (m)` == 3000)
# 0.77880078 / 0.22313016

Now we are in a position to model cycling potential to stations in the case study area for both stations:

# pop_cents$station = NA
# for(grid_num in seq_along(pop_cents$POBLACION)) {
#   sdist = st_distance(
#     pop_cents[grid_num,],
#     stations
#   )
#   station = stations[which.min(sdist),]
#   pop_cents$station[grid_num] = as.character(station$name)
#   l1 = st_linestring(rbind(st_coordinates(station), st_coordinates(pop_cents[grid_num,])))
#   l1 = st_sf(st_sfc(l1), crs = st_crs(pop))
#   # routing
#   l1w = as(st_transform(l1, crs = 4326), "Spatial")
#   lp = line2df(l1w)
# #   r = viaroute(lp$fy, lp$fx, lp$ty, lp$tx, profile = "bicycle")
# #   if(grid_num != 188)
# #      rv = viaroute2sldf(r)
#   rv = route_graphhopper(from = c(lp$fx, lp$fy), to = c(lp$tx, lp$ty), vehicle = "bike")
#   rv$grid_num = grid_num
#   rv$ncycle = pop_cents$POBLACION[grid_num] *
#     ck[1]
#   if(grid_num == 1) {
#     r_all = rv
#   } else {
#     r_all = sp::rbind.SpatialLinesDataFrame(r_all, rv)
#   }
# }
# saveRDS(r_all, "data/r_all_graph.Rds")
r_all = readRDS("../data/r_all_graph.Rds")
ncycle_agg = group_by(r_all@data, grid_num) %>% summarise(ncycle = mean(ncycle))
plot(r_all$ncycle)
nrow(ncycle_agg) == nrow(pop_cents)
pop = pop[pop_cents,]
pop_cents$station
pop_cents$Cycling_potential = ncycle_agg$ncycle
pop_cents = readRDS("../data/pop_cents.Rds")
pop$Cycling_potential = ncycle_agg$ncycle

# sort the name crap out
sn = unlist(pop_cents$station)
snn = as.numeric(sn)
snnn = levels(stations$name)[snn]
snnn[is.na(snnn)] = sn[is.na(snnn)]
pop_cents$station = snnn
# saveRDS(pop_cents, "data/pop_cents.Rds")
# find the stations with routes that go outside the study area

Now we can plot the results:

plot(r_all, lwd = r_all$ncycle / mean(r_all$ncycle))

And on an interactive map:

qtm(r_all, lines.lwd = "ncycle", scale = 20) +
  qtm(stations, scale = 2)

Adding values of overlapping lines:

r_all_orig = r_all
r_all_sf = st_as_sf(r_all)
r_all_sf = st_transform(r_all_sf, st_crs(expo))
r_all = r_all_sf %>% filter(st_within(., expo, sparse = FALSE))
r_all_out = r_all_sf %>% filter(!st_within(., expo, sparse = FALSE))
r_all = as(r_all, "Spatial")
r_all_out = as(r_all_out, "Spatial")
rnet = overline(r_all, attrib = "ncycle")
rnet_out = overline(r_all_out, attrib = "ncycle")
pop <- rename(pop, Population = pob_tot1)
qtm(expo) +
tm_shape(pop) +
  tm_fill(col = "Population", n = 4, breaks = c(0, 10, 100, 1000)) +
tm_shape(rnet) +
  tm_lines(lwd = "ncycle", scale = 20) +
  tm_shape(rnet_out) +
  tm_lines(lwd = "ncycle", scale = 5, col = "blue") +
  qtm(stations) +
  tm_shape(expo) +
  tm_borders() +
  tm_scale_bar()
# total cycling potential
sum(pop$Cycling_potential) / sum(pop$POBLACION)

Switch those with overlapping points from San Juan to Expo...

sel = st_within(r_all_sf, expo, sparse = F)
plot(r_all_sf$geometry[sel])
plot(expo, add = T)
plot(r_all_sf$geometry[!sel], add = T, col = "red")
pop_cents$station[!sel & pop_cents$station == "San Juan Alto"] = "Ciudad Expo"

Summary results per station:

plot(pop_cents$geometry)
pop_cents %>% 
  group_by(station) %>% 
  summarise(`Cycling potential` = sum(Cycling_potential),
            Population = sum(POBLACION),
            N_zones = n()) %>% 
  st_set_geometry(NULL) %>% 
  knitr::kable(digits = 0)

readr::write_csv(st_set_geometry(pop_cents, NULL), "/tmp/pop_cells_results.csv")

Saving the results

stations2 = st_sf(name = stations$name, geometry = stations$geometry)
stations2 = st_transform(stations2, 4326)
# write_sf(stations2, "data/stations.geojson")

# write_sf(expo, "data/expo.geojson")

Robinlovelace/pctSeville documentation built on May 9, 2019, 10:30 a.m.