From 0ae2423aaf2340ca8f63052c178142438ea0a3c3 Mon Sep 17 00:00:00 2001 From: StevePem Date: Mon, 22 Jan 2024 17:23:29 +1100 Subject: [PATCH] conforming other scripts to bash/sql refactoring --- NetworkGenerator.R | 74 ++++++++++++------------- functions/addDestinations.R | 97 +++++++++++---------------------- functions/getDestinationTypes.R | 40 +------------- functions/getPTStops.R | 36 ++++++------ functions/writeOutputs.R | 30 +++++----- network.Rproj | 13 +++++ 6 files changed, 115 insertions(+), 175 deletions(-) create mode 100644 network.Rproj diff --git a/NetworkGenerator.R b/NetworkGenerator.R index fe7f9db..37d52ca 100644 --- a/NetworkGenerator.R +++ b/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,14 +175,19 @@ 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)) { @@ -190,6 +195,8 @@ makeNetwork<-function(city, outputFileName = "test"){ 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")) @@ -197,38 +204,24 @@ makeNetwork<-function(city, outputFileName = "test"){ } } - # 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") diff --git a/functions/addDestinations.R b/functions/addDestinations.R index 2b376a2..af6d847 100644 --- a/functions/addDestinations.R +++ b/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 diff --git a/functions/getDestinationTypes.R b/functions/getDestinationTypes.R index 9570e17..e25ee86 100644 --- a/functions/getDestinationTypes.R +++ b/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) -# -# } - diff --git a/functions/getPTStops.R b/functions/getPTStops.R index e6d4a09..a9f6b0a 100644 --- a/functions/getPTStops.R +++ b/functions/getPTStops.R @@ -2,24 +2,23 @@ # requires tidytransit (loaded in NetworkGenerator.R) -getPTStops <- function(city, gtfs_feed, outputCrs, edges_current) { +getPTStops <- function(city, gtfs_feed, outputCrs, region, regionBufferDist) { # city = "Melbourne" - # gtfs_feed = "./data/gtfs.zip" - # outputCrs = 28355 - # edges_current = networkDensified[[2]] + # gtfs_feed = "../data/processed/gtfs.zip" + # outputCrs = 7899 + # region = "../data/processed/greater_melbourne.sqlite" + # regionBufferDist = 10000 # read in GTFS feed gtfs <- read_gtfs(gtfs_feed) %>% gtfs_as_sf(., crs = 4326) - # extract stops with their locations + # extract stops with their locations, filtered to study area + study.area <- st_buffer(st_read(region), regionBufferDist) 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) + st_transform(outputCrs) %>% + st_set_geometry("geom") %>% + st_filter(study.area, predicate = st_intersects) # table of stops and route types stops.routetypes <- gtfs$stop_times %>% @@ -33,9 +32,12 @@ getPTStops <- function(city, gtfs_feed, outputCrs, edges_current) { # 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 - - + if (city == "exception_city") { + # if a city is known to have exceptional route types, build a condition that + # applies where the city is specified and the route types are in the expected + # list - eg city == "Gotham City" & all(route_types %in% c("13", "14", "15")) - + # and provide an appropriate message and function (similar to 'else') + } 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: @@ -48,7 +50,7 @@ PT stops will not be included in destinations.") 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.") +Adjust 'getPTStops' function if these don't match the codes used in your GTFS feed.") stops.routetypes.coded <- stops.routetypes %>% mutate(pt_stop_type = case_when( route_type == 0 ~ "tram", @@ -78,6 +80,4 @@ Check that these match the codes used in your GTFS feed.") return(c()) # empty vector if no stops can be returned } - - -} \ No newline at end of file +} diff --git a/functions/writeOutputs.R b/functions/writeOutputs.R index e82d7e5..e6ecf40 100644 --- a/functions/writeOutputs.R +++ b/functions/writeOutputs.R @@ -181,21 +181,21 @@ exportXML <- function(networkFinal, outputDir){ if(!("id" %in% colnames(links))) links <- links %>% mutate(id=NA) if(!("fwd_slope_pct" %in% colnames(links))) links <- links %>% mutate(fwd_slope_pct=NA, rvs_slope_pct=NA) - # Adding a reverse links for bi-directionals - bi_links <- links %>% - filter(is_oneway==0) %>% - rename(from_id=to_id, to_id=from_id, toX=fromX, toY=fromY, fromX=toX, - fromY=toY, slope=rvs_slope_pct) %>% - mutate(id=NA) %>% - dplyr::select(from_id, to_id, fromX, fromY, toX, toY, length, freespeed, - permlanes, capacity, is_oneway, cycleway, highway, surface, - slope, is_cycle, is_walk, is_car, modes, id) - - links <- rbind( - {links %>% dplyr::select(from_id, to_id, fromX, fromY, toX, toY, length, freespeed, - permlanes, capacity, is_oneway, cycleway, highway, surface, - slope=fwd_slope_pct, is_cycle, is_walk, is_car, modes, id)}, - bi_links) + # # Adding a reverse links for bi-directionals - not required as makeEdgesOneway has made them all oneway + # bi_links <- links %>% + # filter(is_oneway==0) %>% + # rename(from_id=to_id, to_id=from_id, toX=fromX, toY=fromY, fromX=toX, + # fromY=toY, slope=rvs_slope_pct) %>% + # mutate(id=NA) %>% + # dplyr::select(from_id, to_id, fromX, fromY, toX, toY, length, freespeed, + # permlanes, capacity, is_oneway, cycleway, highway, surface, + # slope, is_cycle, is_walk, is_car, modes, id) + # + # links <- rbind( + # {links %>% dplyr::select(from_id, to_id, fromX, fromY, toX, toY, length, freespeed, + # permlanes, capacity, is_oneway, cycleway, highway, surface, + # slope=fwd_slope_pct, is_cycle, is_walk, is_car, modes, id)}, + # bi_links) # Adding bicycle and extra information links <- fncols(links, c("id","osm_id", "highway", "cycleway","slope", diff --git a/network.Rproj b/network.Rproj new file mode 100644 index 0000000..8e3c2eb --- /dev/null +++ b/network.Rproj @@ -0,0 +1,13 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX