0
# created rasters from .tif file (r1,r2,r3,r4,r5,r6,r7)‎

‎# created raster stack‎
r_stack <- stack(r1,r2,r3,r4,r5,r6,r7)‎
class      : RasterStack ‎
dimensions : 743, 1893, 1406499, 7  (nrow, ncol, ncell, nlayers)‎
resolution : 0.00832774, 0.00832774  (x, y)‎
extent     : 9.391724, 25.15614, 26.9779, 33.16542  (xmin, xmax, ymin, ymax)‎
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ‎
names      : npp_2000, npp_2001, npp_2002, npp_2003, npp_2004, npp_2005, npp_2006 ‎
min values :   0.0085,   0.0009,   0.0000,   0.0000,   0.0061,   0.0000,   0.0031 ‎
max values :   1.2009,   1.1083,   1.2283,   1.1489,   1.1562,   1.1747,   1.1447‎

‎# added date slot to my stack‎
my_s <- setZ(r_stack, dt[,1], "Date")‎
my_s
class      : RasterStack ‎
dimensions : 743, 1893, 1406499, 7  (nrow, ncol, ncell, nlayers)‎
resolution : 0.00832774, 0.00832774  (x, y)‎
extent     : 9.391724, 25.15614, 26.9779, 33.16542  (xmin, xmax, ymin, ymax)‎
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ‎
names      : npp_2000, npp_2001, npp_2002, npp_2003, npp_2004, npp_2005, npp_2006 ‎
min values :   0.0085,   0.0009,   0.0000,   0.0000,   0.0061,   0.0000,   0.0031 ‎
max values :   1.2009,   1.1083,   1.2283,   1.1489,   1.1562,   1.1747,   1.1447 ‎
Date        : 2000-03-01, 2001-03-01, 2002-03-01, 2003-03-01, 2004-03-01, 2005-03-01, 2006-03-01‎
raster
6
  • Unless your data is fairly small, you do not want to convert a raster stack to a zoo object. The best way is to write a function that runs bfast on a vector, including coercion to a zoo object. Then you pass this function to raster::calc or raster::overlay Oct 9, 2020 at 1:11
  • Jeffrey I tried to coerce a vector created from brick data frame. The problem is how do I get the dates because bfast run only ts objects timeseries Oct 9, 2020 at 22:43
  • You need to create an argument in your function that accepts the date vector and set the default to your dates argument or create a dates vector within the function. When I am back on my computer I will post an example. Oct 10, 2020 at 0:01
  • I am trying different approach by building dataframe. Duplicate one of my raster layers ID <- r1. Make index layer ID <- setValues(ID, 1:ncell(ID)). Create a vector for each pixel v <- getValues(ID). Coerce the vector to a dataframe d <- as.data.frame(v). Populate columns for time series d$t1 <- getValues(r1); d$t2 <- getValues(r2); etc. Not sure if create vector for dates using as.date function how to include it in my dataframe and create ts object time series. Really appreciat if you can help. Oct 23, 2020 at 16:01
  • You can simplify your approach considerably by just using as(D, "SpatialPixelsDataFrame") This will result in an sp object with a data.frame in the @data slot. Each row will be the time-series for a given pixel. You can then simply use apply to run a function for each row (treated as a vector). There have been many instance where the raster packge would not run a function that ran effortlessly on an sp pixels object. The cavieat is that the entire raster stack is read into memory but, if you have the memory can also easly multthred it using someting like `parallel::mcmapply. Oct 23, 2020 at 16:58

1 Answer 1

1

Here is a fairly simple approach to running the BFAST model in R on a temporal raster array.

First, you have to install the bfastSpatial package from GitHub (not on CRAN). You will need to first install the devtools package. There is no compile necessary so, you should not need RTools however, this is a very useful install for installing packages that require compiling. You can go download the install or use the installr::install.Rtools function.

Install package from GitHub and then add required libraries. We can use the irregular time-series tura data for example.

devtools::install_github("loicdtx/bfastSpatial")

library(raster)
library(bfastSpatial)
  data(tura)

Here we extract the provided dates object but, you can create your own in a Dates format (not Zoo).

( d <- tura@z[[1]] )

Fill NA values in the time-series using a lowess nonlinear regression

plot(tura[[c(1,100)]])  
  tura <- smooth.time.series(tura, f = 0.8, smooth.data = FALSE)
    plot(tura[[c(1,100)]])   

Here we just pass the start argument from befastmonitor function but, you can also use the Dates (d) object.

bfm <- bfmSpatial(tura, start=c(2009, 1), order=1)

Now that we have an implementation of the BFAST model, I would highly recommend looking at Forkle's work on the phenological analysis of remote sensing time-series data, which is an expansion of the BFAST model. Of his proposed methods, the season-trend model (STM) is an approximation of the BFAST model. There are, however, other options available including Annual Aggregated Time Series (AAT) and Trend Seasonal Adjusted (TSA). The greenbrown package is distributed on R-Forge and can be installed as such.

install.packages("greenbrown", repos="http://R-Forge.R-project.org")

We will use the regular time-series data provided.

library(greenbrown)
  data(ndvimap)

Here is the trend and breakpoint analysis akin to the BFAST model using the STM method. The AAT method is however, preferred.

( bfm.gb <- TrendRaster(ndvimap, start=c(1982, 1), freq=12, method="STM", 
                        breaks=1, funAnnual=max) )
  plot(bfm.gb)

You can also derive metrics of land surface phenology and greenness

pheno <- PhenologyRaster(ndvimap, start=c(1982, 1), freq=12, tsgf="TSGFspline", 
                           approach="Deriv")

  plot(pheno, c(grep("SOS.1982", names(pheno)), 
       grep("EOS.1982", names(pheno)))) 

References

Forkel, M., Carvalhais, N., Verbesselt, J., Mahecha, M., Neigh, C., Reichstein, M., 2013. Trend Change Detection in NDVI Time Series: Effects of Inter-Annual Variability and Methodology. Remote Sensing 5, 2113–2144. doi:10.3390/rs5052113

Forkel, M., Migliavacca, M., Thonicke, K., Reichstein, M., Schaphoff, S., Weber, U., Carvalhais, N., 2015. Codominant water control on global interannual variability and trends in land surface phenology and greenness. Glob Change Biol 21, 3414–3435. doi:10.1111/gcb.12950

3
  • Great thanks. Really appreciate your help Jeffrey. The first part works well. However, The TrendRaster seems doesn't accept my dates. Here is how I created my dates: created a vector for dates d then I used NPP_Brick <- setZ(NPP_Brick, d). When I check the class for "d" returns "Date". I tried to use the "d" instead of start = c(2000,3) but both didn't work Oct 26, 2020 at 15:24
  • I checked and the bfmSpatial is, in fact, not accepting date vectors. This may be tracing back to the bfastmonitor function and may be something to check into. This is something to contact the authors of the package over, or file a bug report on their GitHub repository, just make sure to trace the source of the error to the correct package. Oct 26, 2020 at 15:54
  • I tried to run the function smooth.time.series for another set of raster brick/stack and getting this warning "In FUN(newX[, i], ...) : Fewer than 8 real-value observations, assigning NA" any Idea please? Nov 2, 2020 at 1:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.