Bike accessibility

CyclingImpedances.R

@@ -31,6 +31,8 @@ addImpedances <- function() {
   library(dplyr)
   library(sf)
   library(fs)
+  library(osmextract)
+  library(stringr)
 
   dir_walk(path="./functions/",source, recurse=T, type = "file")
 
@@ -48,7 +50,7 @@ addImpedances <- function() {
                                traffic.links, 
                                traffic.multiplier)
   ## TO DO - maybe traffic can just be joined on link_id?  See whether traffic
-  ## file neatly uses the link_id's from the one-way input
+  ## file neatly uses the link_id's from the one-way input (also in cycling-adoption)
 
   echo("Adding LTS and its impedance\n")
   networkLTS <- addLTS(networkTraffic[[1]], networkTraffic[[2]])
@@ -59,12 +61,17 @@ addImpedances <- function() {
   networkSlope <- addSlopeImped(networkLTS[[1]], networkLTS[[2]])
 
 
+  # Add surface impedance --------------------------------------------------------
+  echo("Adding surface impedance")
+  networkSurf <- addSurfImped(networkSlope[[1]], networkSlope[[2]])
+
+
   # Calculate total weight -----------------------------------------------------
   echo("Calculating cycling weight")
   networkWeighted <- 
-    list(networkSlope[[1]],
-         networkSlope[[2]] %>%
-           mutate(cycle.weight = length + LTS_imped + slope_imped))
+    list(networkSurf[[1]],
+         networkSurf[[2]] %>%
+           mutate(cycle.weight = length + LTS_imped + slope_imped + surf_imped))
 
 
   # write output ---------------------------------------------------------------

NetworkGenerator.R

@@ -6,7 +6,8 @@ makeNetwork<-function(outputFileName="test"){
   # Set city
   city = "Melbourne"
   # city = "Brisbane"
-
+  # city = "Munich"
+
   # City parameters to be set
   # •	outputCrs: desired coordinate system for network
   # •	osmExtract: if 'processOsm=T', OSM extract file in .osm format (.osm.pbf 
@@ -21,9 +22,9 @@ makeNetwork<-function(outputFileName="test"){
   #   in .osm.pbf format
   # •	ndviFile: if 'addNDVI=T', raster file with NDVI values (must be in same
   #   coordinate system as network)
-  # •	gtfs_feed: if 'addGtfs=T', zip file containing GTFS data (must also set
-  #   start and end dates in GTFS section)
-  
+  # •	gtfs_feed: if 'addGtfs=T' or 'addDestinationLayer=T, zip file containing 
+  #   GTFS data (and, if 'addGtfs=T', also set start and end dates in GTFS section)
+
   if (city == "Melbourne") {
     outputCrs = 28355
     osmExtract = "./data/melbourne.osm"
@@ -32,8 +33,8 @@ makeNetwork<-function(outputFileName="test"){
     demFile = "./data/DEM_melbourne.tif"
     osmPbfExtract = "./data/melbourne_australia.osm.pbf"
     ndviFile = "./data/NDVI_1600mBuffer_Melbourne_reprojected.tif"
-    gtfs_feed = "data/gtfs_au_vic_ptv_20191004.zip"
-    
+    gtfs_feed = "./data/gtfs.zip"
+
   } else if (city == "Brisbane") {
     outputCrs = 28356
     osmExtract = ""  # must set 'processOsm=F'
@@ -42,8 +43,17 @@ makeNetwork<-function(outputFileName="test"){
     demFile = "./data/5m_DEM_reprojected.tif" # MIGHT NOT BE FINAL FILE
     osmPbfExtract = "./data/brisbane_australia.osm.pbf"
     ndviFile = ""  # must set 'addNDVI=F'
-    gtfs_feed = ""  # must set 'addGtfs=F'
+    gtfs_feed = "./data/SEQ_GTFS.zip"
 
+  } else if (city == "Munich") {
+    outputCrs = 25832
+    osmExtract = ""  # must set 'processOsm=F'
+    networkSqlite = "./data/munich_network_unconfigured.sqlite"
+    cropAreaPoly = ""  # must set 'crop2TestArea=F'
+    demFile = "" # must set 'addElevation=F'
+    osmPbfExtract = "./data/munich_germany.osm.pbf"
+    ndviFile = ""  # must set 'addNDVI=F'
+    gtfs_feed = "./data/mvv_gtfs.zip" # to test; if not then >> # must set 'addGtfs=F' and 'addDestinationLayer=F'
   }
 
   # INPUT NETWORK 
@@ -278,6 +288,8 @@ makeNetwork<-function(outputFileName="test"){
     destinations <- addDestinations(networkDensified[[1]],
                                     networkDensified[[2]],
                                     osmPbfExtract,
+                                    city,
+                                    gtfs_feed,
                                     outputCrs)
   }
 

functions/addDestinations.R

@@ -7,19 +7,23 @@
 # in 'getDestinationTypes.R'
 
 addDestinations <- function(nodes_current, 
-                             edges_current, 
-                             osmPbfExtract,
-                             outputCrs) {
-
+                            edges_current, 
+                            osmPbfExtract,
+                            city, 
+                            gtfs_feed,
+                            outputCrs) {
+
   # nodes_current = networkDensified[[1]]
   # edges_current = networkDensified[[2]]
   # osmPbfExtract = "./data/melbourne_australia.osm.pbf"
+  # city = "Melbourne"
+  # gtfs_feed = "./data/gtfs.zip"
   # outputCrs = 28355
 
   # # check layers
   # st_layers(osmPbfExtract)
   # # only multipolygons, points and lines are required (not multilinestrings
-  # # or other_relations)
+  # # or other_relations) [and lines not required when using GTFS for PT stops]
 
   # # check keys
   # options(max.print = 2000)
@@ -40,9 +44,9 @@ addDestinations <- function(nodes_current,
   #   some don't exist for some layer types
   extra.tags <- c("access", "amenity", "building", "grades", "healthcare", 
                   "healthcare:speciality", "isced:level", "landuse", "leisure", 
-                  "network", "operator", "operator:type", "public_transport", 
-                  "railway", "school", "shop", "social_facility",  "sport",
-                  "tourism", "train")
+                  "network", "operator", "operator:type", "parking", 
+                  "public_transport", "railway", "school", "shop", 
+                  "social_facility",  "sport", "tourism", "train")
   # oe_vectortranslate(osmPbfExtract, layer = "multipolygons", extra_tags = extra.tags)
   # oe_vectortranslate(osmPbfExtract, layer = "points", extra_tags = extra.tags)
   # oe_vectortranslate(osmPbfExtract, layer = "lines", extra_tags = extra.tags)
@@ -55,8 +59,8 @@ addDestinations <- function(nodes_current,
     st_transform(outputCrs)
   points <- oe_read(osmPbfExtract, layer = "points", extra_tags = extra.tags) %>% 
     st_transform(outputCrs)
-  lines <- oe_read(osmPbfExtract, layer = "lines", extra_tags = extra.tags) %>% 
-    st_transform(outputCrs)
+  # lines <- oe_read(osmPbfExtract, layer = "lines", extra_tags = extra.tags) %>% 
+  #   st_transform(outputCrs)
 
 
   # function to extract specific destination types from point or polygon layers ----
@@ -85,26 +89,43 @@ addDestinations <- function(nodes_current,
         getPost(layer) %>% mutate(dest_type = "post_office"),
         getBank(layer) %>% mutate(dest_type = "bank"),
         getRestaurant(layer) %>% mutate(dest_type = "restaurant"),
-        getCafe(layer) %>% mutate(dest_type = "cafe")
+        getCafe(layer) %>% mutate(dest_type = "cafe"),
+        getParking(layer) %>% mutate(dest_type = "parking")
       ))
   }
 
   # create tables of point and polygon destinations ----
   # ----------------------------------#
   echo("Finding destinations and their nearby nodes\n")
 
-  # create tables for points and polygons, and allocate unique id's (so features 
-  # multiple multiple nodes can be grouped by the id where required)
+  # create tables for points and polygons, allocate unique id's (so features 
+  # multiple multiple nodes can be grouped by the id where required),
+  # and store area and location details
   destination.pt <- 
     bind_rows(destination.layer(points),
-              # add stations (from point, polygons and lines) to point table
-              getStation(points, polygons, lines) %>% 
-                mutate(dest_type = "railway_station")) %>%
-    mutate(dest_id = row_number())
-
+
+              # # add stations (from point, polygons and lines) to point table
+              # getStation(points, polygons, lines) %>% 
+              # mutate(dest_type = "railway_station")) %>%
+
+              # add PT stops (from GTFS feed) to point table
+              getPTStops(city, gtfs_feed, outputCrs, edges_current) %>%
+                mutate(dest_type = "pt_stop")) %>%
+
+    mutate(dest_id = row_number(),
+           area_m2 = 0,
+           centroid_x = st_coordinates(.)[, 1],
+           centroid_y = st_coordinates(.)[, 2])
+
   destination.poly <- 
     destination.layer(polygons) %>%
-    mutate(dest_id = max(destination.pt$dest_id) + row_number())
+    mutate(dest_id = max(destination.pt$dest_id) + row_number(),
+           area_m2 = as.numeric(st_area(.)),
+           centroid_x = st_coordinates(st_centroid(.))[, 1],
+           centroid_y = st_coordinates(st_centroid(.))[, 2])
+
+  # Remove any invalid polygons as they may cause errors
+  destination.poly <- destination.poly[which(st_is_valid(destination.poly$geometry)), ]
 
   # Remove any invalid polygons as they may cause errors
   destination.poly <- destination.poly[which(st_is_valid(destination.poly$geometry)), ]

functions/addNDVI.R

@@ -3,7 +3,7 @@
 
 addNDVI2Links <- function(links, ndviFile, ndviBuffDist) {
 
-  # links = networkRestructured[[2]]
+  # links = networkDensified[[2]]
   # ndviFile = "./data/NDVI_1600mBuffer_Melbourne_reprojected.tif"
   # ndviBuffDist = 30
 
@@ -28,12 +28,22 @@ addNDVI2Links <- function(links, ndviFile, ndviBuffDist) {
     group_by(ID) %>%
     summarise(ndvi = mean(NDVI, na.rm = TRUE))
 
+  # find the mean AND OTHER VALUES of the values for each link
+  ndvi_values_mean <- ndvi_values %>%
+    group_by(ID) %>%
+    summarise(ndvi = mean(NDVI, na.rm = TRUE),
+              ndvi_md = median(NDVI, na.rm = TRUE),
+              ndvi_75 = quantile(NDVI, na.rm = TRUE, probs = 0.75),
+              ndvi_90 = quantile(NDVI, na.rm = TRUE, probs = 0.9))
+
   # join to the links, using the row number and ID
   links.with.ndvi <- links %>%
     mutate(row_no = row_number()) %>%
     left_join(., ndvi_values_mean, by = c("row_no" = "ID")) %>%
     dplyr::select(-row_no)
 
+  # st_write(links.with.ndvi, "./SP_working/links_with_NDVI.sqlite", delete_layer = TRUE)
+
   return(links.with.ndvi)
 
 }

functions/addSurfImped.R

@@ -0,0 +1,71 @@
+# function to add impedance to network links, based on surface
+addSurfImped <- function(nodes_current, edges_current) {
+
+  edges_current <- edges_current %>%
+
+    # add surface level, based on OSM surface categories
+    # [CONSIDER REVISING THIS LIST WHEN USING IN OTHER CITIES - SEE 
+    #  'checkSurfTypes()' BELOW FOR A FUNCTION TO HELP INVESTIGATE]
+    mutate(surf_lvl = case_when(
+      str_detect(tolower(surface), 
+                 "bluestone|boardwalk|brick|cobblestone|composite|fibre|metal|paving|stone|plastic|sett|stone|tile|timber|wood")             
+      ~ "rough",
+      str_detect(tolower(surface), 
+                 "clay|compacted|dirt|earth|grass|gravel|mud|rock|sand|turf|unpacved|unpaved|unsealed")
+      ~ "unpaved",
+      TRUE ~ "smooth"
+    )) %>%
+
+    # add impedance from surface level [ILLUSTRATIVE IMPEDANCES ONLY, TO BE 
+    # REVISED BASED ON FURTHER RESEARCH OR SURVEY RESULTS]
+    mutate(surf_imped = case_when(
+      surf_lvl == "smooth"  ~ 0,
+      surf_lvl == "unpaved" ~ length * 0.1,
+      surf_lvl == "rough"   ~ length * 0.2
+    ))
+
+  return(list(nodes_current, edges_current))
+}
+
+
+# investigation function to check number and type of surfaces
+checkSurfTypes <- function(osmPbfExtract) {
+
+  # osmPbfExtract = "./data/melbourne_australia.osm.pbf"
+
+  echo("Reading in the .osm.pbf extract line layer\n")
+  lines <- oe_read(osmPbfExtract, layer = "lines", extra_tags = "surface")
+
+  highways <- lines %>%
+    filter (highway %in% 
+              c("motorway", "motorway_link", "trunk", "trunk_link",
+                "primary", "primary_link", "secondary", "secondary_link",
+                "tertiary", "tertiary_link", "residential", "road", "unclassified",
+                "living_street", "cycleway", "track", "service", 
+                "pedestrian", "footway", "path", "corridor", "steps"))
+
+  surface.totals <- highways %>%
+    st_drop_geometry() %>%
+    group_by(surface) %>%
+    summarise(n = n()) %>%
+    ungroup()
+
+  # # print outputs as comma-separated string
+  # surface.totals.output <- ""
+  # for (i in 1:nrow(surface.totals)) {
+  #   line <- paste0(surface.totals$surface[i], " (", surface.totals$n[i], ")")
+  #   surface.totals.output <- paste0(surface.totals.output, line)
+  #   if (i < nrow(surface.totals)) {
+  #     surface.totals.output <- paste0(surface.totals.output, ", ")
+  #   }
+  # }
+  # print(surface.totals.output)
+
+  return(surface.totals)
+
+}
+
+# # using checkSurfTypes to compile table of surface types, no. and surf_lvl
+# check <- checkSurfTypes("./data/melbourne_australia.osm.pbf")
+# check.totals <- addSurfImped(NA, check %>% mutate(length = 0))[[2]] %>%
+#   dplyr::select(surface, n, surf_lvl)

functions/getDestinationTypes.R

@@ -125,35 +125,44 @@ getCafe <- function(layer) {
   return(layer %>% filter(amenity == "cafe"))
 }
 
-# 8 railway stations ----
-# Returns list of  stations as points
-# Note the buffer distance of 100m below; closest railway stations in Melbourne are
-# Riversdale & Willison (about 420m)
-getStation <- function(points, polygons, lines) {
-  # general filter to find station objects
-  filterStation <- function(layer) {
-    return(layer %>%
-             filter((public_transport == "station" | public_transport == "stop_position") & 
-                      (railway == "station" | railway == "stop" | train == "yes" | 
-                         grepl("train", tolower(network)) | grepl("train", tolower(operator))) &
-                      (is.na(tourism) |  tourism != "yes") &
-                      (is.na(railway) | railway != "construction")))
-  }
-
-  # find each object, and buffer to 100m
-  buff.dist <- 100
-  station.pt <- filterStation(points) %>% st_buffer(buff.dist)
-  station.poly <- filterStation(polygons) %>% st_buffer(buff.dist)
-  station.line <- filterStation(lines) %>% st_buffer(buff.dist)
-
-  # dissolve, then separate to individual polygons
-  stations <- bind_rows(station.pt, station.poly, station.line) %>%
-    st_union() %>%
-    st_as_sf() %>%
-    st_cast("POLYGON") %>%
-    st_centroid() %>%
-    # label geometry column
-    rename(geometry = x)
-
-}
+# 8 parking ----
+getParking <- function(layer) {
+  return(layer %>% filter(amenity == "parking" &
+                            !access %in% c("no", "private")))
+}
+
+
+# 9 railway stations [not used] ----
+# See getPtStops.R instead
+
+# # Returns list of  stations as points
+# # Note the buffer distance of 100m below; closest railway stations in Melbourne are
+# # Riversdale & Willison (about 420m)
+# getStation <- function(points, polygons, lines) {
+#   # general filter to find station objects
+#   filterStation <- function(layer) {
+#     return(layer %>%
+#              filter((public_transport == "station" | public_transport == "stop_position") & 
+#                       (railway == "station" | railway == "stop" | train == "yes" | 
+#                          grepl("train", tolower(network)) | grepl("train", tolower(operator))) &
+#                       (is.na(tourism) |  tourism != "yes") &
+#                       (is.na(railway) | railway != "construction")))
+#   }
+#   
+#   # find each object, and buffer to 100m
+#   buff.dist <- 100
+#   station.pt <- filterStation(points) %>% st_buffer(buff.dist)
+#   station.poly <- filterStation(polygons) %>% st_buffer(buff.dist)
+#   station.line <- filterStation(lines) %>% st_buffer(buff.dist)
+#   
+#   # dissolve, then separate to individual polygons
+#   stations <- bind_rows(station.pt, station.poly, station.line) %>%
+#     st_union() %>%
+#     st_as_sf() %>%
+#     st_cast("POLYGON") %>%
+#     st_centroid() %>%
+#     # label geometry column
+#     rename(geometry = x)
+#   
+# }