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