2

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?

1 Answer 1

4

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')

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.