I am trying to convert the outputs of a model into a raster based on both a temperature raster and a seasonality raster. I have over 2100 bands so ifelse statements aren't going to cut it.

An example would be that I had my model outputs in a dataframe - with the range of temperature values in column 1, the seasonality in column 2 and the outputs in column 3:


And I have a raster of temperature

t = raster(matrix(1:100,10,10))

And a raster of seasonality

s = raster(matrix(rep(1:10,10),10,10))

I need to write a new raster which is based on the values of s and t, but consists of the model outputs of column 3 of the data frame.

  • Is df supposed to have all the unique combinations of 1:100 and 1:10? Because I don't think it does. df[df$a==35,] shows ten rows, all with b=5 and with 10 different values of c... – Spacedman Feb 4 at 11:42
  • Maybe you mean to start df as expand.grid(a=1:100, b=1:10) and then set the c column? – Spacedman Feb 4 at 11:44
  • Oops yes it is, my bad that was sloppy on my part - should be fixed now – redferry Feb 4 at 14:01

First make a raster of NA in the same shape as t (and s):

> v = raster(t); v[]=NA

Now extract the values from t and s into a data frame, do a merge (like a left join) with the lookup table, extract the c column and replace the NA values in v:

> v[] = merge(data.frame(a=t[], b=s[]), df, all.x=TRUE)$c

and I think v is now what you want. Any values that don't match in the lookup will be NA.

  • Having given this a bash I get the error: 'Warning message: In v[] <- value : number of items to replace is not a multiple of replacement length' – redferry Feb 4 at 15:02
  • It works for your test data as in the question doesn't it? Then something in the data you are trying doesn't match the pattern of your test data. But make sure you have the edit with all.x=TRUE version, because if you have non-matching data that will fix it. – Spacedman Feb 4 at 16:40
  • I note you mention 2100 bands so maybe you are trying this on a raster stack or brick? That's not tested... – Spacedman Feb 4 at 17:04
  • No, no stack or brick just two rasters. Just trying to work through why it isn't working with the actual data. It seems to be going wrong at the merge phase - it's difficult for me to recreate the issue with simulated rasters but going to keep trying. – redferry Feb 4 at 17:14
  • Store the merge somewhere (merged = merge(...)) and see if its the same number of rows as the number of cells in your target raster... If not, see which rows are missing... – Spacedman Feb 4 at 17:17

The above works with an extra step, where you order the dataframe. This needs to be done because merge re-orders your data frame meaning that all the values end up in the wrong raster squares when it is written to the empty raster. So the solution is:

v = raster(t); v[]=NA

f = data.frame(a=t[], b=s[])

vec = c(1:nrow(f))

f[,3] = vec

m = merge(f, df, all.x=TRUE)

colnames(m)[3] = "ord"

m = m[order(m$ord),]

v[] = m$c

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.