Error when trying to calculate mean and SD of environmental dataset with loop from .nc data

46 views Asked by At

I was trying to calculate mean and SD per month of a variable from an environmental dataset (.nc file of Sea surface temp/day during 2 years) and the loop I used gives me the following error

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'mean': recursive indexing failed at level 2

I have no idea where my error could be but if you are curious I was using the following .nc dataset just for SST for 2018-2019 from copernicus sstdata

Here is the script I used so far and the packages I'm using:

# Load required libraries (install the required libraries using the Packages tab, if necessary)
library(raster)
library(ncdf4)

#Opern the .nc file with the environmental data 

ENV = nc_open("SST.nc")
ENV

#create an index of the month for every (daily) capture from 2018 to 2019 (in this dataset)

m_index = c()
for (y in 2018:2019) {
  # if bisestile year (do not apply for this data but in case a larger year set is used)
  if (y%%4==0) { m_index = c(m_index, rep(1:12 , times = c(31,29,31,30,31,30,31,31,30,31,30,31))) }
  # if non-bisestile year
  else { m_index = c(m_index, rep(1:12 , times = c(31,28,31,30,31,30,31,31,30,31,30,31))) }
}

length(m_index) #  expected length (730)
table(m_index) # expected number of records assigned to each of the twelve months

# computing of monthly mean and standard deviation.
# We first create two empty raster stack...
SST_MM = stack() # this stack will contain the twelve average SST (one per month)
SST_MSD = stack() # this stack will contain the twelve SST st. dev. (one per month)

# We run the following loop (this can take a while)
for (m in 1:12) { # for every month
  
  print(m) # print current month to track the progress of the loop...
  
  sstMean = mean(ENV[[which(m_index==m)]], na.rm=T) # calculate the mean SST for all the records of the current month
  sstSd = calc(ENV[[which(m_index==m)]], sd, na.rm=T) # calculate the st. dev. of SST for all the records of the current month
  
  # add the monthly records to the stacks
  
  SST_MM = stack(SST_MM, sstMean)
  SST_MSD = stack(SST_MSD, sstSd)
  
}


And as mentioned, the output of the loop including the error:

SST_MM = stack() # this stack will contain the twelve average SST (one per month)
> SST_MSD = stack() # this stack will contain the twelve SST st. dev. (one per month)
> for (m in 1:12) { # for every month
+   
+   print(m) # print current month to track the progress of the loop...
+   
+   sstMean = mean(ENV[[which(m_index==m)]], na.rm=T) # calculate the mean SST for all the records of the current month
+   sstSd = calc(ENV[[which(m_index==m)]], sd, na.rm=T) # calculate the st. dev. of SST for all the records of the current month
+   
+   # add the monthly records to the stacks
+   
+   SST_MM = stack(SST_MM, sstMean)
+   SST_MSD = stack(SST_MSD, sstSd)
+   
+ }
[1] 1
**Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'mean': recursive indexing failed at level 2**
1

There are 1 answers

2
Robert Hijmans On

It seems that you make things too complicated. I think the easiest way to do this is with terra::tapp like this:

library(terra)
x <- rast("SST.nc")
xmn <- tapp(x, "yearmonths", mean)
xsd <- tapp(x, "yearmonths", sd)

or more manually:

library(terra)
x <- rast("SST.nc")

y <- format(time(x),"%Y")
m <- format(time(x),"%m")
ym <- paste0(y, "_", m)

r <- tapp(x, ym, mean)