diff --git a/.Rbuildignore b/.Rbuildignore index 304dbdb..06615be 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -10,3 +10,4 @@ ^pkgdown$ ^tmp$ ^notes\.md$ +^\.vscode$ diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index fdf6a55..ed74714 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -16,14 +16,11 @@ jobs: fail-fast: false matrix: config: - # - {os: macOS-latest, r: 'devel'} - # - {os: macOS-latest, r: 'release'} - # - {os: windows-latest, r: 'devel'} - # - {os: windows-latest, r: 'release'} - # - {os: windows-latest, r: 'oldrel'} - # - {os: ubuntu-22.04, r: 'devel'} - - {os: ubuntu-22.04, r: 'release'} - # - {os: ubuntu-22.04, r: 'oldrel'} + #- {os: macOS-latest, r: 'release'} + - {os: windows-latest, r: 'release'} + - {os: ubuntu-latest, r: 'devel'} + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'oldrel'} env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true @@ -41,13 +38,8 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: - extra-packages: any::rcmdcheck, any::covr + extra-packages: any::rcmdcheck needs: check - uses: r-lib/actions/check-r-package@v2 - # - name: Test coverage - # if: matrix.config.os == 'ubuntu-22.04' && matrix.config.r == 'release' - # run: | - # covr::codecov(token = "${{secrets.CODECOV_TOKEN}}") - # shell: Rscript {0} diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..2bf2d23 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,49 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + release: + types: [published] + workflow_dispatch: + +name: pkgdown.yaml + +permissions: read-all + +jobs: + pkgdown: + runs-on: ubuntu-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::pkgdown, local::. + needs: website + + - name: Build site + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + + - name: Deploy to GitHub pages πŸš€ + if: github.event_name != 'pull_request' + uses: JamesIves/github-pages-deploy-action@v4.5.0 + with: + clean: false + branch: gh-pages + folder: docs \ No newline at end of file diff --git a/DESCRIPTION b/DESCRIPTION index 13e608d..8cabed1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: SAVM Title: Submerged aquatic vegetation model -Version: 0.0.1 -Date: 2025-03-21 +Version: 0.0.1.9004 +Date: 2025-07-15 Authors@R: c( person(given = "Kevin", family = "Cazelles", @@ -12,11 +12,17 @@ Authors@R: c( family = "Beauchesne", role = c("aut"), email = "david.beauchesne@insileco.io", - comment = c(ORCID = "0000-0002-3590-8161")) + comment = c(ORCID = "0000-0002-3590-8161")), + person( + given = "Paul", + family = "Bozneck", + role = c("aut"), + email = "Paul.Bzonek@dfo-mpo.gc.ca" ) -Description: Here is the description made of sentences. -URL: https://github.com/inSileco/SAVM -BugReports: https://github.com/inSileco/SAVM/issues + ) +Description: Submerged aquatic vegetation model developed by the Fish Ecology Science Lab at DFO. +URL: https://github.com/FishEcologyScience/SAVM +BugReports: https://github.com/FishEcologyScience/SAVM/issues License: GPL-3 Encoding: UTF-8 LazyData: true diff --git a/NAMESPACE b/NAMESPACE index bc44942..a4a2709 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,11 +1,15 @@ # Generated by roxygen2: do not edit by hand +S3method(preview_grid,sav_data) export(compute_fetch) +export(invert_polygon) export(plot_sav_density) export(plot_sav_distribution) export(plot_sav_tmap) +export(preview_grid) export(read_sav) export(read_sav_aoi) export(read_sav_csv) export(read_sav_pts) export(sav_model) +import(randomForest) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..14e1e05 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,10 @@ +# SAVM (devel) + +* Model and plot functions have been adjusted to handle `sf` objects (see #11). +* The element `mean_fetch` returned by `compute_fetch()` is now a `sf` object (see #9). +* `preview_grid()` allows to preview grid (see #8). +* Add more guidance on reading shapefiles in the vignette (see #6). +* `invert_polygon()` has been added to invert polygon (see #6). +* `compute_fetch()` has a new argument `n_bearings` that provides the number of bearings, it replaces `n_quad_seg` (see #4 and #5). +* `compute_fetch()` only compute the mean fetch for all bearings, columns with +suffix `_all` where therefore removed (see #4). \ No newline at end of file diff --git a/R/compute_fetch.R b/R/compute_fetch.R index b1bd671..81d3f63 100644 --- a/R/compute_fetch.R +++ b/R/compute_fetch.R @@ -6,9 +6,9 @@ #' Polygon defining land boundaries used to compute fetch distances. #' @param max_dist {`numeric`}\cr{} #' Maximum fetch distance in kilometers. Fetch beyond this distance is capped. -#' @param n_quad_seg {`integer`}\cr{} -#' Number of segments per quadrant for fetch calculation. -#' Ignored if `wind_weights` is provided. +#' @param n_bearings {`integer`}\cr{} +#' Total number of bearings for fetch calculation (minimal number required is +#' 4, default is 167). Ignored if `wind_weights` is provided. #' @param wind_weights {`data.frame`}\cr{} #' A data frame specifying directional weights for wind exposure. #' Must contain two columns: `direction` (numeric, in degrees) and `weight` @@ -17,12 +17,13 @@ #' Coordinate reference system (CRS) passed to [sf::st_crs()], used to #' transform `points` and `polygon`. #' -#' @details Wind fetch is the unobstructed distance over which wind travels +#' @details +#' Wind fetch is the unobstructed distance over which wind travels #' across a body of water before reaching a specific point. It plays a crucial #' role in wave generation, as longer fetch distances allow wind to transfer #' more energy to the water surface, leading to larger waves. #' -#' For all points in `points`, 4 Γ— `n_quad_seg` radial transects are generated +#' For all points in `points`, `n_bearings` radial transects are generated #' by default. If `wind_weights` is specified, the column `direction`, which #' contains angles in degrees, is used instead to generate the transects. The #' transects are then clipped with the polygon using [`sf::st_intersection()`], @@ -32,17 +33,16 @@ #' element in the returned list and it used to generate the second element: #' `mean_fetch` that included wind fetch averages. #' -#' Ensure that max_dist is specified in meters. An error will be thrown if the +#' Ensure that `max_dist` is specified in meters. An error will be thrown if the #' spatial projection of points and polygon is not in a meter-based coordinate #' system. #' -#' @return A list of two elements: -#' * `mean_fetch`: data frame with 5 columns: +#' @return +#' A list of two elements: +#' * `mean_fetch`: a `sf` object with 3 features: #' * `id_point`: point identifier -#' * `fetch`: mean wind fetch based on the four highest values -#' * `weighted_fetch`: mean weighted wind fetch based on the four highest values -#' * `fetch_all`: mean wind fetch based on all values -#' * `weighted_fetch_all`: mean wind weighted fetch based on all values +#' * `fetch_km`: mean wind fetch based on all bearings. +#' * `weighted_fetch_km`: mean weighted wind fetch based on all bearings. #' * `transect_lines`: a `sf` object containing all radial transect with the #' same columns as `points` and the following additional columns: #' * `id_point`: point identifier @@ -63,22 +63,23 @@ #' #' @examples #' \donttest{ -#' #' le_bound <- system.file("example", "lake_erie.gpkg", package = "SAVM") |> #' sf::st_read() #' le_pt <- system.file("example", "le_points.geojson", package = "SAVM") |> #' sf::st_read(quiet = TRUE) #' res <- compute_fetch(le_pt, le_bound, crs = 32617) -#' # use wind-weight +#' # use wind-weight #' res2 <- compute_fetch( -#' le_pt, le_bound, max_dist = 20, -#' wind_weights = data.frame( -#' direction = seq(0, 360, by = 360 / 16)[-1], -#' weight = rep(c(0, 1), each = 8) -#' ), -#' crs = 32617) +#' le_pt, le_bound, +#' max_dist = 20, +#' wind_weights = data.frame( +#' direction = seq(0, 360, by = 360 / 16)[-1], +#' weight = rep(c(0, 1), each = 8) +#' ), +#' crs = 32617 +#' ) #' -#' # resultat +#' # results #' res$mean_fetch #' res2$mean_fetch #' @@ -86,13 +87,19 @@ #' plot(le_bound |> sf::st_transform(crs = 32617) |> sf::st_geometry()) #' plot(res$transect_lines |> sf::st_geometry(), add = TRUE, col = 2, lwd = 0.5) #' } -compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_weights = NULL, crs = NULL) { +compute_fetch <- function( + points, polygon, max_dist = 15, n_bearings = 16, wind_weights = NULL, crs = NULL) { valid_points(points) points$id_point <- seq_len(nrow(points)) valid_polygon(polygon) - sav_stop_if_not(max_dist > 0) + sav_stop_if_not(max_dist > 0, "`max_dist` must be strictly positive.") max_dist <- 1e3 * max_dist - sav_stop_if_not(n_quad_seg > 0) + sav_stop_if_not(n_bearings >= 4, "`n_bearings` should be equal or greater than 4.") + if (n_bearings > 64) { + sav_msg_warning( + "Large number of bearings detected, computation may take a long time." + ) + } if (!is.null(crs)) { if (!is_proj_unit_meter(crs)) { @@ -127,11 +134,11 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w if (is.null(wind_weights)) { d_direction <- data.frame( - direction = utils::head(seq(0, 360, by = 360 / (n_quad_seg * 4)), -1), + direction = utils::head(seq(0, 360, by = 360 / n_bearings), -1), weight = 1 ) } else { - sav_msg_info("Using `wind_weights`, ignoring `n_quad_seg`") + sav_msg_info("Using `wind_weights`, ignoring `n_bearings`") if (all(c("direction", "weight") %in% names(wind_weights))) { d_direction <- wind_weights[c("direction", "weight")] valid_direction(d_direction$direction) @@ -146,7 +153,7 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w sav_msg_info("Cropping fetch lines") fetch_crop <- suppressWarnings(fetch_lines |> sf::st_intersection(polygon)) geom_type <- sf::st_geometry_type(fetch_crop) - # sf::st_intersection() generates multilinestring with extra lines if there + # sf::st_intersection() generates MULTILINESTRING with extra lines if there # are intersections within the fetch lines transect_lines <- rbind( fetch_crop |> @@ -156,28 +163,33 @@ compute_fetch <- function(points, polygon, max_dist = 15, n_quad_seg = 9, wind_w remove_detached_ends(points) ) |> dplyr::arrange(id_point, direction) - transect_lines <- transect_lines |> dplyr::mutate(transect_length = sf::st_length(transect_lines)) |> dplyr::group_by(id_point) |> - dplyr::mutate(rank = rank(transect_length, ties.method = "min")) + # using -transect so that the longest are ranked 1 + dplyr::mutate(rank = rank(-transect_length, ties.method = "min")) list( - mean_fetch = transect_lines |> - sf::st_drop_geometry() |> - dplyr::group_by(id_point) |> - # dplyr::mutate(rank = rank(transect_length)) |> - dplyr::summarise( - fetch_km = mean(transect_length[rank < 5]), - weighted_fetch_km = mean(transect_length[rank < 5] * weight[rank < 5]), - fetch_km_all = mean(transect_length), - weighted_fetch_km_all = mean(transect_length * weight) + mean_fetch = points |> + dplyr::left_join( + transect_lines |> + sf::st_drop_geometry() |> + dplyr::group_by(id_point) |> + # dplyr::mutate(rank = rank(transect_length)) |> + dplyr::summarise( + fetch_km = mean(transect_length), + weighted_fetch_km = mean(transect_length * weight) + ) |> + dplyr::mutate( + dplyr::across( + !c(id_point), + ~ as.numeric(units::set_units(.x, "km")) + ) + ), + by = "id_point" ) |> - dplyr::mutate( - dplyr::across( - !id_point, - ~ as.numeric(units::set_units(.x, "km")) - ) + dplyr::select( + c("id_point", "fetch_km", "weighted_fetch_km") ), transect_lines = transect_lines ) diff --git a/R/invert_polygon.R b/R/invert_polygon.R new file mode 100644 index 0000000..e5e16bf --- /dev/null +++ b/R/invert_polygon.R @@ -0,0 +1,38 @@ +#' Invert a Polygon +#' +#' @param polygon {`sf`}\cr{} +#' Polygon defining land boundaries that needs to be inverted. +#' @param ratio {`numeric [0-1]`}\cr{} +#' Fraction convex, see [sf::st_concave_hull()]. This may need to be tweaked +#' depending on the form of the polygon to be inverted. +#' +#' @details +#' Utility function that inverts a polygon by drawing a concave hull around +#' `polygon` (see [sf::st_concave_hull()]) and then computes the differences +#' between the polygon and the concave hull (see [sf::st_concave_hull()]). +#' +#' @export +#' +#' @examples +#' \donttest{ +#' erie_land <- system.file( +#' "example", "lake_erie_land", "LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx", +#' package = "SAVM", mustWork = TRUE +#' ) |> sf::st_read() +#' +#' erie_land |> +#' sf::st_geometry() |> +#' plot(col = 1) +#' +#' erie_land |> +#' sf::st_geometry() |> +#' invert_polygon() |> +#' plot(col = 1) +#' } + +invert_polygon <- function(polygon, ratio = 0.5) { + os2 <- sf::sf_use_s2() + on.exit(suppressMessages(sf::sf_use_s2(os2))) + suppressMessages(os2 <- sf::sf_use_s2(FALSE)) + sf::st_difference(polygon |> sf::st_concave_hull(ratio = ratio), polygon) +} \ No newline at end of file diff --git a/R/model.R b/R/model.R index fcc6379..024390d 100644 --- a/R/model.R +++ b/R/model.R @@ -3,7 +3,7 @@ #' Apply Random Forest models to predict SAV cover and presence/absence, with #' optional post-hoc processing. #' -#' @param dat {`data.frame`}\cr{} A `data.frame` containing some or all of the +#' @param dat {`data.frame`|`sf`}\cr{} A `data.frame` or a `sf` object containing some or all of the #' following columns: #' - `depth_m`: Numeric, depth in meters. #' - `fetch_km`: Numeric, fetch in kilometers. @@ -12,28 +12,30 @@ #' limitations. (post_hoc) #' - `limitation`: Binary (0 = absent, 1 = present), indicating user-supplied #' limitations. -#' Additionnal columns will be ignored. +#' Additional columns will be ignored. #' @param type {`character vector`, either `"cover"` or `"pa"`}\cr{} #' Model type(s). -#' @param depth,fetch {`character`} Column specification for the predictors, +#' @param depth,fetch {`character`}\cr{} Column specification for the predictors, #' see *Details*. -#' @param substrate,secchi,limitation Column specification for post_hoc +#' @param substrate,secchi,limitation {`character`}\cr{}Column specification for post_hoc #' variables, see *Details*. #' @param vmax_par {`named list`}\cr{} intercept and slope of the equation from #' Chambers and Kalff (1985) to compute the maximum depth of plant colonization #' (Vmax), see *Details* below. #' #' @return -#' A data frame containing the input columns along with model predictions. +#' A data frame (or a sf object) containing the input columns along with model +#' predictions. #' -#' The prediction column names match the values specified in `type`, and contain -#' the raw model outputs (i.e., without post-hoc adjustment). +#' The prediction column names match the values specified in `type` followed +#' by the suffix `_pred`, and contain the raw model outputs (i.e., without +#' post-hoc adjustment). #' #' Post-hoc adjusted predictions (see *Details*) are included in additional #' columns with the same names as the `type` values, but with the suffix #' `_post_hoc`. #' -#' If a column `secchi` is present, then two additionnal columns are +#' If a column `secchi` is present, then two additional columns are #' returned: `vmax` and `limitation_secchi`, see details for further #' explanation. #' @@ -92,7 +94,7 @@ #' sav_model(data.frame(depth = c(5, 10))) #' sav_model(data.frame(depth = c(5, 10), fetch = c(1, 2)), type = "pa") #' -#' # using post-hoc tretment +#' # using post-hoc treatment #' sav_model( #' data.frame( #' depth = c(5, 10, 5), @@ -106,7 +108,19 @@ sav_model <- function( dat, type = c("cover", "pa"), depth = NULL, fetch = NULL, substrate = NULL, secchi = NULL, limitation = NULL, vmax_par = list(intercept = 1.40, slope = 1.33)) { - sav_stop_if_not(inherits(dat, "data.frame")) + + + # this should rather be handle via S3 + geom <- NULL + if (inherits(dat, "sf")) { + + geom <- dat |> + dplyr::select(geometry) + dat <- dat |> sf::st_drop_geometry() + } else { + sav_stop_if_not(inherits(dat, "data.frame")) + } + type <- unique(type) if (!all(type %in% c("cover", "pa"))) { @@ -202,7 +216,15 @@ sav_model <- function( scrub_if_present("limitation_secchi", "cover_post_hoc") } - out + out <- out |> + rename_if_present("pa", "pa_pred") |> + rename_if_present("cover", "cover_pred") + + if (is.null(geom)) { + out + } else { + cbind(geom, out) + } } diff --git a/R/plot_sav.R b/R/plot_sav.R index 5bda434..10716fe 100644 --- a/R/plot_sav.R +++ b/R/plot_sav.R @@ -1,23 +1,32 @@ #' Plot SAV Data Distribution #' -#' This function generates up to four plots representing the distribution of submerged aquatic vegetation (SAV) -#' data by depth (m) and fetch (km). It visualizes SAV presence/absence (PA) and percent cover (Cover) when -#' the corresponding columns are available in the input data. +#' This function generates up to four plots representing the distribution of +#' submerged aquatic vegetation (SAV) data by depth (m) and fetch (km). It +#' visualizes SAV presence/absence (PA) and percent cover (Cover) when the +#' corresponding columns are available in the input data. #' -#' @param dat A `data.frame` containing some or all of the following columns: +#' @param dat {`data.frame`}\cr{} +#' A data frame containing some or all of the following columns: #' - `depth_m`: Numeric, depth in meters. #' - `fetch_km`: Numeric, fetch in kilometers. #' - `pa`: Binary (0 = absent, 1 = present), indicating SAV presence/absence. #' - `cover`: Numeric, percent cover of SAV. -#' @param type Character vector specifying the type of plots to generate. Options: +#' @param type {`character vector`}\cr{} +#' Character vector specifying the type of plots to generate. Options: #' - `"pa"` (default) for presence/absence plots #' - `"cover"` (default) for cover percentage plots -#' @param predictors Character vector specifying which predictors to use. Options: +#' @param predictors {`character vector`}\cr{} +#' Character vector specifying which predictors to use. Options: #' - `"depth"` (default) for depth-based plots #' - `"fetch"` (default) for fetch-based plots -#' @param post_hoc Logical value indicating whether to use post-hoc analyzed columns (`pa_post_hoc`, `cover_post_hoc`) instead of raw columns (`pa`, `cover`). Default is `TRUE`. -#' @param max_depth Numeric value specifying the maximum depth bin (default: 30 meters). -#' @param max_fetch Numeric value specifying the maximum fetch bin (default: 15 km). +#' @param post_hoc {`logical`}\cr{} +#' Logical value indicating whether to use post-hoc analyzed columns +#' (`pa_post_hoc`, `cover_post_hoc`) instead of raw columns (`pa`, `cover`). +#' Default is `TRUE`. +#' @param max_depth {`numeric`}\cr{} +#' Numeric value specifying the maximum depth bin (default: 30 meters). +#' @param max_fetch {`numeric`}\cr{} +#' Numeric value specifying the maximum fetch bin (default: 15 km). #' #' @return A set of ggplot2 plots displayed in a grid layout. #' @@ -28,8 +37,8 @@ #' dat <- data.frame( #' depth_m = runif(100, 0, 15), #' fetch_km = runif(100, 0, 15), -#' pa = sample(0:1, 100, replace = TRUE), -#' cover = runif(100, 0, 100), +#' pa_pred = sample(0:1, 100, replace = TRUE), +#' cover_pred = runif(100, 0, 100), #' pa_post_hoc = sample(0:1, 100, replace = TRUE), #' cover_post_hoc = runif(100, 0, 100) #' ) @@ -50,6 +59,8 @@ plot_sav_distribution <- function(dat, type = c("pa", "cover"), predictors = c("depth", "fetch"), post_hoc = TRUE, max_depth = 30, max_fetch = 15) { plots <- list() + dat <- as.data.frame(dat) + # Check post-hoc if ("pa" %in% type && post_hoc && !"pa_post_hoc" %in% names(dat)) { rlang::abort("Requested post-hoc presence/absence predictions, but they are missing from the provided data.") @@ -59,45 +70,51 @@ plot_sav_distribution <- function(dat, type = c("pa", "cover"), predictors = c(" } # Determine which columns to use based on post_hoc flag - pa_col <- if (post_hoc) "pa_post_hoc" else "pa" - cover_col <- if (post_hoc) "cover_post_hoc" else "cover" + pa_col <- if (post_hoc) "pa_post_hoc" else "pa_pred" + cover_col <- if (post_hoc) "cover_post_hoc" else "cover_pred" # Check for requested columns, abort if unavailable has_depth <- "depth_m" %in% names(dat) - if (!has_depth && "depth" %in% predictors) rlang::abort("Requested layer `depth` is unavailable in provided data)") + if (!has_depth && "depth" %in% predictors) + rlang::abort("Requested layer `depth` is unavailable in provided data)") has_fetch <- "fetch_km" %in% names(dat) - if (!has_fetch && "fetch" %in% predictors) rlang::abort("Requested layer `fetch` is unavailable in provided data)") + if (!has_fetch && "fetch" %in% predictors) + rlang::abort("Requested layer `fetch` is unavailable in provided data)") has_pa <- pa_col %in% names(dat) - if (!has_pa && "pa" %in% type) rlang::abort("Requested layer `pa` is unavailable in provided data)") + if (!has_pa && "pa" %in% type) + rlang::abort("Requested layer `pa` is unavailable in provided data)") has_cover <- cover_col %in% names(dat) - if (!has_cover && "cover" %in% type) rlang::abort("Requested layer `cover` is unavailable in provided data)") + if (!has_cover && "cover" %in% type) + rlang::abort("Requested layer `cover` is unavailable in provided data)") # Process data - dat <- dplyr::mutate( - dat, - Depth_Bin = if (has_depth) { - cut( - dat$depth_m, - breaks = c(seq(0, max_depth, by = 1), Inf), - include.lowest = TRUE, right = FALSE, - labels = c(paste0(seq(0, max_depth - 1, by = 1), "-", seq(1, max_depth, by = 1)), paste0(max_depth, "+")) - ) - } else { - NULL - }, - Fetch_Bin = if (has_fetch) { - cut( - dat$fetch_km, - breaks = c(seq(0, max_fetch, by = 1), Inf), - include.lowest = TRUE, right = FALSE, - labels = c(paste0(seq(0, max_fetch - 1, by = 1), "-", seq(1, max_fetch, by = 1)), paste0(max_fetch, "+")) - ) - } else { - NULL - }, - PA_Factor = if (has_pa) factor(dat[, pa_col], labels = c("Absent", "Present")) else NULL, - Cover_Bin = if (has_cover) cut(dat[, cover_col], breaks = seq(0, 100, by = 10), include.lowest = TRUE, right = FALSE) else NULL - ) + dat <- dat |> + dplyr::mutate( + Depth_Bin = if (has_depth) { + cut( + dat$depth_m, + breaks = c(seq(0, max_depth, by = 1), Inf), + include.lowest = TRUE, right = FALSE, + labels = c(paste0(seq(0, max_depth - 1, by = 1), "-", seq(1, max_depth, by = 1)), paste0(max_depth, "+")) + ) + } else { + NULL + }, + Fetch_Bin = if (has_fetch) { + cut( + dat$fetch_km, + breaks = c(seq(0, max_fetch, by = 1), Inf), + include.lowest = TRUE, right = FALSE, + labels = c(paste0(seq(0, max_fetch - 1, by = 1), "-", seq(1, max_fetch, by = 1)), paste0(max_fetch, "+")) + ) + } else { + NULL + }, + PA_Factor = if (has_pa) { + factor(c("Absent", "Present")[dat[[pa_col]] + 1]) + } else NULL, + Cover_Bin = if (has_cover) cut(dat[, cover_col], breaks = seq(0, 100, by = 10), include.lowest = TRUE, right = FALSE) else NULL + ) # Define colors cover_palette <- c( @@ -178,7 +195,7 @@ plot_sav_distribution <- function(dat, type = c("pa", "cover"), predictors = c(" #' dat <- data.frame( #' depth_m = runif(100, 0, 15), #' fetch_km = runif(100, 0, 15), -#' pa = sample(0:1, 100, replace = TRUE), +#' pa_pred = sample(0:1, 100, replace = TRUE), #' pa_post_hoc = sample(0:1, 100, replace = TRUE) #' ) #' @@ -196,9 +213,10 @@ plot_sav_distribution <- function(dat, type = c("pa", "cover"), predictors = c(" #' @export plot_sav_density <- function(dat, predictors = c("depth", "fetch"), max_depth = 30, post_hoc = TRUE) { plots <- list() + dat <- as.data.frame(dat) # Determine which column to use based on post_hoc flag - pa_col <- if (post_hoc) "pa_post_hoc" else "pa" + pa_col <- if (post_hoc) "pa_post_hoc" else "pa_pred" # Check if PA column exists if (!pa_col %in% names(dat)) { @@ -207,13 +225,14 @@ plot_sav_density <- function(dat, predictors = c("depth", "fetch"), max_depth = # Check for requested predictors, abort if unavailable has_depth <- "depth_m" %in% names(dat) - if (!has_depth && "depth" %in% predictors) rlang::abort("Requested layer `depth` is unavailable in provided data)") + if (!has_depth && "depth" %in% predictors) + rlang::abort("Requested layer `depth` is unavailable in provided data)") has_fetch <- "fetch_km" %in% names(dat) - if (!has_fetch && "fetch" %in% predictors) rlang::abort("Requested layer `fetch` is unavailable in provided data)") - + if (!has_fetch && "fetch" %in% predictors) + rlang::abort("Requested layer `fetch` is unavailable in provided data)") # Convert PA to factor - dat$PA_Factor <- factor(dat[[pa_col]], labels = c("Absent", "Present")) + dat$PA_Factor <- factor(c("Absent", "Present")[dat[[pa_col]] + 1]) # Colors cols <- c("#56B4E9", "#52854C") @@ -287,8 +306,8 @@ plot_sav_density <- function(dat, predictors = c("depth", "fetch"), max_depth = #' @param layers Character vector specifying the layers to generate. Options: #' - `"pa"` (default) for presence/absence model predictions #' - `"cover"` (default) for cover percentage model predictions -#' - `"depth"` (default) for depth predictor vavlues -#' - `"fetch"` (default) for fetch predictor vavlues +#' - `"depth"` (default) for depth predictor values +#' - `"fetch"` (default) for fetch predictor values #' @param post_hoc Logical value indicating whether to use post-hoc analyzed columns (`pa_post_hoc`, `cover_post_hoc`) instead of raw columns (`pa`, `cover`). Default is `FALSE`. #' @param interactive Logical. If `TRUE` (default), generates an interactive map. If `FALSE`, creates a static map. #' @param export_path Character. If provided, saves the map to the specified file path. Default is `NULL`. @@ -307,8 +326,8 @@ plot_sav_density <- function(dat, predictors = c("depth", "fetch"), max_depth = #' points <- sf::st_sample(polygon, 100) |> #' sf::st_sf() |> #' dplyr::mutate( -#' cover = runif(100, 0, 100), -#' pa = sample(0:1, 100, replace = TRUE) |> +#' cover_pred = runif(100, 0, 100), +#' pa_pred = sample(0:1, 100, replace = TRUE) |> #' factor(levels = c(0, 1), labels = c("Absent", "Present")), #' depth_m = runif(100, 0, 15), #' fetch_km = runif(100, 0, 10), @@ -345,8 +364,8 @@ plot_sav_tmap <- function(study_zone, layers = c("pa", "cover", "depth", "fetch" cols <- c("#56B4E9", "#52854C") # Determine which columns to use based on post_hoc flag - pa_col <- if (post_hoc) "pa_post_hoc" else "pa" - cover_col <- if (post_hoc) "cover_post_hoc" else "cover" + pa_col <- if (post_hoc) "pa_post_hoc" else "pa_pred" + cover_col <- if (post_hoc) "cover_post_hoc" else "cover_pred" # Check for requested columns, abort if unavailable has_depth <- "depth_m" %in% names(study_zone$points) @@ -375,7 +394,7 @@ plot_sav_tmap <- function(study_zone, layers = c("pa", "cover", "depth", "fetch" if ("cover" %in% layers) { map <- map + tmap::tm_shape(study_zone$points) + tmap::tm_dots( - col = "cover", + col = cover_col, group = "Cover", size = .5, fill.scale = tmap::tm_scale("viridis"), @@ -385,14 +404,15 @@ plot_sav_tmap <- function(study_zone, layers = c("pa", "cover", "depth", "fetch" } if ("pa" %in% layers) { - if (!is.factor(study_zone$points$pa)) { - study_zone$points <- study_zone$points |> - dplyr::mutate(pa = factor(pa, levels = c(0, 1), labels = c("Absent", "Present"))) + if (!is.factor(study_zone$points[[pa_col]])) { + study_zone$points[[pa_col]] <- factor( + c("Absent", "Present")[study_zone$points[[pa_col]] + 1] + ) } map <- map + tmap::tm_shape(study_zone$points) + tmap::tm_dots( - col = "pa", + col = pa_col, group = "Presence", size = 0.5, palette = cols diff --git a/R/preview_grid.R b/R/preview_grid.R new file mode 100644 index 0000000..a2af78e --- /dev/null +++ b/R/preview_grid.R @@ -0,0 +1,18 @@ +#' Preview grid +#' +#' @param x an object of class `sav_data`. +#' +#' @export +#' +preview_grid <- function(x) { + UseMethod("preview_grid") +} + +#' @describeIn preview_grid Preview spatail grid. +#' +#' @export +#' +preview_grid.sav_data <- function(x) { + plot(x$polygon |> sf::st_geometry(), border = 1) + plot(x$points |> sf::st_geometry(), col = "grey50", pch = 19, add = TRUE) +} \ No newline at end of file diff --git a/R/data_input.R b/R/read.R similarity index 93% rename from R/data_input.R rename to R/read.R index e47a679..575c5ad 100644 --- a/R/data_input.R +++ b/R/read.R @@ -1,11 +1,17 @@ #' Read SAV data from different formats #' #' @param file_path {`character`}\cr{} Path to the input file. -#' @param spacing {`numeric`}\cr{} Spacing for grid generation (if AOI is used). +#' @param spacing {`numeric`}\cr{} +#' Spacing for grid generation. Used if an area of interest (AOI) is provided, +#' ignored otherwise. #' @param layer {`character`}\cr{} Layer name for multi-layer spatial files (default: NULL). -#' @param crs {`numeric`}\cr{} Coordinate Reference System (CRS) of the output data. Default is 32617. -#' @param crs_input {`numeric`}\cr{} Coordinate Reference System (CRS) of the input data, used for CSV import. Default is 4326. -#' @param export {`character`}\cr{} Optional. Folder path to export outputs as .gpkg. +#' @param crs {`numeric`}\cr{} +#' Coordinate Reference System (CRS) of the output data. Default is 32617. +#' @param crs_input {`numeric`}\cr{} +#' Coordinate Reference System (CRS) of the input data, used for CSV import. +#' Default is 4326. +#' @param export {`character`}\cr{} +#' Optional. Folder path to export outputs as .gpkg. #' #' @return A list with `points` (sf object) and `polygon` (sf object). #' @@ -52,12 +58,15 @@ read_sav <- function(file_path, spacing = 500, layer = NULL, crs = 32617, crs_in file_ext <- tools::file_ext(file_path) if (file_ext == "csv") { + sav_msg_info("csv detected") points_sf <- read_sav_csv(file_path, crs = crs) polygon_sf <- sf::st_sf(geometry = sf::st_union(points_sf) |> sf::st_convex_hull()) } else if (file_ext %in% c("shp", "geojson", "gpkg", "gbd")) { if (file_ext == "gbd") { + sav_msg_info("gdb detected") sf_obj <- sf::st_read(file_path, layer = layer, quiet = TRUE) } else { + sav_msg_info("spatial file detected") sf_obj <- sf::st_read(file_path, quiet = TRUE) } @@ -83,7 +92,12 @@ read_sav <- function(file_path, spacing = 500, layer = NULL, crs = 32617, crs_in sav_msg_success("Exported outputs to {export}.") } - return(list(points = points_sf, polygon = polygon_sf)) + return( + structure( + list(points = points_sf, polygon = polygon_sf), + class = "sav_data" + ) + ) } diff --git a/R/zzz.R b/R/zzz.R index a5c6477..6786e64 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,11 +1,13 @@ #' Helper functions #' +#' @import randomForest +#' #' @noRd #' globalVariables(c( "direction", "fetch", "id_point", "weight","Cover_Bin", "Depth_Bin", "Fetch_Bin", "Mean_Value", "PA_Factor", "depth_m", "fetch_km", - "limitation_secchi", "transect_length", "vmax", "pa" + "limitation_secchi", "transect_length", "vmax", "pa", "geometry" )) diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.cpg b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.cpg new file mode 100644 index 0000000..3ad133c --- /dev/null +++ b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.cpg @@ -0,0 +1 @@ +UTF-8 \ No newline at end of file diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.dbf b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.dbf new file mode 100644 index 0000000..9b027f0 Binary files /dev/null and b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.dbf differ diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.prj b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.prj new file mode 100644 index 0000000..f45cbad --- /dev/null +++ b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.prj @@ -0,0 +1 @@ +GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] \ No newline at end of file diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbn b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbn new file mode 100644 index 0000000..16bca23 Binary files /dev/null and b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbn differ diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbx b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbx new file mode 100644 index 0000000..b709bc6 Binary files /dev/null and b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbx differ diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp new file mode 100644 index 0000000..8f854a3 Binary files /dev/null and b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp differ diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp.xml b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp.xml new file mode 100644 index 0000000..73a8535 --- /dev/null +++ b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp.xml @@ -0,0 +1,3 @@ + +gl_high_resolution_islands_nhd_wbd_nhn_ohnraster digital dataUnknown-93.456528-72.47206851.29027839.282111NoneGreat Lakes, Laurentian Great Lakes, frameworkISO 19115 Topic CategoriesinlandWatersplanningCadastreNoneGreat Lakes, Laurentian Great Lakes, Erie, Huron, Michigan, Ontario, Superior, St. Clair, Niagara River, St. Marys River, Ontario, Minnesota, Wisconsin, Illinois, Indiana, Michigan, Ohio, Pennsylvania, New York, frameworkNoneNoneUniversity of Michigan SNRE, Institute for Fisheries ResearchLacey MasonComputer Research Specialistmailing and physical
400 North Ingalls Building, NIB 400 G250
Ann ArborMI48109-5480US734-663-3554 x.12155lmas@umich.eduUniversity of Michigan SNRECatherine RisengResearch Associatemailing and physical
440 Church St., G170 Dana
Ann ArborMI48109-1041US734-763-9422criseng@umich.eduMichigan DNR, Insitute for Fisheries ResearchDanielle ForsythResource Analystmailing and physical
400 North Ingalls Building, NIB 400 G250
Ann ArborMI48109-5480US734-663-3554 x.11755ForsythD1@michigan.govUniversity of Michigan SNRE, Institute for Fisheries ResearchLacey MasonComputer Research Specialistmailing and physical
400 North Ingalls Building, NIB 400 G250
Ann ArborMI48109-5480US734-663-3554 x.12155lmas@umich.eduMichigan DNR, Insitute for Fisheries ResearchDanielle ForsythResource Analystmailing and physical
400 North Ingalls Building, NIB 400 G250
Ann ArborMI48109-5480US734-663-3554 x.11755ForsythD1@michigan.govThe Great Lakes Aquatic Habitat Framework (GLAHF) project has been funded by the Great Lakes Fishery Trust and lead by the University of Michigan School of Natural Resources and Environment, with partners from Michigan Department of Natural Resources-Institute for Fisheries Research, NOAA Great Lakes Environmental Research Laboratory, International Joint Commission, Michigan State University, The Nature Conservancy, Ontario Ministry of Natural Resources, University of Minnesota-Duluth, U.S. Fish & Wildlife Service, U.S. Geological Survey and many collaborating partners in both the USA and Canada. More information about this project can be found at http://ifrgis.snre.umich.edu/projects/GLAHF/glahf.shtml.Microsoft Windows 7 Version 6.1 (Build 7601) Service Pack 1; Esri ArcGIS 10.1.1.3143RasterGrid CellSimpleFALSE0FALSETRUED North American 1983GRS 19806378137.0298.257222101VAT_metadata_basicsFeature Class0FIDFIDOID400Internal feature number.EsriSequential unique whole numbers that are automatically generated.ShapeShapeGeometry000Feature geometry.EsriCoordinates defining the features.FID_LakeErFID_LakeErInteger10100OBJECTID_1OBJECTID_1Integer10100Internal feature number.EsriSequential unique whole numbers that are automatically generated.GNIS_NameGNIS_NameString6500AreaSqKmAreaSqKmDouble1900Shape_LengShape_LengDouble1900Shape_AreaShape_AreaDouble1900Area of feature in internal units squared.EsriPositive real numbers that are automatically generated.FID_Lake_1FID_Lake_1Integer10100IdIdInteger66020140806University of Michigan SNRE, Institute for Fisheries ResearchLacey MasonComputer Research Specialistmailing and physical
400 North Ingalls Building, NIB 400 G250
Ann ArborMI48109-5480US734-663-3554 x.12155lmas@umich.eduFGDC Content Standard for Digital Geospatial MetadataFGDC-STD-001-1998local time Version 6.2 (Build 9200) ; Esri ArcGIS 10.6.1.9270shoreline v1.10.0000000.0000001-92.542434-73.17064549.19073740.237439In order to create an updated high resolution shoreline layer (including islands) for the Great Lakes, the Great Lakes and connecting channel open water polygons from the NHD 1:24,000 were used to create a new shoreline, with a few exceptions (see below). A corresponding high resolution island layer was created by erasing the finalized shoreline layers from a polygon that covered the extent of the shoreline. This shoreline includes additional EPA bays that were enforced into the GLAHF Hydrology Data V1 data package, and follows the extent of the GLAHF Hydrology Data V1 layers, which extend far into the St. Lawrence Seaway. +<DIV STYLE="text-align:Left;"><DIV><DIV><P><SPAN>The following outlines the deviations in this dataset from the NHD 1:24,000 open water polygons:</SPAN></P><P><SPAN /></P><P><SPAN>1. NHD 1:24,000 data did not exist for the Canadian side of Lake Huron. For this shoreline, the summarized NHN 1:20,0000 shoreline was used from the synchronized Watershed Boundary Dataset (NHN data is available in ~50 work units for the Great Lakes and would have been time consuming to process)</SPAN></P><P><SPAN /></P><P><SPAN>2. The NHN shoreline included in the WBD did not cover the St. Mary’s River, so a shoreline layer provided by the Ontario Ministry of Natural Resources (OMNR) that was created from the Ontario Hydro Network (OHN) was used. </SPAN></P><P><SPAN /></P><P><SPAN>3. Islands for the St. Mary’s River and some larger islands in Lake Huron were added from the NHN</SPAN></P><P><SPAN /></P><P><SPAN>4. A small number of islands that were notably missing were digitized and added to the Lake Michigan layer at the 1:24,000 scale.</SPAN></P><P><SPAN /></P><P><SPAN>5. A set of EPA bays was used to add additional, prominent and bays and some river mouths to the Great Lakes shoreline basin wide.</SPAN></P><P><SPAN /></P><P><SPAN>Year data layer created: 2014</SPAN></P><P><SPAN /></P><P><SPAN>Data Download Date:</SPAN></P><P><SPAN /></P><P><SPAN>National Hydrography Dataset 1:24,000 (1/13/2013)</SPAN></P><P><SPAN /></P><P><SPAN>Watershed Boundary Dataset (acquired 1/7/13 from Kimberly Jones, USGS kjones@usgs.gov)</SPAN></P><P><SPAN /></P><P><SPAN>NHN (downloaded 1/23/2013)</SPAN></P><P><SPAN /></P><P><SPAN>OHN Shoreline (acquired 12/16/2013 from John Gaiot, OMNR john.gaiot@ontario.ca)</SPAN></P></DIV></DIV></DIV>The Great Lakes Aquatic Habitat Framework (GLAHF) project has been funded by the Great Lakes Fishery Trust and lead by the University of Michigan School of Natural Resources and Environment, with partners from Michigan Department of Natural Resources-Institute for Fisheries Research, NOAA Great Lakes Environmental Research Laboratory, International Joint Commission, Michigan State University, The Nature Conservancy, Ontario Ministry of Natural Resources, University of Minnesota-Duluth, U.S. Fish & Wildlife Service, U.S. Geological Survey and many collaborating partners in both the USA and Canada. More information about this project can be found at http://ifrgis.snre.umich.edu/projects/GLAHF/glahf.shtml.Great LakesLaurentian Great LakesframeworkinlandWatersplanningCadastre<DIV STYLE="text-align:Left;"><DIV><DIV><P><SPAN>None</SPAN></P></DIV></DIV></DIV>LakeErie_LandWater_Union002284967.9116001718112.8696002099886.4391002914949.8719001file://\\BUR02HNAS01A.ENT.dfo-mpo.ca\gllfas\Fish Ecology Science\SAV_Tool\SAV Datasets\Fetch_Calc\LakePolygons\LakeErie_LandWater_Union.shpLocal Area Network0.000ProjectedGCS_North_American_1983Linear Unit: Meter (1.000000)<ProjectedCoordinateSystem xsi:type='typens:ProjectedCoordinateSystem' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/10.6'><WKT>PROJCS[&quot;USA_Contiguous_Albers_Equal_Area_Conic_USGS_version&quot;,GEOGCS[&quot;GCS_North_American_1983&quot;,DATUM[&quot;D_North_American_1983&quot;,SPHEROID[&quot;GRS_1980&quot;,6378137.0,298.257222101]],PRIMEM[&quot;Greenwich&quot;,0.0],UNIT[&quot;Degree&quot;,0.0174532925199433]],PROJECTION[&quot;Albers&quot;],PARAMETER[&quot;False_Easting&quot;,0.0],PARAMETER[&quot;False_Northing&quot;,0.0],PARAMETER[&quot;Central_Meridian&quot;,-96.0],PARAMETER[&quot;Standard_Parallel_1&quot;,29.5],PARAMETER[&quot;Standard_Parallel_2&quot;,45.5],PARAMETER[&quot;Latitude_Of_Origin&quot;,23.0],UNIT[&quot;Meter&quot;,1.0],AUTHORITY[&quot;Esri&quot;,102039]]</WKT><XOrigin>-16901100</XOrigin><YOrigin>-6972200</YOrigin><XYScale>266467840.99085236</XYScale><ZOrigin>-100000</ZOrigin><ZScale>10000</ZScale><MOrigin>-100000</MOrigin><MScale>10000</MScale><XYTolerance>0.001</XYTolerance><ZTolerance>0.001</ZTolerance><MTolerance>0.001</MTolerance><HighPrecision>true</HighPrecision><WKID>102039</WKID><LatestWKID>102039</LatestWKID></ProjectedCoordinateSystem>USA_Contiguous_Albers_Equal_Area_Conic_USGS_versionUnion "LakeErie_HEC_Poly_WGS #;LakeErie_Land_forErase #" "I:\Fish Ecology Science\SAV_Tool\SAV Datasets\Fetch_Calc\LakePolygons\LakeErie_LandWater_Union.shp" ALL # GAPSProject LkErie_Land_fromGLAF_Water_Feb2020 "I:\Fish Ecology Science\SAV_Tool\SAV Datasets\GIS\LakePolygons\LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp" GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]] WGS_1984_(ITRF00)_To_NAD_1983 PROJCS['USA_Contiguous_Albers_Equal_Area_Conic_USGS_version',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['Latitude_Of_Origin',23.0],UNIT['Meter',1.0]] NO_PRESERVE_SHAPE # NO_VERTICAL202002041158520020200204115852001500000005000FGDC1.02020021110384600FGDC CSDGM MetadataFALSEShapefile0.000datasetEsri8.1.2020200204Catherine RisengUniveristy of Michigan, Michigan Sea GrantG170 Dana, School of Natural Resources and Environment, 440 Church StAnn ArborMichigan48109-1041criseng@umich.eduUS734-763-9422Lacey MasonUniversity of MichiganInstitute for Fisheries Research, 400 North Ingalls Building, NIB G250Ann ArborMichigan48109-5480lmas@umich.eduUS734.663.3554 ext. 12155Danielle ForsythMichigan Department of Natural ResourcesInstitute for Fisheries Research, 400 North Ingalls Building, NIB G250Ann ArborMichigan48109-5480ForsythD1@Michigan.govUS734.663.3554 ext. 11755Lacey Mason diff --git a/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx new file mode 100644 index 0000000..7b8dfc8 Binary files /dev/null and b/inst/example/lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx differ diff --git a/man/compute_fetch.Rd b/man/compute_fetch.Rd index 6783a78..df923f9 100644 --- a/man/compute_fetch.Rd +++ b/man/compute_fetch.Rd @@ -8,7 +8,7 @@ compute_fetch( points, polygon, max_dist = 15, - n_quad_seg = 9, + n_bearings = 16, wind_weights = NULL, crs = NULL ) @@ -23,9 +23,9 @@ Polygon defining land boundaries used to compute fetch distances.} \item{max_dist}{{\code{numeric}}\cr{} Maximum fetch distance in kilometers. Fetch beyond this distance is capped.} -\item{n_quad_seg}{{\code{integer}}\cr{} -Number of segments per quadrant for fetch calculation. -Ignored if \code{wind_weights} is provided.} +\item{n_bearings}{{\code{integer}}\cr{} +Total number of bearings for fetch calculation (minimal number required is +4, default is 167). Ignored if \code{wind_weights} is provided.} \item{wind_weights}{{\code{data.frame}}\cr{} A data frame specifying directional weights for wind exposure. @@ -39,13 +39,11 @@ transform \code{points} and \code{polygon}.} \value{ A list of two elements: \itemize{ -\item \code{mean_fetch}: data frame with 5 columns: +\item \code{mean_fetch}: a \code{sf} object with 3 features: \itemize{ \item \code{id_point}: point identifier -\item \code{fetch}: mean wind fetch based on the four highest values -\item \code{weighted_fetch}: mean weighted wind fetch based on the four highest values -\item \code{fetch_all}: mean wind fetch based on all values -\item \code{weighted_fetch_all}: mean wind weighted fetch based on all values +\item \code{fetch_km}: mean wind fetch based on all bearings. +\item \code{weighted_fetch_km}: mean weighted wind fetch based on all bearings. } \item \code{transect_lines}: a \code{sf} object containing all radial transect with the same columns as \code{points} and the following additional columns: @@ -67,7 +65,7 @@ across a body of water before reaching a specific point. It plays a crucial role in wave generation, as longer fetch distances allow wind to transfer more energy to the water surface, leading to larger waves. -For all points in \code{points}, 4 Γ— \code{n_quad_seg} radial transects are generated +For all points in \code{points}, \code{n_bearings} radial transects are generated by default. If \code{wind_weights} is specified, the column \code{direction}, which contains angles in degrees, is used instead to generate the transects. The transects are then clipped with the polygon using \code{\link[sf:geos_binary_ops]{sf::st_intersection()}}, @@ -77,28 +75,29 @@ all clipped transects is computed using \code{\link[sf:geos_measures]{sf::st_len element in the returned list and it used to generate the second element: \code{mean_fetch} that included wind fetch averages. -Ensure that max_dist is specified in meters. An error will be thrown if the +Ensure that \code{max_dist} is specified in meters. An error will be thrown if the spatial projection of points and polygon is not in a meter-based coordinate system. } \examples{ \donttest{ - le_bound <- system.file("example", "lake_erie.gpkg", package = "SAVM") |> sf::st_read() le_pt <- system.file("example", "le_points.geojson", package = "SAVM") |> sf::st_read(quiet = TRUE) res <- compute_fetch(le_pt, le_bound, crs = 32617) -# use wind-weight +# use wind-weight res2 <- compute_fetch( - le_pt, le_bound, max_dist = 20, - wind_weights = data.frame( - direction = seq(0, 360, by = 360 / 16)[-1], - weight = rep(c(0, 1), each = 8) - ), - crs = 32617) + le_pt, le_bound, + max_dist = 20, + wind_weights = data.frame( + direction = seq(0, 360, by = 360 / 16)[-1], + weight = rep(c(0, 1), each = 8) + ), + crs = 32617 +) -# resultat +# results res$mean_fetch res2$mean_fetch diff --git a/man/data_input.Rd b/man/data_input.Rd index 2519c77..5fd06ce 100644 --- a/man/data_input.Rd +++ b/man/data_input.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data_input.R +% Please edit documentation in R/read.R \name{read_sav} \alias{read_sav} \alias{read_sav_csv} @@ -34,7 +34,8 @@ of the output data (see \code{\link[sf:st_crs]{sf::st_crs()}}). Default is 32617 \item{crs_input}{{\code{numeric}}\cr{} Coordinate Reference System (CRS) of the input data, used for CSV import. Default is 4326.} -\item{export}{{\code{character}}\cr{} Optional. Folder path to export outputs as .gpkg.} +\item{export}{{\code{character}}\cr{} +Optional. Folder path to export outputs as .gpkg.} \item{...}{Further arguments passed on to \code{\link[utils:read.table]{utils::read.csv()}}.} } diff --git a/man/invert_polygon.Rd b/man/invert_polygon.Rd new file mode 100644 index 0000000..447e47a --- /dev/null +++ b/man/invert_polygon.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/invert_polygon.R +\name{invert_polygon} +\alias{invert_polygon} +\title{Invert a Polygon} +\usage{ +invert_polygon(polygon, ratio = 0.5) +} +\arguments{ +\item{polygon}{{\code{sf}}\cr{} +Polygon defining land boundaries that needs to be inverted.} + +\item{ratio}{{\code{numeric [0-1]}}\cr{} +Fraction convex, see \code{\link[sf:geos_unary]{sf::st_concave_hull()}}. This may need to be tweaked +depending on the form of the polygon to be inverted.} +} +\description{ +Invert a Polygon +} +\details{ +Utility function that inverts a polygon by drawing a concave hull around +\code{polygon} (see \code{\link[sf:geos_unary]{sf::st_concave_hull()}}) and then computes the differences +between the polygon and the concave hull (see \code{\link[sf:geos_unary]{sf::st_concave_hull()}}). +} +\examples{ +\donttest{ +erie_land <- system.file( + "example", "lake_erie_land", "LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx", + package = "SAVM", mustWork = TRUE +) |> sf::st_read() + +erie_land |> + sf::st_geometry() |> + plot(col = 1) + +erie_land |> + sf::st_geometry() |> + invert_polygon() |> + plot(col = 1) +} +} diff --git a/man/plot_sav.Rd b/man/plot_sav.Rd index dc2a255..72c212a 100644 --- a/man/plot_sav.Rd +++ b/man/plot_sav.Rd @@ -38,7 +38,8 @@ plot_sav_tmap( \item \code{pa}: Binary (0 = absent, 1 = present), indicating SAV presence/absence. }} -\item{type}{Character vector specifying the type of plots to generate. Options: +\item{type}{{\verb{character vector}}\cr{} +Character vector specifying the type of plots to generate. Options: \itemize{ \item \code{"pa"} (default) for presence/absence plots \item \code{"cover"} (default) for cover percentage plots @@ -54,7 +55,8 @@ plot_sav_tmap( \item{max_depth}{Numeric value specifying the maximum depth bin (default: 30 meters).} -\item{max_fetch}{Numeric value specifying the maximum fetch bin (default: 15 km).} +\item{max_fetch}{{\code{numeric}}\cr{} +Numeric value specifying the maximum fetch bin (default: 15 km).} \item{study_zone}{A list containing spatial objects: \itemize{ @@ -66,8 +68,8 @@ plot_sav_tmap( \itemize{ \item \code{"pa"} (default) for presence/absence model predictions \item \code{"cover"} (default) for cover percentage model predictions -\item \code{"depth"} (default) for depth predictor vavlues -\item \code{"fetch"} (default) for fetch predictor vavlues +\item \code{"depth"} (default) for depth predictor values +\item \code{"fetch"} (default) for fetch predictor values }} \item{interactive}{Logical. If \code{TRUE} (default), generates an interactive map. If \code{FALSE}, creates a static map.} @@ -82,9 +84,10 @@ One or two ggplot2 density plots visualizing SAV presence/absence by depth and/o A \code{tmap} map object, optionally saved to a file. } \description{ -This function generates up to four plots representing the distribution of submerged aquatic vegetation (SAV) -data by depth (m) and fetch (km). It visualizes SAV presence/absence (PA) and percent cover (Cover) when -the corresponding columns are available in the input data. +This function generates up to four plots representing the distribution of +submerged aquatic vegetation (SAV) data by depth (m) and fetch (km). It +visualizes SAV presence/absence (PA) and percent cover (Cover) when the +corresponding columns are available in the input data. This function generates one or two density plots for submerged aquatic vegetation (SAV) presence (green) and absence (blue) as a function of depth (m) and/or fetch (km). It includes vertical dotted lines @@ -98,8 +101,8 @@ visualizing cover, presence/absence, depth, and fetch values. dat <- data.frame( depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 15), - pa = sample(0:1, 100, replace = TRUE), - cover = runif(100, 0, 100), + pa_pred = sample(0:1, 100, replace = TRUE), + cover_pred = runif(100, 0, 100), pa_post_hoc = sample(0:1, 100, replace = TRUE), cover_post_hoc = runif(100, 0, 100) ) @@ -120,7 +123,7 @@ plot_sav_distribution(dat, post_hoc = TRUE) dat <- data.frame( depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 15), - pa = sample(0:1, 100, replace = TRUE), + pa_pred = sample(0:1, 100, replace = TRUE), pa_post_hoc = sample(0:1, 100, replace = TRUE) ) @@ -144,8 +147,8 @@ set.seed(42) points <- sf::st_sample(polygon, 100) |> sf::st_sf() |> dplyr::mutate( - cover = runif(100, 0, 100), - pa = sample(0:1, 100, replace = TRUE) |> + cover_pred = runif(100, 0, 100), + pa_pred = sample(0:1, 100, replace = TRUE) |> factor(levels = c(0, 1), labels = c("Absent", "Present")), depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 10), diff --git a/man/preview_grid.Rd b/man/preview_grid.Rd new file mode 100644 index 0000000..72b1f9d --- /dev/null +++ b/man/preview_grid.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/preview_grid.R +\name{preview_grid} +\alias{preview_grid} +\alias{preview_grid.sav_data} +\title{Preview grid} +\usage{ +preview_grid(x) + +\method{preview_grid}{sav_data}(x) +} +\arguments{ +\item{x}{an object of class \code{sav_data}.} +} +\description{ +Preview grid +} +\section{Methods (by class)}{ +\itemize{ +\item \code{preview_grid(sav_data)}: Preview spatail grid. + +}} diff --git a/man/sav_model.Rd b/man/sav_model.Rd index 2b51dc5..08a38be 100644 --- a/man/sav_model.Rd +++ b/man/sav_model.Rd @@ -16,7 +16,7 @@ sav_model( ) } \arguments{ -\item{dat}{{\code{data.frame}}\cr{} A \code{data.frame} containing some or all of the +\item{dat}{{\code{data.frame}|\code{sf}}\cr{} A \code{data.frame} or a \code{sf} object containing some or all of the following columns: \itemize{ \item \code{depth_m}: Numeric, depth in meters. @@ -26,16 +26,16 @@ following columns: limitations. (post_hoc) \item \code{limitation}: Binary (0 = absent, 1 = present), indicating user-supplied limitations. -Additionnal columns will be ignored. +Additional columns will be ignored. }} \item{type}{{\verb{character vector}, either \code{"cover"} or \code{"pa"}}\cr{} Model type(s).} -\item{depth, fetch}{{\code{character}} Column specification for the predictors, +\item{depth, fetch}{{\code{character}}\cr{} Column specification for the predictors, see \emph{Details}.} -\item{substrate, secchi, limitation}{Column specification for post_hoc +\item{substrate, secchi, limitation}{{\code{character}}\cr{}Column specification for post_hoc variables, see \emph{Details}.} \item{vmax_par}{{\verb{named list}}\cr{} intercept and slope of the equation from @@ -43,16 +43,18 @@ Chambers and Kalff (1985) to compute the maximum depth of plant colonization (Vmax), see \emph{Details} below.} } \value{ -A data frame containing the input columns along with model predictions. +A data frame (or a sf object) containing the input columns along with model +predictions. -The prediction column names match the values specified in \code{type}, and contain -the raw model outputs (i.e., without post-hoc adjustment). +The prediction column names match the values specified in \code{type} followed +by the suffix \verb{_pred}, and contain the raw model outputs (i.e., without +post-hoc adjustment). Post-hoc adjusted predictions (see \emph{Details}) are included in additional columns with the same names as the \code{type} values, but with the suffix \verb{_post_hoc}. -If a column \code{secchi} is present, then two additionnal columns are +If a column \code{secchi} is present, then two additional columns are returned: \code{vmax} and \code{limitation_secchi}, see details for further explanation. } @@ -105,7 +107,7 @@ corresponding prediction is set to 0 in the post-hoc adjusted output. sav_model(data.frame(depth = c(5, 10))) sav_model(data.frame(depth = c(5, 10), fetch = c(1, 2)), type = "pa") -# using post-hoc tretment +# using post-hoc treatment sav_model( data.frame( depth = c(5, 10, 5), diff --git a/tests/testthat/_snaps/plot/plot-sav-tmap.svg b/tests/testthat/_snaps/plot/plot-sav-tmap.svg index bb5afe6..2d783b4 100644 --- a/tests/testthat/_snaps/plot/plot-sav-tmap.svg +++ b/tests/testthat/_snaps/plot/plot-sav-tmap.svg @@ -18,143 +18,143 @@ - - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - -cover -0 to 10 -10 to 20 -20 to 30 -30 to 40 -40 to 50 -50 to 60 -60 to 70 -70 to 80 -80 to 90 -90 to 100 -Missing - - - - - - - - - - - + + +cover_post_hoc +0 to 10 +10 to 20 +20 to 30 +30 to 40 +40 to 50 +50 to 60 +60 to 70 +70 to 80 +80 to 90 +90 to 100 +Missing + + + + + + + + + + + + diff --git a/tests/testthat/test-compute-fetch.R b/tests/testthat/test-compute-fetch.R index dfe9a59..6a61a94 100644 --- a/tests/testthat/test-compute-fetch.R +++ b/tests/testthat/test-compute-fetch.R @@ -45,13 +45,21 @@ test_that("helpers work", { test_that("helpers work", { withr::with_options( list(savm.verbose = "quiet"), - { + { + expect_error( + compute_fetch(le_pt, le_bound, n_bearings = 3), + "`n_bearings` should be equal or greater than 4." + ) + expect_error( + compute_fetch(le_pt, le_bound, max_dist = 0), + "`max_dist` must be strictly positive." + ) expect_error( - compute_fetch(le_pt, le_bound, 10, max_dist = 15000), + compute_fetch(le_pt, le_bound, max_dist = 15000, n_bearings = 4), "Projection units must be meters." ) expect_error( - compute_fetch(le_pt_out, le_bound_merc, 10, max_dist = 15000), + compute_fetch(le_pt_out, le_bound_merc), "`polygon` must include `points`" ) } @@ -63,15 +71,15 @@ test_that("compute_fetch() work", { list(savm.verbose = "quiet"), { res <- compute_fetch(le_pt, le_bound_merc) + expect_s3_class(res$mean_fetch, "sf") expect_identical(names(res), c("mean_fetch", "transect_lines")) expect_identical( names(res$mean_fetch), c( - "id_point", "fetch_km", "weighted_fetch_km", "fetch_km_all", - "weighted_fetch_km_all" + "id_point", "fetch_km", "weighted_fetch_km", "geometry" ) ) - expect_true(inherits(res$transect_lines, "sf")) + expect_s3_class(res$transect_lines, "sf") expect_true( all( c("direction", "weight", "transect_length", "rank") %in% @@ -82,8 +90,6 @@ test_that("compute_fetch() work", { res2 <- compute_fetch(le_pt_mid, le_bound_merc, max_dist = 15) expect_equal(res2$mean_fetch$fetch_km, 15) expect_equal(res2$mean_fetch$weighted_fetch_km, 15) - expect_equal(res2$mean_fetch$fetch_km_all, 15) - expect_equal(res2$mean_fetch$weighted_fetch_km_all, 15) } ) }) @@ -100,8 +106,6 @@ test_that("compute_fetch() with wind_weight", { ) expect_equal(res$mean_fetch$fetch_km, 15) expect_equal(res$mean_fetch$weighted_fetch_km, 15) - expect_equal(res$mean_fetch$fetch_km_all, 15) - expect_equal(res$mean_fetch$weighted_fetch_km_all, 15) # res2 <- compute_fetch(le_pt[1, ], le_bound_merc, wind_weights = data.frame( @@ -109,10 +113,8 @@ test_that("compute_fetch() with wind_weight", { weight = rep(c(0, 1), each = 8) ) ) - expect_equal(round(res2$mean_fetch$fetch_km, 4), 1.0327) - expect_equal(round(res2$mean_fetch$weighted_fetch_km, 4), 0.2533) - expect_equal(round(res2$mean_fetch$fetch_km_all, 4), 2.3682) - expect_equal(round(res2$mean_fetch$weighted_fetch_km_all, 4), 1.6387) + expect_equal(round(res2$mean_fetch$fetch_km, 4), 2.3682) + expect_equal(round(res2$mean_fetch$weighted_fetch_km, 4), 1.6387) } ) }) diff --git a/tests/testthat/test-plot.R b/tests/testthat/test-plot.R index 9b6a957..1d4bb12 100644 --- a/tests/testthat/test-plot.R +++ b/tests/testthat/test-plot.R @@ -4,8 +4,8 @@ withr::with_seed(123, { test_data <- data.frame( depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 15), - pa = sample(0:1, 100, replace = TRUE), - cover = runif(100, 0, 100), + pa_pred = sample(0:1, 100, replace = TRUE), + cover_pred = runif(100, 0, 100), pa_post_hoc = sample(0:1, 100, replace = TRUE), cover_post_hoc = runif(100, 0, 100) ) @@ -83,14 +83,14 @@ test_that("Function errors if post_hoc columns are missing", { }) test_that("Function works with only pa and fetch_km", { - minimal_data <- test_data[, c("pa", "fetch_km")] + minimal_data <- test_data[, c("pa_pred", "fetch_km")] pdf(NULL) expect_silent(plot_sav_distribution(minimal_data, type = "pa", predictors = "fetch", post_hoc = FALSE)) dev.off() }) test_that("Function works with only cover and depth", { - data_subset <- test_data[, c("depth_m", "cover", "cover_post_hoc")] + data_subset <- test_data[, c("depth_m", "cover_pred", "cover_post_hoc")] pdf(NULL) expect_silent(plot_sav_distribution(data_subset, type = "cover", predictors = "depth")) dev.off() @@ -124,7 +124,7 @@ withr::with_seed(123, { test_data <- data.frame( depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 15), - pa = sample(0:1, 100, replace = TRUE), + pa_pred = sample(0:1, 100, replace = TRUE), pa_post_hoc = sample(0:1, 100, replace = TRUE) ) }) @@ -191,8 +191,8 @@ withr::with_seed(123, { points <- sf::st_sample(polygon, 100) |> sf::st_sf() |> dplyr::mutate( - cover = runif(100, 0, 100), - pa = sample(0:1, 100, replace = TRUE), + cover_pred = runif(100, 0, 100), + pa_pred = sample(0:1, 100, replace = TRUE), depth_m = runif(100, 0, 15), fetch_km = runif(100, 0, 10), pa_post_hoc = sample(0:1, 100, replace = TRUE), diff --git a/tests/testthat/test_model.R b/tests/testthat/test_model.R index c7152ad..9e9f001 100644 --- a/tests/testthat/test_model.R +++ b/tests/testthat/test_model.R @@ -12,7 +12,7 @@ df_ok_1 <- data.frame(depth = c(5, 10)) df_ok_2 <- data.frame(FETCH_km = c(1, 2)) df_ok_3 <- cbind(df_ok_1, df_ok_2) -res1 <- structure(list(depth_m = c(5, 10), pa = 1:0, cover = c( +res1 <- structure(list(depth_m = c(5, 10), pa_pred = 1:0, cover_pred = c( 85.4316666666667, 21.144 ), pa_post_hoc = 1:0, cover_post_hoc = c( @@ -20,21 +20,21 @@ res1 <- structure(list(depth_m = c(5, 10), pa = 1:0, cover = c( 0 )), row.names = c(NA, -2L), class = "data.frame") -res2 <- structure(list(fetch_km = c(1, 2), pa = 0:1, cover = c( +res2 <- structure(list(fetch_km = c(1, 2), pa_pred = 0:1, cover_pred = c( 78.8616666666666, 98.5113333333333 ), pa_post_hoc = 0:1, cover_post_hoc = c(0, 98.5113333333333)), row.names = c(NA, -2L), class = "data.frame") -res3 <- structure(list(depth_m = c(5, 10), fetch_km = c(1, 2), pa = c( +res3 <- structure(list(depth_m = c(5, 10), fetch_km = c(1, 2), pa_pred = c( 0L, 0L -), cover = c(77.1123333333333, 44.6653333333333), pa_post_hoc = c( +), cover_pred = c(77.1123333333333, 44.6653333333333), pa_post_hoc = c( 0L, 0L ), cover_post_hoc = c(0, 0)), row.names = c(NA, -2L), class = "data.frame") -res1_cover <- structure(list(depth_m = c(5, 10), cover = c( +res1_pred <- structure(list(depth_m = c(5, 10), cover_pred = c( 85.4316666666667, 21.144 ), cover_post_hoc = c(85.4316666666667, 21.144)), row.names = c( @@ -43,7 +43,7 @@ res1_cover <- structure(list(depth_m = c(5, 10), cover = c( ), class = "data.frame") -res1_pa <- structure(list(depth_m = c(5, 10), pa = 1:0, pa_post_hoc = 1:0), row.names = c( +res1_pa <- structure(list(depth_m = c(5, 10), pa_pred = 1:0, pa_post_hoc = 1:0), row.names = c( NA, -2L ), class = "data.frame") @@ -57,7 +57,7 @@ test_that("sav_model() works", { expect_equal(sav_model(df_ok_2), res2) expect_equal(sav_model(df_ok_3), res3) # - expect_equal(sav_model(df_ok_1, type = "cover"), res1_cover) + expect_equal(sav_model(df_ok_1, type = "cover"), res1_pred) expect_equal(sav_model(df_ok_1, type = "pa"), res1_pa) # } @@ -99,7 +99,7 @@ res_ph1a <- structure(list(depth_m = c(2, 2, 5), fetch_km = c(1, 1, 1), substrat ), secchi = c(20, 1, 20), limitation_secchi = c( TRUE, FALSE, TRUE -), vmax = c(28.242038767358, 1.7689, 28.242038767358), pa = c(1L, 1L, 0L), cover = c( +), vmax = c(28.242038767358, 1.7689, 28.242038767358), pa_pred = c(1L, 1L, 0L), cover_pred = c( 62.7603333333334, 62.7603333333334, 77.1123333333333 ), pa_post_hoc = c(1L, 0L, 0L), cover_post_hoc = c( diff --git a/vignettes/get_started.Rmd b/vignettes/get_started.Rmd index eff34f1..cfde169 100644 --- a/vignettes/get_started.Rmd +++ b/vignettes/get_started.Rmd @@ -20,29 +20,107 @@ library(sf) library(stars) ``` + + # Introduction -This vignette provides an overview of the workflow offered by the **Submerged Aquatic Vegetation (SAV) Model R package**. The package allows users to import spatial and tabular data related to SAV presence and habitat conditions in aquatic ecosystems. The presented in this vignette provide an overview of the functionalities built into the package and presents a workflow from user input to model predictions. +This vignette provides an overview of the workflow offered by the **Submerged Aquatic Vegetation Model (SAVM) R package**. The package allows users to import spatial and tabular data related to SAV presence and habitat conditions in aquatic ecosystems. The presented in this vignette provide an overview of the functionalities built into the package and presents a workflow from user input to model predictions. # User workflow -## Reading files +## R as a GIS -As an example, we consider a zone near Buffalo in the Lake Erie. Note that all -the data required data are included on in the package SAVM. +To use SAVM, basic knowledge about using R to work with spatial objects. +There are several packages that turn R into a powerful Geospatial Information System (GIS). +One of the most popular is [`sf`](https://r-spatial.github.io/sf/), +which can read a wide variety of spatial file formats, including ESRI shapefiles. +Assuming you have the following shapefile files: -```{r} -# Lake Erie boundaries polygon +``` +lake_erie_land +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.cpg +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.dbf +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.prj +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbn +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.sbx +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp +β”œβ”€β”€ LkErie_Land_fromGLAF_Water_WGS_Feb2020.shp.xml +└── LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx +``` + +then `st_read(lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx)` is used to read the shapefile into a `sf` object used to do spatial operations with R. For instance, the shapefile can be plot as follows: + + +```R +erie_land <- st_read("lake_erie_land/LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx") +plot(erie_land) +``` + +```{r echo = FALSE} +erie_land <- system.file( + "example", "lake_erie_land", "LkErie_Land_fromGLAF_Water_WGS_Feb2020.shx", + package = "SAVM", mustWork = TRUE +) |> + st_read() + +erie_land |> + st_geometry() |> # to plot only the geometry + plot(col = 1) +``` + +The package `SAVM` includes a util function to invert a polygon. + +```{r invert} +erie_lake <- erie_land |> + invert_polygon() + +erie_lake |> + st_geometry() |> + plot(col = 1) +``` + +Using `sf` capacities, the file can then be exported in various format, for instance in [GeoPackage](https://www.geopackage.org/): + +```R +erie_lake |> sf::st_write("lake_erie.gpkg") +``` + +Note that `sf` focuses on vector files, to read and manipulate raster file +we will use the package [`stars`](https://r-spatial.github.io/stars/) that works seamlessly with `sf`. + + + +## Reading file + +Now we consider a zone near Buffalo in the Lake Erie. Note that all +the data required are included on in the package SAVM. We first load the geopackage of Lake Erie. + +```{r reading1} +# Lake Erie boundaries polygon (reading from SAVM internal files) le_bound <- system.file("example", "lake_erie.gpkg", package = "SAVM") |> sf::st_read() |> sf::st_transform(crs = 3857) +``` + +We then load the shapefile of the area of interest that is stored as a GeoJSON file. Note that for the area (or points) of interest, we use `read_sav()` that returns an object of class `sav_data`. +```{r reading2} # Lake Erie study zone: read study_zone <- system.file("example", "study_zone.geojson", package = "SAVM") |> read_sav(spacing = 2000) +``` -# Depth +Note that the corresponfing grid can be visualized using the `preview_grid()` function. + +```{r preview_grid} +preview_grid(study_zone) +``` + +We now load the depth data. + +```{r reading3} +# Depth (reading with stars) study_depth <- stars::read_stars(system.file("example", "le_bathy.tiff", package = "SAVM")) ``` @@ -53,17 +131,26 @@ study_depth <- stars::read_stars(system.file("example", "le_bathy.tiff", package ### Compute fetch -To compute wind fetch (in kilometers, we called `compute_fetch()`), by default, the maximum distance is set to 15 km (see argument `max_dist`) and the default number of radial transect by quadrant is 9 (see argument `n_quad_seg`). +To compute wind fetch (in kilometers), we called `compute_fetch()`. By default, the maximum distance is set to 15 km (see argument `max_dist`) and the default number of bearings is 16 (see argument `n_bearings`). ```{r fetch} fetch <- compute_fetch(study_zone$points, le_bound) fetch ``` +`fetch` is a list of two elements, `mean_fetch`, an `sf` object with all +summarized computations and `transect_lines`, a `sf` object that includes all +transect lines. As `sf` objects, both elements can easily be exported, for +instance, the following code writes `fetch$mean_fetch` as a shapefile: + +```{r fetch_export, eval = FALSE} +fetch$mean_fetch |> sf::st_write("fetch.shp") +``` + ### Extract depth -To extract depth data from an raster file for our study points, we use `st_extract()` from the package `stars`. Note that the points and the raster must use the same projection. +To extract depth data from a raster file for our study points, we use `st_extract()` from the package `stars`. Note that the points and the raster must use the same projection. ```{r extract_depth} depth <- st_extract( @@ -80,43 +167,47 @@ This return an sf object ad the `le_bathy.tiff` contain the depth values. ## Model -With the `depth` and `fetch`, we create an data frame that contains the predictor value for all points. +With the `depth` and `fetch`, we create a `sf` object that contains predictor values for all points. -```{r } -dat <- data.frame( - depth_m = depth$le_bathy.tiff, - fetch_km = fetch$mean_fetch$fetch_km -) +```{r create_df} +sp_dat <- fetch$mean_fetch +sp_dat$depth_m <- depth$le_bathy.tiff ``` -You can generate predictions using `sav_model(dat)`. By default, the function uses the column names in the input data to identify the relevant variables. If the column names are not recognized, you can manually specify them using the `depth` and `fetch` parameters. +You can generate predictions using `sav_model(sp_dat)`. By default, the function uses the column names in the input data to identify the relevant variables. If the column names are not recognized, you can manually specify them using the `depth` and `fetch` parameters. -```{r} +```{r model1} # get input -res <- sav_model(dat) +res <- sav_model(sp_dat) res ``` -Even without additional columns, the model will return columns containing post hoc results. For SAV cover, if presence-absence predictions are available, the cover_post_hoc column will contain cover values set to 0 where absence is predicted. +Even without additional columns, the model will return columns containing post hoc results. For SAV cover, if presence-absence predictions are available, the `cover_post_hoc` column will contain `cover_pred` values set to 0 where absence is predicted. There are three additionnal limitaions that can be used as posthoc treatement. The `secchi` parameter should contain Secchi depth values and is used to compute Vmax following the equation from Chambers and Kalff (1985). If `depth` is also provided, it is compared to Vmax: when the depth exceeds Vmax, both cover and presence-absence are set to 0. The `substrate` and `limitation` parameters are binary variables that can be used to set cover and presence-absence to 0 based on external knowledge or site-specific constraints. Below we used an example with random Secchi depth values: -```{r} +```{r model2} res_secchi <- sav_model( - cbind(dat, secchi = runif(nrow(dat), 4, 12)) + cbind(sp_dat, secchi = runif(nrow(sp_dat), 4, 12)) ) res_secchi ``` +Note that `res_secchi` is a `sf` object and can therefore be exported as a shapefile: + +```{r final_export, eval = FALSE} +res_secchi |> sf::st_write("results.shp") +``` + ## Vizualize The package includes plotting helper functions. `plot_sav_distribution(res)` allows users to visualize the distribution of SAV cover and presence/absence across predictor bins. Alternatively, `plot_sav_density(res)` provides a density-based visualization rather than using binned values. -```{r} +```{r plot1} # Visualize results plot_sav_distribution(res_secchi) plot_sav_density(res_secchi) @@ -124,7 +215,7 @@ plot_sav_density(res_secchi) Note that these functions include an option to force the use of raw predictions without applying the post hoc treatment. -```{r} +```{r, plot2} plot_sav_distribution(res_secchi, post_hoc = FALSE) plot_sav_density(res_secchi, post_hoc = FALSE) ```