0

I have two distinct tasks to complete- find the index of the maximum value in each cell of a stack of 12 rasters. There are 12 rasters because each raster is the mean monthly precipitation over 30 years. Then I want to reorder another raster stack of 360 layers, where each layer is the monthly precip over 30 years such that the resulting raster stack has cells that start at the index found to be the of the 12 mean monthly precip values and then continues with the rest of the monthly precip values.

So for example, I have monthly precip values Jan -Dec 1981, 1982 etc. However, I have found that the wettest month is July i.e. 7th month), so I want to reorder such for each pixel the stack starts from July 1981 until Dec 2020.

If I were to do this with a list of vector and not rasters, I would do

twleve_values<-c(11.630832  19.372524  16.788150   4.051670  14.857163 186.173774 420.738594 304.162227 161.079527  23.955036   8.945979   3.019485)
max_tweleve_value<-which.max(twelve_values) #this would result in 7 because the max value 420.738 is the 7th value

if(max_twelve_value > 1) other_360_values[c(max_twelve_value:(length(other_360_values) - (13-max_twelve_value)))] else other_360_values 

I am stuck at the first task of finding the index of the wettest month in the stack of 12 rasters. I have tried

x<-stackApply(12_raster_stack, indices = rep(1,12), fun=which.max, na.rm = TRUE)

Error in FUN(newX[, i], ...) : unused argument (na.rm = TRUE) and

x<-calc(mean_monthly_ppt, fun= whiches.max)

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function

How do I move ahead?

1
  • As a nudge take a loos at the base function sort.int with the index.return = TRUE argument. It allows you to sort a vector AND returns the index of the sort. As such you can use this index to referent the order in the raster stack. It does not make much sense to reorder your pixels but, you could calculate a global max for each raster, store in a vector and then pass this vector to sort.int. The raster stack could then be reordered using the resulting index values in a double bracket. Something like r[[sort.int(max.vals, index.return = TRUE)$ix]] Commented Sep 20, 2021 at 17:54

2 Answers 2

0

For the first task use calc:

library(raster)

r <- raster()

set.seed(123)

l <- list()

for(i in 1:12){
  l[[i]] <- setValues(r,sample(c(NA,1:100),size = 64800,replace = T))
}

s <- stack(l)

max_ <- calc(s,which.max)

I used NA values for showing you the function works with that kind of values:

plot(max_)

enter image description here

I'm not really sure about your second task. If you're trying to reorder layers by pixel, I don't think you can accomplish it. You can reoder pixels among layers, but you aim isn't clear for me.

4
  • Your suggestion does not work. I have tried calc(stack, which.max) and the error I get now is Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘writeValues’ for signature ‘"RasterBrick", "list"’
    – tg110
    Commented Mar 1, 2021 at 20:57
  • I tried just which.max(s) and it worked. I do not understand why this is working with which.max() and not with calc() when logically it should work with calc()
    – tg110
    Commented Mar 1, 2021 at 21:22
  • The error says the issue, did you try loading rasters with stack or brick?
    – aldo_tapia
    Commented Mar 2, 2021 at 1:02
  • I loaded all the rasters and made it a stack. As per my comment, I tried which.max(stack) and it worked. I am moving on to the second task now
    – tg110
    Commented Mar 2, 2021 at 8:00
0

I know that this is an old post but, it is still quite relevant (now using terra). Honestly, I have found that this type of timeseries analysis is most efficiently performed at the function level rather that manipulating the source data, in the case the raster stack. When using functions such as terra:app (formally raster::calc) you are functionally operating on individual vectors representing the timeseries. This means that you can write a function that performs any type of reordering, in support of your specific analysis. I often have a vector, independent of my raster stack, representing a dates class object, so that I can call it from a function and specify a ts or zoo (timeseries) object in support formal temporal statistical models (eg., partialled out residuals using an ARIMA term).

So, rather than thinking about the raster stack, simplify the problem and just think about the timeseries vector and how you would reorder the vector, especially since you want pixel-level results. However, it would be best to not write out a new raster stack but rather perform reordering in a function that implements the analysis that is necessitating it in the first place.

I did actually directly address the question in a comment; "As a nudge take a look at the base function sort.int with the index.return = TRUE argument. It allows you to sort a vector AND returns the index of the sort. As such you can use this index to reference the order in the raster stack. It does not make much sense to reorder your pixels but, you could calculate a global max for each raster, store in a vector and then pass this vector to sort.int. The raster stack could then be reordered using the resulting index values in a double bracket. Something like r[[sort.int(max.vals, index.return = TRUE)$ix]] "

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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