How to index individual layers from a SpatRaster object by time?-CodePudding

I found out the other day, more or less by chance, that it is possible to query layers from SpatRaster objects based on the time attribute in general (c.f. here), e.g based on years (r["2017"]) and yearmonths (r["2017-10"]).

Now I wanted to deep-dive a little bit into this because of the great flexibility you receive when working with spatio-temporal data. Unfortunately, I seem to be failing from the beginning and can't reproduce the behaviour from before because of the following error: [subset] no (valid) layer selected

library(terra)
#> terra 1.6.3

r <- rast(ncols = 10, nrows = 10, nlyr = 365)


time(r) <- seq(from = as.POSIXlt("2001-01-01", tz = "UTC"),
               to = as.POSIXlt("2001-12-31", tz = "UTC"),
               by = "day")

r
#> class       : SpatRaster 
#> dimensions  : 10, 10, 365  (nrow, ncol, nlyr)
#> resolution  : 36, 18  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> time        : 2001-01-01 to 2001-12-31 UTC

time(r) |> class()
#> [1] "POSIXct" "POSIXt"

r["2001"]
#> Error: [subset] no (valid) layer selected


time(r) <- as.Date("2002-01-01")   0:364

r
#> class       : SpatRaster 
#> dimensions  : 10, 10, 365  (nrow, ncol, nlyr)
#> resolution  : 36, 18  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> time (days) : 2002-01-01 to 2002-12-31

time(r) |> class()
#> [1] "Date"

r["2002"]
#> Error: [subset] no (valid) layer selected

I had a look at ?time as well as at the docs at rspatial.org but was not able to find relevant documentation related to indexing.

What am I missing here?

CodePudding user response：

You are mixing up layer names (that may look like a time-stamp) with an actual time-stamp.

Here is how you can use time-stamps (using Date here, as it is a bit less verbose)

library(terra)
r <- rast(ncols = 10, nrows = 10, nlyr = 365)
time(r) <- as.Date("2001-01-01")   0:364
head(names(r))
# "lyr.1" "lyr.2" "lyr.3" "lyr.4" "lyr.5" "lyr.6"

You can subset layers by name like this

r[["lyr.1"]]

Note the double brackets for sub-setting layers, although you can use single brackets when using a name as opposed to a numerical index.

To subset by time, you can do

r[[time(r) == as.Date("2001-06-01")]]
#class       : SpatRaster 
#dimensions  : 10, 10, 1  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#time (days) : 2001-06-01

With dates, you do not even need to use "as.Date"

r[[time(r) >= "2001-12-01"]]
#class       : SpatRaster 
#dimensions  : 10, 10, 31  (nrow, ncol, nlyr)
#resolution  : 36, 18  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#time (days) : 2001-12-01 to 2001-12-31