4

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?

dimfalk
  • 853
  • 1
  • 5
  • 15

1 Answers1

4

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
Robert Hijmans
  • 40,301
  • 4
  • 55
  • 63
  • 1
    Ah, I see what I did there unconsciously. With `names(r)` being `c("2017-09", "2017-10", "2017-11", ...)` I was able to query them as intended because of partial matching. Tricky! Thank you very much! Having worked with `xts` a lot in the past months I was somehow expecting the z-dimension of the stack to be the time automatically, but that's not true of course for multi-layered raster data. Ok, think I'll consider setting names per default in the future to enable partial matching - and maybe also write some kind of a `subset()` to accept indexing via `"2001-01-01/2001-03-31"` like xts does. – dimfalk Aug 06 '22 at 15:12