R Terra: Write multiple time series variables to NetCDF

Question

I have multiple variables in a timeseries which I want to save as NetCDF. But when I combine all variables with terra::sds() the time information disappears. I tried the following:

library("terra")


ts <- c("20210101", "20210301", "20210501", "20210701", "20210901", "20211101")
ts_dt <- as.Date(ts, "%Y%m%d")

r1 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r1) <- 1:ncell(r1)
r2 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r2) <- 1:ncell(r2)
r3 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r3) <- 1:ncell(r3)
r4 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r4) <- 1:ncell(r4)
r5 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r5) <- 1:ncell(r5)
r6 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r6) <- 1:ncell(r6)

r_p <- rast(sds(r1, r2, r3, r4, r5, r6))
names(r_p) <- ts_dt
terra::time(r_p) <- ts_dt

Everything is fine when I check the SpatRaster information:

class : SpatRaster
dimensions : 18, 36, 6 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
sources : memory
memory
memory
... and 3 more source(s)
names : 2021-01-01, 2021-03-01, 2021-05-01, 2021-07-01, 2021-09-01, 2021-11-01
min values : 1, 1, 1, 1, 1, 1
max values : 648, 648, 648, 648, 648, 648
time : 2021-01-01 to 2021-11-01

So I create another variable:

s1 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s1) <- 1:ncell(s1)
s2 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s2) <- 1:ncell(s2)
s3 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s3) <- 1:ncell(s3)
s4 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s4) <- 1:ncell(s4)
s5 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s5) <- 1:ncell(s5)
s6 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(s6) <- 1:ncell(s6)

s_p <- rast(sds(s1, s2, s3, s4, s5, s6))
names(s_p) <- ts_dt
terra::time(s_p) <- ts_dt

In a last step I combine both variables with sds():

p <- sds(r_p, s_p)

But p hasn't any time variable:

class : SpatRasterDataset
subdatasets : 2
dimensions : 18, 36 (nrow, ncol)
nlyr : 6, 6
resolution : 10, 10 (x, y)
extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84
source(s) : memory
names : Mean Temp, Mean Prec

When I write p to NetCDF, the resulting .nc file, as expected, hasn't any time information. So my question is as the headline says: How can I write multiple timeseries variables to one NetCDF file while preserving information about time?

Answer 1

A SpatRasterDataset does not have a time dimension. The SpatRasters it contains can have that dimension, and it may be different for each SpatRaster.

Here is a simplified script that suggests that all is good (you did not show what you get, so I cannot compare)

library("terra") 
#terra 1.5.36
ts <- c("20210101", "20210301", "20210501", "20210701", "20210901", "20211101")
ts_dt <- as.Date(ts, "%Y%m%d")

r1 <- rast(xmin=-180, xmax=180,ymin=-90, ymax=90, nrows=18, ncols=36)
values(r1) <- 1:ncell(r1)
# use "c()" instead of "rast(sds())"
r_p <- c(r1, r1, r1, r1, r1, r1)
terra::time(r_p) <- ts_dt
# not setting layer names as that concept does not exist in NetCDF 
s_p <- c(r1, r1, r1, r1, r1, r1)
terra::time(s_p) <- ts_dt

p <- sds(r_p, s_p)
names(p) <- c("first", "second")

p
#class       : SpatRasterDataset 
#subdatasets : 2 
#dimensions  : 18, 36 (nrow, ncol)
#nlyr        : 6, 6 
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source(s)   : memory 
#names       : first, second 

You can see that time is preserved in the SpatRasterDataset

p[1] #or  p["first"]
#class       : SpatRaster 
#dimensions  : 18, 36, 6  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#sources     : memory  
#              memory  
#              memory  
#              ... and 3 more source(s)
#names       : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1, lyr.1 
#min values  :          1,          1,          1,          1,          1,          1 
#max values  :        648,        648,        648,        648,        648,        648 
#time        : 2021-01-01 to 2021-11-01 
 

The the NetCDF file also preserved the time stamps:

writeCDF(p, "test.nc", overwrite=TRUE)

# Read back 
pp <- sds("test.nc")
pp[1]
#class       : SpatRaster 
#dimensions  : 18, 36, 6  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : test.nc:first 
#varname     : first 
#names       : first_1, first_2, first_3, first_4, first_5, first_6 
#time        : 2021-01-01 to 2021-11-01 

# Or read a single sub-dataset

rast("test.nc", 1)
#class       : SpatRaster 
#dimensions  : 18, 36, 6  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#source      : test.nc:first 
#varname     : first 
#names       : first_1, first_2, first_3, first_4, first_5, first_6 
#time        : 2021-01-01 to 2021-11-01 
 
# Or read both sub-datasets as a single SpatRaster 

rast("test.nc")
#class       : SpatRaster 
#dimensions  : 18, 36, 12  (nrow, ncol, nlyr)
#resolution  : 10, 10  (x, y)
#extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#coord. ref. : lon/lat WGS 84 
#sources     : test.nc:first  (6 layers) 
#              test.nc:second  (6 layers) 
#varnames    : first 
#              second 
#names       : first_1, first_2, first_3, first_4, first_5, first_6, ... 
#time        : 2021-01-01 to 2021-11-01 

If you get different results, please use the development version of terra instead: install.packages('terra', repos='https://rspatial.r-universe.dev')