diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml
index c25daf3..d46a617 100644
--- a/.github/workflows/R-CMD-check.yml
+++ b/.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")'
diff --git a/DESCRIPTION b/DESCRIPTION
index c07f1d9..5e85826 100644
--- a/DESCRIPTION
+++ b/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 = "Jan.Verbesselt@wur.nl"),
              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
diff --git a/NAMESPACE b/NAMESPACE
index e29fed7..617bfc0 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/NEWS b/NEWS
index 6c9eadf..550bf16 100644
--- a/NEWS
+++ b/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().
diff --git a/R/bfast-package.R b/R/bfast-package.R
index b59bf42..1393dcf 100644
--- a/R/bfast-package.R
+++ b/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
diff --git a/R/bfast.R b/R/bfast.R
index 5bbe81e..0bd2d32 100644
--- a/R/bfast.R
+++ b/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")
diff --git a/R/bfastts.R b/R/bfastts.R
index a0ab465..956e85d 100644
--- a/R/bfastts.R
+++ b/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)
   }
   
diff --git a/R/history_roc.R b/R/history_roc.R
index 5a5e46a..2a7faa2 100644
--- a/R/history_roc.R
+++ b/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) {
diff --git a/examples/bfastlite.r b/examples/bfastlite.r
index 33a9373..75224c7 100644
--- a/examples/bfastlite.r
+++ b/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)
+}
diff --git a/examples/bfastmonitor.r b/examples/bfastmonitor.r
index 3cdb8bd..aa1b50a 100644
--- a/examples/bfastmonitor.r
+++ b/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
+}
diff --git a/examples/bfastts.r b/examples/bfastts.r
index 9b678bb..593131f 100644
--- a/examples/bfastts.r
+++ b/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)
+}
diff --git a/inst/extdata/modisraster.grd b/inst/extdata/modisraster.grd
deleted file mode 100644
index e5b48eb..0000000
--- a/inst/extdata/modisraster.grd
+++ /dev/null
@@ -1,25 +0,0 @@
-[general]
-creator=R package 'raster'
-created=2021-04-27 11:20:48
-[georeference]
-nrows=5
-ncols=5
-xmin=41.9
-ymin=-0.15
-xmax=42.15
-ymax=0.1
-projection=+proj=longlat +datum=NAD27 +no_defs
-[data]
-datatype=FLT4S
-byteorder=little
-nbands=275
-bandorder=BIL
-minvalue=4052:4190:3788:3326:3789:4481:6721:3790:5071:3950:3533:2709:3583:2749:3170:2652:3121:5993:6464:5077:4536:4281:3913:3849:3774:3267:5644:6368:5804:4695:4622:2371:3512:4350:2887:3790:2923:3175:2330:5940:7088:5118:6601:5287:4315:4015:3862:4042:2068:2847:2696:6961:5674:3152:3885:4164:3534:3941:3232:5154:5702:6348:7113:7546:6942:7329:5445:4792:4267:4188:3771:3794:3198:6271:6393:7328:6045:4718:5666:5768:4456:3786:3319:3297:2667:6771:6325:7796:6610:5640:4688:4129:3747:3753:3953:2837:4587:4449:4711:4742:4300:5344:5411:5091:3777:4016:3353:5913:6910:7596:6150:6377:5655:4490:3899:3909:3642:3631:4126:5535:6031:5017:4503:2032:3893:4266:3884:3897:3032:2737:4060:3636:6703:6675:4699:4217:4013:3836:3518:3625:3387:5077:5576:5990:3825:4285:4040:2180:4134:4040:4327:3689:3462:4675:6674:6043:7061:7489:7085:5622:5247:4425:4140:3325:3635:6540:6314:6222:6592:6191:5821:5469:4359:3563:3187:5946:5538:7272:4786:7013:5472:4830:4332:3664:3405:3136:3417:4691:6470:4518:5731:4215:4364:4646:1914:3736:3783:3305:2984:4076:3275:5915:5927:4327:4388:3533:3638:3177:3013:3104:3887:5387:5667:6574:4284:4126:4479:4114:4441:3732:3114:3492:3097:7140:6672:6198:5459:4838:4192:3469:3323:3255:3766:5476:4841:6309:5703:5011:3407:2879:4428:3585:2960:2801:3671:2583:1895:4554:3806:3507:3429:3224:3059:2992:2800:2957:2749:5189:5471:5080:4258:2558:2148:3573:2606:2171:2963:2296:5586:5626:4462:7314:6624:5827:5095
-maxvalue=5042:5376:4767:7207:8028:8042:8287:8051:6989:5550:5505:4480:4878:4468:5583:5139:6748:7405:8054:7187:6075:5563:4810:4757:4491:4692:7923:8597:7909:8565:6818:6909:5502:5356:4929:4733:4566:4342:4347:8091:8384:8576:8074:7004:5180:4885:4488:5317:3677:6433:7097:8327:7962:6536:6200:5638:5428:5060:5301:7506:8083:8315:8250:8335:8250:8223:7195:6094:5215:4725:4655:5786:6455:8345:8509:8588:8677:7661:7983:7109:6068:5934:4809:4106:4308:8632:8433:8538:7974:7370:5986:5523:4981:4528:4784:6861:7772:7450:6768:7057:8075:7835:7637:6459:5477:5596:5216:8163:8637:8808:8025:7580:6930:5533:5487:4624:4209:4216:6938:7408:7881:8007:6735:5857:7409:6110:5965:5193:4627:4037:6200:7626:7920:7945:6298:6176:5560:5086:4331:4604:4906:6906:7881:7654:6612:7515:6566:6856:5665:5361:6323:6165:5402:7547:8349:7974:8371:8509:8448:7399:6612:5696:5364:5169:6166:8532:8471:7984:8303:7706:7018:6350:6823:5553:6172:8255:8127:8001:8064:8213:7949:7424:6977:4910:4543:4207:6652:7321:8301:7674:7963:7843:7215:6231:5416:4914:6489:6773:6358:8081:7355:7562:7333:5956:6088:4829:4828:4288:4188:4653:5330:7988:7937:7723:7161:6863:7611:6849:5568:4936:4160:4538:6920:8055:7638:7935:6906:6513:5715:5201:4599:4720:6239:7173:7160:8067:7720:7428:6569:6038:6003:5397:4380:4296:4877:4792:7457:7753:6808:6056:5741:4813:4432:4499:4604:4112:4829:7657:8149:6817:6052:5335:4712:5037:4066:3813:4265:3665:7673:7907:6916:9020:7860:7084:6480
-nodatavalue=-3.4e+38
-[legend]
-legendtype=
-values=
-color=
-[description]
-layername=X2000.02.18:X2000.03.05:X2000.03.21:X2000.04.06:X2000.04.22:X2000.05.08:X2000.05.24:X2000.06.09:X2000.06.25:X2000.07.11:X2000.07.27:X2000.08.12:X2000.08.28:X2000.09.13:X2000.09.29:X2000.10.15:X2000.10.31:X2000.11.16:X2000.12.02:X2000.12.18:X2001.01.01:X2001.01.17:X2001.02.02:X2001.02.18:X2001.03.06:X2001.03.22:X2001.04.07:X2001.04.23:X2001.05.09:X2001.05.25:X2001.06.10:X2001.06.26:X2001.07.12:X2001.07.28:X2001.08.13:X2001.08.29:X2001.09.14:X2001.09.30:X2001.10.16:X2001.11.01:X2001.11.17:X2001.12.03:X2001.12.19:X2002.01.01:X2002.01.17:X2002.02.02:X2002.02.18:X2002.03.06:X2002.03.22:X2002.04.07:X2002.04.23:X2002.05.09:X2002.05.25:X2002.06.10:X2002.06.26:X2002.07.12:X2002.07.28:X2002.08.13:X2002.08.29:X2002.09.14:X2002.09.30:X2002.10.16:X2002.11.01:X2002.11.17:X2002.12.03:X2002.12.19:X2003.01.01:X2003.01.17:X2003.02.02:X2003.02.18:X2003.03.06:X2003.03.22:X2003.04.07:X2003.04.23:X2003.05.09:X2003.05.25:X2003.06.10:X2003.06.26:X2003.07.12:X2003.07.28:X2003.08.13:X2003.08.29:X2003.09.14:X2003.09.30:X2003.10.16:X2003.11.01:X2003.11.17:X2003.12.03:X2003.12.19:X2004.01.01:X2004.01.17:X2004.02.02:X2004.02.18:X2004.03.05:X2004.03.21:X2004.04.06:X2004.04.22:X2004.05.08:X2004.05.24:X2004.06.09:X2004.06.25:X2004.07.11:X2004.07.27:X2004.08.12:X2004.08.28:X2004.09.13:X2004.09.29:X2004.10.15:X2004.10.31:X2004.11.16:X2004.12.02:X2004.12.18:X2005.01.01:X2005.01.17:X2005.02.02:X2005.02.18:X2005.03.06:X2005.03.22:X2005.04.07:X2005.04.23:X2005.05.09:X2005.05.25:X2005.06.10:X2005.06.26:X2005.07.12:X2005.07.28:X2005.08.13:X2005.08.29:X2005.09.14:X2005.09.30:X2005.10.16:X2005.11.01:X2005.11.17:X2005.12.03:X2005.12.19:X2006.01.01:X2006.01.17:X2006.02.02:X2006.02.18:X2006.03.06:X2006.03.22:X2006.04.07:X2006.04.23:X2006.05.09:X2006.05.25:X2006.06.10:X2006.06.26:X2006.07.12:X2006.07.28:X2006.08.13:X2006.08.29:X2006.09.14:X2006.09.30:X2006.10.16:X2006.11.01:X2006.11.17:X2006.12.03:X2006.12.19:X2007.01.01:X2007.01.17:X2007.02.02:X2007.02.18:X2007.03.06:X2007.03.22:X2007.04.07:X2007.04.23:X2007.05.09:X2007.05.25:X2007.06.10:X2007.06.26:X2007.07.12:X2007.07.28:X2007.08.13:X2007.08.29:X2007.09.14:X2007.09.30:X2007.10.16:X2007.11.01:X2007.11.17:X2007.12.03:X2007.12.19:X2008.01.01:X2008.01.17:X2008.02.02:X2008.02.18:X2008.03.05:X2008.03.21:X2008.04.06:X2008.04.22:X2008.05.08:X2008.05.24:X2008.06.09:X2008.06.25:X2008.07.11:X2008.07.27:X2008.08.12:X2008.08.28:X2008.09.13:X2008.09.29:X2008.10.15:X2008.10.31:X2008.11.16:X2008.12.02:X2008.12.18:X2009.01.01:X2009.01.17:X2009.02.02:X2009.02.18:X2009.03.06:X2009.03.22:X2009.04.07:X2009.04.23:X2009.05.09:X2009.05.25:X2009.06.10:X2009.06.26:X2009.07.12:X2009.07.28:X2009.08.13:X2009.08.29:X2009.09.14:X2009.09.30:X2009.10.16:X2009.11.01:X2009.11.17:X2009.12.03:X2009.12.19:X2010.01.01:X2010.01.17:X2010.02.02:X2010.02.18:X2010.03.06:X2010.03.22:X2010.04.07:X2010.04.23:X2010.05.09:X2010.05.25:X2010.06.10:X2010.06.26:X2010.07.12:X2010.07.28:X2010.08.13:X2010.08.29:X2010.09.14:X2010.09.30:X2010.10.16:X2010.11.01:X2010.11.17:X2010.12.03:X2010.12.19:X2011.01.01:X2011.01.17:X2011.02.02:X2011.02.18:X2011.03.06:X2011.03.22:X2011.04.07:X2011.04.23:X2011.05.09:X2011.05.25:X2011.06.10:X2011.06.26:X2011.07.12:X2011.07.28:X2011.08.13:X2011.08.29:X2011.09.14:X2011.09.30:X2011.10.16:X2011.11.01:X2011.11.17:X2011.12.03:X2011.12.19:X2012.01.01:X2012.01.17
diff --git a/inst/extdata/modisraster.gri b/inst/extdata/modisraster.gri
deleted file mode 100644
index 9a6e953..0000000
Binary files a/inst/extdata/modisraster.gri and /dev/null differ
diff --git a/inst/extdata/modisraster.tif b/inst/extdata/modisraster.tif
new file mode 100644
index 0000000..4a37c3b
Binary files /dev/null and b/inst/extdata/modisraster.tif differ
diff --git a/man/bfast.Rd b/man/bfast.Rd
index 20be669..e662c94 100644
--- a/man/bfast.Rd
+++ b/man/bfast.Rd
@@ -39,8 +39,12 @@ time series is 1) "none" can be selected to avoid fitting a seasonal model.}
 \item{max.iter}{maximum amount of iterations allowed for estimation of
 breakpoints in seasonal and trend component.}
 
-\item{breaks}{integer specifying the maximal number of breaks to be
-calculated. By default the maximal number allowed by h is used.}
+\item{breaks}{either an integer specifying the number of breaks to compute,
+or a string specifying a statistic with which to compute
+the optimal number of breakpoints (see \code{\link[strucchangeRcpp:breakpoints]{strucchangeRcpp::breakpoints()}} for
+more information). If one value is given, it is used for both the trend and
+seasonal components, and if two are given, the first one is considered to be
+for the trend and the second for the seasonal component.}
 
 \item{hpc}{A character specifying the high performance computing support.
 Default is "none", can be set to "foreach". Install the "foreach" package
diff --git a/man/bfastlite.Rd b/man/bfastlite.Rd
index e306194..a2ee462 100644
--- a/man/bfastlite.Rd
+++ b/man/bfastlite.Rd
@@ -135,6 +135,36 @@ 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)
+}
 }
 \references{
 \insertRef{dainiusbfastlite}{bfast}
diff --git a/man/bfastmonitor.Rd b/man/bfastmonitor.Rd
index ca7ae5a..93aa531 100644
--- a/man/bfastmonitor.Rd
+++ b/man/bfastmonitor.Rd
@@ -175,19 +175,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
@@ -197,19 +197,17 @@ 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
+}
 }
 \references{
 \insertRef{janbfastmonitor}{bfast}
diff --git a/man/bfastts.Rd b/man/bfastts.Rd
index 74a3e78..5a13357 100644
--- a/man/bfastts.Rd
+++ b/man/bfastts.Rd
@@ -45,23 +45,19 @@ bfastts(timedf$y, timedf$dates, type = "16-day")
 # Irregular
 head(bfastts(timedf$y, timedf$dates, type = "irregular"), 50)
 
-
-\dontrun{
+\donttest{
 # 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
+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(as.data.frame(modisbrick))[,1:4]
+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)
 }
-
 }
 \seealso{
 \code{\link[strucchangeRcpp]{monitor}}, \code{\link[strucchangeRcpp]{mefp}},
diff --git a/man/dates.Rd b/man/dates.Rd
index f111aaf..880c296 100644
--- a/man/dates.Rd
+++ b/man/dates.Rd
@@ -4,7 +4,7 @@
 \name{dates}
 \alias{dates}
 \title{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)}
 \source{
 \insertRef{janbfastmonitor}{bfast}
 }
diff --git a/man/modisraster.Rd b/man/modisraster.Rd
index 1d1163f..30d2fd8 100644
--- a/man/modisraster.Rd
+++ b/man/modisraster.Rd
@@ -3,13 +3,13 @@
 \docType{data}
 \name{modisraster}
 \alias{modisraster}
-\title{A raster brick of 16-day satellite image NDVI time series for a small subset
+\title{A raster datacube of 16-day satellite image NDVI time series for a small subset
 in south eastern Somalia.}
 \source{
 \insertRef{janbfastmonitor}{bfast}
 }
 \description{
-A raster brick containing 16-day NDVI satellite images (MOD13C1 product).
+A Cloud-Optimised GeoTIFF containing 16-day NDVI satellite images (MOD13C1 product).
 }
 \examples{
 
diff --git a/man/plot.bfast.Rd b/man/plot.bfast.Rd
index 6b4353d..2bbbe5d 100644
--- a/man/plot.bfast.Rd
+++ b/man/plot.bfast.Rd
@@ -46,7 +46,7 @@ The type of plot shown depends on the value of \code{type}.
 \item{components} Shows the final estimated components with breakpoints.
 \item{all} Plots the estimated components and breakpoints from all iterations.
 \item{data} Just plots the original time series data.
-\item{seasonal} Shows the trend component including breakpoints.
+\item{seasonal} Shows the seasonal component including breakpoints.
 \item{trend} Shows the trend component including breakpoints.
 \item{noise} Plots the noise component along with its acf and pacf.
 }
diff --git a/tests/Examples/bfast-Ex.Rout.save b/tests/Examples/bfast-Ex.Rout.save
index 9bf6fa7..9876dc5 100644
--- a/tests/Examples/bfast-Ex.Rout.save
+++ b/tests/Examples/bfast-Ex.Rout.save
@@ -1,7 +1,7 @@
 
-R version 4.1.0 (2021-05-18) -- "Camp Pontanezen"
-Copyright (C) 2021 The R Foundation for Statistical Computing
-Platform: x86_64-suse-linux-gnu (64-bit)
+R version 4.4.1 (2024-06-14) -- "Race for Your Life"
+Copyright (C) 2024 The R Foundation for Statistical Computing
+Platform: x86_64-suse-linux-gnu
 
 R is free software and comes with ABSOLUTELY NO WARRANTY.
 You are welcome to redistribute it under certain conditions.
@@ -111,7 +111,7 @@ NULL
 	 of optimal 4-segment partition: 
 
 Call:
-confint.breakpointsfull(object = bp.Vt, het.err = FALSE)
+confint.breakpointsfull(object = bp.Vt, breaks = breaks.Vt, het.err = FALSE)
 
 Breakpoints at observation number:
   2.5 % breakpoints 97.5 %
@@ -171,7 +171,7 @@ NULL
 	 of optimal 5-segment partition: 
 
 Call:
-confint.breakpointsfull(object = bp.Vt, het.err = FALSE)
+confint.breakpointsfull(object = bp.Vt, breaks = breaks.Vt, het.err = FALSE)
 
 Breakpoints at observation number:
   2.5 % breakpoints 97.5 %
@@ -361,6 +361,7 @@ RE = 6.3887, p-value < 2.2e-16
 > 
 > 
 > 
+> 
 > cleanEx()
 > nameEx("bfastmonitor")
 > ### * bfastmonitor
@@ -581,52 +582,10 @@ F-statistic: 14.82 on 10 and 227 DF,  p-value: < 2.2e-16
 > AIC(mon$model)
 [1] -428.9029
 > 
-> ## Example for processing raster bricks (satellite image time series of 16-day NDVI images)
-> f <- system.file("extdata/modisraster.grd", package = "bfast")
-> library("raster")
-Loading required package: sp
-> modisbrick <- raster::brick(f)
-> data <- as.vector(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)
-> 
-> ## helper function to be used with the calc() function
-> xbfastmonitor <- function(x, timestamps = dates) {
-+ 	ndvi <- bfastts(x, timestamps, type = c("16-day"))
-+ 	ndvi <- window(ndvi, end = c(2011, 14))/10000
-+ 	## delete end of the time to obtain a dataset similar to RSE paper (Verbesselt et al.,2012)
-+ 	bfm <- bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))
-+ 	return(c(breakpoint = bfm$breakpoint, magnitude = bfm$magnitude))
-+ }
-> 
-> ## 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
-[1] -0.08451347
-> plot(bfm)
-> xbfastmonitor(modisbrick[1], dates) ## helper function applied on one pixel
-   breakpoint     magnitude 
-2011.34782609   -0.08709291 
-> 
-> ## apply the bfastmonitor function onto a raster brick
-> timeofbreak <- raster::calc(modisbrick, fun=xbfastmonitor)
-> 
-> raster::plot(timeofbreak) ## time of break and magnitude of change
-> raster::plot(timeofbreak,2) ## magnitude of change
 > 
 > 
 > 
 > cleanEx()
-
-detaching ‘package:raster’, ‘package:sp’
-
 > nameEx("bfastpp")
 > ### * bfastpp
 > 
@@ -776,23 +735,6 @@ Frequency = 365
 [41]     NA     NA     NA     NA     NA     NA     NA 0.3216     NA     NA
 > 
 > 
-> ## Not run: 
-> ##D # Example of use with a raster
-> ##D 
-> ##D library("raster")
-> ##D f <- system.file("extdata/modisraster.grd", package="bfast")
-> ##D modisbrick <- brick(f)
-> ##D ndvi <- bfastts(as.vector(modisbrick[1]), dates, type = c("16-day")) ## data of pixel 1
-> ##D plot(ndvi/10000) 
-> ##D 
-> ##D # Time series of 4 pixels
-> ##D modis_ts = t(as.data.frame(modisbrick))[,1:4]
-> ##D # Data with multiple columns, 2-4 are external regressors
-> ##D ndvi <- bfastts(modis_ts, dates, type = c("16-day"))
-> ##D plot(ndvi/10000)
-> ## End(Not run)
-> 
-> 
 > 
 > 
 > cleanEx()
@@ -803,7 +745,8 @@ Frequency = 365
 > 
 > ### Name: dates
 > ### Title: A vector with date information (a Datum type) to be linked with
-> ###   each NDVI layer within the modis raster brick (modisraster data set)
+> ###   each NDVI layer within the modis raster datacube (modisraster data
+> ###   set)
 > ### Aliases: dates
 > ### Keywords: datasets
 > 
@@ -841,8 +784,8 @@ Frequency = 365
 > flush(stderr()); flush(stdout())
 > 
 > ### Name: modisraster
-> ### Title: A raster brick of 16-day satellite image NDVI time series for a
-> ###   small subset in south eastern Somalia.
+> ### Title: A raster datacube of 16-day satellite image NDVI time series for
+> ###   a small subset in south eastern Somalia.
 > ### Aliases: modisraster
 > ### Keywords: datasets ts
 > 
@@ -973,7 +916,7 @@ Frequency = 365
 > cleanEx()
 > options(digits = 7L)
 > base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
-Time elapsed:  3.721 0.084 3.84 0.011 0.024 
+Time elapsed:  1.761 0.067 1.827 0 0 
 > grDevices::dev.off()
 null device 
           1