conforming other scripts to bash/sql refactoring

NetworkGenerator.R

@@ -33,22 +33,22 @@ makeNetwork<-function(city, outputFileName = "test"){
     outputCrs = 7899
     osmGpkg = "../data/processed/bendigo_osm.gpkg"
     unconfiguredSqlite = "../data/processed/bendigo_network_unconfigured.sqlite"
-    # cropAreaPoly = ""  # must set 'crop2Area=F'
+    cropAreaPoly = ""  # must set 'crop2Area=F'
     # demFile = "./data/5m_DEM_reprojected.tif" # MIGHT NOT BE FINAL FILE
     # osmPbfExtract = "./data/brisbane_australia.osm.pbf"
     # ndviFile = ""  # must set 'addNDVI=F'
-    # gtfs_feed = "./data/SEQ_GTFS.zip"
+    gtfs_feed = "../data/processed/gtfs.zip"
 
   } else if (city == "Melbourne") {
     region = "../data/processed/greater_melbourne.sqlite"
     outputCrs = 7899
     osmGpkg = "../data/processed/melbourne_osm.gpkg"
     unconfiguredSqlite = "../data/processed/melbourne_network_unconfigured.sqlite"
-    # cropAreaPoly = "city-of-melbourne_victoria"
+    cropAreaPoly = "city-of-melbourne_victoria"
     # demFile = "./data/DEM_melbourne.tif"
     # osmPbfExtract = "./data/melbourne_australia.osm.pbf"
     # ndviFile = "./data/NDVI_1600mBuffer_Melbourne_reprojected.tif"
-    # gtfs_feed = "./data/gtfs.zip"
+    gtfs_feed = "../data/processed/gtfs.zip"
 
   } else {
     echo(paste("City parameters for", city, "have not been set; unable to proceed\n"))
@@ -65,7 +65,7 @@ makeNetwork<-function(city, outputFileName = "test"){
   # NETWORK FROM OSM 
   # A flag for whether to build unconfigured network from osm extract (if not,
   # must already have unconfigured sqlite)
-  networkFromOsm=T
+  networkFromOsm=F
   saveUnconfigured=T
 
   # SIMPLIFICATION
@@ -74,7 +74,7 @@ makeNetwork<-function(city, outputFileName = "test"){
   crop2Area=F
 
   # DENSIFICATION
-  desnificationMaxLengh=500
+  densificationMaxLength=500
   densifyBikeways=T
 
   # CAPACITY ADJUSTMENT
@@ -84,16 +84,16 @@ makeNetwork<-function(city, outputFileName = "test"){
 
   # ELEVATION
   # A flag for whether to add elevation or not
-  addElevation=T
+  addElevation=F
   ElevationMultiplier=1
 
   # DESTINATIONS
-  # A flag for whether to add a destinations layer (drawn from OSM) or not
+  # A flag for whether to add a destinations layer (drawn from OSM, and GTFS for PT) or not
   addDestinationLayer=T
 
   # NDVI
   # A flag for whether to add NDVI or not
-  addNDVI=T
+  addNDVI=F
   # Buffer distance for finding average NDVI for links
   ndviBuffDist=30
 
@@ -175,60 +175,53 @@ makeNetwork<-function(city, outputFileName = "test"){
   # Processing OSM, or loading existing layers if not required
   if(networkFromOsm) {
     echo(paste0("Starting to process osm extract file, ", osmGpkg, "\n"))
-    networkUnconfigured <- processOsm(osmGpkg, outputCrs)
+    networkUnconfiguredOutputs <- processOsm(osmGpkg, outputCrs)
 
     if (saveUnconfigured) {
-      st_write(networkUnconfigured[[1]], unconfiguredSqlite, layer = "nodes", delete_layer = T)
-      st_write(networkUnconfigured[[2]], unconfiguredSqlite, layer = "edges", delete_layer = T)
-      st_write(networkUnconfigured[[3]], unconfiguredSqlite, layer = "osm_metadata", delete_layer = T)
+      if (file_exists(unconfiguredSqlite)) st_delete(unconfiguredSqlite)
+      st_write(networkUnconfiguredOutputs[[1]], unconfiguredSqlite, layer = "nodes")
+      st_write(networkUnconfiguredOutputs[[2]], unconfiguredSqlite, layer = "edges")
+      st_write(networkUnconfiguredOutputs[[3]], unconfiguredSqlite, layer = "osm_metadata")
     }
 
+    networkUnconfigured <- list(networkUnconfiguredOutputs[[1]],
+                                networkUnconfiguredOutputs[[2]])
+    osm_metadata <- networkUnconfiguredOutputs[[3]]
+
   } else {
 
     if (file_exists(unconfiguredSqlite)) {
       echo(paste("Reading in existing unconfigured network,", unconfiguredSqlite, "\n"))
       networkUnconfigured <- 
         list(st_read(unconfiguredSqlite, layer = "nodes") %>% st_set_geometry("geom"),
              st_read(unconfiguredSqlite, layer = "edges") %>% st_set_geometry("geom"))
+      osm_metadata <- st_read(unconfiguredSqlite, layer = "osm_metadata") %>%
+        filter(osm_id %in% networkUnconfigured[[2]]$osm_id)
 
     } else {
       echo(paste("Unconfigured network file", unconfiguredSqlite, "not found; unable to proceed\n"))
       return()
     }
   }
 
-  # TO DO FROM HERE, and fix up 'networkInput' (now 'networkUnconfigured', and probably 
-  # read in osm_metadata above, after 'networkUnconfigured <-')
   # crop to test area if required
-  if(crop2Area)system.time(networkInput <- crop2Poly(networkInput,
-                                                     cropAreaPoly,
-                                                     outputCrs))
+  if(crop2Area)system.time(networkUnconfigured <- crop2Poly(networkUnconfigured,
+                                                            cropAreaPoly,
+                                                            outputCrs))
+  # process OSM metadata
   echo("processing OSM meta data\n")
-  osm_metadata <- st_read(unconfiguredSqlite,layer="osm_metadata",quiet=T) %>%
-    filter(osm_id%in%networkInput[[2]]$osm_id)
   echo("Building default OSM attribute tables\n")
   defaults_df <- buildDefaultsDF()
   highway_lookup <- defaults_df %>% dplyr::select(highway, highway_order)
   echo("Processing OSM tags and joining with defaults\n")
   system.time( osmAttributes <- processOsmTags(osm_metadata,defaults_df))
 
-  edgesAttributed <- networkInput[[2]] %>%
+  edgesAttributed <- networkUnconfigured[[2]] %>%
     inner_join(osmAttributes, by="osm_id") %>%
-    # dplyr::select(-osm_id,highway,highway_order)
-    dplyr::select(-highway,highway_order)
-
-  cat(paste0("edgesAttributed:\n"))
-  str(edgesAttributed)
-  cat(paste0("\n"))
+    dplyr::select(-highway, highway_order)
 
   # keep only the largest connected component
-  largestComponent <- largestConnectedComponent(networkInput[[1]],edgesAttributed)
-
-  cat(paste0("largestComponent, nodes:\n"))
-  str(largestComponent[[1]])
-  cat(paste0("\nlargestComponent, edges:\n"))
-  str(largestComponent[[2]])
-  cat(paste0("\n"))
+  largestComponent <- largestConnectedComponent(networkUnconfigured[[1]], edgesAttributed)
 
   # simplify intersections while preserving attributes and original geometry.
   system.time(intersectionsSimplified <- simplifyIntersections(largestComponent[[1]],
@@ -287,8 +280,8 @@ makeNetwork<-function(city, outputFileName = "test"){
   networkConnected <- largestNetworkSubgraph(networkNonDisconnected,'walk')
 
   # densify the network so that no residential streets are longer than 500m
-  if (addElevation==T & densifyBikeways==F) message("Consider changing densifyBikeways to true when addElevation is true to ge a more accurate slope esimation for bikeways")
-  networkDensified <- densifyNetwork(networkConnected,desnificationMaxLengh,
+  if (addElevation==T & densifyBikeways==F) message("Consider changing densifyBikeways to true when addElevation is true to get a more accurate slope esimation for bikeways")
+  networkDensified <- densifyNetwork(networkConnected,densificationMaxLength,
                                      densifyBikeways)
 
   # Adding NDVI to links
@@ -302,10 +295,12 @@ makeNetwork<-function(city, outputFileName = "test"){
   if (addDestinationLayer) {
     destinations <- addDestinations(networkDensified[[1]],
                                     networkDensified[[2]],
-                                    osmPbfExtract,
+                                    osmGpkg,
                                     city,
                                     gtfs_feed,
-                                    outputCrs)
+                                    outputCrs,
+                                    region,
+                                    regionBufferDist)
   }
 
   # simplify geometry so all edges are straight lines
@@ -381,4 +376,5 @@ makeNetwork<-function(city, outputFileName = "test"){
 }
 
 ## JUST FOR TESTING
-output <- makeNetwork(city = "Bendigo")
+makeNetwork(city = "Bendigo")
+makeNetwork(city = "Melbourne")

functions/addDestinations.R

@@ -1,76 +1,53 @@
 # function to create a destination layer to add to output network
 
-# assumes input file (OSMextract) is in .osm.pbf format, for example,
-# as downloaded from  https://www.interline.io/osm/extracts/ 
-
 # uses functions for various destination types with tag combinations set out
 # in 'getDestinationTypes.R'
 
-# NOTE - WILL REQUIRE REFACTORING IF THE .GPKG, RATHER THAN THE .OSM.PBF, IS
-# SAVED AS THE BASE FILE - NEED TO EXTRACT THE OTHER TAGS FROM THE .GPKG
-
 addDestinations <- function(nodes_current, 
                             edges_current, 
-                            osmPbfExtract,
+                            osmGpkg,
                             city, 
                             gtfs_feed,
-                            outputCrs) {
-
+                            outputCrs,
+                            region,
+                            regionBufferDist) {
+
   # nodes_current = networkDensified[[1]]
   # edges_current = networkDensified[[2]]
-  # osmPbfExtract = "./data/melbourne_australia.osm.pbf"
+  # osmGpkg = "../data/processed/melbourne_osm.gpkg"
   # city = "Melbourne"
-  # gtfs_feed = "./data/gtfs.zip"
+  # gtfs_feed = "../data/processed/gtfs.zip"
   # outputCrs = 28355
 
-  # # check layers
-  # st_layers(osmPbfExtract)
-  # # only multipolygons, points and lines are required (not multilinestrings
-  # # or other_relations) [and lines not required when using GTFS for PT stops]
-
   # # check keys
   # options(max.print = 2000)
-  # polygon.tags <- oe_get_keys(osmPbfExtract, layer = "multipolygons") %>% sort()
-  # point.tags <- oe_get_keys(osmPbfExtract, layer = "points") %>% sort()
-  # line.tags <- oe_get_keys(osmPbfExtract, layer = "lines") %>% sort()
-
+  # point.tags <- oe_get_keys(osmGpkg, layer = "points") %>% sort()
+  # polygon.tags <- oe_get_keys(osmGpkg, layer = "multipolygons") %>% sort()
+
   # reading layers ----
   # ----------------------------------#
-  echo("Reading in the .osm.pbf extract layers\n")
+  echo("Reading in the OSM layers\n")
 
-  # create gpkg file in same directory as osmPbfExtract, using the 'extra_tags'
-  # Note:
-  # - the gpkg does not need to be retained permanently, but its creation is part
-  #   of the process of reading the layers; if already created, the reading 
-  #   process will be quicker)
-  # - for simplicity, the same extra tags are added for all layers, though
-  #   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", "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)
-  # 
-  # # read in the .gpkg file (same directory and name as .osm.pbf file, but .gpkg extension)
-  # gpkg <- paste0(path_dir(osmPbfExtract), "/", 
-  #                gsub(".osm.pbf", ".gpkg", path_file(osmPbfExtract)))
-  # read in the layers
-  polygons <- oe_read(osmPbfExtract, layer = "multipolygons", extra_tags = extra.tags) %>% 
-    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)
+  extra.tag.string <- "SELECT *, 
+                      hstore_get_value(other_tags, 'access') AS access,
+                      hstore_get_value(other_tags, 'amenity') AS amenity,
+                      hstore_get_value(other_tags, 'grades') AS grades,
+                      hstore_get_value(other_tags, 'healthcare') AS healthcare,
+                      hstore_get_value(other_tags, 'isced:level') AS isced_level,
+                      hstore_get_value(other_tags, 'leisure') AS leisure,
+                      hstore_get_value(other_tags, 'parking') AS parking,
+                      hstore_get_value(other_tags, 'school') AS school,
+                      hstore_get_value(other_tags, 'shop') AS shop,
+                      hstore_get_value(other_tags, 'sport') AS sport"
 
+  # read in the layers
+  points <- oe_read(osmGpkg, query = paste(extra.tag.string, "FROM points"), quiet = TRUE)
+  polygons <- oe_read(osmGpkg, query = paste(extra.tag.string, "FROM multipolygons"), quiet = TRUE)
 
   # function to extract specific destination types from point or polygon layers ----
   # ----------------------------------#
   # all the tag combination functions in 'getDestinationTypes.R' apply to both
-  # points and polygons, except 'railway station', which are a combination of
-  # point, polygon and line features
+  # points and polygons
 
   destination.layer <- function(layer) {
     return(
@@ -106,44 +83,36 @@ 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 PT stops (from GTFS feed) to point table
-              getPTStops(city, gtfs_feed, outputCrs, edges_current) %>%
+              getPTStops(city, gtfs_feed, outputCrs, region, regionBufferDist) %>%
                 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) %>%
+    filter(st_is_valid(geom)) %>%
     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)), ]
-
+
   # # check numbers of each destination type
   # chk <- full_join(destination.poly %>%
   #                    st_drop_geometry() %>%
   #                    group_by(dest_type) %>%
-  #                    summarise(poly = n()), 
+  #                    summarise(poly = n()),
   #                  destination.pt %>%
   #                    st_drop_geometry() %>%
   #                    group_by(dest_type) %>%
-  #                    summarise(pt = n()), 
+  #                    summarise(pt = n()),
   #                  by = "dest_type")
-  
+
 
   # find relevant nodes ----
   # For all destinations except parks and schools ('small features'), relevant 

functions/getDestinationTypes.R

@@ -1,8 +1,6 @@
 # functions to locate specific types of destinations
 
-## All tag combinations below can be applied to both points and polygons, except
-## railway stations which are a combination of points, polygons and lines, and
-## require aggregation within a boundary distance
+## All tag combinations below can be applied to both points and polygons
 
 # 1 open space ----
 getPlayground <- function(layer) {
@@ -130,39 +128,3 @@ 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)
-#   
-# }
-