CODEWITHSUNDEEP
rr-starsdata-cube
Rocìo Beatriz Cortès
Rocìo Beatriz Cortès

Reputation: 13

stars R package: add a new dimension in stars object

I've been trying to solve a problem for days. I downloaded and imported 4 bioclimatic variables from CHELSA for 3 different months: January, February, and March. I created three stacks, one for each month, containing the variables. Here is how the files look:

> jan_stack
class       : SpatRaster 
dimensions  : 146, 212, 4  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 13.01653, 14.78319, 41.68319, 42.89986  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : CHELSA_bio_01_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_07_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_13_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_14_1981-2010_V2.1_clipped.tif  
names       : CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1 
min values  :           -0.95,            16.2,            59.8,            18.9 
max values  :           16.65,            29.5,           241.7,            92.8 

> feb_stack
class       : SpatRaster 
dimensions  : 146, 212, 4  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 13.01653, 14.78319, 41.68319, 42.89986  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : CHELSA_bio_01_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_07_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_13_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_14_1981-2010_V2.1_clipped.tif  
names       : CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1 
min values  :           -0.95,            16.2,            59.8,            18.9 
max values  :           16.65,            29.5,           241.7,            92.8 

> mar_stack
class       : SpatRaster 
dimensions  : 146, 212, 4  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 13.01653, 14.78319, 41.68319, 42.89986  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
sources     : CHELSA_bio_01_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_07_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_13_1981-2010_V2.1_clipped.tif  
              CHELSA_bio_14_1981-2010_V2.1_clipped.tif  
names       : CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1, CHELSA_~10_V2.1 
min values  :           -0.95,            16.2,            59.8,            18.9 
max values  :           16.65,            29.5,           241.7,            92.

My goal is to create data cubes, so I need to add the temporal dimension. In this case, since these are climate series, each layer represents a monthly average over a 30-year period. Therefore, I thought of simply assigning a representative date of the month to each stack.

But if I proceed to create a 'stars' object, I get an object with a single attribute and three dimensions.

> jan_stars <- st_as_stars(jan_stack)
> jan_stars
stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min. 1st Qu. Median     Mean 3rd Qu.  Max.
CHELSA_bio_01_1981-2010_V2....  -0.95   16.55   28.1 50.30252    71.1 241.7
dimension(s):
           from  to offset     delta refsys point                                                        values x/y
x             1 212  13.02  0.008333 WGS 84 FALSE                                                          NULL [x]
y             1 146   42.9 -0.008333 WGS 84 FALSE                                                          NULL [y]
attributes    1   4     NA        NA     NA    NA CHELSA_bio_01_1981-2010_V2.1,...,CHELSA_bio_14_1981-2010_V2.1

'split' option

At this point, I tried to convert the 'attributes' dimension into separate attributes using 'split', but I lost one dimension and couldn't find any way to add it back while keeping the four attributes in place.

> split(jan_stars)
stars object with 2 dimensions and 4 attributes
attribute(s):
                               Min. 1st Qu. Median      Mean 3rd Qu.   Max.
CHELSA_bio_01_1981-2010_V2.1  -0.95    9.35  12.35  12.03118   15.45  16.65
CHELSA_bio_07_1981-2010_V2.1  16.20   25.30  27.20  25.69690   27.80  29.50
CHELSA_bio_13_1981-2010_V2.1  59.80   84.90 111.20 119.10602  147.40 241.70
CHELSA_bio_14_1981-2010_V2.1  18.90   33.20  43.90  44.37600   53.60  92.80
dimension(s):
  from  to offset     delta refsys point x/y
x    1 212  13.02  0.008333 WGS 84 FALSE [x]
y    1 146   42.9 -0.008333 WGS 84 FALSE [y]

manually add date

However, if I manually add the date, I'm left with only one attribute and I can no longer recover the other three. I need to keep this information as it will be needed for creating an SDM.

# Random January date
my_date <- as.Date("2009-01-01")

# Since I have 4 bands, I repeated the same date
repeated_dates <- rep(my_date, 4)

jan_time <- stars::st_set_dimensions(jan_stars,
      3,
      values = repeated_dates,
      names = "date")
jan_time

stars object with 3 dimensions and 1 attribute
attribute(s):
                                 Min. 1st Qu. Median     Mean 3rd Qu.  Max.
CHELSA_bio_01_1981-2010_V2....  -0.95   16.55   28.1 50.30252    71.1 241.7
dimension(s):
     from  to offset     delta refsys point                    values x/y
x       1 212  13.02  0.008333 WGS 84 FALSE                      NULL [x]
y       1 146   42.9 -0.008333 WGS 84 FALSE                      NULL [y]
date    1   4     NA        NA   Date    NA 2009-01-01,...,2009-01-01

I need to obtain something like this:

# 3 dimensions instead of 2, 4 attributes
stars object with 3 dimensions and 4 attributes
attribute(s):
                               Min. 1st Qu. Median      Mean 3rd Qu.   Max.
CHELSA_bio_01_1981-2010_V2.1  -0.95    9.35  12.35  12.03118   15.45  16.65
CHELSA_bio_07_1981-2010_V2.1  16.20   25.30  27.20  25.69690   27.80  29.50
CHELSA_bio_13_1981-2010_V2.1  59.80   84.90 111.20 119.10602  147.40 241.70
CHELSA_bio_14_1981-2010_V2.1  18.90   33.20  43.90  44.37600   53.60  92.80
dimension(s):
  from  to offset     delta refsys point x/y
x    1 212  13.02  0.008333 WGS 84 FALSE [x]
y    1 146   42.9 -0.008333 WGS 84 FALSE [y]
# ... and the time here

I have the impression that the solution is simple, but I just can't seem to find it. I really hope for some help. Thank you so much in advance!

Upvotes: 1

Views: 70

Answers (0)

Related Questions