* Switched raster/vector operations to terra

.lintr

@@ -1,5 +1,4 @@
-linters: with_defaults(
-    line_length_linter(100),
-    assignment_linter = NULL,
-    object_name_linter = NULL
-    )
+linters: linters_with_defaults(
+	line_length_linter = line_length_linter(120),
+	object_name_linter = NULL
+	)

DESCRIPTION

@@ -32,17 +32,17 @@ Imports:
   jsonlite,
   lazyeval,
   RColorBrewer,
-  sp,
+  magrittr,
   stats,
-  stars
+  terra,
+  tidyterra
 Suggests:
   lattice,
   leaflet,
   leafsync,
   lidR,
   plot3D,
   rasterVis,
-  sf,
   viridis
 RoxygenNote: 7.2.3
 Roxygen: list(markdown = TRUE)

NAMESPACE

@@ -16,6 +16,9 @@ export(clipLevel2BPAVDProfile)
 export(clipLevel2BPAVDProfileGeometry)
 export(clipLevel2BVPM)
 export(clipLevel2BVPMGeometry)
+export(clipLevel3)
+export(clipLevel3Geometry)
+export(crop)
 export(gedi.fullwaveform)
 export(gedi.level1b)
 export(gedi.level2a)
@@ -39,24 +42,26 @@ export(polyStatsLevel2BVPM)
 export(readLevel1B)
 export(readLevel2A)
 export(readLevel2B)
+export(readLevel3)
 exportClasses(gedi.fullwaveform)
 exportClasses(gedi.level1b)
 exportClasses(gedi.level2a)
 exportClasses(gedi.level2b)
 exportClasses(gedi.level3)
 exportMethods(close)
+exportMethods(crop)
+exportMethods(geom_spatraster)
 exportMethods(plot)
 import(curl)
 import(data.table)
-import(fs)
 import(ggplot2)
+import(graphics)
 import(hdf5r)
 import(jsonlite)
 import(lazyeval)
 import(methods)
-import(sp)
 import(stats)
-import(utils)
+importClassesFrom(terra,SpatRaster)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,element_rect)
@@ -67,13 +72,6 @@ importFrom(ggplot2,scale_fill_gradientn)
 importFrom(ggplot2,theme)
 importFrom(ggplot2,xlab)
 importFrom(ggplot2,ylab)
-importFrom(graphics,abline)
-importFrom(graphics,arrows)
-importFrom(graphics,axis)
-importFrom(graphics,mtext)
-importFrom(graphics,par)
-importFrom(graphics,text)
 importFrom(hdf5r,H5File)
-importFrom(stats,na.omit)
-importFrom(stats,quantile)
-importFrom(stats,setNames)
+importFrom(magrittr,`%>%`)
+importFrom(tidyterra,geom_spatraster)

R/class.gedi.R

@@ -1,7 +1,5 @@
-# setClass("gedi.level1b", representation(h5="H5File",level1b.spdf='SpatialPointsDataFrame'))
-#' @importFrom hdf5r H5File
-setRefClass("H5File")
-requireNamespace("data.table")
+#' @import methods
+methods::setRefClass("H5File")
 
 #' Class for GEDI level1B
 #'
@@ -11,9 +9,8 @@ requireNamespace("data.table")
 #' @seealso [`H5File`][hdf5r::H5File-class] in the `hdf5r` package and
 #' \url{https://lpdaac.usgs.gov/products/gedi01_bv002/}
 #'
-#' @import methods
 #' @export
-gedi.level1b <- setClass(
+gedi.level1b <- methods::setClass(
   Class = "gedi.level1b",
   slots = list(h5 = "H5File")
 )
@@ -26,9 +23,8 @@ gedi.level1b <- setClass(
 #' @seealso [`H5File`][hdf5r::H5File-class] in the `hdf5r` package and
 #' \url{https://lpdaac.usgs.gov/products/gedi02_av002/}
 #'
-#' @import methods
 #' @export
-gedi.level2a <- setClass(
+gedi.level2a <- methods::setClass(
   Class = "gedi.level2a",
   slots = list(h5 = "H5File")
 )
@@ -42,9 +38,8 @@ gedi.level2a <- setClass(
 #' @seealso [`H5File`][hdf5r::H5File-class] in the `hdf5r` package and
 #' \url{https://lpdaac.usgs.gov/products/gedi02_bv002/}
 #'
-#' @import methods
 #' @export
-gedi.level2b <- setClass(
+gedi.level2b <- methods::setClass(
   Class = "gedi.level2b",
   slots = list(h5 = "H5File")
 )
@@ -55,39 +50,38 @@ gedi.level2b <- setClass(
 #' @slot dt Object of class data.table from \emph{data.table} package containing
 #' the extracted GEDI full-waveform elevation and amplitude.
 #'
-#' @import methods
 #' @export
-gedi.fullwaveform <- setClass(
+gedi.fullwaveform <- methods::setClass(
   Class = "gedi.fullwaveform",
   slots = list(dt = "data.table")
 )
 
 
 #' Class for GEDI level3
 #'
-#' @slot raster Object of class [`SpatRaster`][terra::SpatRaster] from `terra` package containing the
-#' GEDI level3 products: rasterized metrics
+#' @slot raster Object of class [`SpatRaster`][terra::SpatRaster]
+#' from `terra` package, read by `terra::rast` containing the GEDI level3
+#' products: rasterized metrics
 #'
 #' @seealso
 #' \url{https://daac.ornl.gov/GEDI/guides/GEDI_L3_LandSurface_Metrics_V2.html}
 #'
-#' @import methods
+#' @importClassesFrom terra SpatRaster
 #' @export
-gedi.level3 <- setClass(
+gedi.level3 <- methods::setClass(
   Class = "gedi.level3",
   slots = list(raster = "SpatRaster")
 )
 
-
 #' Plot GEDI* object
-#' 
+#'
 #' @description For [`gedi.fullwaveform-class`]: will plot the full waveform
-#' 
+#'
 #' @param x An object of class [`gedi.fullwaveform-class`] (output of [getLevel1BWF()] function)
 #' @param relative if TRUE, the Waveform Amplitude will be showed in percentage (%)
 #' @param polygon if TRUE, the polygon will be added to the plot
 #' @param ... will be passed to the main plot
-#' 
+#'
 #' @return No return value
 #'
 #' @examples
@@ -143,11 +137,18 @@ setMethod(
       }
 
       if (polygon == TRUE) {
-        xstart <- x[which(z == min(z, na.rm = T))]
-        xend <- x[which(z == max(z, na.rm = T))]
+        xstart <- x[which(z == min(z, na.rm = TRUE))]
+        xend <- x[which(z == max(z, na.rm = TRUE))]
 
         xl <- c(min(x), min(x), xstart, rev(x), xend, min(x))
-        yl <- c(max(z, na.rm = T), min(z, na.rm = T), min(z, na.rm = T), rev(z), max(z, na.rm = T), max(z, na.rm = T))
+        yl <- c(
+          max(z, na.rm = TRUE),
+          min(z, na.rm = TRUE),
+          min(z, na.rm = TRUE),
+          rev(z),
+          max(z, na.rm = TRUE),
+          max(z, na.rm = TRUE)
+        )
 
         suppressWarnings({
           plot(xl, yl, ...)
@@ -177,7 +178,7 @@ h5closeall <- function(con, ...) {
 #'
 #' @param con An object of class `gedi.level1b`
 #' @param ... Inherited from base
-#' 
+#'
 #' @export
 #' @rdname close
 #' @method close gedi.level1b
@@ -188,7 +189,7 @@ setMethod("close", signature = c("gedi.level1b"), h5closeall)
 #'
 #' @description
 #' Closing files will avoid locking HDF5 GEDI files.
-#' 
+#'
 #' @param con An object of class `gedi.level2a`
 #' @param ... Inherited from base
 #' @method close gedi.level2a

R/clipLevel1B.R

@@ -50,7 +50,7 @@
 #' close(level1b)
 #' close(level1b_clip)
 #' }
-#' @import hdf5r fs
+#' @import hdf5r
 #' @export
 clipLevel1B <- function(level1b, xmin, xmax, ymin, ymax, output = "") {
   output <- checkOutput(output)
@@ -81,7 +81,7 @@ clipLevel1B <- function(level1b, xmin, xmax, ymin, ymax, output = "") {
 #'
 #' @param level1b A [`gedi.level1b-class`] object (output of [readLevel1B()] function).
 #' An S4 object of class "gedi.level1b".
-#' @param polygon Polygon or Multipolygon. An object opened with `sf::st_read`,
+#' @param polygon SpatVect. An object opened with `terra::vect`,
 #' @param output Optional character path where to save the new
 #' [`hdf5r::H5File`][hdf5r::H5File-class]. The default stores a temporary file only.
 #' @param split_by Polygon id. If defined, GEDI data will be clipped by each polygon using
@@ -109,9 +109,9 @@ clipLevel1B <- function(level1b, xmin, xmax, ymin, ymax, output = "") {
 #' # Specifying the path to shapefile
 #' polygon_filepath <- system.file("extdata", "stands_cerrado.shp", package = "rGEDI")
 #'
-#' # Reading shapefile as sf object
-#' library(sf)
-#' polygon <- sf::st_read(polygon_filepath)
+#' # Reading shapefile as SpatVect object
+#' library(terra)
+#' polygon <- terra::vect(polygon_filepath)
 #'
 #' # Spepecifing output file and path
 #' output <- file.path(outdir, "GEDI01_B_2019108080338_O01964_T05337_02_003_01_clip")
@@ -134,7 +134,7 @@ clipLevel1BGeometry <- function(level1b, polygon, output = "", split_by = NULL)
 
   spData <- getSpatialData1B(level1b)
 
-  bbox <- sf::st_bbox(polygon)
+  bbox <- terra::ext(polygon)
   xmin <- bbox$xmin
   xmax <- bbox$xmax
   ymin <- bbox$ymin

R/clipLevel1BGeo.R

@@ -59,19 +59,18 @@
 #' @export
 clipLevel1BGeo <- function(level1BGeo, xmin, xmax, ymin, ymax) {
   # xmin ymin xmax ymax
-  mask <-
-    level1BGeo$longitude_bin0 >= xmin &
-      level1BGeo$longitude_bin0 <= xmax &
-      level1BGeo$latitude_bin0 >= ymin &
-      level1BGeo$latitude_bin0 <= ymax &
-      level1BGeo$longitude_lastbin >= xmin &
-      level1BGeo$longitude_lastbin <= xmax &
-      level1BGeo$latitude_lastbin >= ymin &
-      level1BGeo$latitude_lastbin <= ymax
+  mask <- level1BGeo$longitude_bin0 >= xmin &
+    level1BGeo$longitude_bin0 <= xmax &
+    level1BGeo$latitude_bin0 >= ymin &
+    level1BGeo$latitude_bin0 <= ymax &
+    level1BGeo$longitude_lastbin >= xmin &
+    level1BGeo$longitude_lastbin <= xmax &
+    level1BGeo$latitude_lastbin >= ymin &
+    level1BGeo$latitude_lastbin <= ymax
   mask[!stats::complete.cases(mask)] <- FALSE
   mask <- (seq_along(level1BGeo$longitude_bin0))[mask]
   newFile <- level1BGeo[mask, ]
-  # newFile<- new("gedi.level1b.dt", dt = x[mask,])
+
   if (nrow(newFile) == 0) {
     print("The polygon does not overlap the GEDI data")
   } else {
@@ -85,8 +84,8 @@ clipLevel1BGeo <- function(level1BGeo, xmin, xmax, ymin, ymax) {
 #' data within a given geometry
 #'
 #' @param level1BGeo A [`data.table::data.table-class`] resulting from [getLevel1BGeo()] function.
-#' @param polygon Polygon. An object of class [`sf::sf`],
-#' which can be loaded as an ESRI shapefile using [sf::st_read] function in the \emph{sf} package.
+#' @param polygon SpatVect. An object of class `SpatVect`,
+#' which can be loaded as an ESRI shapefile using [terra::vect] function in the \emph{terra} package.
 #' @param split_by Polygon id. If defined, GEDI data will be clipped by each polygon using the
 #' polygon id from table of attribute defined by the user.
 #'
@@ -115,9 +114,9 @@ clipLevel1BGeo <- function(level1BGeo, xmin, xmax, ymin, ymax) {
 #' # Specifying the path to shapefile
 #' polygon_filepath <- system.file("extdata", "stands_cerrado.shp", package = "rGEDI")
 #'
-#' # Reading shapefile as sf object
-#' library(sf)
-#' polygon <- sf::st_read(polygon_filepath)
+#' # Reading shapefile as SpatVect object
+#' library(terra)
+#' polygon <- terra::vect(polygon_filepath)
 #'
 #' # Clipping GEDI Full Waveform Geolocations by Geometry
 #' level1BGeo_clip <- clipLevel1BGeoGeometry(level1BGeo, polygon, split_by = "id")
@@ -143,7 +142,7 @@ clipLevel1BGeo <- function(level1BGeo, xmin, xmax, ymin, ymax) {
 #' close(level1b)
 #' @export
 clipLevel1BGeoGeometry <- function(level1BGeo, polygon, split_by = "id") {
-  exshp <- sf::st_bbox(polygon)
+  exshp <- terra::ext(polygon)
   level1BGeo <- clipLevel1BGeo(
     level1BGeo,
     xmin = exshp$xmin,
@@ -155,27 +154,24 @@ clipLevel1BGeoGeometry <- function(level1BGeo, polygon, split_by = "id") {
   if (nrow(level1BGeo) == 0) {
     print("The polygon does not overlap the GEDI data")
   } else {
-    points <- sf::st_as_sf(
+    points <- terra::vect(
       level1BGeo,
-      coords = c("longitude_bin0", "latitude_bin0"),
-      crs = sf::st_crs(polygon)
+      geom = c("longitude_bin0", "latitude_bin0"),
+      crs = terra::crs(polygon)
     )
 
-    names(points) <- gsub("^(?!geometry)", "x_\\1", names(points), perl = TRUE)
-    pts <- sf::st_intersection(sf::st_make_valid(points), sf::st_make_valid(polygon))
+    points$rowNumber <- as.integer(seq_along(points))
+    pts <- terra::intersect(terra::makeValid(points), terra::makeValid(polygon))
 
+    mask <- pts$rowNumber
+    newFile <- level1BGeo[mask, ]
     if (!is.null(split_by)) {
       if (any(names(polygon) == split_by)) {
-        mask <- as.integer(pts$x_id)
-        newFile <- level1BGeo[mask, ]
         newFile$poly_id <- pts[[split_by]]
       } else {
         stop(paste("The", split_by, "is not included in the attribute table.
                        Please check the names in the attribute table"))
       }
-    } else {
-      mask <- as.integer(rownames(pts))
-      newFile <- level1BGeo[mask, ]
     }
 
     return(newFile)

R/clipLevel2A.R

@@ -80,8 +80,8 @@ clipLevel2A <- function(level2a, xmin, xmax, ymin, ymax, output = "") {
 #'
 #' @param level2a A GEDI Level2A object (output of [readLevel2A()] function).
 #' An S4 object of class "gedi.level2a".
-#' @param polygon Polygon. An object of class [`sf::sf`],
-#' which can be loaded as an ESRI shapefile using [sf::st_read()] function in the \emph{sf} package.
+#' @param polygon Polygon. An object of class [`SpatVector`],
+#' which can be loaded as an ESRI shapefile using [terra::vect()] function in the \emph{terra} package.
 #' @param output optional character path where to save the new h5file. Default "" (temporary file).
 #' @param split_by Polygon id. If defined, GEDI data will be clipped by each polygon using the
 #' attribute specified by `split_by` from the attribute table.
@@ -109,9 +109,9 @@ clipLevel2A <- function(level2a, xmin, xmax, ymin, ymax, output = "") {
 #' # Specifying the path to shapefile
 #' polygon_filepath <- system.file("extdata", "stands_cerrado.shp", package = "rGEDI")
 #'
-#' # Reading shapefile as sf object
-#' library(sf)
-#' polygon <- sf::st_read(polygon_filepath)
+#' # Reading shapefile as SpatVect object
+#' library(terra)
+#' polygon <- terra::vect(polygon_filepath)
 #'
 #' # Specifying output file and path
 #' output <- file.path(outdir, "GEDI02_A_2019108080338_O01964_T05337_02_001_01_clip")
@@ -132,7 +132,7 @@ clipLevel2AGeometry <- function(level2a, polygon, output = "", split_by = NULL)
 
   spData <- getSpatialData2A(level2a)
 
-  bbox <- sf::st_bbox(polygon)
+  bbox <- terra::ext(polygon)
   xmin <- bbox$xmin
   xmax <- bbox$xmax
   ymin <- bbox$ymin