0

Consider a vector dataset of 228 observations

precip<-c(3.287516e+00,3.012129e+00,7.450114e+01,5.232168e+00,8.328000e+00,9.426842e+01,3.952095e+02,
 3.116342e+02,2.276140e+02,1.608153e+01,3.258420e+01,0.000000e+00,5.485738e+00,9.192645e+01,
 2.943402e-01,2.280000e-01,1.752476e+01,1.351948e+02,4.482827e+02,3.531432e+02,6.959639e+02,
 1.379384e+02,9.000000e-01,6.720001e-01,0.000000e+00,1.974000e+00,0.000000e+00,2.982000e+00,
 2.726213e+01,1.399524e+02,6.766090e+02,1.935723e+02,4.636400e+01,3.309727e-01,8.099999e-01,
 2.952000e+00,5.683278e+00,7.083724e-03,2.668146e+00,1.628295e+01,3.548195e+01,2.934440e+02,
 3.024660e+02, 1.382988e+02, 3.707115e+01, 5.961257e+01, 0.000000e+00, 3.600000e-02, 0.000000e+00,
 3.055869e+01,1.139865e+01,3.723898e+00,1.036747e+01,1.062829e+02,1.005832e+02,5.299438e+02,
 7.339809e+01, 4.536283e+00, 8.922304e+00, 0.000000e+00, 1.098449e+00,8.125899e+01,6.968601e-03,
 9.240000e-01, 0.000000e+00, 1.238425e+02, 2.596544e+02, 2.374932e+02, 3.476996e+02, 0.000000e+00,
 3.909159e-01, 0.000000e+00, 5.711232e+01, 9.000000e-02, 0.000000e+00, 2.136000e+00, 2.434159e+01, 
 3.015469e+02, 2.186706e+02, 2.877532e+02,6.169636e+01, 2.438971e+01, 2.373664e+00, 1.315564e+00,
 7.095953e+00, 1.740000e-01, 3.288097e+01, 2.742000e+00, 1.638000e+00, 1.247376e+02, 9.973801e+02,
 1.199484e+02, 1.901760e+02, 7.984202e-01, 0.000000e+00, 0.000000e+00, 0.000000e+00, 4.324105e-01,
 7.113221e+01, 4.140000e+00, 5.404327e+01, 8.203735e+01, 4.207430e+02, 4.862388e+02, 1.329529e+02,
 7.671657e+00, 0.000000e+00, 0.000000e+00, 6.360000e-01, 2.518633e+01,0.000000e+00, 2.124000e+00,
 2.677887e+01, 8.775770e+01, 2.946802e+02, 1.800176e+02, 1.420970e+02, 0.000000e+00, 2.400000e-02,
 0.000000e+00,0.000000e+00, 0.000000e+00, 4.649093e+00, 4.426157e+00, 4.831618e+00, 3.774218e+02,
 2.771018e+02, 2.202528e+02, 1.217157e+02, 7.763163e+00, 9.600000e-01, 6.000000e-03, 1.256319e+01,
 0.000000e+00, 6.000000e-02, 1.224000e+00, 1.692969e+01, 1.054706e+02, 3.525440e+02, 2.389812e+02, 
 2.442882e+02, 7.527153e+01, 6.778238e+01, 1.358918e+01, 1.102953e+01, 3.551746e+00, 0.000000e+00,
 1.140000e-01, 4.440000e-01, 7.078428e+01, 3.542093e+02, 2.461660e+02, 1.025227e+02, 3.087413e+01,
 5.262804e+01, 6.300000e-01, 0.000000e+00, 1.676579e+00,1.398000e+00, 1.664469e+00, 2.012430e+01,
 4.735457e+02, 3.973598e+02, 3.091013e+02, 1.914775e+02, 0.000000e+00, 0.000000e+00,0.000000e+00,
 1.831472e+01, 0.000000e+00, 1.524202e-02, 7.800000e-01, 5.477224e+00, 7.821223e+01, 3.506086e+02,
 3.987886e+02,8.562125e+01, 0.000000e+00, 0.000000e+00, 0.000000e+00,0.000000e+00, 2.404258e+01,
 2.571477e+01, 6.611002e+00, 0.000000e+00, 5.205880e+02, 5.405767e+02, 5.016281e+02 ,5.492355e+01,
 5.706278e+01,2.250000e+00, 0.000000e+00, 6.612912e+01, 9.587150e+01, 3.960000e+00, 1.074052e+01,
 1.614000e+00, 9.181252e+01, 4.643691e+02, 1.932034e+02, 1.541364e+02, 9.265393e+00, 0.000000e+00,
 3.816947e+01, 2.979000e+01, 6.305477e+00, 6.317532e+01, 8.266569e+00, 5.865227e+00, 1.847435e+02,
 3.620533e+02, 2.988630e+02 ,5.844741e+01, 0.000000e+00, 3.480975e-01, 0.000000e+00, 2.760000e+00,
 2.010000e+00, 2.712000e+01, 2.640000e+00, 2.123399e+01, 1.456584e+02, 7.809321e+02, 5.340542e+02,
 9.234533e+01, 2.354917e+01, 0.000000e+00,0.000000e+00)

These are monthly precipitation at a location over a number of years. I find the mean of all the Jans, all the Febs, all the Mar etc

x<-rowMeans(matrix(precip, nrow = 12, byrow = F))
x
11.630832  19.372524  16.788150   4.051670  14.857163 186.173774 420.738594 304.162227 161.079527 23.955036   8.945979   3.019485

Then I find the index of the maximum value of x

wetMo<-which(x==max(x)) 
wetMo
7 # 420.738 is the maximum value in x

Then I reorder the precip vector such that the 7th value is the start of the series. I do this by

if(wetMo > 1) precip[c(wetMo:(length(precip) - (13-wetMo)))] else precip

Above reorders precip, such that the first value is 3.952 (which is the 7th value in the original vector).

I want to do the above for a stack of raster over a region. I have 360 rasters of monthly precip for 30 years. I have found the mean of all the Jans, Febs etc and created a raster by

mean_monthly_ppt<-stackApply(pptStack, c(1:12), fun = mean) #12 rasters created which I stack again

Then I have found the maximum precip value in the above stack by

x<-which.max(mean_monthly_ppt_stack)

Now I want to do the reordering with a similar if statement as I have written above for the vector, but in a function that I will use with calc

reorder<-function(r_s){
 if (any(is.na(r_s))) return(r_s*NA)
 r_s_reorder <- if(x > 1) r_s[c(x:((length(r_s) - (13-x)))] else r_s
 return(r_s_reorder)
}
reorder_ppt<-calc(pptStack, fun=reorder)

However, I get the error Error in if (x > 1) r_s[c(x:((length(r_s) - 1) - (13 - x)))] else r_s : argument is not interpretable as logical

I cannot figure this out because I have tried the function on a vector and it works but it does not with the calc().

1
  • I would just create a vector of the new ordering and not attempt to apply a function to calc. You can then simply pass that vector as an index to the stack eg., reorder_ppt <- pptStack[[order.idx]] I do this type of operation to ensure that after reading in a timeseries of separate rasters, that the resulting stack is then ordered correctly. Mar 9 at 18:57
1

The error is because x as you have defined it is a raster layer, not a single value. So when you run the logical test (x > 1), it evaluates to another raster layer and not the TRUE/FALSE value expected by the if statement, hence the error "not interpretable as logical".

You need to rewrite the function to include the calculation of x at a pixel-level, i.e. for a vector rather than RasterStack. Since I don't think there's an equivalent of stackApply() for a vector, you'll need to do the calculation step-by-step. One option is to use the modulo function %% to group the values for like months (all Januarys, Februarys etc.) across the vector and use sapply to calculate the mean for each.

The first line will also give you an error here, as you're returning an NA vector of the same length as r_s, while the later code returns a different length vector. Similarly, I don't think you actually want the if statement here, since if x = 1 you would still want to cut off the final 12 months to avoid returning a vector of a different length to other cells.

Here's a modified version that seems to work on a test rasterstack:

reorder <- function(r_s) {
  if(any(is.nan(r_s)) | any(is.na(r_s[1]))) {
    return(rep(NA, length.out = (length(r_s) - 12)))
  } else {
    months <- 1:12

    # New function to calculate mean for each month across years
    monthly_mean <- function(month_number, vector) {
      month_indices <- (1:length(vector) - 1) %% 12 + 1 == month_number
      month_mean <- mean(vector[month_indices])
      month_mean
    }
    
    # Apply function for each month
    monthly_mean_ppt <- sapply(months, monthly_mean, vector = r_s)

    x <- which.max(monthly_mean_ppt)

    r_s_reorder <- r_s[x:(length(r_s) - (13 - x))]
    return(r_s_reorder)
  }
}

reorder_ppt <- calc(pptStack, fun = reorder)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.