Merge pull request #114 from bfast2/111-terra-examples

Update from raster to terra package

.github/workflows/R-CMD-check.yml

@@ -1,14 +1,14 @@
+# 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
+    branches: [main, master]
   pull_request:
-    branches:
-      - main
-      - master
+    branches: [main, master]
 
-name: R-CMD-check
+name: R-CMD-check.yaml
+
+permissions: read-all
 
 jobs:
   R-CMD-check:
@@ -20,66 +20,33 @@ jobs:
       fail-fast: false
       matrix:
         config:
+          - {os: macos-latest,   r: 'release'}
           - {os: windows-latest, r: 'release'}
-          - {os: macOS-latest, r: 'release'}
-          - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
-          - {os: ubuntu-20.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+          - {os: ubuntu-latest,   r: 'devel', http-user-agent: 'release'}
+          - {os: ubuntu-latest,   r: 'release'}
+          - {os: ubuntu-latest,   r: 'oldrel-1'}
 
     env:
-      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
-      RSPM: ${{ matrix.config.rspm }}
       GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+      R_KEEP_PKG_SOURCE: yes
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
+
+      - uses: r-lib/actions/setup-pandoc@v2
 
-      - uses: r-lib/actions/setup-r@v1
+      - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
+          http-user-agent: ${{ matrix.config.http-user-agent }}
+          use-public-rspm: true
 
-      - uses: r-lib/actions/setup-pandoc@v1
-
-      - name: Query dependencies
-        run: |
-          install.packages('remotes')
-          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
-          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
-        shell: Rscript {0}
-
-      - name: Restore R package cache
-        if: runner.os != 'Windows'
-        uses: actions/cache@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
         with:
-          path: ${{ env.R_LIBS_USER }}
-          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
-          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
-
-      - name: Install system dependencies
-        if: runner.os == 'Linux'
-        run: |
-          while read -r cmd
-          do
-            eval sudo $cmd
-          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
+          extra-packages: any::rcmdcheck
+          needs: check
 
-      - name: Install dependencies
-        run: |
-          remotes::install_deps(dependencies = TRUE)
-          remotes::install_cran("rcmdcheck")
-        shell: Rscript {0}
-
-      - name: Check
-        env:
-          _R_CHECK_CRAN_INCOMING_REMOTE_: false
-        run: |
-          options(crayon.enabled = TRUE)
-          rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check")
-        shell: Rscript {0}
-
-      - name: Upload check results
-        if: failure()
-        uses: actions/upload-artifact@main
+      - uses: r-lib/actions/check-r-package@v2
         with:
-          name: ${{ runner.os }}-r${{ matrix.config.r }}-results
-          path: check
-
+          upload-snapshots: true
+          build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")'

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: bfast
-Version: 1.6.1.9999
+Version: 1.7.0
 Title: Breaks for Additive Season and Trend
 Authors@R: c(person(given = "Jan", family = "Verbesselt", role = c("aut"), email = "[email protected]"),
              person(given = "Dainius", family = "Masili\u016Bnas", role = c("aut", "cre"),
@@ -32,7 +32,7 @@ Description: Decomposition of time series into
     time series data and detect structural changes (breaks).
 Depends:
     R (>= 3.0.0),
-    strucchangeRcpp
+    strucchangeRcpp (>= 1.5-4)
 Imports:
     graphics,
     stats,
@@ -44,14 +44,14 @@ Suggests:
     MASS,
     sfsmisc,
     stlplus,
-    raster
+    terra
 License: GPL (>= 2)
 URL: https://bfast2.github.io/
 BugReports: https://github.com/bfast2/bfast/issues
 LazyLoad: yes
 LazyData: yes
 LinkingTo: Rcpp
-RoxygenNote: 7.1.1
+RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)
 RdMacros: Rdpack
 Encoding: UTF-8

NAMESPACE

@@ -8,6 +8,8 @@ S3method(breakpoints,bfast01)
 S3method(coef,bfast01)
 S3method(deviance,bfast01)
 S3method(fitted,bfast01)
+S3method(history_roc,formula)
+S3method(history_roc,matrix)
 S3method(logLik,bfast01)
 S3method(model.frame,bfast01)
 S3method(model.matrix,bfast01)

NEWS

@@ -12,6 +12,11 @@ Changes in Version 1.7.0
 
   o Added a reference to the paper describing the BFAST Lite algorithm.
 
+  o Updated the examples, migrating from the 'raster' to the 'terra' package.
+
+  o The 'history_roc()' generic is now consistent in parameters with the
+  'history_roc.matrix()' and 'history_roc.formula()' methods.
+
 Changes in Version 1.6.1
 
   o Fixed a heap overflow issue in the C++ version of bfastts().

R/bfast-package.R

@@ -73,7 +73,7 @@ NULL
 
 
 #' A vector with date information (a Datum type) to be linked with each NDVI
-#' layer within the modis raster brick (modisraster data set)
+#' layer within the modis raster datacube (modisraster data set)
 #' 
 #' \code{dates} is an object of class "Date" and contains the "Date"
 #' information to create a 16-day time series object.
@@ -113,10 +113,10 @@ NULL
 
 
 
-#' A raster brick of 16-day satellite image NDVI time series for a small subset
+#' A raster datacube of 16-day satellite image NDVI time series for a small subset
 #' in south eastern Somalia.
 #' 
-#' A raster brick containing 16-day NDVI satellite images (MOD13C1 product).
+#' A Cloud-Optimised GeoTIFF containing 16-day NDVI satellite images (MOD13C1 product).
 #' 
 #' 
 #' @name modisraster

R/bfast.R

@@ -108,7 +108,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"),
   level  <- rep(level, length.out = 2)
   ti <- time(Yt)
   f <- frequency(Yt) # on cycle every f time points (seasonal cycle)
-  if (class(Yt) != "ts")
+  if (!inherits(Yt, "ts"))
     stop("Not a time series object")
   if (f < 1)
     stop("Time series frequency needs to be more than 1")

R/bfastts.R

@@ -35,24 +35,6 @@
 #' @keywords ts
 #' 
 #' @example examples/bfastts.r
-#' @examples
-#'
-#'
-#' \dontrun{
-#' # Example of use with a raster
-#'
-#' library("raster")
-#' f <- system.file("extdata/modisraster.grd", package="bfast")
-#' modisbrick <- brick(f)
-#' ndvi <- bfastts(as.vector(modisbrick[1]), dates, type = c("16-day")) ## data of pixel 1
-#' plot(ndvi/10000) 
-#' 
-#' # Time series of 4 pixels
-#' modis_ts = t(as.data.frame(modisbrick))[,1:4]
-#' # Data with multiple columns, 2-4 are external regressors
-#' ndvi <- bfastts(modis_ts, dates, type = c("16-day"))
-#' plot(ndvi/10000)
-#' }
 #' 
 #' @export
 bfastts <- function(data, dates, type = c("irregular", "16-day", "10-day")) {
@@ -171,7 +153,7 @@ bfastts <- function(data, dates, type = c("irregular", "16-day", "10-day")) {
 # this is an alternative (faster) implementation that tries to work more precisely 
 # with leap years and arbitrary frequencies but needs more testing. 
 .bfast_construct_ts <- function(data, dates, freq=23) {
-  if (class(dates) != "Date") {
+  if (!inherits(dates, "Date")) {
     dates = as.Date(dates)
   }
 

R/history_roc.R

@@ -7,12 +7,12 @@ history_roc <- function (x, ...) {
   UseMethod("history_roc", x)
 }
 
-
-history_roc.matrix <- function(X, y,time, level = 0.05) {
-  n <- nrow(X)
+#' @export
+history_roc.matrix <- function(x, y,time, level = 0.05, ...) {
+  n <- nrow(x)
   #data_rev <- data[n:1,]
   #data_rev$response <- ts(data_rev$response)
-  X_rev  <- X[n:1,]
+  X_rev  <- x[n:1,]
   y_rev <-  y[n:1]
   time_rev <-  time[n:1]
 
@@ -32,11 +32,12 @@ history_roc.matrix <- function(X, y,time, level = 0.05) {
 ## A technique to verify whether or not the historical period is stable or not
 ## reversely order sample and perform
 ## recursive CUSUM test
-history_roc.formula <- function(formula, data, level = 0.05) {
+#' @export
+history_roc.formula <- function(x, data, level = 0.05, ...) {
   n <- nrow(data)
   data_rev <- data[n:1,]
   data_rev$response <- ts(data_rev$response)
-  y_rcus <- efp(formula, data = data_rev, type = "Rec-CUSUM")
+  y_rcus <- efp(x, data = data_rev, type = "Rec-CUSUM")
 
   pval = sctest(y_rcus)$p.value
   y_start <- if(!is.na(pval) && pval < level) {

examples/bfastlite.r

@@ -26,3 +26,33 @@ plot(bp)
 
 # Details of the structural change test with the type RE
 bfastlite(datats, level=0.05, type="RE")$sctest
+
+\donttest{
+## Run bfastlite() on a raster
+f <- system.file("extdata/modisraster.tif", package="bfast")
+modisbrick <- terra::rast(f)
+
+# Run on pixel 10
+data <- unlist(modisbrick[10])
+ndvi <- bfastts(data, dates, type = c("16-day"))
+bfl <- bfastlite(ndvi, breaks = "BIC")
+# Get max magnitude by RMSD
+max(magnitude(bfl[["breakpoints"]])$Mag[,"RMSD"])
+
+# Wrapper function
+bflSpatial <- function(pixels)
+{
+    ts <- bfastts(pixels, dates, type = c("16-day"))
+    bfl <- bfastlite(ts, breaks="BIC")
+    bp <- bfl[["breakpoints"]]
+    # Return number of breakpoints and max breakpoint magnitude
+    if (length(bp[["breakpoints"]]) == 1 && is.na(bp[["breakpoints"]]))
+        return(c(0, 0))
+
+    return(c(length(bp[["breakpoints"]]), max(magnitude(bp)$Mag[,"RMSD"])))
+}
+
+# Run function on each raster pixel
+rastbfl <- terra::app(modisbrick, bflSpatial)
+terra::plot(rastbfl)
+}

examples/bfastmonitor.r

@@ -38,19 +38,19 @@ mon <- bfastmonitor(NDVIb, formula = response ~ season,
 summary(mon$model)
 AIC(mon$model)
 
+\donttest{
 ## Example for processing raster bricks (satellite image time series of 16-day NDVI images)
-f <- system.file("extdata/modisraster.grd", package = "bfast")
-library("raster")
-modisbrick <- raster::brick(f)
-data <- as.vector(modisbrick[1])
+f <- system.file("extdata/modisraster.tif", package="bfast")
+modisbrick <- terra::rast(f)
+data <- unlist(modisbrick[1])
 ndvi <- bfastts(data, dates, type = c("16-day"))
 plot(ndvi/10000)
 
 ## derive median NDVI of a NDVI raster brick
-medianNDVI <- raster::calc(modisbrick, fun = function(x) median(x, na.rm = TRUE))
-raster::plot(medianNDVI)
+medianNDVI <- terra::app(modisbrick, fun = "median")
+terra::plot(medianNDVI)
 
-## helper function to be used with the calc() function
+## helper function to be used with the app() function
 xbfastmonitor <- function(x, timestamps = dates) {
 	ndvi <- bfastts(x, timestamps, type = c("16-day"))
 	ndvi <- window(ndvi, end = c(2011, 14))/10000
@@ -60,16 +60,14 @@ xbfastmonitor <- function(x, timestamps = dates) {
 }
 
 ## apply on one pixel for testing
-ndvi <- bfastts(as.numeric(modisbrick[1])/10000, dates, type = c("16-day"))
-plot(ndvi)
-
 bfm <- bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))
 bfm$magnitude
 plot(bfm)
-xbfastmonitor(modisbrick[1], dates) ## helper function applied on one pixel
+xbfastmonitor(data, dates) ## helper function applied on one pixel
 
 ## apply the bfastmonitor function onto a raster brick
-timeofbreak <- raster::calc(modisbrick, fun=xbfastmonitor)
+timeofbreak <- terra::app(modisbrick, fun=xbfastmonitor)
 
-raster::plot(timeofbreak) ## time of break and magnitude of change
-raster::plot(timeofbreak,2) ## magnitude of change
+terra::plot(timeofbreak) ## time of break and magnitude of change
+terra::plot(timeofbreak,2) ## magnitude of change
+}

examples/bfastts.r

@@ -4,3 +4,17 @@ bfastts(timedf$y, timedf$dates, type = "16-day")
 
 # Irregular
 head(bfastts(timedf$y, timedf$dates, type = "irregular"), 50)
+
+\donttest{
+# Example of use with a raster
+f <- system.file("extdata/modisraster.tif", package="bfast")
+modisbrick <- terra::rast(f)
+ndvi <- bfastts(unlist(modisbrick[1]), dates, type = c("16-day")) ## data of pixel 1
+plot(ndvi/10000) 
+
+# Time series of 4 pixels
+modis_ts = t(modisbrick[1:4])
+# Data with multiple columns, 2-4 are external regressors
+ndvi <- bfastts(modis_ts, dates, type = c("16-day"))
+plot(ndvi/10000)
+}