From 1dbba5f01d5db5131e76af91b694a66b4b0c20b7 Mon Sep 17 00:00:00 2001 From: wangzhao0217 <74598734+wangzhao0217@users.noreply.github.com> Date: Mon, 4 Sep 2023 21:54:09 +0100 Subject: [PATCH 1/4] update funtions get_vector and calculate_angle --- R/rnet_join.R | 37 +++++++++++++++++++++++++++++++++---- 1 file changed, 33 insertions(+), 4 deletions(-) diff --git a/R/rnet_join.R b/R/rnet_join.R index d97b02e1..e7ca0118 100644 --- a/R/rnet_join.R +++ b/R/rnet_join.R @@ -106,6 +106,7 @@ 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) + rnet_y_centroids$corr_line_geometry = rnet_y$geometry rnetj = sf::st_join(rnet_x_buffer[key_column], rnet_y_centroids) } @@ -250,15 +251,43 @@ rnet_merge <- function(rnet_x, rnet_y, dist = 5, funs = NULL, sum_flows = TRUE, res_sf = dplyr::left_join(rnet_x, res_df) if (sum_flows) { res_sf$length_x = as.numeric(sf::st_length(res_sf)) + # Calculate the angle between 'corr_line_geometry' and 'geometry' for each row + res_sf$angle = sapply(1:nrow(res_sf), function(i) { + calculate_angle(get_vector(res_sf$corr_line_geometry[i,]), get_vector(res_sf$geometry[i,])) + }) for(i in sum_cols) { # TODO: deduplicate length_y = as.numeric(sf::st_length(rnet_y)) + mask <- (res_sf$angle < 15) | (res_sf$angle > 160) # i = sum_cols[1] - res_sf[[i]] = res_sf[[i]] / res_sf$length_x - over_estimate = sum(res_sf[[i]] * res_sf$length_x, na.rm = TRUE) / - sum(rnet_y[[i]] * length_y, na.rm = TRUE) - res_sf[[i]] = res_sf[[i]] / over_estimate + res_sf[mask, i] <- res_sf[mask, i] / res_sf[mask, "length_x"] + over_estimate <- sum(res_sf[mask, i] * res_sf[mask, "length_x"], na.rm = TRUE) / + sum(rnet_y[mask, i] * length_y, na.rm = TRUE) + + res_sf[mask, i] <- res_sf[mask, i] / over_estimate } } res_sf } + +get_vector <- function(line) { + if (class(line) == "LINESTRING") { + coords <- st_coordinates(line) + start <- coords[1, 1:2] + end <- coords[2, 1:2] + } else { # For MultiLineStrings, just use the first line + first_line <- st_cast(line, "LINESTRING")[[1]] + coords <- st_coordinates(first_line) + 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) +} From 5b2b40d320173a572f10d3f5f18962140b9a53c4 Mon Sep 17 00:00:00 2001 From: wangzhao0217 <74598734+wangzhao0217@users.noreply.github.com> Date: Mon, 4 Sep 2023 22:20:25 +0100 Subject: [PATCH 2/4] corr_line_geometry is not exist --- R/rnet_join.R | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/R/rnet_join.R b/R/rnet_join.R index e7ca0118..9acb5b2e 100644 --- a/R/rnet_join.R +++ b/R/rnet_join.R @@ -224,7 +224,7 @@ rnet_merge <- function(rnet_x, rnet_y, dist = 5, funs = NULL, sum_flows = TRUE, sum_cols = sapply(funs, function(f) identical(f, sum)) sum_cols = names(funs)[which(sum_cols)] rnetj = rnet_join(rnet_x, rnet_y, dist = dist, ...) - names(rnetj) + print(names(rnetj)) rnetj_df = sf::st_drop_geometry(rnetj) # Apply functions to columns with lapply: res_list = lapply(seq(length(funs)), function(i) { @@ -250,11 +250,13 @@ rnet_merge <- function(rnet_x, rnet_y, dist = 5, funs = NULL, sum_flows = TRUE, res_df = dplyr::bind_cols(res_list) res_sf = dplyr::left_join(rnet_x, res_df) if (sum_flows) { + print('-------------start----------') res_sf$length_x = as.numeric(sf::st_length(res_sf)) # Calculate the angle between 'corr_line_geometry' and 'geometry' for each row res_sf$angle = sapply(1:nrow(res_sf), function(i) { calculate_angle(get_vector(res_sf$corr_line_geometry[i,]), get_vector(res_sf$geometry[i,])) }) + print(names(res_sf)) for(i in sum_cols) { # TODO: deduplicate length_y = as.numeric(sf::st_length(rnet_y)) From 0335a94980ca0288f39a6112d48d2ff08b9ae0cb Mon Sep 17 00:00:00 2001 From: wangzhao0217 <74598734+wangzhao0217@users.noreply.github.com> Date: Wed, 6 Sep 2023 23:45:45 +0100 Subject: [PATCH 3/4] the side road not getting value from main road now --- R/rnet_join.R | 91 ++++++++++++++++++---------- vignettes/merging-route-networks.Rmd | 10 +-- 2 files changed, 65 insertions(+), 36 deletions(-) diff --git a/R/rnet_join.R b/R/rnet_join.R index 9acb5b2e..9b7ee2b5 100644 --- a/R/rnet_join.R +++ b/R/rnet_join.R @@ -82,12 +82,13 @@ #' @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 = "FLAT", contains = FALSE, ...) { if (is.null(dist_subset)) { dist_subset = dist + 1 } if (subset_x) { rnet_x = rnet_subset(rnet_x, rnet_y, dist = dist_subset, ...) + print('rnet_subset done') } rnet_x_buffer = geo_buffer(rnet_x, dist = dist, nQuadSegs = 2, endCapStyle = endCapStyle) if (segment_length > 0) { @@ -109,7 +110,6 @@ rnet_join = function(rnet_x, rnet_y, dist = 5, length_y = TRUE, key_column = 1, rnet_y_centroids$corr_line_geometry = rnet_y$geometry rnetj = sf::st_join(rnet_x_buffer[key_column], rnet_y_centroids) } - rnetj } @@ -223,69 +223,96 @@ rnet_merge <- function(rnet_x, rnet_y, dist = 5, funs = NULL, sum_flows = TRUE, } sum_cols = sapply(funs, function(f) identical(f, sum)) sum_cols = names(funs)[which(sum_cols)] - rnetj = rnet_join(rnet_x, rnet_y, dist = dist, ...) - print(names(rnetj)) + # First, calculate the angle and store it in the `rnetj` + rnetj = rnet_join(rnet_x, rnet_y, dist = dist) + rnetj$angle = sapply(1:nrow(rnetj), function(i) { + calculate_angle(get_vector(rnetj$corr_line_geometry[[i]]), get_vector(rnetj$geometry[[i]])) + }) rnetj_df = sf::st_drop_geometry(rnetj) + # Apply functions to columns with lapply: res_list = lapply(seq(length(funs)), function(i) { # i = 1 nm = names(funs[i]) fn = funs[[i]] - + + # Keep the first non-NA value of 'corr_line_geometry' and 'angle' within each group + intermediate_df = rnetj_df %>% + dplyr::group_by_at(1) %>% + dplyr::mutate(corr_line_geometry = first(corr_line_geometry), + angle = first(angle)) %>% + dplyr::ungroup() + if (identical(fn, sum) && sum_flows) { - res = rnetj_df %>% + res = intermediate_df %>% dplyr::group_by_at(1) %>% - dplyr::summarise(dplyr::across(dplyr::matches(nm), function(x) sum(x * length_y))) + dplyr::summarise(dplyr::across(dplyr::matches(nm), function(x) sum(x * length_y)), .groups = 'drop') } else { - res = rnetj_df %>% + res = intermediate_df %>% dplyr::group_by_at(1) %>% - dplyr::summarise(dplyr::across(dplyr::matches(nm), fn)) + dplyr::summarise(dplyr::across(dplyr::matches(nm), fn), .groups = 'drop') } - names(res)[2] = nm + + # Add back the 'corr_line_geometry' and 'angle' columns + res = dplyr::left_join(res, unique(intermediate_df %>% dplyr::select(identifier, corr_line_geometry, angle)), by = "identifier") + + # Rename the summarised column + names(res)[which(names(res) == nm)[1]] = nm + if(i > 1) { - res = res[-1] + # Drop the 'identifier' column to avoid duplication + res = res %>% dplyr::select(-matches("^identifier")) } + res }) + res_df = dplyr::bind_cols(res_list) + res_df <- res_df %>% + dplyr::select(identifier, value, Quietness, corr_line_geometry = corr_line_geometry...3, angle = angle...4) + names(res_df) + mask <- (res_df$angle < 60) | (res_df$angle > 110) + # mask[is.na(mask)] <- FALSE + filtered_res_df <- res_df[mask, ] + res_sf = dplyr::left_join(rnet_x, res_df) + if (sum_flows) { - print('-------------start----------') res_sf$length_x = as.numeric(sf::st_length(res_sf)) - # Calculate the angle between 'corr_line_geometry' and 'geometry' for each row - res_sf$angle = sapply(1:nrow(res_sf), function(i) { - calculate_angle(get_vector(res_sf$corr_line_geometry[i,]), get_vector(res_sf$geometry[i,])) - }) - print(names(res_sf)) for(i in sum_cols) { # TODO: deduplicate length_y = as.numeric(sf::st_length(rnet_y)) - mask <- (res_sf$angle < 15) | (res_sf$angle > 160) # i = sum_cols[1] - res_sf[mask, i] <- res_sf[mask, i] / res_sf[mask, "length_x"] - over_estimate <- sum(res_sf[mask, i] * res_sf[mask, "length_x"], na.rm = TRUE) / - sum(rnet_y[mask, i] * length_y, na.rm = TRUE) - - res_sf[mask, i] <- res_sf[mask, i] / over_estimate + res_sf[[i]] = res_sf[[i]] / res_sf$length_x + over_estimate = sum(res_sf[[i]] * res_sf$length_x, na.rm = TRUE) / + sum(rnet_y[[i]] * length_y, na.rm = TRUE) + res_sf[[i]] = res_sf[[i]] / over_estimate } } res_sf } get_vector <- function(line) { - if (class(line) == "LINESTRING") { - coords <- st_coordinates(line) - start <- coords[1, 1:2] - end <- coords[2, 1:2] - } else { # For MultiLineStrings, just use the first line - first_line <- st_cast(line, "LINESTRING")[[1]] - coords <- st_coordinates(first_line) - start <- coords[1, 1:2] - end <- coords[2, 1:2] + 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)) diff --git a/vignettes/merging-route-networks.Rmd b/vignettes/merging-route-networks.Rmd index acfedbff..cdfb13bf 100644 --- a/vignettes/merging-route-networks.Rmd +++ b/vignettes/merging-route-networks.Rmd @@ -17,7 +17,7 @@ knitr::opts_chunk$set( message = FALSE, warning = FALSE ) -# devtools::load_all() +devtools::load_all() sf::sf_use_s2(FALSE) ``` @@ -43,15 +43,17 @@ tmap_mode("view") # nrow(rnet_x) # summary(sf::st_length(rnet_x)) plot(sf::st_geometry(rnet_x)) +dim(rnet_x) rnet_x = rnet_subset(rnet_x, rnet_y, dist = 20) -# nrow(rnet_x) +dim(rnet_x)# nrow(rnet_x) # plot(sf::st_geometry(rnet_x)) rnet_x = rnet_subset(rnet_x, rnet_y, dist = 20, min_length = 5) +dim(rnet_x) # summary(sf::st_length(rnet_x)) # nrow(rnet_x) # plot(sf::st_geometry(rnet_x)) rnet_x = rnet_subset(rnet_x, rnet_y, dist = 20, rm_disconnected = TRUE) -# nrow(rnet_x) +dim(rnet_x)# nrow(rnet_x) plot(sf::st_geometry(rnet_x)) ``` @@ -71,7 +73,7 @@ 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) +rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 10, funs = funs) 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) From be57336d7723328a926fa7d4cc9cd6c6d85c3a1a Mon Sep 17 00:00:00 2001 From: wangzhao0217 <74598734+wangzhao0217@users.noreply.github.com> Date: Wed, 6 Sep 2023 23:52:08 +0100 Subject: [PATCH 4/4] update merging-route-networks.Rmd --- vignettes/merging-route-networks.Rmd | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/vignettes/merging-route-networks.Rmd b/vignettes/merging-route-networks.Rmd index cdfb13bf..e25443b0 100644 --- a/vignettes/merging-route-networks.Rmd +++ b/vignettes/merging-route-networks.Rmd @@ -73,7 +73,11 @@ tmap_arrange(m1, m2, sync = TRUE) We can more reduce the minimum segment length to ensure fewer NA values in the outputs: ```{r} +tmap_mode("view") +funs = list(value = sum, Quietness = mean) +brks = c(0, 100, 500, 1000, 5000) rnet_merged = rnet_merge(rnet_x, rnet_y, dist = 10, segment_length = 10, funs = funs) +rnet_merged$identifier <- NULL 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)