diff --git a/NetworkGenerator.R b/NetworkGenerator.R index 2a52e95..fe7f9db 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -1,17 +1,21 @@ -makeNetwork<-function(outputFileName="test"){ - # outputFileName="network" +makeNetwork<-function(city, outputFileName = "test"){ + # city = "Bendigo" + # city = "Melbourne" + + # outputFileName = "network" + # Parameters -------------------------------------------------------------- - # CITY AND ITS PARAMETERS - # Set city - city = "Bendigo" - # city = "Melbourne" - + # CITY PARAMETERS # City parameters to be set + # • region: if 'downloadOsm=T', file delineating the boundary of the area for + # which Osm extract is to be downloaded (assumed to be in sqlite format + # with a single layer) # • outputCrs: desired coordinate system for network - # • osmExtract: if 'processOsm=T', OSM extract file in .osm format (.osm.pbf - # not supported for this step) - # • networkSqlite: if 'processOsm=F', network sqlite file + # • osmGpkg: location where downloaded OSM extract for region is to be stored + # (if 'downloadOsm=T') and/or read from (if 'processOsm=T') + # • unconfiguredSqlite: location where processed OSM file is to be stored + # (if 'processOsm=T') or read from (if 'processOsm=F') # • cropAreaPoly: if 'crop2TestArea=T' cropArea location from # https://github.com/JamesChevalier/cities/tree/master/australia/victoria # (only supported for Victoria at this stage) @@ -28,8 +32,8 @@ makeNetwork<-function(outputFileName="test"){ region = "../data/processed/greater_bendigo.sqlite" outputCrs = 7899 osmGpkg = "../data/processed/bendigo_osm.gpkg" - # networkSqlite = "./data/brisbane_network_unconfigured.sqlite" - # cropAreaPoly = "" # must set 'crop2TestArea=F' + unconfiguredSqlite = "../data/processed/bendigo_network_unconfigured.sqlite" + # 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' @@ -39,13 +43,16 @@ makeNetwork<-function(outputFileName="test"){ region = "../data/processed/greater_melbourne.sqlite" outputCrs = 7899 osmGpkg = "../data/processed/melbourne_osm.gpkg" - # networkSqlite = "./data/melbourne_network_unconfigured.sqlite" + unconfiguredSqlite = "../data/processed/melbourne_network_unconfigured.sqlite" # 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" + } else { + echo(paste("City parameters for", city, "have not been set; unable to proceed\n")) + return() } # DOWNLOAD OSM EXTRACT @@ -55,11 +62,11 @@ makeNetwork<-function(outputFileName="test"){ regionBufferDist=10000 # Distance to buffer region when getting osm extract retainDownload=F # Whether to retain downloaded file after region extracted - # NETWORK FROM OSM # A flag for whether to build unconfigured network from osm extract (if not, - # must already have network sqlite) + # must already have unconfigured sqlite) networkFromOsm=T + saveUnconfigured=T # SIMPLIFICATION shortLinkLength=20 @@ -109,6 +116,7 @@ makeNetwork<-function(outputFileName="test"){ library(sf) library(fs) library(dplyr) + library(tidyr) library(data.table) library(stringr) library(igraph) @@ -129,12 +137,11 @@ makeNetwork<-function(outputFileName="test"){ # Building the output folder structure ------------------------------------ - ## COMMENTING THIS OUT FOR NOW BECAUSE IT'S ANNOYING; ADD BACK LATER - # outputDir <- paste0("output/",outputFileName) - # if(dir.exists(outputDir)) dir_delete(outputDir) - # dir_create(paste0('./',outputDir)) - # sink(paste0('./',outputDir,'/makeMatsimNetwork.log'), append=FALSE, split=TRUE) - # if (addGtfs) dir_create(paste0(outputDir,"/gtfs")) + outputDir <- paste0("output/",outputFileName) + if(dir.exists(outputDir)) dir_delete(outputDir) + dir_create(paste0('./',outputDir)) + sink(paste0('./',outputDir,'/makeMatsimNetwork.log'), append=FALSE, split=TRUE) + if (addGtfs) dir_create(paste0(outputDir,"/gtfs")) # Functions -------------------------------------------------------------- @@ -165,47 +172,39 @@ makeNetwork<-function(outputFileName="test"){ getOsmExtract(region, outputCrs, regionBufferDist, osmGpkg) } - # Processing OSM + # Processing OSM, or loading existing layers if not required if(networkFromOsm) { echo(paste0("Starting to process osm extract file, ", osmGpkg, "\n")) - # networkSqlite="./data/network.sqlite" - # if(file_exists(osmExtract)){ - # system(paste("./processOSM.sh ", osmExtract, outputCrs, networkSqlite)) - # }else{ - # warning("OSM extract not found, skipping this step") - # } - # FINALISE THIS CODE BLOCK DEPENDING ON FINAL STRUCTURE OF processOsm.R - # AND ARRANGMEENTS FOR SAVING - networkSqlite <- processOsm(osmGpkg, outputCrs) + networkUnconfigured <- 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) + } + + } 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")) + + } else { + echo(paste("Unconfigured network file", unconfiguredSqlite, "not found; unable to proceed\n")) + return() + } } - return(networkSqlite) ## JUST FOR TESTING - DELETE! - - # Note: writing logical fields to sqlite is a bad idea, so switching to integers - networkInput <- list(st_read(networkSqlite,layer="nodes",quiet=T), - st_read(networkSqlite,layer="edges",quiet=T)) - - # We run into trouble if the geometry column is 'geom' instead of 'GEOMETRY' - if('GEOMETRY'%in%colnames(networkInput[[1]])) { - networkInput[[1]]<-networkInput[[1]]%>%rename(geom=GEOMETRY) - } - if('GEOMETRY'%in%colnames(networkInput[[2]])) { - networkInput[[2]]<-networkInput[[2]]%>%rename(geom=GEOMETRY) - } - - cat(paste0("Network input, nodes:\n")) - str(networkInput[[1]]) - # print.data.frame(head(networkInput[[1]])) - cat(paste0("\nNetwork input, edges:\n")) - str(networkInput[[2]]) - cat(paste0("\n")) - + # 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)) echo("processing OSM meta data\n") - osm_metadata <- st_read(networkSqlite,layer="osm_metadata",quiet=T) %>% + 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() @@ -382,4 +381,4 @@ makeNetwork<-function(outputFileName="test"){ } ## JUST FOR TESTING -# output <- makeNetwork() \ No newline at end of file +output <- makeNetwork(city = "Bendigo") diff --git a/functions/processOsm.R b/functions/processOsm.R index ff41a1d..b1ea997 100644 --- a/functions/processOsm.R +++ b/functions/processOsm.R @@ -1,24 +1,10 @@ # function to convert OSM .gpkg file into network of nodes and edges - -library(sf) -library(fs) -library(dplyr) -library(lwgeom) -library(nngeo) -library(stringr) -library(doSNOW) -library(parallel) -library(foreach) -library(ggplot2) -source("./functions/etc/logging.R") -source("./functions/splitPathsAtPoints.R") -outputCrs = 7899 - processOsm <- function(osmGpkg, outputCrs) { - osmGpkg = "../data/processed/bendigo_osm.gpkg" + # osmGpkg = "../data/processed/bendigo_osm.gpkg" # osmGpkg = "../data/processed/melbourne_osm.gpkg" + # outputCrs = 7899 # read in OSM data # -----------------------------------# @@ -87,7 +73,6 @@ processOsm <- function(osmGpkg, outputCrs) { ungroup() ) - # temp dev notes (SP): # (1) compared to network.sql, this only places intersections where both are # bridges, both are tunnels, or both are neither (whereas network.sql, @@ -105,7 +90,7 @@ processOsm <- function(osmGpkg, outputCrs) { split.path.list <- splitPathsAtPoints(paths, intersections, 0.001, "osm_id") # convert to dataframe, snap to grid, remove empty geometries, add unique id - echo("Combining the split paths into a single dataframe") + echo("Combining the split paths into a single dataframe\n") system.time( split.paths <- bind_rows(split.path.list) %>% st_snap_to_grid(1) %>% @@ -116,19 +101,20 @@ processOsm <- function(osmGpkg, outputCrs) { ) # find endpoints of each split path - echo("Finding endpoints of the split paths") + echo("Finding endpoints of the split paths\n") endpoints <- rbind(lwgeom::st_startpoint(split.paths) %>% st_sf(), lwgeom::st_endpoint(split.paths) %>% st_sf()) %>% # remove duplicates (produces multipoint) summarise() %>% # convert multipoint to point st_cast("POINT") %>% - # add unique id - mutate(endpoint_id = row_number()) + # add unique id and set geometry column name + mutate(endpoint_id = row_number()) %>% + st_set_geometry("geom") # find split.paths that have more than 2 endpoints within 0.1m, in order # to re-split at ajdacent endpoints - echo("Finding paths that need to be re-split") + echo("Finding paths that need to be re-split\n") paths.with.nearby.endpoints <- split.paths %>% # joint endpoints within 0.1m st_join(endpoints %>% st_buffer(0.1), join = st_intersects) @@ -150,7 +136,7 @@ processOsm <- function(osmGpkg, outputCrs) { # do a second round of splitting: re-split the paths that have adjacent endpoints, # using 0.1 distance this time, but only where adjacent endpoint is an endpoint # for a path that has the same bridge_tunnel status as the path to be resplit - echo(paste("Re-splitting", nrow(paths.to.resplit), "paths at adjacent endpoints")) + echo(paste("Re-splitting", nrow(paths.to.resplit), "paths at adjacent endpoints\n")) endpoints.for.resplit <- paths.with.nearby.endpoints %>% # just keep paths that need to be resplit, with their bridge_tunnel status @@ -174,7 +160,7 @@ processOsm <- function(osmGpkg, outputCrs) { splitPathsAtPoints(paths.to.resplit, endpoints.for.resplit, 0.1, "path_id") # convert to dataframe, snap to grid, remove empty geometries - echo("Combining the resplit paths into a single dataframe") + echo("Combining the resplit paths into a single dataframe\n") system.time( resplit.paths <- bind_rows(resplit.path.list) %>% st_snap_to_grid(1) %>% @@ -185,7 +171,6 @@ processOsm <- function(osmGpkg, outputCrs) { combined.paths <- split.paths %>% filter(!path_id %in% paths.to.resplit$path_id) %>% rbind(resplit.paths) %>% - dplyr::select(osm_id) %>% # add a new id field, for joining to from and to id's mutate(combined_path_id = row_number()) @@ -199,7 +184,7 @@ processOsm <- function(osmGpkg, outputCrs) { # or 'level' tag status - # finalise paths + # finalise paths with metadata # -----------------------------------# # find from and to id's from endpoints @@ -220,26 +205,41 @@ processOsm <- function(osmGpkg, outputCrs) { st_drop_geometry() # assemble final paths with length, from_id and to_id - final.paths <- combined.paths %>% + final.paths.with.metadata <- combined.paths %>% # add length column mutate(length = as.integer(st_length(geom))) %>% # join from_id and to_id left_join(from_ids, by = "combined_path_id") %>% - left_join(to_ids, by = "combined_path_id") %>% - - # select final fields - dplyr::select(osm_id, length, from_id, to_id) + left_join(to_ids, by = "combined_path_id") - - - - # nodes..... + # nodes # -----------------------------------# + + echo("Finding roundabout and traffic signal status of nodes\n") + # extract from_id and to_id of edges that are roundabouts + roundabout.edges <- final.paths.with.metadata %>% + st_drop_geometry() %>% + # left_join(osm.metadata, by = "osm_id") %>% # <<< don't now need this + # is_roundabout: 1 if 'roundabout' in 'other_tags'; 0 if not, or if 'other_tags' is NA + mutate(is_roundabout = case_when( + is.na(other_tags) ~ 0, + str_detect(other_tags, "roundabout") ~ 1, + TRUE ~ 0)) %>% + dplyr::select(from_id, to_id, is_roundabout) + + # find node ids that connect to roundabouts + roundabout.nodes <- roundabout.edges %>% + filter(is_roundabout == 1) %>% + # combine from_id and to_id columns into new 'id' column + tidyr::gather(key = "key", value = "id", from_id, to_id) %>% + # keep distinct ids + dplyr::select(id) %>% + distinct() - # extract the traffic signals [MOVE THIS DOWN TO THE BIT DEALING WITH INTERSECTIONS] + # extract the traffic signals traffic.signals <- osm.points %>% # filter to traffic signals filter(str_detect(highway, "traffic_signals")) %>% @@ -247,61 +247,34 @@ processOsm <- function(osmGpkg, outputCrs) { st_snap_to_grid(1) %>% dplyr::select(osm_id, highway, other_tags) + # find endpoints within 20m of traffic signals + endpoints.near.signals <- endpoints %>% + st_filter(st_buffer(traffic.signals, 20), .predicate = st_intersects) %>% + .$endpoint_id + # finalise nodes: attribute with roundabout and signal status + final.nodes <- endpoints %>% + mutate(is_roundabout = ifelse(endpoint_id %in% roundabout.nodes$id, 1, 0), + is_signal = ifelse(endpoint_id %in% endpoints.near.signals, 1, 0)) - - - - #------ working section - - # # read in temporary tables - # sql.tables <- "./SP_working/temp_melb_sql_tables.sqlite" - # roads <- st_read(sql.tables, layer = "roads") - # - # write out temporary Bendigo outputs - bend.out <- "./SP_working/temp_bendigo.sqlite" - st_write(paths, bend.out, layer = "paths", delete_layer = T) - st_write(intersections, bend.out, layer = "intersections", delete_layer = T) - st_write(split.paths, bend.out, layer = "split_paths", delete_layer = T) - st_write(endpoints, bend.out, layer = "endpoints", delete_layer = T) - st_write(final.paths, bend.out, layer = "final_paths", delete_layer = T) - - st_write(split.paths, bend.out, layer = "split_paths_diff_loop", delete_layer = T) - st_write(paths.to.resplit, bend.out, layer = "paths_to_resplit", delete_layer = T) - st_delete(bend.out, layer = "split_paths_unbuffered_point_diff") - - paths <- st_read(bend.out, layer = "paths") - intersections <- st_read(bend.out, layer = "intersections") - split.paths <- st_read(bend.out, layer = "split_paths") - endpoints <- st_read(bend.out, layer = "endpoints") - final.paths <- st_read(bend.out, layer = "final_paths") - - # write out / read in temporary Melbourne outputs - melb.out <- "./SP_working/temp_melbourne.sqlite" - st_write(paths, melb.out, layer = "paths", delete_layer = T) - st_write(intersections, melb.out, layer = "intersections", delete_layer = T) - st_write(split.paths, melb.out, layer = "split_paths", delete_layer = T) - st_write(endpoints, melb.out, layer = "endpoints", delete_layer = T) - st_write(final.paths, melb.out, layer = "final_paths", delete_layer = T) - - st_write(paths.to.resplit, melb.out, layer = "paths_to_resplit", delete_layer = T) - st_write(problem.paths, melb.out, layer = "problem_paths", delete_layer = T) - - st_delete(melb.out, layer = "problem_paths") - - paths <- st_read(melb.out, layer = "paths") - intersections <- st_read(melb.out, layer = "intersections") - split.paths <- st_read(melb.out, layer = "split_paths") - endpoints <- st_read(melb.out, layer = "endpoints") - final.paths <- st_read(melb.out, layer = "final_paths") + # separate final paths and osm metadata + # -----------------------------------# - #--------end working section + # remove metadata from final paths + final.paths <- final.paths.with.metadata %>% + dplyr::select(osm_id, length, from_id, to_id) - # temporary return statement for testing - return(list(paths, intersections, endpoints, final.paths)) + # extract the non-spatial data for the paths + osm.metadata <- paths %>% + st_drop_geometry() %>% + dplyr::select(osm_id, highway, other_tags) %>% + filter(osm_id %in% final.paths$osm_id) + # return outputs + # -----------------------------------# + return(list(final.nodes, final.paths, osm.metadata)) }