diff --git a/functions/addDestinations.R b/functions/addDestinations.R
index 40bd573..79f4a93 100644
--- a/functions/addDestinations.R
+++ b/functions/addDestinations.R
@@ -81,6 +81,7 @@ addDestinations <- function(nodes_current,
   # create tables for points and polygons, allocate unique id's (so features 
   # multiple multiple nodes can be grouped by the id where required),
   # and store area and location details
+  echo("Destination point features\n")
   destination.pt <- 
     bind_rows(destination.layer(points),
               
@@ -93,6 +94,7 @@ addDestinations <- function(nodes_current,
            centroid_x = st_coordinates(.)[, 1],
            centroid_y = st_coordinates(.)[, 2])
   
+  echo("Destination polygon features\n")
   destination.poly <- 
     destination.layer(polygons) %>%
     filter(st_is_valid(geom)) %>%
@@ -117,70 +119,219 @@ addDestinations <- function(nodes_current,
   # 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)
+  # Find nodes available for each mode
+  links.cycle <- edges_current %>% filter(is_cycle == 1)
+  nodes.cycle <- nodes_current %>% filter(id %in% links.cycle$from_id | id %in% links.cycle$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) 
+  links.walk <- edges_current %>% filter(is_walk == 1)
+  nodes.walk <- nodes_current %>% filter(id %in% links.walk$from_id | id %in% links.walk$to_id)
   
+  links.car <- edges_current %>% filter(is_car == 1)
+  nodes.car <- nodes_current %>% filter(id %in% links.car$from_id | id %in% links.car$to_id)
+  
+  # 'small' features
+  dest.small <- bind_rows(destination.pt, destination.poly) %>%
+    filter(!(dest_type %in% c("park", "primary_school", "secondary_school")))
+  # nearest node of each mode
+  node_cycle <- nodes.cycle$id[st_nearest_feature(dest.small %>% st_centroid(), nodes.cycle)]
+  node_walk <- nodes.walk$id[st_nearest_feature(dest.small %>% st_centroid(), nodes.walk)]
+  node_car <- nodes.car$id[st_nearest_feature(dest.small %>% st_centroid(), nodes.car)]
+  # join nearest nodes to features
+  dest.small.with.nodes <- cbind(dest.small, node_cycle, node_walk, node_car) 
   
-  # 'large' destinations
-  dest.large <- bind_rows(destination.pt %>% st_buffer(50),
-                          destination.poly) %>%
+  # 'large' features - points
+  dest.large.pt <- destination.pt %>%
+    filter(dest_type %in% c("park", "primary_school", "secondary_school"))
+  # nearest node of each mode
+  node_cycle <- nodes.cycle$id[st_nearest_feature(dest.large.pt %>% st_centroid(), nodes.cycle)]
+  node_walk <- nodes.walk$id[st_nearest_feature(dest.large.pt %>% st_centroid(), nodes.walk)]
+  node_car <- nodes.car$id[st_nearest_feature(dest.large.pt %>% st_centroid(), nodes.car)]
+  # join nearest nodes to features
+  dest.large.pt.with.nodes <- cbind(dest.large.pt, node_cycle, node_walk, node_car) 
+
+  # 'large' features - polygons
+  dest.large.poly <- 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)) %>%
+  echo(paste("Finding entry nodes for", nrow(dest.large.poly), "parks and schools\n"))
+  
+  # internal nodes
+  internal.nodes.cycle <- nodes.cycle %>%
+    st_intersection(dest.large.poly) %>%
     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) 
+    dplyr::select(id, dest_id)
+  internal.nodes.walk <- nodes.walk %>%
+    st_intersection(dest.large.poly) %>%
+    st_drop_geometry() %>%
+    dplyr::select(id, dest_id)
+  internal.nodes.car <- nodes.car %>%
+    st_intersection(dest.large.poly) %>%
+    st_drop_geometry() %>%
+    dplyr::select(id, dest_id)
+  
+  # pseudo entry points
+  # buffered links
+  links.cycle.buffered <- st_buffer(links.cycle, 30)
+  links.walk.buffered <- st_buffer(links.walk, 30)
+  links.car.buffered <- st_buffer(links.car, 30)
+  
+  # setup for parallel processing and progress reporting
+  cores <- detectCores()
+  cluster <- parallel::makeCluster(cores)
+  doSNOW::registerDoSNOW(cluster)
+  pb <- txtProgressBar(max = nrow(dest.large.poly), style = 3)
+  progress <- function(n) setTxtProgressBar(pb, n)
+  opts <- list(progress = progress)
+  
+  # report
+  echo(paste("Finding boundary points for", nrow(dest.large.poly), 
+             "parks and schools; parallel processing with", cores, "cores\n"))
+  
+  # loop to find list of boundary points
+  boundary.points <-
+    foreach(i = 1:nrow(dest.large.poly),
+            # foreach(i = 1:8,
+            # .combine = rbind,
+            .packages = c("dplyr", "sf"),
+            .options.snow = opts) %dopar% {
+              
+              dest <- dest.large.poly[i,]
+              
+              dest.boundary.points <- dest %>% 
+                # convert destination polygons boundaries to linestring
+                st_cast(to = "MULTILINESTRING") %>%
+                st_cast(to = "LINESTRING") %>%
+                # locate points at 20m along boundaries
+                st_line_sample(., density = units::set_units(20, m))
+              
+              return(dest.boundary.points)
+            }
+  
+  # convert the list of boundary points into a dataframe
+  # extract coordinates and ids (row numbers, corresponding to dest.large.poly rows)
+  coordinates <- lapply(boundary.points, function(x) st_coordinates(x[[1]]))
+  ids <- seq_along(boundary.points)
+  # make dataframe of ids and x and y coordinates
+  boundary.df <- data.frame(
+    id = rep(ids, sapply(coordinates, nrow)),
+    x = unlist(lapply(coordinates, function(x) x[,1])),
+    y = unlist(lapply(coordinates, function(x) x[,2])))
+  # convert to sf object
+  boundary.sf <- st_as_sf(boundary.df, coords = c("x", "y")) %>%
+    st_set_crs(outputCrs)
+  
+  # pseudo entry points with nodes
+  pseudo.entry.nodes.cycle <- boundary.sf %>%
+    # keep only those within 30m of cyclable links
+    st_filter(links.cycle.buffered, .predicate = st_intersects) %>%
+    # nearest nodes to the entry points
+    mutate(node_cycle = nodes.cycle$id[st_nearest_feature(., nodes.cycle)]) %>%
+    st_drop_geometry()
+  
+  pseudo.entry.nodes.walk <- boundary.sf %>%
+    # keep only those within 30m of walkable links
+    st_filter(links.walk.buffered, .predicate = st_intersects) %>%
+    # nearest nodes to the entry points
+    mutate(node_walk = nodes.walk$id[st_nearest_feature(., nodes.walk)]) %>%
+    st_drop_geometry()
+  
+  pseudo.entry.nodes.car <- boundary.sf %>%
+    # keep only those within 30m of driveable links
+    st_filter(links.car.buffered, .predicate = st_intersects) %>%
+    # nearest nodes to the entry points
+    mutate(node_car = nodes.car$id[st_nearest_feature(., nodes.car)]) %>%
+    st_drop_geometry()
+  
+  # loop through large destination polygons, find their internal and pseudo entry
+  # nodes (and fallback nearest node), and combine
+  echo(paste("Assembling entry points for", nrow(dest.large.poly), 
+             "parks and schools; parallel processing with", cores, "cores\n"))
+  
+  # loop to assemble entry nodes
+  dest.large.poly.with.nodes <-
+    foreach(i = 1:nrow(dest.large.poly),
+            # foreach(i = 1:8,
+            .combine = rbind,
+            .packages = c("dplyr"),
+            .options.snow = opts) %dopar% {
+              
+              dest <- dest.large.poly[i,]
+              
+              # nodes within features (identifier is 'dest_id')
+              dest.internal.nodes.cycle <- internal.nodes.cycle %>%
+                filter(dest_id == dest$dest_id) %>%
+                .$id
+              dest.internal.nodes.walk <- internal.nodes.walk %>%
+                filter(dest_id == dest$dest_id) %>%
+                .$id
+              dest.internal.nodes.car <- internal.nodes.car %>%
+                filter(dest_id == dest$dest_id) %>%
+                .$id
+               
+              # pseudo entry nodes (identifier is 'id', which is row number)
+              dest.pseudo.entry.nodes.cycle <- pseudo.entry.nodes.cycle %>%
+                filter(id == i) %>%
+                .$node_cycle
+              dest.pseudo.entry.nodes.walk <- pseudo.entry.nodes.walk %>%
+                filter(id == i) %>%
+                .$node_walk
+              dest.pseudo.entry.nodes.car <- pseudo.entry.nodes.car %>%
+                filter(id == i) %>%
+                .$node_car
+              
+              # combined nodes
+              entry.nodes.cycle <- unique(c(dest.internal.nodes.cycle, dest.pseudo.entry.nodes.cycle))
+              entry.nodes.walk <- unique(c(dest.internal.nodes.walk, dest.pseudo.entry.nodes.walk))
+              entry.nodes.car <- unique(c(dest.internal.nodes.car, dest.pseudo.entry.nodes.car))
+              
+              # fallback if no internal or pseudo entry nodes
+              if (length(entry.nodes.cycle) == 0) {
+                entry.nodes.cycle <- nodes.cycle$id[st_nearest_feature(dest, nodes.cycle)]
+              }
+              if (length(entry.nodes.walk) == 0) {
+                entry.nodes.walk <- nodes.walk$id[st_nearest_feature(dest, nodes.walk)]
+              }
+              if (length(entry.nodes.car) == 0) {
+                entry.nodes.car <- nodes.car$id[st_nearest_feature(dest, nodes.car)]
+              }
+              
+              # convert the entry nodes to strings, and add to the dest row
+              dest$node_cycle <- toString(entry.nodes.cycle)
+              dest$node_walk <- toString(entry.nodes.walk)
+              dest$node_car <- toString(entry.nodes.car)
+              
+              return(dest)
+            }
   
-  # combine the large destinations
-  dest.large.with.nodes <- bind_rows(dest.large.found.nodes,
-                                     dest.large.other.nodes)
+  # close the progress bar and cluster
+  close(pb)
+  stopCluster(cluster)
   
   
   # combine all destinations for output ----
-  dest.with.nodes <- bind_rows(dest.small.with.nodes,
-                               dest.large.with.nodes) %>%
+  dest.with.nodes <- bind_rows(dest.small.with.nodes %>%
+                                 mutate(node_cycle = as.character(node_cycle),
+                                        node_walk = as.character(node_walk),
+                                        node_car = as.character(node_car)),
+                               dest.large.pt.with.nodes %>%
+                                 mutate(node_cycle = as.character(node_cycle),
+                                        node_walk = as.character(node_walk),
+                                        node_car = as.character(node_car)),
+                               dest.large.poly.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"))
-    
+    relocate(node_cycle, .after = dest_type) %>%
+    relocate(node_walk, .after = node_cycle) %>%
+    relocate(node_car, .after = node_walk) %>%
+    relocate(other_tags, .after = last_col())
+  
   return(dest.with.nodes)
   
 }
diff --git a/functions/getPTStops.R b/functions/getPTStops.R
index a9f6b0a..9a57a32 100644
--- a/functions/getPTStops.R
+++ b/functions/getPTStops.R
@@ -9,6 +9,8 @@ getPTStops <- function(city, gtfs_feed, outputCrs, region, regionBufferDist) {
   # region = "../data/processed/greater_melbourne.sqlite"
   # regionBufferDist = 10000
   
+  echo("Reading in GTFS data to find public transport stop locations\n")
+  
   # read in GTFS feed
   gtfs <- read_gtfs(gtfs_feed) %>%
     gtfs_as_sf(., crs = 4326)
@@ -41,7 +43,7 @@ getPTStops <- function(city, gtfs_feed, outputCrs, region, regionBufferDist) {
   } else if (!all(route_types %in% c("0", "1", "2", "3", "4", "5", "6", "7", "11", "12"))) {
     message("GTFS Feed contains the following route type codes: ", paste(route_types, collapse = ", "), ". Unable to process these using 
 the standard route codes from https://developers.google.com/transit/gtfs/reference, which are:
-   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12 monorail. 
+   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12-monorail. 
 Edit getPTStops.R to specify the meanings of the codes used in the GTFS Feed.
 PT stops will not be included in destinations.")
     stops.found = FALSE
@@ -49,7 +51,7 @@ PT stops will not be included in destinations.")
   } else {
     message("GTFS Feed contains the following route type codes: ", paste(route_types, collapse = ", "), ".
 Using standard route_type codes from https://developers.google.com/transit/gtfs/reference:
-   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12 monorail. 
+   0-tram, 1-metro, 2-train, 3-bus, 4-ferry, 5-cable tram, 6-cable car, 7-funicular, 11-trolleybus, 12-monorail. 
 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(
diff --git a/functions/writeOutputs.R b/functions/writeOutputs.R
index e6ecf40..60165cd 100644
--- a/functions/writeOutputs.R
+++ b/functions/writeOutputs.R
@@ -70,7 +70,13 @@ exportShp <- function(networkFinal, outputDir, outputCrs){
            driver = "ESRI Shapefile", layer_options = 'GEOMETRY=AS_XY', 
            delete_layer = T)
   if (length(networkFinal) > 2) {
-    st_write(networkFinal[[3]], paste0(shpDir,'/destinations.shp'),
+    message("When writing destinations to shapefile, long 'other_tags' may be truncated; consider using sqlite instead.")
+    dest.pt <- st_collection_extract(networkFinal[[3]], type = "POINT")
+    dest.poly <- st_collection_extract(networkFinal[[3]], type = "POLYGON")
+    st_write(dest.pt, paste0(shpDir,'/destinations_point.shp'),
+             driver = "ESRI Shapefile", layer_options = 'GEOMETRY=AS_XY', 
+             delete_layer = T)
+    st_write(dest.poly, paste0(shpDir,'/destinations_polygon.shp'),
              driver = "ESRI Shapefile", layer_options = 'GEOMETRY=AS_XY', 
              delete_layer = T)
   }