From 232d64d26a7d9941d1a8e6293daf335bc3a26688 Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 09:20:22 +1000 Subject: [PATCH 01/11] Trunk roads walkable/cyclable by default --- functions/buildDefaultsDF.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/functions/buildDefaultsDF.R b/functions/buildDefaultsDF.R index 65b8216..63c32d5 100644 --- a/functions/buildDefaultsDF.R +++ b/functions/buildDefaultsDF.R @@ -4,8 +4,8 @@ buildDefaultsDF <- function(){ ~highway , ~permlanes, ~freespeed, ~is_oneway, ~laneCapacity, ~is_cycle, ~is_walk, ~is_car, ~highway_order, "motorway" , 4 , (110/3.6), 0 , 2000 , 0 , 0 , 1 , 1 , "motorway_link" , 2 , (80/3.6) , 0 , 1500 , 0 , 0 , 1 , 8 , - "trunk" , 3 , (100/3.6), 0 , 2000 , 0 , 0 , 1 , 2 , - "trunk_link" , 2 , (80/3.6) , 0 , 1500 , 0 , 0 , 1 , 9 , + "trunk" , 3 , (100/3.6), 0 , 2000 , 1 , 1 , 1 , 2 , + "trunk_link" , 2 , (80/3.6) , 0 , 1500 , 1 , 1 , 1 , 9 , "primary" , 2 , (80/3.6) , 0 , 1500 , 1 , 1 , 1 , 3 , "primary_link" , 1 , (60/3.6) , 0 , 1500 , 1 , 1 , 1 , 10 , From 1743e138661e53e964421aeb52530510d1aeb540 Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 09:52:46 +1000 Subject: [PATCH 02/11] Remove specific Melbourne corrections --- NetworkGenerator.R | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/NetworkGenerator.R b/NetworkGenerator.R index 1282384..2d165a9 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -24,11 +24,7 @@ makeNetwork<-function(outputFileName="test"){ desnificationMaxLengh=500 densifyBikeways=F - # CORRECTION - # To add/remove specified links - see osmCorrection.R - # Change to TRUE if running on Greater Melbourne OSM, Otherwise, keep FALSE - # Also you can use the same function to correct networks for your region if needed - correctNetwork=F + # CAPACITY ADJUSTMENT # A flag for whether to multiply capacity of links shorter than 100m by 2 or not # In some cases such as when building network for simulation of small samples (e.g. <1%) it might be desired adjustCapacity=F @@ -145,15 +141,8 @@ makeNetwork<-function(outputFileName="test"){ echo("Processing OSM tags and joining with defaults\n") system.time( osmAttributes <- processOsmTags(osm_metadata,defaults_df)) - # There are some roads in OSM that are not correctly attributed - # Use the function below to manually add their attributes based osm id - osmAttributesCorrected <- osmMetaCorrection(osmAttributes) - edgesOsm <- networkInput[[2]] - # Some network link corrections (+/-) specifically for Greater Melbourne OSM - if(correctNetwork) edgesOsm <- osmNetworkCorrection(networkInput) - - edgesAttributed <- edgesOsm %>% - inner_join(osmAttributesCorrected, by="osm_id") %>% + edgesAttributed <- networkInput[[2]] %>% + inner_join(osmAttributes, by="osm_id") %>% # dplyr::select(-osm_id,highway,highway_order) dplyr::select(-highway,highway_order) From fba3434ec7f47a6440169b7ca4fcccaa05f8b2a0 Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 16:13:41 +1000 Subject: [PATCH 03/11] CRS as parameter not hard-coded --- NetworkGenerator.R | 25 ++++++++++++------- functions/combineRedundantEdges.R | 4 +-- functions/combineUndirectedAndDirectedEdges.R | 4 +-- functions/gtfs2PtNetwork.R | 2 +- functions/makeEdgesDirect.R | 4 +-- functions/removeRedundantUndirectedEdges.R | 4 +-- functions/simplifyIntersections.R | 6 ++--- functions/simplifyNetwork.R | 6 ++--- functions/writeOutputs.R | 12 ++++----- 9 files changed, 37 insertions(+), 30 deletions(-) diff --git a/NetworkGenerator.R b/NetworkGenerator.R index 2d165a9..d3ed555 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -162,7 +162,8 @@ makeNetwork<-function(outputFileName="test"){ # simplify intersections while preserving attributes and original geometry. system.time(intersectionsSimplified <- simplifyIntersections(largestComponent[[1]], largestComponent[[2]], - shortLinkLength)) + shortLinkLength, + outputCrs)) # Merge edges going between the same two nodes, picking the shortest geometry. # * One-way edges going in the same direction will be merged @@ -171,14 +172,16 @@ makeNetwork<-function(outputFileName="test"){ # * One-way edges will NOT be merged with two-way edges. # * Non-car edges do NOT count towards the merged lane count (permlanes) system.time(edgesCombined <- combineRedundantEdges(intersectionsSimplified[[1]], - intersectionsSimplified[[2]])) + intersectionsSimplified[[2]], + outputCrs)) # Merge one-way and two-way edges going between the same two nodes. In these # cases, the merged attributes will be two-way. # This guarantees that there will only be a single edge between any two nodes. system.time(combinedUndirectedAndDirected <- combineUndirectedAndDirectedEdges(edgesCombined[[1]], - edgesCombined[[2]])) + edgesCombined[[2]], + outputCrs)) # If there is a chain of edges between intersections, merge them together system.time(edgesSimplified <- simplifyLines(combinedUndirectedAndDirected[[1]], @@ -190,15 +193,18 @@ makeNetwork<-function(outputFileName="test"){ # Do a second round of simplification. system.time(edgesCombined2 <- combineRedundantEdges(noDangles[[1]], - noDangles[[2]])) + noDangles[[2]], + outputCrs)) system.time(combinedUndirectedAndDirected2 <- combineUndirectedAndDirectedEdges(edgesCombined2[[1]], - edgesCombined2[[2]])) + edgesCombined2[[2]], + outputCrs)) system.time(edgesSimplified2 <- simplifyLines(combinedUndirectedAndDirected2[[1]], combinedUndirectedAndDirected2[[2]])) system.time(edgesCombined3 <- combineRedundantEdges(edgesSimplified2[[1]], - edgesSimplified2[[2]])) + edgesSimplified2[[2]], + outputCrs)) networkMode <- addMode(edgesCombined3) @@ -217,7 +223,8 @@ makeNetwork<-function(outputFileName="test"){ # simplify geometry so all edges are straight lines system.time(networkDirect <- makeEdgesDirect(networkDensified[[1]], - networkDensified[[2]])) + networkDensified[[2]], + outputCrs)) # add mode to edges, add type to nodes, change cycleway from numbers to text networkRestructured <- restructureData(networkDirect, highway_lookup, @@ -270,8 +277,8 @@ makeNetwork<-function(outputFileName="test"){ echo("| **Launching Output Writing** |\n") echo("--------------------------------------------------------\n") - if(writeSqlite) system.time(exportSQlite(networkFinal, outputDir)) - if(writeShp) system.time(exportShp(networkFinal, outputDir)) + if(writeSqlite) system.time(exportSQlite(networkFinal, outputDir, outputCrs)) + if(writeShp) system.time(exportShp(networkFinal, outputDir, outputCrs)) if(writeXml) system.time(exportXML(networkFinal, outputDir)) } diff --git a/functions/combineRedundantEdges.R b/functions/combineRedundantEdges.R index 3294708..f6abb36 100644 --- a/functions/combineRedundantEdges.R +++ b/functions/combineRedundantEdges.R @@ -1,7 +1,7 @@ # nodes_current<-intersectionsSimplified[[1]] # edges_current<-intersectionsSimplified[[2]] -combineRedundantEdges <- function(nodes_current,edges_current){ +combineRedundantEdges <- function(nodes_current,edges_current,outputCrs){ # assuming a dataframe with a 'current_group' column, merge edges together groupingFunction <- function(grouped_edges) { @@ -130,7 +130,7 @@ combineRedundantEdges <- function(nodes_current,edges_current){ inner_join(edges_all, by="uid") %>% dplyr::select(-uid) %>% st_sf() %>% - st_set_crs(28355) + st_set_crs(outputCrs) return(list(nodes_current,edges_all_geom)) } diff --git a/functions/combineUndirectedAndDirectedEdges.R b/functions/combineUndirectedAndDirectedEdges.R index 566dfeb..17748e9 100644 --- a/functions/combineUndirectedAndDirectedEdges.R +++ b/functions/combineUndirectedAndDirectedEdges.R @@ -1,7 +1,7 @@ # nodes_current<-edgesCombined[[1]] # edges_current<-edgesCombined[[2]] -combineUndirectedAndDirectedEdges <- function(nodes_current,edges_current){ +combineUndirectedAndDirectedEdges <- function(nodes_current,edges_current,outputCrs){ edges_current <- edges_current %>% mutate(uid=row_number()) %>% @@ -80,7 +80,7 @@ combineUndirectedAndDirectedEdges <- function(nodes_current,edges_current){ inner_join(edges_grouped2, by="uid") %>% dplyr::select(-uid) %>% st_sf() %>% - st_set_crs(28355) + st_set_crs(outputCrs) return(list(nodes_current,edges_all_geom)) } diff --git a/functions/gtfs2PtNetwork.R b/functions/gtfs2PtNetwork.R index 5892739..209556f 100644 --- a/functions/gtfs2PtNetwork.R +++ b/functions/gtfs2PtNetwork.R @@ -5,7 +5,7 @@ addGtfsLinks <- function(outputLocation="./test/", analysis_start = as.Date("2019-10-11","%Y-%m-%d"), analysis_end = as.Date("2019-10-17","%Y-%m-%d"), studyRegion=NA, - outputCrs=28355){ + outputCrs=outputCrs){ # outputLocation="./gtfs/" # nodes=networkRestructured[[1]] # links=networkRestructured[[2]] diff --git a/functions/makeEdgesDirect.R b/functions/makeEdgesDirect.R index 5afbaf6..12bea06 100644 --- a/functions/makeEdgesDirect.R +++ b/functions/makeEdgesDirect.R @@ -1,6 +1,6 @@ # nodes_current <-combinedUndirectedAndDirected2[[1]] # edges_current <-combinedUndirectedAndDirected2[[2]] -makeEdgesDirect <- function(nodes_current,edges_current){ +makeEdgesDirect <- function(nodes_current,edges_current,outputCrs){ # nodes_coords <- nodes_current %>% # st_drop_geometry() %>% @@ -15,7 +15,7 @@ makeEdgesDirect <- function(nodes_current,edges_current){ left_join(st_drop_geometry(nodes_current),by=c("to_id"="id")) %>% rename(toX=X,toY=Y) %>% mutate(geom=paste0("LINESTRING(",fromX," ",fromY,",",toX," ",toY,")")) %>% - st_as_sf(wkt = "geom", crs = 28355) + st_as_sf(wkt = "geom", crs = outputCrs) return(list(nodes_current,edges_current)) } diff --git a/functions/removeRedundantUndirectedEdges.R b/functions/removeRedundantUndirectedEdges.R index 7959110..0ae1ea4 100644 --- a/functions/removeRedundantUndirectedEdges.R +++ b/functions/removeRedundantUndirectedEdges.R @@ -1,7 +1,7 @@ # nodes_current<-networkSimplified[[1]] # edges_current<-networkSimplified[[2]] -removeRedundantUndirectedEdges <- function(nodes_current,edges_current,road_types){ +removeRedundantUndirectedEdges <- function(nodes_current,edges_current,road_types,outputCrs){ isOneWay <- road_types %>% dplyr::select(road_type,oneway) @@ -29,7 +29,7 @@ removeRedundantUndirectedEdges <- function(nodes_current,edges_current,road_type ungroup() %>% data.frame() %>% st_sf() %>% - st_set_crs(28355) + st_set_crs(outputCrs) return(list(nodes_current,edgesNoRedundancies)) } diff --git a/functions/simplifyIntersections.R b/functions/simplifyIntersections.R index af5fc19..9f45dc3 100644 --- a/functions/simplifyIntersections.R +++ b/functions/simplifyIntersections.R @@ -1,4 +1,4 @@ -simplifyIntersections <- function(n_df, l_df, shortLinkLength=10){ +simplifyIntersections <- function(n_df, l_df, shortLinkLength=10, outputCrs){ # shortLinkLength = 20 # l_df=largestComponent[[2]] # n_df=largestComponent[[1]] @@ -48,7 +48,7 @@ simplifyIntersections <- function(n_df, l_df, shortLinkLength=10){ filter(!id %in% comp_df$node_id) %>% rbind(dplyr::select(comp_df_centroid,id=cluster_id,is_roundabout,is_signal,X,Y)) %>% mutate(geom=paste0("POINT(",X," ",Y,")")) %>% - st_as_sf(wkt = "geom", crs = 28355) + st_as_sf(wkt = "geom", crs = outputCrs) # this function adds endpoints to the geometries addEndpoints <- function(fromX,fromY,toX,toY,geom) { @@ -84,7 +84,7 @@ simplifyIntersections <- function(n_df, l_df, shortLinkLength=10){ st_drop_geometry() %>% mutate(geom=geomExtended) %>% st_sf() %>% - st_set_crs(28355) %>% + st_set_crs(outputCrs) %>% # remove any loops filter(from_id != to_id) %>% dplyr::select(-fromX,-fromY,-toX,-toY) %>% diff --git a/functions/simplifyNetwork.R b/functions/simplifyNetwork.R index e05b98e..d6c5106 100644 --- a/functions/simplifyNetwork.R +++ b/functions/simplifyNetwork.R @@ -1,4 +1,4 @@ -simplifyNetwork <- function(n_df,l_df,osm_metadata, shortLinkLength = 20){ +simplifyNetwork <- function(n_df,l_df,osm_metadata, shortLinkLength = 20, outputCrs){ # l_df=lines_np # n_df=nodes_np # shortLinkLength = 20 @@ -55,7 +55,7 @@ simplifyNetwork <- function(n_df,l_df,osm_metadata, shortLinkLength = 20){ filter(!id %in% comp_df$node_id) %>% rbind(dplyr::select(comp_df_centroid,id=cluster_id,is_roundabout,is_signal,X,Y)) %>% mutate(geom=paste0("POINT(",X," ",Y,")")) %>% - st_as_sf(wkt = "geom", crs = 28355) + st_as_sf(wkt = "geom", crs = outputCrs) # keeping only the long edges, we alter any endpoints that are part of a # cluster, replacing them with the cluster id and the centroid coordinates @@ -72,7 +72,7 @@ simplifyNetwork <- function(n_df,l_df,osm_metadata, shortLinkLength = 20){ mutate(toY=ifelse(is.na(cluster_id),toY,Y)) %>% dplyr::select(road_type,length,from_id,to_id,fromX,fromY,toX,toY) %>% mutate(geom=paste0("LINESTRING(",fromX," ",fromY,",",toX," ",toY,")")) %>% - st_as_sf(wkt = "geom", crs = 28355) + st_as_sf(wkt = "geom", crs = outputCrs) # remove the disconnected bits diff --git a/functions/writeOutputs.R b/functions/writeOutputs.R index 5568ad1..c2eba00 100644 --- a/functions/writeOutputs.R +++ b/functions/writeOutputs.R @@ -1,6 +1,6 @@ # SQlite ------------------------------------------------------------------ -exportSQlite <- function(networkFinal, outputDir){ +exportSQlite <- function(networkFinal, outputDir, outputCrs){ cat('\n') echo(paste0('Writing the sqlite output: ', nrow(networkFinal[[2]]), @@ -11,14 +11,14 @@ exportSQlite <- function(networkFinal, outputDir){ if(class(networkFinal[[1]])[1]!="sf"){ networkFinal[[1]] <- networkFinal[[1]] %>% mutate(GEOMETRY=paste0("POINT(",x," ",y,")")) %>% - st_as_sf(wkt = "GEOMETRY", crs = 28355) %>% + st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% as.data.frame() %>% st_sf() } if(class(networkFinal[[2]])[1]!="sf"){ networkFinal[[2]] <- networkFinal[[2]] %>% mutate(GEOMETRY=paste0("LINESTRING(",fromX," ",fromY,",",toX," ",toY,")")) %>% - st_as_sf(wkt = "GEOMETRY", crs = 28355) %>% + st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% as.data.frame() %>% st_sf() } @@ -35,7 +35,7 @@ exportSQlite <- function(networkFinal, outputDir){ } # ShapeFile --------------------------------------------------------------- -exportShp <- function(networkFinal, outputDir){ +exportShp <- function(networkFinal, outputDir, outputCrs){ cat('\n') echo(paste0('Writing the ShapeFile output: ', nrow(networkFinal[[2]]), @@ -44,14 +44,14 @@ exportShp <- function(networkFinal, outputDir){ if(class(networkFinal[[1]])!="sf"){ networkFinal[[1]] <- networkFinal[[1]] %>% mutate(GEOMETRY=paste0("POINT(",x," ",y,")")) %>% - st_as_sf(wkt = "GEOMETRY", crs = 28355) %>% + st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% as.data.frame() %>% st_sf() } if(class(networkFinal[[2]])!="sf"){ networkFinal[[2]] <- networkFinal[[2]] %>% mutate(GEOMETRY=paste0("LINESTRING(",fromX," ",fromY,",",toX," ",toY,")")) %>% - st_as_sf(wkt = "GEOMETRY", crs = 28355) %>% + st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% as.data.frame() %>% st_sf() } From 0d83dbf298ab68b14ef044ba89e1fda50a662f7c Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 16:14:14 +1000 Subject: [PATCH 04/11] Fix missing CRS --- functions/densifyNetwork.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/functions/densifyNetwork.R b/functions/densifyNetwork.R index dc6bc2a..1fb3f93 100644 --- a/functions/densifyNetwork.R +++ b/functions/densifyNetwork.R @@ -86,14 +86,14 @@ densifyNetwork <- function(networkList, minimum_length=400, densifyBikeways=F){ links_combined <- bind_rows( links_unsegmented, - links_segmented + links_segmented %>% st_set_crs(st_crs(links_unsegmented)) ) %>% dplyr::select(-tmp_id) %>% st_sf() nodes_combined <- bind_rows( nodes_df, - nodes_segmented + nodes_segmented %>% st_set_crs(st_crs(nodes_df)) ) %>% st_sf() From 0c3e3faeb3cbf1eb2c1b0de992d0771a9416d537 Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 16:15:06 +1000 Subject: [PATCH 05/11] Fix error when can't find geometry for NA --- functions/gtfs2PtNetwork.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/functions/gtfs2PtNetwork.R b/functions/gtfs2PtNetwork.R index 209556f..65ea206 100644 --- a/functions/gtfs2PtNetwork.R +++ b/functions/gtfs2PtNetwork.R @@ -123,7 +123,7 @@ processGtfs <- function(outputLocation="./test/", st_snap_to_grid(1) # only want stops within the study region - if(!is.na(st_geometry(studyRegion))){ + if(!is.na(studyRegion)){ message("Cropping to study region") validStops <- validStops %>% filter(lengths(st_intersects(., studyRegion)) > 0) From fdc1fce9adf643d4c2a2a4cf108fa4d42f1d0dd0 Mon Sep 17 00:00:00 2001 From: StevePem Date: Fri, 4 Aug 2023 23:01:37 +1000 Subject: [PATCH 06/11] handle densification for lines digitised in either direction --- functions/densifyNetwork.R | 53 +++++++++++++++++++++++++++++++++----- 1 file changed, 47 insertions(+), 6 deletions(-) diff --git a/functions/densifyNetwork.R b/functions/densifyNetwork.R index 1fb3f93..8994606 100644 --- a/functions/densifyNetwork.R +++ b/functions/densifyNetwork.R @@ -61,6 +61,22 @@ densifyNetwork <- function(networkList, minimum_length=400, densifyBikeways=F){ if(i==length(links_list)) cat(paste0(i,"/",length(links_list)," rows segmented\n")) } + # add flag for direction in which geometry of links_to_segmentize is recorded + links_to_segmentize <- links_to_segmentize %>% + # join X & Y coordinates for from_id + left_join(nodes_df %>% + st_drop_geometry() %>% + dplyr::select(id, fromx = X, fromy = Y), + by = c("from_id" = "id")) %>% + # check whether startpoint of geometry matches from_id ("forward" if yes, "reverse" if no) + mutate(startpoint = st_coordinates(st_startpoint(geom))) %>% + rowwise() %>% + mutate(direction = + ifelse(startpoint[[1]] == fromx & startpoint[[2]] == fromy, + "forward", + "reverse")) %>% + ungroup() + links_segmented <- links_to_segmentize %>% st_set_geometry(st_sfc(links_list_segmented)) %>% mutate(group_id=row_number()) %>% @@ -69,21 +85,46 @@ densifyNetwork <- function(networkList, minimum_length=400, densifyBikeways=F){ st_cast(to="LINESTRING") %>% mutate(new_node_id=row_number()+max(nodes_df$id,na.rm=T)) %>% group_by(group_id) %>% - mutate(from_id=ifelse(row_number()!=1,new_node_id-1,from_id)) %>% - mutate(to_id=ifelse(row_number()!=max(row_number()),new_node_id,to_id)) %>% + # mutate(from_id=ifelse(row_number()!=1,new_node_id-1,from_id)) %>% + mutate(from_id = case_when( + direction == "forward" & row_number() == 1 ~ from_id, + direction == "forward" & row_number() != 1 ~ new_node_id-1, + direction == "reverse" & row_number() != max(row_number()) ~ new_node_id, + direction == "reverse" & row_number() == max(row_number()) ~ from_id + )) %>% + # mutate(to_id=ifelse(row_number()!=max(row_number()),new_node_id,to_id)) %>% + mutate(to_id = case_when( + direction == "forward" & row_number() != max(row_number()) ~ new_node_id, + direction == "forward" & row_number() == max(row_number()) ~ to_id, + direction == "reverse" & row_number() == 1 ~ to_id, + direction == "reverse" & row_number() != 1 ~ new_node_id-1, + )) %>% mutate(length=round(as.numeric(st_length(geom)),3)) %>% - dplyr::select(-group_id,-new_node_id) + dplyr::select(-group_id, -new_node_id, -fromx, -fromy, -startpoint) - nodes_segmented <- links_segmented %>% + nodes_segmented_forward <- links_segmented %>% + filter(direction == "forward") %>% dplyr::select(id=to_id) %>% filter(id>max(nodes_df$id,na.rm=T)) %>% - st_set_geometry(st_endpoint(.)) %>% - mutate(is_roundabout=0,is_signal=0) + st_set_geometry(st_endpoint(.)) + nodes_segmented_reverse <- links_segmented %>% + filter(direction == "reverse") %>% + dplyr::select(id=to_id) %>% + filter(id>max(nodes_df$id,na.rm=T)) %>% + st_set_geometry(st_startpoint(.)) + + nodes_segmented <- rbind(nodes_segmented_forward, + nodes_segmented_reverse) %>% + mutate(is_roundabout=0,is_signal=0) + nodes_segmented <- bind_cols(nodes_segmented, data.frame(st_coordinates(nodes_segmented))) %>% dplyr::select(id,is_roundabout,is_signal,X,Y) + links_segmented <- links_segmented %>% + dplyr::select(-direction) + links_combined <- bind_rows( links_unsegmented, links_segmented %>% st_set_crs(st_crs(links_unsegmented)) From c99bea949c635e8662fd0096ab3ce651912f2fd7 Mon Sep 17 00:00:00 2001 From: StevePem Date: Sat, 5 Aug 2023 09:19:09 +1000 Subject: [PATCH 07/11] correctly recognise class of network --- functions/writeOutputs.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/functions/writeOutputs.R b/functions/writeOutputs.R index c2eba00..672a8d1 100644 --- a/functions/writeOutputs.R +++ b/functions/writeOutputs.R @@ -41,14 +41,14 @@ exportShp <- function(networkFinal, outputDir, outputCrs){ echo(paste0('Writing the ShapeFile output: ', nrow(networkFinal[[2]]), ' links and ', nrow(networkFinal[[1]]),' nodes\n')) - if(class(networkFinal[[1]])!="sf"){ + if(class(networkFinal[[1]])[1]!="sf"){ networkFinal[[1]] <- networkFinal[[1]] %>% mutate(GEOMETRY=paste0("POINT(",x," ",y,")")) %>% st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% as.data.frame() %>% st_sf() } - if(class(networkFinal[[2]])!="sf"){ + if(class(networkFinal[[2]])[1]!="sf"){ networkFinal[[2]] <- networkFinal[[2]] %>% mutate(GEOMETRY=paste0("LINESTRING(",fromX," ",fromY,",",toX," ",toY,")")) %>% st_as_sf(wkt = "GEOMETRY", crs = outputCrs) %>% From d5a50055bfd1741a7089758fb529fa0b16a8a635 Mon Sep 17 00:00:00 2001 From: StevePem Date: Sat, 5 Aug 2023 11:11:38 +1000 Subject: [PATCH 08/11] fix writing vehicleTripMatching --- functions/gtfs2PtNetwork.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/functions/gtfs2PtNetwork.R b/functions/gtfs2PtNetwork.R index 65ea206..23be99c 100644 --- a/functions/gtfs2PtNetwork.R +++ b/functions/gtfs2PtNetwork.R @@ -503,12 +503,12 @@ exportGtfsSchedule <- function(links, str<-paste0(str," \n") } - if (i%%writeInterval==0 || i==nrow(vehicleTripMatching)) { + if (i%%writeInterval==0 || i==length(transitRoutes)) { cat(str,file=outxml,append=TRUE) str<-"" # clear the buffer after writing it out } # report progress - if (i%%50==0 || i==nrow(vehicleTripMatching)) printProgress(i,nrow(vehicleTripMatching),' vehicleTripMatching') + if (i%%50==0 || i==length(transitRoutes)) printProgress(i,length(transitRoutes),' vehicleTripMatching') } cat(paste0(" \n"),file=outxml,append=TRUE) cat(paste0("\n"),file=outxml,append=TRUE) From b661471d8d7fb79d187b768955d7ca8d00b9c05f Mon Sep 17 00:00:00 2001 From: StevePem Date: Thu, 10 Aug 2023 08:41:31 +1000 Subject: [PATCH 09/11] cycling file locations and flags --- NetworkGenerator.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/NetworkGenerator.R b/NetworkGenerator.R index d3ed555..4395e04 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -11,7 +11,7 @@ makeNetwork<-function(outputFileName="test"){ # Note that osm.pbf format is not yet supported osmExtract='./data/melbourne.osm' # If procesOsm=F, set the following to the network sqlite file - networkSqlite="data/network.sqlite" + networkSqlite="./data/melbourne_network_unconfigured.sqlite" # SIMPLIFICATION shortLinkLength=20 @@ -22,7 +22,7 @@ makeNetwork<-function(outputFileName="test"){ # DENSIFICATION desnificationMaxLengh=500 - densifyBikeways=F + densifyBikeways=T # CAPACITY ADJUSTMENT # A flag for whether to multiply capacity of links shorter than 100m by 2 or not @@ -31,11 +31,11 @@ makeNetwork<-function(outputFileName="test"){ # ELEVATION # A flag for whether to add elevation or not - addElevation=F + addElevation=T # Digital elevation model file - make sure it is in the same coordinate system as your network - demFile= 'data/DEMx10EPSG28355.tif' + demFile= "./data/DEM_melbourne.tif" # DEM's multiplier- set to 1 if DEM contains actual elevation - ElevationMultiplier=10 + ElevationMultiplier=1 # GTFS addGtfs=F From 51c8fc8777ab168a0794dc74c67c8487189e35e3 Mon Sep 17 00:00:00 2001 From: StevePem Date: Tue, 15 Aug 2023 11:45:11 +1000 Subject: [PATCH 10/11] LTS, LTS impedance, slope impedance --- CyclingImpedances.R | 74 ++++++++++++++ Destinations.R | 169 ++++++++++++++++++++++++++++++++ NetworkGenerator.R | 8 +- functions/addLTS.R | 152 ++++++++++++++++++++++++++++ functions/addSlopeImped.R | 26 +++++ functions/addTraffic.R | 46 +++++++++ functions/getDestinationTypes.R | 159 ++++++++++++++++++++++++++++++ functions/makeEdgesOneway.R | 70 +++++++++++++ 8 files changed, 703 insertions(+), 1 deletion(-) create mode 100644 CyclingImpedances.R create mode 100644 Destinations.R create mode 100644 functions/addLTS.R create mode 100644 functions/addSlopeImped.R create mode 100644 functions/addTraffic.R create mode 100644 functions/getDestinationTypes.R create mode 100644 functions/makeEdgesOneway.R diff --git a/CyclingImpedances.R b/CyclingImpedances.R new file mode 100644 index 0000000..e38e27d --- /dev/null +++ b/CyclingImpedances.R @@ -0,0 +1,74 @@ +# Using a network created by 'NetworkGenerator.R', add impedances for +# calculating cycling accessibility + +# Assumes input network is a one-way network that includes elevation, +# and that one-way daily traffic volumes are available + +addImpedances <- function() { + + # Parameters ----------------------------------------------------------------- + # Input network, to which cycling impedances are to be added, with layer names + # input.network <- "./output/test/network.sqlite" + input.network <- "./output/test/melbourne_network.sqlite" #<<< OLD VERSION FOR TESTING + input.node.layer <- "nodes" + input.link.layer <- "links" + + # Traffic file, with links layer name - file must match input.network, + # and with a 'total_vol' column containing 1-way daily traffic + # traffic.file <- "./output/test/network_traffic.sqlite" + traffic.file <- "./output/test/links_with_traffic.sqlite" #<<< OLD VERSION FOR TESTING + traffic.link.layer <- "cars_aht" + + # Traffic multiplier - where volumes are for a sample of traffic only (eg + # multiplier of 20 if the volumes are a 5% sample; use 1 if full volumes) + traffic.multiplier <- 10 + + # Output location - same directory as input.network + output.location <- paste0(path_dir(input.network), "/networkWeighted.sqlite") + + + # Packages and functions ----------------------------------------------------- + library(dplyr) + library(sf) + library(fs) + + dir_walk(path="./functions/",source, recurse=T, type = "file") + + + # Load input network and traffic file ---------------------------------------- + input.nodes <- st_read(input.network, layer = input.node.layer) + input.links <- st_read(input.network, layer = input.link.layer) + traffic.links <- st_read(traffic.file, layer = traffic.link.layer) + + + # Add LTS and its impedance -------------------------------------------------- + echo("Adding daily traffic volumes\n") + networkTraffic <- addTraffic(input.nodes, + input.links, + 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 + + echo("Adding LTS and its impedance\n") + networkLTS <- addLTS(networkTraffic[[1]], networkTraffic[[2]]) + + + # Add slope impedance -------------------------------------------------------- + echo("Adding slope impedance") + networkSlope <- addSlopeImped(networkLTS[[1]], networkLTS[[2]]) + + + # Calculate total weight ----------------------------------------------------- + echo("Calculating cycling weight") + networkWeighted <- + list(networkSlope[[1]], + networkSlope[[2]] %>% + mutate(cycle.weight = length + LTS_imped + slope_imped)) + + + # write output --------------------------------------------------------------- + st_write(networkWeighted[[1]], "./output/test/networkWeighted.sqlite", layer = "nodes") + st_write(networkWeighted[[2]], "./output/test/networkWeighted.sqlite", layer = "links") + +} \ No newline at end of file diff --git a/Destinations.R b/Destinations.R new file mode 100644 index 0000000..d60c9d8 --- /dev/null +++ b/Destinations.R @@ -0,0 +1,169 @@ +# Make file of destinations required for accessibilty routing + +library(dplyr) +library(sf) +library(osmextract) + +# 1 Download OSM extract ---- +# -----------------------------------------------------------------------------# +# Download from https://www.interline.io/osm/extracts/ + +## Downloaded for Melbourne - melbourne_australia.osm.pbf + + +# 2 Converting to .gpkg format ---- +# -----------------------------------------------------------------------------# +# input file name and project CRS +INPUTFILE <- "./data/melbourne_australia.osm.pbf" + +PROJECT.CRS = 28355 + +# check layers +st_layers(INPUTFILE) + +# check keys +options(max.print = 2000) +polygon.tags <- oe_get_keys(INPUTFILE, layer = "multipolygons") %>% sort() +point.tags <- oe_get_keys(INPUTFILE, layer = "points") %>% sort() +line.tags <- oe_get_keys(INPUTFILE, layer = "lines") %>% sort() + +# create gpkg file in same directory as INPUTFILE, using the 'extra_tags' +# argument for specific extra tags required for various destination types +oe_vectortranslate(INPUTFILE, layer = "multipolygons", + extra_tags = c("access", "building", "grades", "healthcare", + "healthcare:speciality","isced:level", + "network", "operator", + "operator:type", "public_transport", "railway", + "school", "social_facility", "sport", + "tourism", "train")) +oe_vectortranslate(INPUTFILE, layer = "points", + 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")) +oe_vectortranslate(INPUTFILE, layer = "lines", + 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", + "smoothness", "surface")) +oe_vectortranslate(INPUTFILE, layer = "multilinestrings") +oe_vectortranslate(INPUTFILE, layer = "other_relations") + + +# 3 Read in the .gpkg file ---- +# -----------------------------------------------------------------------------# +GPKG <- "./data/melbourne_australia.gpkg" + +polygons <- st_read(GPKG, layer = "multipolygons") %>% st_transform(PROJECT.CRS) +points <- st_read(GPKG, layer = "points") %>% st_transform(PROJECT.CRS) +lines <- st_read(GPKG, layer = "lines") %>% st_transform(PROJECT.CRS) +multilines <- st_read(GPKG, layer = "multilinestrings") %>% st_transform(PROJECT.CRS) +other_relations <- st_read(GPKG, layer = "other_relations") %>% st_transform(PROJECT.CRS) + + +# 4 Extract required destinations ---- +# -----------------------------------------------------------------------------# + +## 4.1 Tag combinations for feature types and network ---- +## ----------------------------------------------------------------------------# +# load functions for locating specific feature types +source("./functions/getDestinationTypes.R") + +# load network +NETWORK <- "./output/test/network.sqlite" #<<< CHECK FINAL NAME +NODE.LAYER <- "nodes" +LINK.LAYER <- "links" + +network.nodes <- st_read(NETWORK, layer = NODE.LAYER) +network.links <- st_read(NETWORK, layer = LINK.LAYER) + + +## 4.2 Compile point and polygon destinations ---- +## ----------------------------------------------------------------------------# +destination.layer <- function(layer) { + return( + bind_rows( + getPlayground(layer) %>% mutate(dest_type = "playground"), + getPark(layer) %>% mutate(dest_type = "park"), + getSport(layer) %>% mutate(dest_type = "sport"), + getKindergarten(layer) %>% mutate(dest_type = "kindergarten"), + getCommunity(layer) %>% mutate(dest_type = "community_centre"), + getLibrary(layer) %>% mutate(dest_type = "library"), + getPrimary(layer) %>% mutate(dest_type = "primary_school"), + getSecondary(layer) %>% mutate(dest_type = "secondary_school"), + getClinic(layer) %>% mutate(dest_type = "health_clinic"), + getDentist(layer) %>% mutate(dest_type = "dentist"), + getPharmacy(layer) %>% mutate(dest_type = "pharmacy"), + getConvenience(layer) %>% mutate(dest_type = "convenience_store"), + getSupermarket(layer) %>% mutate(dest_type = "supermarket"), + getShop(layer) %>% mutate(dest_type = "shop"), + 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") + )) +} + +# create tables of destinations, and allocate unique id's (so features with +# multiple nodes can be grouped by the id where required) +destination.pt <- + bind_rows(destination.layer(points), + getStation() %>% mutate(dest_type = "railway_station")) %>% + mutate(dest_id = row_number()) + +destination.poly <- + destination.layer(polygons) %>% + mutate(dest_id = max(destination.pt$dest_id) + row_number()) + + +## 4.3 Find relevant nodes ---- +## ----------------------------------------------------------------------------# +# TO CONFIRM:- +# For all destinations except parks and schools ('small features'), relevant +# node is nearest node to point or to polygon centroid +# For parks and schools ('large features'): +# - points are buffered to 50m to create a polygon feature, +# - for buffered points and polygons, relevant nodes are all nodes within the +# feature and terminal nodes of links within 30m of boundary, or if none, +# then nearest node to boundary + +# Maybe this should be all nodes within 30m of buffered feature, and if link is within +# 30m of boundary but doesn't have a node within the buffer, then also its closest terminal +# node ??? + +dest.small <- bind_rows(destination.pt, + destination.poly %>% st_centroid()) %>% + filter(!(dest_type %in% c("park", "primary_school", "secondary_school"))) +near_node <- network.nodes$id[st_nearest_feature(dest.small, network.nodes)] +dest.small.with.nodes <- cbind(dest.small %>% st_drop_geometry(), near_node) + + +## NOTE - the code below is a simplified version which just finds nodes within +## feature or its 30m buffer, or nearest node if none - doesn't extend to terminal +## nodes of nearby features +dest.large <- bind_rows(destination.pt %>% st_buffer(50), + destination.poly) %>% + filter(dest_type %in% c("park", "primary_school", "secondary_school")) + +dest.large.found.nodes <- dest.large %>% + st_buffer(30) %>% + st_intersection(., network.nodes %>% dplyr::select(near_node = id)) + +dest.large.need.nodes <- dest.large %>% + filter(!(dest_id %in% dest.large.found.nodes$dest_id)) +near_node <- network.nodes$id[st_nearest_feature(dest.large.need.nodes, network.nodes)] + +dest.large.with.nodes <- bind_rows(dest.large.found.nodes %>% st_drop_geometry(), + cbind(dest.large.need.nodes %>% st_drop_geometry(), + near_node)) + +dest.with.nodes <- bind_rows(dest.small.with.nodes, + dest.large.with.nodes) + diff --git a/NetworkGenerator.R b/NetworkGenerator.R index 4395e04..0f25329 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -270,7 +270,13 @@ makeNetwork<-function(outputFileName="test"){ outputCrs=outputCrs)) } - networkFinal <- networkRestructured + # Make network oneway (required because cycling impedances such as level of + # traffic stress and slope may be different in each direction) + echo("Making all links one way\n") + networkOneway <- makeEdgesOneway(networkRestructured[[1]], + networkRestructured[[2]]) + + networkFinal <- networkOneway # writing outputs echo("========================================================\n") diff --git a/functions/addLTS.R b/functions/addLTS.R new file mode 100644 index 0000000..9b1baa3 --- /dev/null +++ b/functions/addLTS.R @@ -0,0 +1,152 @@ +# function to add level of traffic stress to network, based on highway class, +# traffic (ADT) and speed, and related impedance + +# traffic volumes (eg ADT 10000) are for two-way traffic, so are halved +# (eg 10000 / 2) in order to apply to the one-way links in edges_current + + + +addLTS <- function(nodes_current, edges_current) { + + # testing + # nodes_current <- networkTraffic[[1]] + # edges_current <- networkTraffic[[2]] + + # assign LTS to edges + # '1' to '4' are categories of increasing stress, as per table below] + + # road groups + local <- c("residential", "road", "unclassified", "living_street", "service") + tertiary <- c("tertiary", "tertiary_link") + secondary <- c("secondary", "secondary_link") + primary <- c("primary", "primary_link") + + + edges_current <- edges_current %>% + # make speed field (rounded, to avoid floating point issues) + mutate(speed = round(freespeed * 3.6)) %>% + # add LTS + mutate(lvl_traf_stress = case_when( + + # LTS 1 - off-road paths + cycleway %in% c("bikepath", "shared_path") ~ 1, + highway %in% c("cycleway", "track", "pedestrian", + "footway", "path", "corridor", "steps") ~ 1, + + # LTS 1 - separated cycle lanes + cycleway == "separated_lane" & speed <= 50 ~ 1, + + # LTS 1 - on-road cycle lanes + cycleway == "simple_lane" & + highway %in% c(local, tertiary, secondary) & + ADT <= 10000 /2 & speed <= 30 ~ 1, + + # LTS 1 - mixed traffic + highway %in% local & ADT <= 2000 / 2 & speed <= 30 ~ 1, + + # LTS 2 - separated cycle lanes + cycleway == "separated_lane" & speed <= 60 ~ 2, + + # LTS 2 - on-road cycle lanes + cycleway == "simple_lane" & + highway %in% c(local, tertiary, secondary) & + ADT <= 10000 / 2 & speed <= 50 ~ 2, + cycleway == "simple_lane" & + (highway %in% primary | + (highway %in% c(local, tertiary, secondary) & + ADT > 10000 / 2)) & + speed <= 40 ~ 2, + + # LTS 2 - mixed traffic + highway %in% local & ADT <= 750 / 2 & speed <= 50 ~ 2, + highway %in% local & ADT <= 2000 / 2 & speed <= 40 ~ 2, + highway %in% c(local, tertiary) & ADT <= 3000 / 2 & speed <= 30 ~ 2, + + # LTS 3 - on-road cycle lanes + cycleway == "simple_lane" & speed <= 60 ~ 3, + + # LTS 3 - mixed traffic + highway %in% local & ADT <= 750 / 2 & speed <= 60 ~ 3, + highway %in% c(local, tertiary) & ADT <= 3000 / 2 & speed <= 50 ~ 3, + (highway %in% c(secondary, primary) | + (highway %in% c(local, tertiary) & ADT > 3000 / 2)) & + speed <= 40 ~ 3, + + # LTS 4 - everything not covered above + TRUE ~ 4 + )) + + # check to test how many in each category + # LTS_table <- edges_current %>% + # st_drop_geometry() %>% + # group_by(highway, lvl_traf_stress) %>% + # summarise(n = n()) + + # assign LTS to nodes, based on highest + # begin with all nodes (from and to) and the LTS level of the associated link + node_max_lookup <- rbind(edges_current %>% + st_drop_geometry() %>% + dplyr::select(id = from_id, LTS = lvl_traf_stress), + edges_current %>% + st_drop_geometry() %>% + dplyr::select(id = to_id, LTS = lvl_traf_stress)) %>% + group_by(id) %>% + # find highest level of LTS for links connecting with the node + summarise(max_LTS = max(LTS)) %>% + ungroup() + + # Calculate impedance for intersection, and total impedance + + # Impedance for intersection applies to the to-node (the intersection + # that the link arrives at), and only if it's unsignalised + # penalty is calculated as: + # penaltya = (Buffb – Buffa) * (IFb – 1)for + # where + # a is the link for which the penalty (penaltya) is being calculated + # b is the highest-ranked other link at the relevant intersection + # Buffa and Buffb are the buffer distances for a and b, where the buffer + # distance is 0, 5, 10 or 25m for a link of LTS 1, 2, 3 or 4 respectively + # IFb is the impedance factor for b, where the impedance factor is 1.00, + # 1.05, 1.10 or 1.15 for a link of LTS 1, 2, 3 or 4 respectively + + # LTS impedance, which is to be added to the length of the link and any other + # impedances (outside this function) to create the weight for the link, + # is the length-based impedance for the link plus its intersection impedance. + # - Length-based impedance is the link multiplied by its impedance factor + # minus 1 (that is, subtracting 1 so it os only the additional impedance, + # not the length itself: + # total_imped = length * (IFa - 1) + # where IFa is the impedance factor for a, where the impedance factor is 1.00, + # 1.05, 1.10 or 1.15 for a link of LTS 1, 2, 3 or 4 respectively + # - Intersection impedance is calcualted as above + + buff_imped_df <- data.frame(cbind(LTS = c(1, 2, 3, 4), + buffer = c(0, 5, 10, 25), + imped = c(1, 1.05, 1.10, 1.15))) + + edges_current <- edges_current %>% + # join node intersection details for the to-node + left_join(., nodes_current %>% + st_drop_geometry() %>% + dplyr::select(id, type), + by = c("to_id" = "id")) %>% + # join the node max LTS buffer & impedance details for the to-node + left_join(., node_max_lookup, by = c("to_id" = "id")) %>% + left_join(., buff_imped_df, by = c("lvl_traf_stress" = "LTS")) %>% + # and the buff_imped_df details for the max LTS + left_join(., buff_imped_df, by = c("max_LTS" = "LTS"), suffix = c(".a", ".b")) %>% + + # calculate intersection impedance, using formula above (unsignalised only) + mutate(intersec_imped = ifelse(type %in% c("simple_intersection", + "simple_roundabout"), + (buffer.b - buffer.a) * (imped.b - 1), + 0)) %>% + # calculate total LTS impedance (to be added to length along with other impedances) + mutate(LTS_imped = (length * (imped.a - 1)) + intersec_imped) %>% + + # remove unwanted fields + dplyr::select(-speed, -type, -max_LTS, -buffer.a, -buffer.b, + -imped.a, -imped.b, -intersec_imped) + + return(list(nodes_current, edges_current)) +} diff --git a/functions/addSlopeImped.R b/functions/addSlopeImped.R new file mode 100644 index 0000000..c5ee42e --- /dev/null +++ b/functions/addSlopeImped.R @@ -0,0 +1,26 @@ +# function to add impedance to network links, based on slope + +# impedance is 50m of length per 1m of climb (Ziemke, D, Metzler, S & Nagel, K 2019, +# 'Bicycle traffic and its interaction with motorized traffic in an agent-based +# transport simulation framework', Future Generation Computer Systems, vol. 97, pp. 30-40.) + +addSlopeImped <- function(nodes_current, edges_current) { + + # testing + # nodes_current <- networkLTS[[1]] + # edges_current <- networkLTS[[2]] + + edges_current <- edges_current %>% + + # some coastal links are missing elevation; make slope 0 + mutate(slope_pct = ifelse(is.na(slope_pct), 0, slope_pct)) %>% + + # 50m of length per 1m of climb + # climb = run * slope_pct / 100 + mutate(slope_imped = case_when( + slope_pct <= 0 ~ 0, + slope_pct > 0 ~ (length * slope_pct / 100) * 50 + )) + + return(list(nodes_current, edges_current)) +} diff --git a/functions/addTraffic.R b/functions/addTraffic.R new file mode 100644 index 0000000..22c4cc5 --- /dev/null +++ b/functions/addTraffic.R @@ -0,0 +1,46 @@ +# function to add traffic to network links + +# assumes a network links file, and a second file which is the same (or, at +# least, contains the same links with the same 'from_id', 'to_id' and +# highway type), and with a 'total_vol' column containing 1-way daily traffic + +# traffic volumes are for each link in its operating direction (so one way +# traffic only); when comparing to two-way traffic volumes in 'addLTS.R', the +# two-way volumes need to be halved + +# the multiplier is used where volumes are for a sample of traffic only +# (eg use multiplier of 20 if the volumes are a 5% sample) + +addTraffic <- function(network.nodes, network.links, traffic.links, multiplier = 1) { + + # testing + # network.nodes = st_read("./output/test/melbourne_network.sqlite", layer = "nodes) + # network.links = st_read("./output/test/melbourne_network.sqlite", layer = "links") + # traffic.links = st_read("./output/test/links_with_traffic.sqlite", layer = "cars_aht") + # multiplier = 10 + + # select traffic fields from traffic.links + traffic <- traffic.links %>% + dplyr::select(from_id, to_id, highway, length, total_vol) %>% + st_drop_geometry() + + # join traffic to links, in each direction using 'highway' and 'length' as + # additional join ids: simplification can result in two links sharing nodes (eg loop) + links.with.traffic <- network.links %>% + left_join(traffic, by = c("from_id", "to_id", "highway", "length")) %>% + # convert NAs to zeros, and apply multiplier + mutate(total_vol = if_else(is.na(total_vol), 0, total_vol)) %>% + mutate(total_vol = total_vol * multiplier) %>% + rename(ADT = total_vol) %>% + + # remove duplicates, if any, which have arisen because of original duplicate links, + # summing their ATDs: group, sum the ATDs, then remove the duplicate with 'distinct' + group_by(link_id) %>% + mutate(ADT = sum(ADT)) %>% + distinct() %>% + ungroup() + + # return network with traffic volumes added + return(list(network.nodes, links.with.traffic)) + +} diff --git a/functions/getDestinationTypes.R b/functions/getDestinationTypes.R new file mode 100644 index 0000000..9010e82 --- /dev/null +++ b/functions/getDestinationTypes.R @@ -0,0 +1,159 @@ +# 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 + +# 1 open space ---- +getPlayground <- function(layer) { + return(layer %>% filter(leisure == "playground")) +} + +getPark <- function(layer) { + return(layer %>% filter(leisure == "park")) +} + + +# 2 sport ---- +getSport <- function(layer) { + return(layer %>% filter(!is.na(sport))) +} + + +# 3 lifelong learning ---- +getKindergarten <- function(layer) { + return(layer %>% filter(amenity == "kindergarten" | school == "kindergarten")) +} + +getCommunity <- function(layer) { + return(layer %>% filter(amenity == "community_centre")) +} + +getLibrary <- function(layer) { + return(layer %>% filter(amenity == "library")) +} + + +# 4 schools ---- +getPrimary <- function(layer) { + return( + layer %>% + rowwise() %>% + mutate(lowest_grade = unlist(strsplit(grades, "-"))[1]) %>% + ungroup() %>% + filter(amenity == "school" & + (as.numeric(lowest_grade) < 7 | + lowest_grade %in% c("P", "K") | + school %in% c("primary", "primary;secondary") | + isced_level %in% c("0", "0-1", "0-2", "0-3", "0;1", "0;2", "0;3", + "1", "1-2", "1-3", "1;2", "1;3") | + grepl("Primary", name)) & + # omit other types such as special_education_needs, prison + (school %in% c("primary", "primary;secondary") | is.na(school))) + ) +} + +getSecondary <- function(layer) { + return( + layer %>% + rowwise() %>% + # some are eg "5-8; 10-12' or '0-4;9' - first, find the grade after the + # last hyphen, then find the grade after the last semi=colon + mutate(highest_grade = dplyr::last(unlist(strsplit(grades, "-")))) %>% + mutate(highest_grade = dplyr::last(unlist(strsplit(highest_grade, ";")))) %>% + ungroup() %>% + filter(amenity == "school" & + (as.numeric(highest_grade) >= 7 | + school %in% c("secondary", "primary;secondary") | + isced_level %in% c("0-2", "0-3", "0;2", "0;3", + "1-2", "1-3", "1;2", "1;3", + "2", "3", "2-3", "2;3") | + (grepl("Secondary", name) | grepl("High ", name))) & # space after "High" to avoid eg "Highview Primary" + # omit other types such as special_education_needs, prison + (school %in% c("secondary", "primary;secondary") | is.na(school))) + ) +} + + +# 5 health ---- +getClinic <- function(layer) { + return(layer %>% + filter(amenity %in% c("clinic", "doctor", "doctors") | + healthcare %in% c("clinic", "doctor", "doctors"))) +} + +getDentist <- function(layer) { + return(layer %>% filter(amenity == "dentist" | healthcare == "dentist")) +} + +getPharmacy <- function(layer) { + return(layer %>% + filter(amenity %in% c("chemist", "pharmacy") | + healthcare %in% c("chemist", "pharmacy") | + shop %in% c("chemist", "pharmacy"))) +} + + +# 6 shopping ---- +getConvenience <- function(layer) { + return(layer %>% filter(shop == "convenience")) +} + +getSupermarket <- function(layer) { + return(layer %>% filter(shop == "supermarket")) +} + +getShop <- function(layer) { + return(layer %>% filter(!(is.na(shop)))) +} + +getPost <- function(layer) { + return(layer %>% filter(amenity == "post_office")) +} + +getBank <- function(layer) { + return(layer %>% filter(amenity == "bank")) +} + + +# 7 eating ---- +getRestaurant <- function(layer) { + return(layer %>% filter(amenity == "restaurant")) +} + +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() { + # 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/makeEdgesOneway.R b/functions/makeEdgesOneway.R new file mode 100644 index 0000000..66a2878 --- /dev/null +++ b/functions/makeEdgesOneway.R @@ -0,0 +1,70 @@ +# function to convert two-way edges to one-way + +makeEdgesOneway <- function(nodes_current, edges_current) { + + # testing + # nodes_current <- input.nodes + # edges_current <- input.links + + # ensure fromx, fromy, tox and toy column names are lower case (eg not 'fromX') + names.to.change <- c("fromX", "fromY", "toX", "toY") + edges_current <- rename_with(edges_current, tolower, any_of(names.to.change)) + + # select only two-way edges + edges_twoway <- edges_current %>% + filter(is_oneway == 0) + + # swap from/to details + edges_twoway_reversed <- edges_twoway %>% + # store original from/to details + mutate(orig_from_id = from_id, + orig_to_id = to_id, + orig_fromx = fromx, + orig_fromy = fromy, + orig_tox = tox, + orig_toy = toy) %>% + # swap from/to + mutate(from_id = orig_to_id, + to_id = orig_from_id, + fromx = orig_tox, + fromy = orig_toy, + tox = orig_fromx, + toy = orig_fromy) + + # if elevation is present, use the reverse slope + if("rvs_slope_pct" %in% colnames(edges_twoway_reversed)) { + edges_twoway_reversed <- edges_twoway_reversed %>% + mutate(slope_pct = rvs_slope_pct) + } + + # select required fields (excluding 'is_oneway') [note that "id" is not + # retained here - it is replaced by link_id] + required_fields <- c("from_id", "to_id", "fromx", "fromy", "tox", "toy", + "length", "freespeed", "permlanes", "capacity", "highway", + "cycleway", "surface", "is_cycle", "is_walk", "is_car", + "modes") + if ("slope_pct" %in% colnames(edges_twoway_reversed)) { + required_fields <- c(required_fields, "slope_pct") + } + edges_twoway_reversed <- edges_twoway_reversed %>% + dplyr::select(all_of(required_fields)) + + # modify original edges to rename fwd_slope_pct if present + if ("fwd_slope_pct" %in% colnames(edges_current)) { + edges_current <- edges_current %>% + rename(slope_pct = fwd_slope_pct) + } + + # exclude 'is_oneway' from original edges, and bind with reversed two-way edges + edges_current <- edges_current %>% + dplyr::select(all_of(required_fields)) %>% + rbind(., edges_twoway_reversed) + + # add link_id, based on rownumber (at the end, not beginning, because igraph + # requires from_id and to_id to be the first two columns) + edges_current <- edges_current %>% + mutate(link_id = row_number()) + + + return(list(nodes_current, edges_current)) +} From be526e916f2ec2697db9456cd1a7f17a5665ae57 Mon Sep 17 00:00:00 2001 From: StevePem Date: Wed, 16 Aug 2023 14:55:35 +1000 Subject: [PATCH 11/11] Add destinations to network --- Destinations.R | 169 -------------------------------- NetworkGenerator.R | 20 ++++ functions/addDestinations.R | 190 ++++++++++++++++++++++++++++++++++++ functions/writeOutputs.R | 10 ++ 4 files changed, 220 insertions(+), 169 deletions(-) delete mode 100644 Destinations.R create mode 100644 functions/addDestinations.R diff --git a/Destinations.R b/Destinations.R deleted file mode 100644 index d60c9d8..0000000 --- a/Destinations.R +++ /dev/null @@ -1,169 +0,0 @@ -# Make file of destinations required for accessibilty routing - -library(dplyr) -library(sf) -library(osmextract) - -# 1 Download OSM extract ---- -# -----------------------------------------------------------------------------# -# Download from https://www.interline.io/osm/extracts/ - -## Downloaded for Melbourne - melbourne_australia.osm.pbf - - -# 2 Converting to .gpkg format ---- -# -----------------------------------------------------------------------------# -# input file name and project CRS -INPUTFILE <- "./data/melbourne_australia.osm.pbf" - -PROJECT.CRS = 28355 - -# check layers -st_layers(INPUTFILE) - -# check keys -options(max.print = 2000) -polygon.tags <- oe_get_keys(INPUTFILE, layer = "multipolygons") %>% sort() -point.tags <- oe_get_keys(INPUTFILE, layer = "points") %>% sort() -line.tags <- oe_get_keys(INPUTFILE, layer = "lines") %>% sort() - -# create gpkg file in same directory as INPUTFILE, using the 'extra_tags' -# argument for specific extra tags required for various destination types -oe_vectortranslate(INPUTFILE, layer = "multipolygons", - extra_tags = c("access", "building", "grades", "healthcare", - "healthcare:speciality","isced:level", - "network", "operator", - "operator:type", "public_transport", "railway", - "school", "social_facility", "sport", - "tourism", "train")) -oe_vectortranslate(INPUTFILE, layer = "points", - 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")) -oe_vectortranslate(INPUTFILE, layer = "lines", - 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", - "smoothness", "surface")) -oe_vectortranslate(INPUTFILE, layer = "multilinestrings") -oe_vectortranslate(INPUTFILE, layer = "other_relations") - - -# 3 Read in the .gpkg file ---- -# -----------------------------------------------------------------------------# -GPKG <- "./data/melbourne_australia.gpkg" - -polygons <- st_read(GPKG, layer = "multipolygons") %>% st_transform(PROJECT.CRS) -points <- st_read(GPKG, layer = "points") %>% st_transform(PROJECT.CRS) -lines <- st_read(GPKG, layer = "lines") %>% st_transform(PROJECT.CRS) -multilines <- st_read(GPKG, layer = "multilinestrings") %>% st_transform(PROJECT.CRS) -other_relations <- st_read(GPKG, layer = "other_relations") %>% st_transform(PROJECT.CRS) - - -# 4 Extract required destinations ---- -# -----------------------------------------------------------------------------# - -## 4.1 Tag combinations for feature types and network ---- -## ----------------------------------------------------------------------------# -# load functions for locating specific feature types -source("./functions/getDestinationTypes.R") - -# load network -NETWORK <- "./output/test/network.sqlite" #<<< CHECK FINAL NAME -NODE.LAYER <- "nodes" -LINK.LAYER <- "links" - -network.nodes <- st_read(NETWORK, layer = NODE.LAYER) -network.links <- st_read(NETWORK, layer = LINK.LAYER) - - -## 4.2 Compile point and polygon destinations ---- -## ----------------------------------------------------------------------------# -destination.layer <- function(layer) { - return( - bind_rows( - getPlayground(layer) %>% mutate(dest_type = "playground"), - getPark(layer) %>% mutate(dest_type = "park"), - getSport(layer) %>% mutate(dest_type = "sport"), - getKindergarten(layer) %>% mutate(dest_type = "kindergarten"), - getCommunity(layer) %>% mutate(dest_type = "community_centre"), - getLibrary(layer) %>% mutate(dest_type = "library"), - getPrimary(layer) %>% mutate(dest_type = "primary_school"), - getSecondary(layer) %>% mutate(dest_type = "secondary_school"), - getClinic(layer) %>% mutate(dest_type = "health_clinic"), - getDentist(layer) %>% mutate(dest_type = "dentist"), - getPharmacy(layer) %>% mutate(dest_type = "pharmacy"), - getConvenience(layer) %>% mutate(dest_type = "convenience_store"), - getSupermarket(layer) %>% mutate(dest_type = "supermarket"), - getShop(layer) %>% mutate(dest_type = "shop"), - 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") - )) -} - -# create tables of destinations, and allocate unique id's (so features with -# multiple nodes can be grouped by the id where required) -destination.pt <- - bind_rows(destination.layer(points), - getStation() %>% mutate(dest_type = "railway_station")) %>% - mutate(dest_id = row_number()) - -destination.poly <- - destination.layer(polygons) %>% - mutate(dest_id = max(destination.pt$dest_id) + row_number()) - - -## 4.3 Find relevant nodes ---- -## ----------------------------------------------------------------------------# -# TO CONFIRM:- -# For all destinations except parks and schools ('small features'), relevant -# node is nearest node to point or to polygon centroid -# For parks and schools ('large features'): -# - points are buffered to 50m to create a polygon feature, -# - for buffered points and polygons, relevant nodes are all nodes within the -# feature and terminal nodes of links within 30m of boundary, or if none, -# then nearest node to boundary - -# Maybe this should be all nodes within 30m of buffered feature, and if link is within -# 30m of boundary but doesn't have a node within the buffer, then also its closest terminal -# node ??? - -dest.small <- bind_rows(destination.pt, - destination.poly %>% st_centroid()) %>% - filter(!(dest_type %in% c("park", "primary_school", "secondary_school"))) -near_node <- network.nodes$id[st_nearest_feature(dest.small, network.nodes)] -dest.small.with.nodes <- cbind(dest.small %>% st_drop_geometry(), near_node) - - -## NOTE - the code below is a simplified version which just finds nodes within -## feature or its 30m buffer, or nearest node if none - doesn't extend to terminal -## nodes of nearby features -dest.large <- bind_rows(destination.pt %>% st_buffer(50), - destination.poly) %>% - filter(dest_type %in% c("park", "primary_school", "secondary_school")) - -dest.large.found.nodes <- dest.large %>% - st_buffer(30) %>% - st_intersection(., network.nodes %>% dplyr::select(near_node = id)) - -dest.large.need.nodes <- dest.large %>% - filter(!(dest_id %in% dest.large.found.nodes$dest_id)) -near_node <- network.nodes$id[st_nearest_feature(dest.large.need.nodes, network.nodes)] - -dest.large.with.nodes <- bind_rows(dest.large.found.nodes %>% st_drop_geometry(), - cbind(dest.large.need.nodes %>% st_drop_geometry(), - near_node)) - -dest.with.nodes <- bind_rows(dest.small.with.nodes, - dest.large.with.nodes) - diff --git a/NetworkGenerator.R b/NetworkGenerator.R index 0f25329..f8be832 100644 --- a/NetworkGenerator.R +++ b/NetworkGenerator.R @@ -36,6 +36,12 @@ makeNetwork<-function(outputFileName="test"){ demFile= "./data/DEM_melbourne.tif" # DEM's multiplier- set to 1 if DEM contains actual elevation ElevationMultiplier=1 + + # DESTINATIONS + # A flag for whether to add a destinations layer (drawn from OSM) or not + addDestinationLayer=T + # OSM extract for destinations, in .osm.pbf format + osmPbfExtract="./data/melbourne_australia.osm.pbf" # GTFS addGtfs=F @@ -67,6 +73,8 @@ makeNetwork<-function(outputFileName="test"){ library(tidytransit) library(hablar) library(hms) + library(osmextract) + library(tidyr) # Building the output folder structure ------------------------------------ @@ -220,6 +228,14 @@ makeNetwork<-function(outputFileName="test"){ networkDensified <- densifyNetwork(networkConnected,desnificationMaxLengh, densifyBikeways) + # adding destinations layer + if (addDestinationLayer) { + destinations <- addDestinations(networkDensified[[1]], + networkDensified[[2]], + osmPbfExtract, + outputCrs) + } + # simplify geometry so all edges are straight lines system.time(networkDirect <- makeEdgesDirect(networkDensified[[1]], @@ -278,6 +294,10 @@ makeNetwork<-function(outputFileName="test"){ networkFinal <- networkOneway + if (addDestinationLayer) { + networkFinal[[3]] <- destinations + } + # writing outputs echo("========================================================\n") echo("| **Launching Output Writing** |\n") diff --git a/functions/addDestinations.R b/functions/addDestinations.R new file mode 100644 index 0000000..a0abf38 --- /dev/null +++ b/functions/addDestinations.R @@ -0,0 +1,190 @@ +# 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' + +addDestinations <- function(nodes_current, + edges_current, + osmPbfExtract, + outputCRS) { + + # nodes_current = networkDensified[[1]] + # edges_current = networkDensified[[2]] + # osmPbfExtract = "./data/melbourne_australia.osm.pbf" + # outputCrs = 28355 + + # # check layers + # st_layers(osmPbfExtract) + # # only multipolygons, points and lines are required (not multilinestrings + # # or other_relations) + + # # 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() + + # reading layers ---- + # ----------------------------------# + echo("Reading in the .osm.pbf extract 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", "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) + + + # 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 + + destination.layer <- function(layer) { + return( + bind_rows( + getPlayground(layer) %>% mutate(dest_type = "playground"), + getPark(layer) %>% mutate(dest_type = "park"), + getSport(layer) %>% mutate(dest_type = "sport"), + getKindergarten(layer) %>% mutate(dest_type = "kindergarten"), + getCommunity(layer) %>% mutate(dest_type = "community_centre"), + getLibrary(layer) %>% mutate(dest_type = "library"), + getPrimary(layer) %>% mutate(dest_type = "primary_school"), + getSecondary(layer) %>% mutate(dest_type = "secondary_school"), + getClinic(layer) %>% mutate(dest_type = "health_clinic"), + getDentist(layer) %>% mutate(dest_type = "dentist"), + getPharmacy(layer) %>% mutate(dest_type = "pharmacy"), + getConvenience(layer) %>% mutate(dest_type = "convenience_store"), + getSupermarket(layer) %>% mutate(dest_type = "supermarket"), + getShop(layer) %>% mutate(dest_type = "shop"), + 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") + )) + } + + # 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) + destination.pt <- + bind_rows(destination.layer(points), + # add stations (from point, polygons and lines) to point table + getStation() %>% mutate(dest_type = "railway_station")) %>% + mutate(dest_id = row_number()) + + destination.poly <- + destination.layer(polygons) %>% + mutate(dest_id = max(destination.pt$dest_id) + row_number()) + + # # check numbers of each destination type + # chk <- full_join(destination.poly %>% + # st_drop_geometry() %>% + # group_by(dest_type) %>% + # summarise(poly = n()), + # destination.pt %>% + # st_drop_geometry() %>% + # group_by(dest_type) %>% + # summarise(pt = n()), + # by = "dest_type") + + + # find relevant nodes ---- + # For all destinations except parks and schools ('small features'), relevant + # node is nearest node to point or to polygon centroid + # For parks and schools ('large features'): + # - points are buffered to 50m to create a polygon feature, + # - for buffered points and polygons, relevant nodes are all nodes within the + # feature and terminal nodes of links within 30m of boundary, or if none, + # then nearest node to boundary + # In each case, nodes/links must be cyclable + + cyclable.links <- edges_current %>% + filter(is_cycle == 1) + cyclable.nodes <- nodes_current %>% + filter(id %in% cyclable.links$from_id | id %in% cyclable.links$to_id) + + # 'small' destinations + dest.small <- bind_rows(destination.pt, + destination.poly %>% st_centroid()) %>% + filter(!(dest_type %in% c("park", "primary_school", "secondary_school"))) + near_node <- cyclable.nodes$id[st_nearest_feature(dest.small, cyclable.nodes)] + dest.small.with.nodes <- cbind(dest.small %>% st_drop_geometry(), near_node) + + + # 'large' destinations + dest.large <- bind_rows(destination.pt %>% st_buffer(50), + destination.poly) %>% + filter(dest_type %in% c("park", "primary_school", "secondary_school")) + + # # - nodes within the feature + # dest.large.nodes.within <- dest.large %>% + # st_intersection(., cyclable.nodes %>% dplyr::select(near_node = id)) %>% + # st_drop_geometry() + + # - terminal nodes of links within feature buffered to 30m (will include any + # nodes within feature itself, as their links will fall within the buffered + # feature) + dest.large.found.nodes <- dest.large %>% + st_buffer(30) %>% + st_intersection(., cyclable.links %>% dplyr::select(from_id, to_id)) %>% + st_drop_geometry() %>% + pivot_longer(cols = c("from_id", "to_id"), + names_to = NULL, + values_to = "near_node") %>% + distinct() + + # - nearest node if none within and no links within 30m + dest.large.other <- dest.large %>% + filter(!(dest_id %in% dest.large.found.nodes$dest_id)) + near_node <- cyclable.nodes$id[st_nearest_feature(dest.large.other, cyclable.nodes)] + dest.large.other.nodes <- cbind(dest.large.other %>% st_drop_geometry(), near_node) + + # combine the large destinations + dest.large.with.nodes <- bind_rows(dest.large.found.nodes, + dest.large.other.nodes) + + + # combine all destinations for output ---- + dest.with.nodes <- bind_rows(dest.small.with.nodes, + dest.large.with.nodes) %>% + relocate(dest_id) %>% + relocate(dest_type, .after = dest_id) %>% + relocate(near_node, .after = dest_type) %>% + relocate(other_tags, .after = last_col()) %>% + + # and join nodes for locations + left_join(., nodes_current %>% dplyr::select(id), by = c("near_node" = "id")) + + return(dest.with.nodes) + +} + diff --git a/functions/writeOutputs.R b/functions/writeOutputs.R index 672a8d1..e82d7e5 100644 --- a/functions/writeOutputs.R +++ b/functions/writeOutputs.R @@ -30,6 +30,11 @@ exportSQlite <- function(networkFinal, outputDir, outputCrs){ st_write(networkFinal[[1]], paste0(outputDir,'/network.sqlite'), layer = 'nodes', driver = 'SQLite', layer_options = 'GEOMETRY=AS_XY', delete_layer = T) + if (length(networkFinal) > 2) { + st_write(networkFinal[[3]], paste0(outputDir,'/network.sqlite'), + layer = 'destinations', driver = 'SQLite', + layer_options = 'GEOMETRY=AS_XY', delete_layer = T) + } echo(paste0('Finished generating the sqlite output\n')) } @@ -64,6 +69,11 @@ exportShp <- function(networkFinal, outputDir, outputCrs){ st_write(networkFinal[[1]], paste0(shpDir,'/nodes.shp'), driver = "ESRI Shapefile", layer_options = 'GEOMETRY=AS_XY', delete_layer = T) + if (length(networkFinal) > 2) { + st_write(networkFinal[[3]], paste0(shpDir,'/destinations.shp'), + driver = "ESRI Shapefile", layer_options = 'GEOMETRY=AS_XY', + delete_layer = T) + } echo(paste0('Finished generating the ShapeFile output\n'))