diff --git a/R/rnet_join.R b/R/rnet_join.R
index c4e25b2b..1bbed492 100644
--- a/R/rnet_join.R
+++ b/R/rnet_join.R
@@ -82,7 +82,8 @@
 #' @export
 rnet_join = function(rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1,
                      subset_x = TRUE, dist_subset = NULL, segment_length = 0,
-                     endCapStyle = "FLAT", contains = TRUE, ...) {
+                     endCapStyle = "SQUARE", contains = FALSE, 
+                     mask_function = function(angle) (angle < 40) | (angle > 140),...) {
   if (is.null(dist_subset)) {
     dist_subset = dist + 1
   }
@@ -90,6 +91,9 @@ rnet_join = function(rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1,
     rnet_x = rnet_subset(rnet_x, rnet_y, dist = dist_subset, ...)
   }
   rnet_x_buffer = geo_buffer(rnet_x, dist = dist, nQuadSegs = 2, endCapStyle = endCapStyle)
+  # Store the original geometry of 'rnet_x' in the buffer object
+  rnet_x_buffer$corr_line_geometry_buffer = rnet_x$geometry
+
   if (segment_length > 0) {
     rnet_y = line_segment(rnet_y, segment_length = segment_length)
   }
@@ -98,7 +102,7 @@ rnet_join = function(rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1,
   }
   # browser()
   if (contains) {
-    rnetj = sf::st_join(rnet_x_buffer[key_column], rnet_y, join = sf::st_contains)
+    rnetj = sf::st_join(rnet_x_buffer, rnet_y, join = sf::st_contains)
     # # For debugging:
     # library(tmap)
     # tmap_mode("view")
@@ -106,9 +110,40 @@ rnet_join = function(rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1,
     #   qtm(osm_net_example)
   } else {
     rnet_y_centroids = sf::st_centroid(rnet_y)
-    rnetj = sf::st_join(rnet_x_buffer[key_column], rnet_y_centroids)
-  }
+    # Store the original geometry of 'rnet_y' in the centroid object
+    rnet_y_centroids$corr_line_geometry_point = rnet_y$geometry
+    rnetj = sf::st_join(rnet_x_buffer, rnet_y_centroids)
+    # Calculate angles between the buffer geometry and the point geometry for each row
+    rnetj$angle = sapply(1:nrow(rnetj), function(i) {
+    calculate_angle(get_vector(rnetj$corr_line_geometry_buffer[[i]]), get_vector(rnetj$corr_line_geometry_point[[i]]))
+    })  
+    # Using angle_diff fun to replace calculate_angle/get_vector
+    # rnetj$angle_2 = sapply(1:nrow(rnetj), function(i) {
+    #   # Check if either of the geometries is empty
+    #   if (st_is_empty(rnetj$corr_line_geometry_buffer[i]) || 
+    #       st_is_empty(rnetj$corr_line_geometry_point[i])) {
+    #     return(NA)
+    #   }
+      
+    #   # Extract the two lines into a temporary 'sf' object
+    #   temp_sf <- rbind(
+    #     data.frame(geometry = rnetj$corr_line_geometry_buffer[i]), 
+    #     data.frame(geometry = rnetj$corr_line_geometry_point[i])
+    #   )
+    #   temp_sf <- st_as_sf(temp_sf)
+      
+    #   # Get the angles difference using the temporary 'sf' object
+    #   angles = angle_diff(temp_sf, angle = 0)
+      
+    #   # Return the absolute difference between the two angles
+    #   abs(angles[1] - angles[2])
+    # })
 
+    # Filter rows based on the angle values
+    mask = mask_function(rnetj$angle)
+    filtered_rnetj <- rnetj[mask, ]
+    rnetj = filtered_rnetj
+  }
   rnetj
 }
 
@@ -262,3 +297,27 @@ rnet_merge <- function(rnet_x, rnet_y, dist = 5, funs = NULL, sum_flows = TRUE,
   }
   res_sf
 }
+
+
+get_vector <- function(line) {
+  if (sf::st_is_empty(line)) {
+    # warning("Encountered an empty LINESTRING. Returning NULL.")
+    return(NULL)
+  }
+  coords <- sf::st_coordinates(line)
+  # Check if coords is empty or has insufficient dimensions
+  if (is.null(coords) || nrow(coords) < 2 || ncol(coords) < 2) {
+    stop("Insufficient coordinate data")
+  }
+  start <- coords[1, 1:2]
+  end <- coords[2, 1:2]
+  return(c(end[1] - start[1], end[2] - start[2]))
+}
+calculate_angle <- function(vector1, vector2) {
+  dot_product <- sum(vector1 * vector2)
+  magnitude_product <- sqrt(sum(vector1^2)) * sqrt(sum(vector2^2))
+  cos_angle <- dot_product / magnitude_product
+  angle <- acos(cos_angle) * (180 / pi)
+  return(angle)
+}
+ 		
diff --git a/man/rnet_join.Rd b/man/rnet_join.Rd
index 8c55ac19..f375c34d 100644
--- a/man/rnet_join.Rd
+++ b/man/rnet_join.Rd
@@ -13,8 +13,9 @@ rnet_join(
   subset_x = TRUE,
   dist_subset = NULL,
   segment_length = 0,
-  endCapStyle = "FLAT",
-  contains = TRUE,
+  endCapStyle = "SQUARE",
+  contains = FALSE,
+  mask_function = function(angle) (angle < 40) | (angle > 140),
   ...
 )
 }
diff --git a/vignettes/merging-route-networks 3km.Rmd b/vignettes/merging-route-networks 3km.Rmd
new file mode 100644
index 00000000..b2ed9d14
--- /dev/null
+++ b/vignettes/merging-route-networks 3km.Rmd	
@@ -0,0 +1,120 @@
+---
+title: "Merging route networks"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Merging route networks}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  # # Uncomment to speed-up build
+  # eval = FALSE,
+  comment = "#>",
+  echo = FALSE,
+  message = FALSE,
+  warning = FALSE
+)
+# devtools::load_all()
+sf::sf_use_s2(FALSE)
+```
+
+
+```{r setup}
+library(stplanr)
+library(dplyr)
+library(tmap)
+library(sf)
+
+rnet_x = sf::read_sf("https://github.com/nptscot/networkmerge/releases/download/v0.1/OS_3km.geojson")
+rnet_y = sf::read_sf("https://github.com/nptscot/npt/releases/download/rnet_3km_buffer/rnet_3km_buffer.geojson")
+st_crs(rnet_x) <- st_crs(rnet_y)
+crs_y <- st_crs(rnet_y)
+rnet_x <- st_transform(rnet_x, crs_y)
+
+# Check CRS for rnet_x
+crs_x <- st_crs(rnet_x)
+print(crs_x)
+
+# Check CRS for rnet_y
+crs_y <- st_crs(rnet_y)
+print(crs_y)
+
+```
+
+# Target network preprocessing
+
+We pre-processed the input simple geometry to make it even simpler as shown below.
+```{r}
+dim(rnet_x)
+rnet_x = rnet_subset(rnet_x, rnet_y, dist = 20)
+dim(rnet_x)
+rnet_x = rnet_subset(rnet_x, rnet_y, dist = 20, min_length = 5)
+dim(rnet_y)
+```
+
+```{r}
+# Extract column names from the rnet_x data frame
+name_list <- names(rnet_y)
+name_list
+# Initialize an empty list
+funs <- list()
+
+# Loop through each name and assign it a function based on specific conditions
+for (name in name_list) {
+  if (name == "geometry") {
+    next  # Skip the current iteration
+  } else if (name %in% c("Gradient", "Quietness")) {
+    funs[[name]] <- mean
+  } else {
+    funs[[name]] <- sum
+  }
+}
+
+tmap_mode("view")
+brks = c(0, 100, 500, 1000, 5000,10000)
+# funs = list(all_fastest_bicycle_go_dutch = sum, Quietness = mean)
+rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 10, funs = funs)
+
+
+dim(rnet_y)
+summary(rnet_y$all_fastest_bicycle)
+dim(rnet_merged)
+summary(rnet_merged$all_fastest_bicycle)
+
+st_write(rnet_merged, "G:/rnet_merged.geojson", driver = "GeoJSON")
+#rnet_merged$identifier <- NULL
+# rnet_merged$value[is.na(rnet_merged$value)] <- 0
+m1 = tm_shape(rnet_y) + tm_lines(col ="all_fastest_bicycle_go_dutch", palette = "viridis", lwd = 5, breaks = brks)
+m2 = tm_shape(rnet_merged) + tm_lines(col ="all_fastest_bicycle_go_dutch", palette = "viridis", lwd = 5, breaks = brks)
+tmap_arrange(m1, m2, sync = TRUE)
+```
+
+We can more reduce the minimum segment length to ensure fewer NA values in the outputs:
+****
+```{r}
+rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 20, segment_length = 10, funs = funs)
+m1 = tm_shape(rnet_y) + tm_lines("all_fastest_bicycle", palette = "viridis", lwd = 5, breaks = brks)
+m2 = tm_shape(rnet_merged) + tm_lines("all_fastest_bicycle", palette = "viridis", lwd = 5, breaks = brks)
+tmap_arrange(m1, m2, sync = TRUE)
+```
+
+Reducing the max length of the complex route network led to the following result.
+
+```{r}
+rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 20, segment_length = 5, funs = funs, dist_subset = 30)
+# Drop the 'identifier' column
+# rnet_merged <- rnet_merged[ , 'value', drop = FALSE]
+m1 = tm_shape(rnet_y) + tm_lines("all_fastest_bicycle", palette = "viridis", lwd = 5, breaks = brks)
+m2 = tm_shape(rnet_merged) + tm_lines("all_fastest_bicycle", palette = "viridis", lwd = 5, breaks = brks)
+tmap_arrange(m1, m2, sync = TRUE)
+
+```
+
+```{r, echo=FALSE, eval=FALSE}
+sum(rnet_merged$all_fastest_bicycle_go_dutch * sf::st_length(rnet_merged), na.rm = TRUE)
+sum(rnet_y$all_fastest_bicycle_go_dutch * sf::st_length(rnet_y), na.rm = TRUE)
+```
+
diff --git a/vignettes/merging-route-networks.Rmd b/vignettes/merging-route-networks.Rmd
index acfedbff..e698fb20 100644
--- a/vignettes/merging-route-networks.Rmd
+++ b/vignettes/merging-route-networks.Rmd
@@ -34,6 +34,20 @@ rnet_y = sf::read_sf("https://github.com/ropensci/stplanr/releases/download/v1.0
 # sf::write_sf(rnet_x, "~/github/ropensci/stplanr/rnet_x_ed.geojson", delete_dsn = TRUE)
 ```
 
+```{r}
+# Checking reults with adding funtions to calcualte the angle
+tmap_mode("view")
+funs = list(value = sum, Quietness = mean)
+brks = c(0, 100, 500, 1000, 5000,10000)
+rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 10, funs = funs)
+#rnet_merged$identifier <- NULL
+# rnet_merged$value[is.na(rnet_merged$value)] <- 0
+m1 = tm_shape(rnet_y) + tm_lines("value", palette = "viridis", lwd = 5, breaks = brks)
+m2 = tm_shape(rnet_merged) + tm_lines("value", palette = "viridis", lwd = 5, breaks = brks)
+tmap_arrange(m1, m2, sync = TRUE)
+
+```
+
 # Target network preprocessing
 
 We pre-processed the input simple geometry to make it even simpler as shown below.
@@ -90,4 +104,3 @@ tmap_arrange(m1, m2, sync = TRUE)
 sum(rnet_merged$value * sf::st_length(rnet_merged), na.rm = TRUE)
 sum(rnet_y$value * sf::st_length(rnet_y), na.rm = TRUE)
 ```
-