rgisrasternetcdftemporal

How to split a raster stack by seasons in R


Issue

I have downloaded multiple Lm3 4 km netcdf files (n = 4640) for sea surface temperature from the Ocean Colour Project, NASA from 2012 to 2024.

I have put all my files into a raster stack:

class       : SpatRaster 
dimensions  : 766, 709, 1  (nrow, ncol, nlyr)
resolution  : 0.04165021, 0.04165796  (x, y)
extent      : 24.61, 54.14, 5.8, 37.71  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 
source      : AQUA_MODIS.20120101.L3m.DAY.SST.x_sst.nc:sst 
varname     : sst (Sea Surface Temperature) 
name        :      sst 
unit        : degree_C 
> sst_stack
class      : RasterStack 
dimensions : 766, 709, 543094, 4640  (nrow, ncol, ncell, nlayers)
resolution : 0.04165021, 0.04165796  (x, y)
extent     : 24.61, 54.14, 5.8, 37.71  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : Sea.Surface.Temperature.1, Sea.Surface.Temperature.2, Sea.Surface.Temperature.3, Sea.Surface.Temperature.4, Sea.Surface.Temperature.5, Sea.Surface.Temperature.6, Sea.Surface.Temperature.7, Sea.Surface.Temperature.8, Sea.Surface.Temperature.9, Sea.Surface.Temperature.10, Sea.Surface.Temperature.11, Sea.Surface.Temperature.12, Sea.Surface.Temperature.13, Sea.Surface.Temperature.14, Sea.Surface.Temperature.15, ... 
   

During fieldwork, we collected data for counts of dolphins at coral reefs per month per year. My dataframe spans from 2012 to 2024, and winter is 7 months and summer is 5 months.

I want to split the raster stack into two seasons (summer and winter):

I have read many sources to see how to do this and this solution looks like the closest. I have no idea how to perform this action as I'm a novice

I'm sorry, I can't share the data due to ownership issues

Many thanks if anyone can point me in the right direction.

R-Code

#Download Sea Surface Temperature 
folder_SST<"~/Documents/GIS_Data/SST."
files_SST <-list.files(folder_SST, pattern='*.nc', full.names ="TRUE")

for (file in files_SST) {
  nc_SST=open.nc(files_SST)
  print.nc(nc_SST)
}

#Stack the AQUA-MODIS files with raster stack
sst_stack<-raster::stack(files_SST)

This is what my dataframe looks like

Survey_Number       Dates Longitude_E_DD Latitude_N_DD    SST
              1   2012-08-01       33.89083      27.26778 23.635
              2   2012-06-02       33.86782      27.40854 23.640
              3   2012-02-07       33.86230      27.44623 23.690
              4   2012-02-12       33.88653      27.26957 23.635
              5   2012-02-13       33.88766      27.26848 23.635
              6   2012-02-14       33.85000      27.36111 23.780
              7   2012-02-15       33.86177      27.41302 23.640

Solution

  • Assuming your SpatRaster layer names:

    Sea.Surface.Temperature.1
    Sea.Surface.Temperature.2
    ...

    correspond to:

    AQUA_MODIS.20120101.L3m.DAY.SST.x_sst.nc
    AQUA_MODIS.20120102.L3m.DAY.SST.x_sst.nc
    ...

    or whatever your actual file dates are, you could first rename your layers so they are 'intuitively' named with date strings corresponding to their .nc source files, then use grep() and regular expressions to split the stack. Ultimately, it would be better to have date strings in the original SpatRaster to remove any uncertainty, but if my assumption is correct, this will work.

    First, an example SpatRaster stack with two years of data:

    library(terra)
    
    r_list <- paste0("sst.", 1:731)
    
    set.seed(42)
    
    rr <- lapply(r_list, function(x) {
      
      r <- rast(nrows = 100, ncols = 100, ext(-180, 180, -90, 90))
      values(r) <- round(runif(10000, 5, 30),2)
      names(r) <- as.character(x)
      return(r)
      
    })
    
    sst_stack <- rast(rr)
    
    sst_stack
    # class       : SpatRaster 
    # dimensions  : 100, 100, 731  (nrow, ncol, nlyr)
    # resolution  : 3.6, 1.8  (x, y)
    # extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
    # source(s)   : memory
    # names       : sst.1, sst.2, sst.3, sst.4, sst.5, sst.6, ... 
    # min values  :  5.01,     5,     5,  5.01,     5,     5, ... 
    # max values  : 30.00,    30,    30, 29.99,    30,    30, ... 
    

    Now create a vector of new layer names with dates the same length as nlyr(sst_stack) and rename your SpatRaster layers. In this example, I have consecutive dates. If your dataframe contains all 4640 dates and these dates correspond, in order, to all of the layers in your SpatRaster stack, for you it will be something like r_list1 <- paste0("sst.", df$Dates) (assuming df is your dataframe).

    # Create vector of equal length to nlyr(sst_stack) to rename SpatRaster layers with dates
    r_list1 <- paste0("sst.", seq(as.Date("2012-01-01"), as.Date("2013-12-31"), by = "day"))
    
    # Rename SpatRaster layes
    names(sst_stack) <- r_list1
    
    sst_stack
    # class       : SpatRaster 
    # dimensions  : 100, 100, 731  (nrow, ncol, nlyr)
    # resolution  : 3.6, 1.8  (x, y)
    # extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
    # source(s)   : memory
    # names       : sst.2~01-01, sst.2~01-02, sst.2~01-03, sst.2~01-04, sst.2~01-05, sst.2~01-06, ... 
    # min values  :        5.01,           5,           5,        5.01,           5,           5, ... 
    # max values  :       30.00,          30,          30,       29.99,          30,          30, ... 
    
    # Note the names are truncated in the print(), but they are like:
    head(names(sst_stack))
    # "sst.2012-01-01" "sst.2012-01-02" "sst.2012-01-03" "sst.2012-01-04" "sst.2012-01-05" "sst.2012-01-06"
    

    Now use two regular expression to create two separate SpatRaster stacks based on month string matches:

    # Create two separate SpatRasters using regex to match month values in layers
    sst_summer <- sst_stack[[grep("-0[5-9]-", r_list1, value = TRUE)]]
    sst_winter <- sst_stack[[grep("-(0[1-4]|1[0-2])-", r_list1, value = TRUE)]]
    
    sst_summer
    # class       : SpatRaster 
    # dimensions  : 100, 100, 306  (nrow, ncol, nlyr)
    # resolution  : 3.6, 1.8  (x, y)
    # extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
    # source(s)   : memory
    # names       : sst.2~05-01, sst.2~05-02, sst.2~05-03, sst.2~05-04, sst.2~05-05, sst.2~05-06, ... 
    # min values  :           5,        5.01,           5,           5,        5.00,           5, ... 
    # max values  :          30,       30.00,          30,          30,       29.99,          30, ... 
    
    sst_winter
    # class       : SpatRaster 
    # dimensions  : 100, 100, 425  (nrow, ncol, nlyr)
    # resolution  : 3.6, 1.8  (x, y)
    # extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
    # coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
    # source(s)   : memory
    # names       : sst.2~01-01, sst.2~01-02, sst.2~01-03, sst.2~01-04, sst.2~01-05, sst.2~01-06, ... 
    # min values  :        5.01,           5,           5,        5.01,           5,           5, ... 
    # max values  :       30.00,          30,          30,       29.99,          30,          30, ...
    

    One advantage of renaming each layer is that you can then check the correct layers have been carried using names(sst_summer) and names(sst_winter).