Raster calculation over multiple layers using a function or do-loop

Question

I have a bunch of raster calculations that are working, but there must be a more efficient way to do this with a do-loop or function.

Essentially I have a bunch of combinations of: rates (county level, zip level, etc), which correspond with age group, population datasets for each age group, two sets of concentration data sets, and multiple betas and lower and upper confidence intervals for the betas. I am running the below attributable fraction calculation.

y = rate * population * (1-exp(-beta*concentration))

It's not important to know exactly which combinations I have to run, just that there are a lot of them. I have tried if statements that haven't worked, functions using (i in length) which haven't worked, and multiplying stacked rasters (stacks of the rates, stacks of the pollutants, vectors of the betas, etc).

Doing this for every combination, lower and upper confidence intervals, etc, is not efficient and I am bound to make mistakes.

Is there a straightforward way to get started?

INPUTS - all are on the same resolution, extent, cell size.

#Population rasters
adult.pop
all.pop
eld.pop 
#Rates
zip.all # Zip, all ages, all-cause, rate per 10,000
cbg.25 # CBG, ages 25-99 years, all-cause, rate per 10,000
cbg.65 # CBG, ages 65-99 years, all-cause, rate per 10,000
co.all # County, all ages, all-cause, rate per 10,000
co.25 # County, ages 25-99 years, all-cause, rate per 10,000
co.65 # County, ages 65-99 years, all-cause, rate per 10,000
#Pollutants
no2.v1 # V1 concentration
no2.v2 # V2 concentration
#Stack rates by age group
rates.adult <- stack(cbg.25, co.25) 
rates.all <- stack(zip.all, co.all)
rates.eld <- stack(cbg.65, co.65)
#Betas
beta.adult <- 0.002078254
adult.lower <- 0.000598207
adult.upper <- 0.003536714
elderly.lower <- 0.001105454
elderly.upper <- 0.001365306

CALCULATION

# (1) V1 concentration
elderly.1 <- eld.pop*rates.eld*(10^-4)*(1-exp(-beta.eldery*no2.1))
# (2) V2 concentration
elderly.2 <- eld.pop*rates.eld*(10^-4)*(1-exp(-beta.elderly*no2.2))

• Spacedman: "It's not important to know exactly which combinations I have to run" seems to be the exact opposite of the problem. You need to put all the various components of each category into lists and then make a data frame where each row says which element goes into the calculation. Your calc then looks like M[[N]] = pops[[J]] * rates[[K]] * (10^-4)* (1-exp(-betas[[L]] * concs[[M]])) where J,K,L,M are the elements of row N of the specification data frame...

Answer 1

Here's what I'd do.

Make named lists of the various data items:

pops = list(
 adult.pop=adult.pop,
 all.pop=all.pop,
 eld.pop=eld.pop
)
rates = list(
 zip.all=zip.all, # Zip, all ages, all-cause, rate per 10,000
 cbg.25=cbg.25, # CBG, ages 25-99 years, all-cause, rate per 10,000
 cbg.65=cbg.65, # CBG, ages 65-99 years, all-cause, rate per 10,000
 co.all=co.all # County, all ages, all-cause, rate per 10,000
)
pollutants = list(
 no2.v1 = no2.v1, # V1 concentration
 no2.v2 = no2.v2 # V2 concentration
)

Then make a data frame of the combinations you want to do:

> models
 pop rate pollutant
1 adult.pop co.all no2.v1
2 adult.pop co.all no2.v2
3 eld.pop co.all no2.v1
4 eld.pop co.all no2.v2

then for N from 1 to 4 you can compute something like this:

 pops[[models[N, "pop"] ]] * 
 rates[[models[N,"rate"] ]] + 
 pollutants[[ models[N,"pollutant"] ]]

And in a loop it looks like this:

lapply(1:4, function(N){pops[[models[N, "pop"] ]] * rates[[models[N,"rate"] ]] + pollutants[[ models[N,"pollutant"] ]]})

which for me produces this, because I've created those lists with simple numeric elements rather than rasters or whatever you've got:

[[1]]
[1] 1005
[[2]]
[1] 1006
[[3]]
[1] 1009
[[4]]
[1] 1010

• vatsouth: Wonderful! This seems like a great solution, thanks Spacedman!
• vatsouth: Hi Spacedman - this has worked for keeping it in a dataframe - however, I would like to keep in gridded raster format as well. Do you have a suggestion how to pull from the data frame of combinations into the raster calculation? I have tried several options using overlay but with no success.