swath profile added

.Rbuildignore

@@ -9,4 +9,3 @@
 ^\.github$
 ^doc$
 ^Meta$
-^R/swath\.R$

DESCRIPTION

@@ -4,7 +4,8 @@ Version: 0.0.0.9000
 Authors@R: 
     person("Tobias", "Stephan", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-9290-014X"))
-Description: Measures distances along lines defined by a direction.
+Description: Toolset to measures distances along lines defined by a direction or
+    create swath profiles.
 License: MIT + file LICENSE
 URL: https://tobiste.github.io/geoprofiler/
 BugReports: https://github.com/tobiste/geoprofiler/issues
@@ -15,6 +16,7 @@ Imports:
     sf,
     structr,
     tectonicr,
+    terra,
     units
 Suggests: 
     knitr,

NAMESPACE

@@ -1,15 +1,26 @@
 # Generated by roxygen2: do not edit by hand
 
+export(line_ends)
 export(point_distance)
 export(profile_azimuth)
 export(profile_coords)
 export(profile_length)
 export(profile_line)
 export(profile_points)
+export(swath_extract)
+export(swath_profile)
 importFrom(dplyr,as_tibble)
+importFrom(dplyr,c_across)
+importFrom(dplyr,everything)
+importFrom(dplyr,matches)
 importFrom(dplyr,mutate)
 importFrom(dplyr,mutate_all)
+importFrom(dplyr,rename)
+importFrom(dplyr,rowwise)
+importFrom(dplyr,select)
+importFrom(dplyr,starts_with)
 importFrom(dplyr,tibble)
+importFrom(dplyr,ungroup)
 importFrom(sf,st_as_sf)
 importFrom(sf,st_cast)
 importFrom(sf,st_combine)
@@ -18,11 +29,19 @@ importFrom(sf,st_crs)
 importFrom(sf,st_geometry_type)
 importFrom(sf,st_is_longlat)
 importFrom(sf,st_length)
+importFrom(sf,st_linestring)
+importFrom(sf,st_point)
+importFrom(sf,st_set_crs)
 importFrom(sf,st_transform)
 importFrom(structr,vrotate)
 importFrom(tectonicr,deg2rad)
 importFrom(tectonicr,dist_greatcircle)
 importFrom(tectonicr,get_azimuth)
 importFrom(tectonicr,rad2deg)
+importFrom(terra,as.lines)
+importFrom(terra,extract)
+importFrom(terra,vect)
+importFrom(terra,xmin)
+importFrom(terra,ymin)
 importFrom(units,drop_units)
 importFrom(units,set_units)

R/profile_points.R

@@ -171,3 +171,28 @@ point_distance <- function(a, b, ...) {
 #'   distances.km <- scales::rescale(x, to = c(0, distance.km.total))
 #'   units::set_units(distances.km, "km")
 #' }
+
+
+#' Extract end points of a line
+#'
+#' @param x `sf` line object
+#'
+#' @returns `sf` point object
+#' @export
+#' @importFrom sf st_cast
+#'
+#' @examples
+#' p1 <- data.frame(lon = -90.8, lat = 48.6) |>
+#'   sf::st_as_sf(coords = c("lon", "lat"), crs = "WGS84")
+#' profile_points(p1,
+#'   profile.azimuth = 135, profile.length = 10000,
+#'   crs = sf::st_crs("EPSG:26915")
+#' ) |>
+#'   profile_line() |>
+#'   line_ends()
+line_ends <- function(x){
+  x_pts <- sf::st_cast(x, "POINT")
+  start <- x_pts[1]
+  end <- x_pts[length(x_pts)]
+  c(start, end)
+}

R/swath.R

@@ -7,39 +7,70 @@
 #' extracted from a given raster file, see argument \code{raster}. CRS of raster
 #' and points have to be the same.
 #'
-#' @param coords matrix(ncol=2, nrow=2) with x and y coordinates of beginning and
-#' end point of the baseline; each point in one row
+#' @param coords either a `sf` object or a matrix(ncol=2, nrow=2) with x and
+#' y coordinates of beginning and end point of the baseline; each point in one row
 #' \describe{
-#'   \item{column 1}{xcoordinates}
-#'   \item{column 2}{ycoordinates}
+#'   \item{column 1}{x coordinates}
+#'   \item{column 2}{y coordinates}
 #' }
-#' @param raster raster file (loaded with [raster::raster()])
-#' @param k integer; number of lines on each side of the baseline
-#' @param dist numeric; distance between lines
-#' @param crs string; CRS
-#' @param method string; method for extraction of raw data, see
-#' [raster::extract()]: default value: "bilinear"
+#' @param raster Raster file (loaded with [terra::rast()])
+#' @param k integer. number of lines on each side of the baseline
+#' @param dist numeric. distance between lines
+#' @param crs character. coordinate reference system.
+#' @param method character. method for extraction of raw data, see
+#' [terra::extract()]: default value: "bilinear"
 #'
-#' @importFrom sp SpatialPoints CRS
-#' @importFrom raster extract spLines
+#' @returns list.
+#' \describe{
+#'  \item{`swath`}{matrix. Statistics of the raster measured along the lines}
+#'  \item{`data`}{list of numeric vector containing the data extracted from the raster along each line}
+#'  \item{`lines`}{list of of the lines as SpatVectors}
+#'  }
+#'
+#' @importFrom sf st_crs st_transform st_point st_set_crs st_linestring
+#' @importFrom terra extract vect ymin xmin as.lines
 #'
-#' @author V. Haburaj
+#' @source The algorithm is a modified version of "swathR"
+#' by Vincent Haburaj (https://github.com/jjvhab/swathR).
 #'
 #' @export
-swathR <- function(coords, raster, k, dist, crs, method = 'bilinear') {
-  #message("Initializing ...")
+#'
+#' @seealso [swath_profile()]
+#'
+#' @examples
+#' # Create a random raster
+#' r <- terra::rast(ncol = 10, nrow = 10, xmin = -150, xmax = -80, ymin = 20, ymax = 60)
+#' values(r) <- runif(terra::ncell(r))
+#'
+#' # Create a random profile
+#' profile <- data.frame(lon = c(-140, -90), lat = c(55, 25)) |>
+#'   sf::st_as_sf(coords = c("lon", "lat"), crs = "WGS84")
+#' swath_extract(profile, r, k = 2, dist = 1)
+swath_extract <- function(coords, raster, k = 1, dist, crs = "EPSG:4326", method = c("bilinear", "simple")) {
+  # message("Initializing ...")
+  method <- match.arg(method)
+
+  if (inherits(coords, "sf") & all(sf::st_geometry_type(coords) == "LINESTRING")) {
+    coords <- line_ends(coords)
+  }
 
   # create SpatialPoints from coords:
-  spt <- SpatialPoints(coords, proj4string = CRS(crs))
+  # spt <- SpatialPoints(coords, proj4string = CRS(crs))
+  if (!inherits(coords, "sf") & is.matrix(coords)) {
+    coords <- sf::st_point(coords) |> sf::st_set_crs(crs)
+  }
+  coords_mat <- sf::st_coordinates(coords)
+
+  spt <- sf::st_transform(coords, crs = sf::st_crs(crs)) |> terra::vect()
 
   # get slope of baseline:
-  m <- (ymin(spt[1]) - ymin(spt[2])) / (xmin(spt[1]) - xmin(spt[2]))
+  m <- (terra::ymin(spt[1]) - terra::ymin(spt[2])) / (terra::xmin(spt[1]) - terra::xmin(spt[2]))
   # get slope of normal function:
   m1 <- -(1 / m)
   # get slope-angle from slope:
   alpha <- atan(m)
   alpha1 <- atan(m1)
-  # get deltax and deltay from pythagoras:
+  # get deltax and deltay from Pythagoras:
   if ((alpha * 180) / pi < 90 & (alpha * 180) / pi > 270) {
     deltax <- cos(alpha1) * dist
   } else {
@@ -59,39 +90,45 @@ swathR <- function(coords, raster, k, dist, crs, method = 'bilinear') {
   # list for spatial lines:
   allLines <- list()
   # add baseline:
-  allLines[[k + 1]] <- spLines(coords, crs = crs)
+  allLines[[k + 1]] <- terra::vect(coords) |> terra::as.lines()
   # set distance for baseline:
   swath[k + 1, 1] <- 0
   # generate k lines parallel to baseline:
   for (n in 1:k) {
     # BELOW BASELINE:
     # new points
     cn <- matrix(nrow = 2, ncol = 2)
-    cn[1, ] <- cbind(coords[1, 1] - (deltax * n), coords[1, 2] - (deltay * n))
-    cn[2, ] <- cbind(coords[2, 1] - (deltax * n), coords[2, 2] - (deltay * n))
+    cn[1, ] <- cbind(coords_mat[1, 1] - (deltax * n), coords_mat[1, 2] - (deltay * n))
+    cn[2, ] <- cbind(coords_mat[2, 1] - (deltax * n), coords_mat[2, 2] - (deltay * n))
     # line between points:
-    allLines[[k + 1 - n]] <- spLines(cn, crs = crs)
+    allLines[[k + 1 - n]] <- terra::vect(cn, crs = crs) |> terra::as.lines()
+
     # distance value:
     swath[k + 1 - n, 1] <- -1 * n * dist
     # ABOVE BASELINE:
     # new points
     cn <- matrix(nrow = 2, ncol = 2)
-    cn[1, ] <- cbind(coords[1, 1] + (deltax * n), coords[1, 2] + (deltay * n))
-    cn[2, ] <- cbind(coords[2, 1] + (deltax * n), coords[2, 2] + (deltay * n))
+    cn[1, ] <- cbind(coords_mat[1, 1] + (deltax * n), coords_mat[1, 2] + (deltay * n))
+    cn[2, ] <- cbind(coords_mat[2, 1] + (deltax * n), coords_mat[2, 2] + (deltay * n))
     # line between points:
-    allLines[[k + n + 1]] <- spLines(cn, crs = crs)
+    allLines[[k + n + 1]] <- terra::vect(cn, crs = crs) |> terra::as.lines()
     # distance value:
     swath[k + n + 1, 1] <- n * dist
   }
-  gc(verbose = FALSE)
+
+  # Expand raster to make sure all lines fall into the extent of the raster
+  lines_extent <- terra::vect(allLines) |> terra::ext()
+  raster_expanded <- terra::extend(raster, lines_extent)
+
+  # gc(verbose = FALSE)
   # get raw data:
-  message("Extracting raw data (this may take some time) ...")
+  # message("Extracting raw data (this may take some time) ...")
   raw.data <- sapply(allLines, FUN = function(x) {
-    extract(raster, x, method = method)
+    terra::extract(raster_expanded, x, method = method, ID = FALSE)
   })
-  gc(verbose = FALSE)
+  # gc(verbose = FALSE)
   # generalise data:
-  message("Generalising data ...")
+  # message("Generalising data ...")
   swath[, 2] <- sapply(raw.data, function(x) {
     mean(x, na.rm = T)
   })
@@ -115,45 +152,65 @@ swathR <- function(coords, raster, k, dist, crs, method = 'bilinear') {
   })
   # return results:
   results <- list(swath = swath, data = raw.data, lines = allLines)
-  message("Operation finished successfully!")
-  message('Structure of results (list): "swath": swath profile data (matrix, numeric), "data": raw data (list, numeric), "lines": generated lines (list, spLines)')
-  gc(verbose = FALSE)
+  # message("Operation finished successfully!")
+  # message('Structure of results (list): "swath": swath profile data (matrix, numeric), "data": raw data (list, numeric), "lines": generated lines (list, spLines)')
+  # gc(verbose = FALSE)
   return(results)
 }
 
 
 
 #' Elevation profile
 #'
-#' Extracts the minimum and maximum elevation data along a swathR profile.
+#' Statistics of the elevation data across a swath profile.
 #'
-#' @param x list. an return object from \code{swathR}
+#' @param x list. The return object of [swath_extract()]
+#' @param profile.length numeric or `units` object. If `NULL` the fractional
+#' distance is returned, i.e. 0 at start and 1 at the end of the profile.
 #' @return tibble
-#' @importFrom dplyr c_across rowwise ungroup mutate tibble as_tibble
+#' @importFrom dplyr c_across rowwise ungroup mutate tibble as_tibble everything matches select starts_with rename
 #' @export
-swath_profile <- function(x) {
-  elevs <- numeric()
-  center <- as.character(median(seq_along(x$data)))
-  for (i in seq_along(x$data)) {
-    elevs <- cbind(elevs, x$data[[i]])
-  }
-  elevs.df <- as_tibble(elevs)
-  names(elevs.df) <- seq_along(x$data)
+#' @examples
+#' # Create a radnom raster
+#' r <- terra::rast(ncol = 10, nrow = 10, xmin = -150, xmax = -80, ymin = 20, ymax = 60)
+#' values(r) <- runif(terra::ncell(r))
+#'
+#' # Create a random profile
+#' profile <- data.frame(lon = c(-140, -90), lat = c(55, 25)) |>
+#'   sf::st_as_sf(coords = c("lon", "lat"), crs = "WGS84")
+#' swath <- swath_extract(profile, r, k = 5, dist = 10)
+#'
+#' swath_profile(swath, profile.length = profile_length(profile_line(profile)))
+swath_profile <- function(x, profile.length = NULL) {
+  if (is.null(profile.length)) profile.length <- 1
+  center <- paste0("X", as.character(median(seq_along(x$data))))
+
+  max_length <- max(sapply(x$data, length))
 
-  elevs.df <- elevs.df |>
+  padded_vectors <- lapply(x$data, function(x) {
+    length(x) <- max_length
+    return(x)
+  })
+
+  # Combine the padded vectors into a matrix
+  elevs <- do.call(cbind, padded_vectors)
+  colnames(elevs) <- as.character(seq_along(x$data))
+
+  data.frame(elevs) |>
     rowwise() |>
     mutate(
-      min = min(c_across()),
-      max = max(c_across())
+      min = min(c_across(everything()), na.rm = TRUE),
+      mean = mean(c_across(everything()), na.rm = TRUE),
+      sd = sd(c_across(everything()), na.rm = TRUE),
+      quantile25 = quantile(c_across(everything()), na.rm = TRUE)[2],
+      median = median(c_across(everything()), na.rm = TRUE),
+      quantile75 = quantile(c_across(everything()), na.rm = TRUE)[4],
+      max = max(c_across(everything()), na.rm = TRUE)
+    ) |>
+    ungroup() |>
+    mutate(
+      "distance" = ((seq_along(elevs[, 1]) - 1) / (length(elevs[, 1]) - 1)) * profile.length,
     ) |>
-    ungroup()
-  elevs.df <-
-    tibble(
-      seq_along(elevs.df$min),
-      elevs.df[center],
-      elevs.df$min,
-      elevs.df$max
-    )
-  names(elevs.df) <- c("distance", "elevation", "min", "max")
-  return(elevs.df)
+    rename("elevation" = matches(center)) |>
+    select(distance, everything(), -starts_with("X"))
 }

README.Rmd

@@ -21,18 +21,18 @@ knitr::opts_chunk$set(
 
 <!-- badges: end -->
 
-The goal of geoprofiler is to get distances along and across user-defined 
+The goal of `{geoprofiler}` is to get distances along and across user-defined 
 profile lines or transects. This is useful when variables depend on distances.
 
 ![](man/figures/fig.png)
 
-The concept of geoprofiler is a coordinate transformation of your 
+The concept of `{geoprofiler}` is a coordinate transformation of your 
 geo-coordinates into "profile coordinates". These coordinates are the 
 distances along and across your profile.
 
 ## Installation
 
-You can install the development version of geoprofiler from [GitHub](https://github.com/) with:
+You can install the development version of `{geoprofiler}` from [GitHub](https://github.com/) with:
 
 ``` r
 # install.packages("devtools")

README.md

@@ -9,19 +9,19 @@
 
 <!-- badges: end -->
 
-The goal of geoprofiler is to get distances along and across
+The goal of `{geoprofiler}` is to get distances along and across
 user-defined profile lines or transects. This is useful when variables
 depend on distances.
 
 ![](man/figures/fig.png)
 
-The concept of geoprofiler is a coordinate transformation of your
+The concept of `{geoprofiler}` is a coordinate transformation of your
 geo-coordinates into “profile coordinates”. These coordinates are the
 distances along and across your profile.
 
 ## Installation
 
-You can install the development version of geoprofiler from
+You can install the development version of `{geoprofiler}` from
 [GitHub](https://github.com/) with:
 
 ``` r