Destinations - PT from GTFS and other improvements

NetworkGenerator.R

@@ -6,6 +6,7 @@ makeNetwork<-function(outputFileName="test"){
   # Set city
   city = "Melbourne"
   # city = "Brisbane"
+  # city = "Munich"
 
   # City parameters to be set
   # •	outputCrs: desired coordinate system for network
@@ -21,8 +22,8 @@ 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
@@ -32,7 +33,7 @@ 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
@@ -42,8 +43,18 @@ 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 +289,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)
@@ -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 ----
@@ -99,14 +103,20 @@ addDestinations <- function(nodes_current,
   # 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")) %>%
+
+              # # 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(),

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,11 +28,21 @@ 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/getDestinationTypes.R

@@ -133,35 +133,37 @@ getParking <- function(layer) {
 }
 
 
-# 9 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)
-
-}
+# 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)
+#   
+# }
 

functions/getPTStops.R

@@ -0,0 +1,83 @@
+# read in PT stops from GTFS
+
+# requires tidytransit (loaded in NetworkGenerator.R)
+
+getPTStops <- function(city, gtfs_feed, outputCrs, nodes_current) {
+  # city = "Melbourne"
+  # gtfs_feed = "./data/gtfs.zip"
+  # outputCrs = 28355
+  # edges_current = networkDensified[[2]]
+
+  # read in GTFS feed
+  gtfs <- read_gtfs(gtfs_feed) %>%
+    gtfs_as_sf(., crs = 4326)
+
+  # extract stops with their locations
+  stops <- gtfs$stops %>%
+    st_transform(outputCrs)
+
+  # limit to stops within the study area (convex hull of edges)
+  stops <- stops %>%
+    st_filter(., st_convex_hull(st_union(edges_current)), 
+              predicate = st_intersects)
+
+  # table of stops and route types
+  stops.routetypes <- gtfs$stop_times %>%
+    left_join(gtfs$trips, by = "trip_id") %>%
+    left_join(gtfs$routes, by = "route_id") %>%
+
+    # keep only distinct stop_id and route_type combinations
+    dplyr::select(stop_id, route_type) %>%
+    distinct()
+
+  # apply route types
+  route_types = stops.routetypes$route_type %>% unique() %>% sort()
+
+  if (city == "ProvisionForMunich") { # test should be city is Munich AND stops are in the expected list
+
+
+  } else if (!all(route_types %in% c("0", "1", "2", "3", "4", "5", "6", "7", "11", "12"))) {
+    message("GTFS Feed contains the following route type codes: ", paste(route_types, collapse = ", "), ". Unable to process these using 
+the standard route codes from https://developers.google.com/transit/gtfs/reference, which are:
+   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12 monorail. 
+Edit getPTStops.R to specify the meanings of the codes used in the GTFS Feed.
+PT stops will not be included in destinations.")
+    stops.found = FALSE
+
+  } else {
+    message("GTFS Feed contains the following route type codes: ", paste(route_types, collapse = ", "), ".
+Using standard route_type codes from https://developers.google.com/transit/gtfs/reference:
+   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12 monorail. 
+Check that these match the codes used in your GTFS feed.")
+    stops.routetypes.coded <- stops.routetypes %>%
+      mutate(pt_stop_type = case_when(
+        route_type == 0  ~ "tram",
+        route_type == 1  ~ "metro",
+        route_type == 2  ~ "train",
+        route_type == 3  ~ "bus",
+        route_type == 4  ~ "ferry",
+        route_type == 5  ~ "cable tram",
+        route_type == 6  ~ "cable car",
+        route_type == 7  ~ "funicular",
+        route_type == 11 ~ "trolleybus",
+        route_type == 12 ~ "monorail"
+      )) %>%
+      dplyr::select(stop_id, pt_stop_type)
+
+    stops.found = TRUE
+
+  }
+
+  if(stops.found) {
+    stops.with.types <- stops %>%
+      left_join(stops.routetypes.coded, by = "stop_id")
+
+    return(stops.with.types)
+
+  } else {
+    return(c())  # empty vector if no stops can be returned
+  }
+
+
+
+}