66# ' Polygon defining land boundaries used to compute fetch distances.
77# ' @param max_dist {`numeric`}\cr{}
88# ' Maximum fetch distance in kilometers. Fetch beyond this distance is capped.
9- # ' @param n_quad_seg {`integer`}\cr{}
10- # ' Number of segments per quadrant for fetch calculation.
11- # ' Ignored if `wind_weights` is provided.
9+ # ' @param n_bearings {`integer`}\cr{}
10+ # ' Total number of bearings for fetch calculation (minimal number required is
11+ # ' 4, default is 167). Ignored if `wind_weights` is provided.
1212# ' @param wind_weights {`data.frame`}\cr{}
1313# ' A data frame specifying directional weights for wind exposure.
1414# ' Must contain two columns: `direction` (numeric, in degrees) and `weight`
1717# ' Coordinate reference system (CRS) passed to [sf::st_crs()], used to
1818# ' transform `points` and `polygon`.
1919# '
20- # ' @details Wind fetch is the unobstructed distance over which wind travels
20+ # ' @details
21+ # ' Wind fetch is the unobstructed distance over which wind travels
2122# ' across a body of water before reaching a specific point. It plays a crucial
2223# ' role in wave generation, as longer fetch distances allow wind to transfer
2324# ' more energy to the water surface, leading to larger waves.
2425# '
25- # ' For all points in `points`, 4 × `n_quad_seg ` radial transects are generated
26+ # ' For all points in `points`, `n_bearings ` radial transects are generated
2627# ' by default. If `wind_weights` is specified, the column `direction`, which
2728# ' contains angles in degrees, is used instead to generate the transects. The
2829# ' transects are then clipped with the polygon using [`sf::st_intersection()`],
3233# ' element in the returned list and it used to generate the second element:
3334# ' `mean_fetch` that included wind fetch averages.
3435# '
35- # ' Ensure that max_dist is specified in meters. An error will be thrown if the
36+ # ' Ensure that ` max_dist` is specified in meters. An error will be thrown if the
3637# ' spatial projection of points and polygon is not in a meter-based coordinate
3738# ' system.
3839# '
39- # ' @return A list of two elements:
40- # ' * `mean_fetch`: data frame with 5 columns:
40+ # ' @return
41+ # ' A list of two elements:
42+ # ' * `mean_fetch`: a `sf` object with 3 features:
4143# ' * `id_point`: point identifier
42- # ' * `fetch`: mean wind fetch based on the four highest values
43- # ' * `weighted_fetch`: mean weighted wind fetch based on the four highest values
44- # ' * `fetch_all`: mean wind fetch based on all values
45- # ' * `weighted_fetch_all`: mean wind weighted fetch based on all values
44+ # ' * `fetch_km`: mean wind fetch based on all bearings.
45+ # ' * `weighted_fetch_km`: mean weighted wind fetch based on all bearings.
4646# ' * `transect_lines`: a `sf` object containing all radial transect with the
4747# ' same columns as `points` and the following additional columns:
4848# ' * `id_point`: point identifier
6363# '
6464# ' @examples
6565# ' \donttest{
66- # '
6766# ' le_bound <- system.file("example", "lake_erie.gpkg", package = "SAVM") |>
6867# ' sf::st_read()
6968# ' le_pt <- system.file("example", "le_points.geojson", package = "SAVM") |>
7069# ' sf::st_read(quiet = TRUE)
7170# ' res <- compute_fetch(le_pt, le_bound, crs = 32617)
72- # ' # use wind-weight
71+ # ' # use wind-weight
7372# ' res2 <- compute_fetch(
74- # ' le_pt, le_bound, max_dist = 20,
75- # ' wind_weights = data.frame(
76- # ' direction = seq(0, 360, by = 360 / 16)[-1],
77- # ' weight = rep(c(0, 1), each = 8)
78- # ' ),
79- # ' crs = 32617)
73+ # ' le_pt, le_bound,
74+ # ' max_dist = 20,
75+ # ' wind_weights = data.frame(
76+ # ' direction = seq(0, 360, by = 360 / 16)[-1],
77+ # ' weight = rep(c(0, 1), each = 8)
78+ # ' ),
79+ # ' crs = 32617
80+ # ' )
8081# '
81- # ' # resultat
82+ # ' # results
8283# ' res$mean_fetch
8384# ' res2$mean_fetch
8485# '
8586# ' # visualizing fetch lines
8687# ' plot(le_bound |> sf::st_transform(crs = 32617) |> sf::st_geometry())
8788# ' plot(res$transect_lines |> sf::st_geometry(), add = TRUE, col = 2, lwd = 0.5)
8889# ' }
89- compute_fetch <- function (points , polygon , max_dist = 15 , n_quad_seg = 9 , wind_weights = NULL , crs = NULL ) {
90+ compute_fetch <- function (
91+ points , polygon , max_dist = 15 , n_bearings = 16 , wind_weights = NULL , crs = NULL ) {
9092 valid_points(points )
9193 points $ id_point <- seq_len(nrow(points ))
9294 valid_polygon(polygon )
93- sav_stop_if_not(max_dist > 0 )
95+ sav_stop_if_not(max_dist > 0 , " `max_dist` must be strictly positive. "
9496 max_dist <- 1e3 * max_dist
95- sav_stop_if_not(n_quad_seg > 0 )
97+ sav_stop_if_not(n_bearings > = 4 , " `n_bearings` should be equal or greater than 4."
98+ if (n_bearings > 64 ) {
99+ sav_msg_warning(
100+ " Large number of bearings detected, computation may take a long time."
101+ )
102+ }
96103
97104 if (! is.null(crs )) {
98105 if (! is_proj_unit_meter(crs )) {
@@ -127,11 +134,11 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w
127134
128135 if (is.null(wind_weights )) {
129136 d_direction <- data.frame (
130- direction = utils :: head(seq(0 , 360 , by = 360 / ( n_quad_seg * 4 ) ), - 1 ),
137+ direction = utils :: head(seq(0 , 360 , by = 360 / n_bearings ), - 1 ),
131138 weight = 1
132139 )
133140 } else {
134- sav_msg_info(" Using `wind_weights`, ignoring `n_quad_seg `"
141+ sav_msg_info(" Using `wind_weights`, ignoring `n_bearings `"
135142 if (all(c(" direction" " weight" %in% names(wind_weights ))) {
136143 d_direction <- wind_weights [c(" direction" " weight"
137144 valid_direction(d_direction $ direction )
@@ -146,7 +153,7 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w
146153 sav_msg_info(" Cropping fetch lines"
147154 fetch_crop <- suppressWarnings(fetch_lines | > sf :: st_intersection(polygon ))
148155 geom_type <- sf :: st_geometry_type(fetch_crop )
149- # sf::st_intersection() generates multilinestring with extra lines if there
156+ # sf::st_intersection() generates MULTILINESTRING with extra lines if there
150157 # are intersections within the fetch lines
151158 transect_lines <- rbind(
152159 fetch_crop | >
@@ -156,28 +163,33 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w
156163 remove_detached_ends(points )
157164 ) | >
158165 dplyr :: arrange(id_point , direction )
159-
160166 transect_lines <- transect_lines | >
161167 dplyr :: mutate(transect_length = sf :: st_length(transect_lines )) | >
162168 dplyr :: group_by(id_point ) | >
163- dplyr :: mutate(rank = rank(transect_length , ties.method = " min"
169+ # using -transect so that the longest are ranked 1
170+ dplyr :: mutate(rank = rank(- transect_length , ties.method = " min"
164171
165172 list (
166- mean_fetch = transect_lines | >
167- sf :: st_drop_geometry() | >
168- dplyr :: group_by(id_point ) | >
169- # dplyr::mutate(rank = rank(transect_length)) |>
170- dplyr :: summarise(
171- fetch_km = mean(transect_length [rank < 5 ]),
172- weighted_fetch_km = mean(transect_length [rank < 5 ] * weight [rank < 5 ]),
173- fetch_km_all = mean(transect_length ),
174- weighted_fetch_km_all = mean(transect_length * weight )
173+ mean_fetch = points | >
174+ dplyr :: left_join(
175+ transect_lines | >
176+ sf :: st_drop_geometry() | >
177+ dplyr :: group_by(id_point ) | >
178+ # dplyr::mutate(rank = rank(transect_length)) |>
179+ dplyr :: summarise(
180+ fetch_km = mean(transect_length ),
181+ weighted_fetch_km = mean(transect_length * weight )
182+ ) | >
183+ dplyr :: mutate(
184+ dplyr :: across(
185+ ! c(id_point ),
186+ ~ as.numeric(units :: set_units(.x , " km"
187+ )
188+ ),
189+ by = " id_point"
175190 ) | >
176- dplyr :: mutate(
177- dplyr :: across(
178- ! id_point ,
179- ~ as.numeric(units :: set_units(.x , " km"
180- )
191+ dplyr :: select(
192+ c(" id_point" " fetch_km" " weighted_fetch_km"
181193 ),
182194 transect_lines = transect_lines
183195 )
0 commit comments