2

This is the function that I want to create:

Erosivity_Factor = 1.74 * log[sum(montly²/annual)]+1.29

I have solved it using a loop:

# Reproducible RasterStack with 12 layers
Monlty<- stack(lapply(1:12, function(i) raster(ncol=20, nrow=20, vals=rnorm(400,30)))) 

# Reproducible Raster
Annual <- raster(ncol=20, nrow=20, vals=rnorm(400,230))

R_month<- list()

for (i in 1:12) {
  R_month[[i]]<- (Monlty[[i]]*Monlty[[i]]/Annual)
          }

R_month_raster <- do.call(stack,R_month)
R_erosivity_factor<- (1.74*log10(calc(R_month_raster,sum))+1.29)

Is there a way to make it with a function?

1 Answer 1

2

First, a reproducible example:

library(raster)

r <- raster()

l <- list()

for(i in 1:12){
  l[[i]] <- setValues(r,rep(i,ncell(r)))
}

s <- stack(l)

You need to create your function and then apply it with calc function. Is annual the sum of all month values? or is annual an average? Well, I considered the sum of all monthly values, you can change it later:

Erosivity_Factor <- function(x){1.74*log10((x^2)/sum(x))+1.29}

R_erosivity_factor <- calc(s,Erosivity_Factor)

Since the value of the raster increases every month, results have sense:ç

plot(R_erosivity_factor,zlim=c(-2.1,1.8))

enter image description here


For only two raster objects the approach is direct:

R_erosivity_factor <- 1.74*log10((Monlty^2)/Annual)+1.29
3
  • It is a good approach, neverless I need to compute the function with two Raster objects, since Annual != sum of the monthly values, and calc() admits only one Raster object. Jan 12 at 1:18
  • I have fixed it but it takes three times longer than the loop. I had thought that a function would speed up the process: x <- stack(Monlty,Annual) Erosivity_Factor <- function(x){1.74*log10(sum(x[c(1:12)]^2/(x[c(13)])))+1.29} R_erosivity_factor <- calc(x,Erosivity_Factor) Jan 12 at 3:15
  • @franciscocorvalan I added an approach with two rasters objects. Although I don't really understand the composition of both rasters. Is monthly raster a 12-band raster while annual is 1-band?
    – aldo_tapia
    Jan 12 at 11:23

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.