0

I am having trouble figuring out the best way to create a new raster from values in two separate rasters. Some solutions involve creating stacks or bricks, others do not.

Here is an example:

First raster is "propSummPCPmod4"

Determine frequency of values

freq(propSummPCPmod4)

 value    count

[1,] 1 1417084

[2,] 3 9103585

[3,] NA 14535331

Second raster is "propWintPCPmod4"

Determine frequency of values

freq(propWintPCPmod4)

 value    count

[1,] 2 241512

[2,] 3 10279157

[3,] NA 14535331

create a brick and plot both rasters side by side

brick1 <- brick(propSummPCPmod4, propWintPCPmod4)
plot(brick1, col=terrain.colors(255))

plot of raster brick

The plot shows that values of 1 from "propSummPCPmod4" [layer 1] and values of 2 from "propWintPCPmod4" [layer 2] do not overlap.

I want to a create a new raster and to employ the "1" values from propSummPCPmod4, the "2" values from propWintPCPmod4, and to have all other values represented by "3" or NA. (NA values are outside the boundaries of the original source rasters)

I tried merge, cvr, overlay, and mask but couldn't get a satisfactory result.

Improve this question

3 Answers 3

Reset to default
1

Your example was not reproduceble so I tryed to create a small set imitating your data. If I understood right, ou can do:

`library(raster)
propSummPCPmod4 <- raster(matrix(c(1,3,3,1,3,3,1,1,NA,NA,NA,NA,NA,NA,NA,NA),4,4))
propWintPCPmod4 <- raster(matrix(c(NA,NA,NA,NA,NA,NA,NA, NA,2,3,3,2,2,2,2,3),4,4))
newraster <- raster(matrix(c(NA,NA,NA,NA,NA,NA,NA, NA,NA,NA,NA,NA,NA,NA,NA, NA),4,4))

summer_indices <- which(propSummPCPmod4[] == 1)
winter_indices <- which(propWintPCPmod4[] == 2)
newraster[summer_indices] <- 1
newraster[winter_indices] <- 2
newraster[is.na(newraster[])] <- 3
plot(newraster)`

propSummPCPmod4 and propWintPCPmod4

enter image description hereenter image description here

newraster

enter image description here

Improve this answer
4
  • I added a reproducible example to clarify my question. Does your suggestion still apply? I haven't used the "which" function from the raster package. Thanks
    – Sean Basquill
    Sep 24, 2019 at 1:08
  • I ran your code using a new blank raster (same structure as the originals outlined above). The newraster doesn't plot as expected. photos.app.goo.gl/dgUPvCQzBzG6D64f7
    – Sean Basquill
    Sep 24, 2019 at 14:14
  • check the edit I made on my answer and see if this is what you need
    – MarujoRe
    Sep 24, 2019 at 17:08
  • I tried again using your code @MarujoRe and it worked! There had been a small formatting mistake in my new new raster code The result looks way better now. However, values of "3" and "NA" are not differentiated. See photos.app.goo.gl/oFBggKnaqmPvpF5L6 I tried this modified version of your code but it didn't work: summer_indicesV3 <- which(propSummPCPmod4[] == 1,3) winter_indicesV3 <- which(propWintPCPmod4[] == 2,3) newrasterV3[summer_indices] <- 1 newrasterV3[winter_indices] <- 2 newrasterV3[winter_indices] <- 3 Sorry to pester you, it’s very close now!
    – Sean Basquill
    Sep 25, 2019 at 15:07
0

For this case, I would use min value between both rasters using overlay() function:

(reproducible example)

library(raster)

r <- raster()

r1 <- setValues(r,3)

r2 <- setValues(r,3)

r1[1:20000] <- 1
r2[44800:64800] <- 2

plot(stack(r1,r2), col=terrain.colors(255))

enter image description here

Minimum between both rasters:

r3 <- overlay(r1,r2,fun=min)

plot(r3)

enter image description here

If you want to keep the other values as NA:

r3[r3==3] <- NA

plot(r3)

enter image description here

Improve this answer
1
  • Thanks @aldo_tapia! I will try this solution next week and see if works for me. Initially I had ran some trials with the overlay function but couldn't get them to provide the result I need.
    – Sean Basquill
    Sep 27, 2019 at 15:44
0

Try to separate which are not NA and which are value 3

summer_indices <- which(propSummPCPmod4[] == 1,3)

summer_notNA <- which(!is.na(propSummPCPmod4[]))

newraster[summer_notNA] <- 3

winter_indices <- which(propWintPCPmod4[] == 2,3) 

winter_notNA <- which(!is.na(propWintPCPmod4[]))

newraster[winter_notNA] <- 3

newraster[summer_indices] <- 1 

newraster[winter_indices] <- 2
Improve this answer
4
  • That should have worked @MarujoRe. There seems to be something wrong with the "3" values, however. I tried trouble shooting and it happens right after the summer indices are added to the new raster. This plot has both summer and winter indices. photos.app.goo.gl/yidQmSUZo9XT6fPm8 You can see that value "3" is only at the top of the map. All of North America should be value 3, except cells that are values 1 or 2. NA values are only those outside the continent but within the raster extent (corresponds with the oceans). Like this photos.app.goo.gl/rG3L1J2REwrjvuM57
    – Sean Basquill
    Sep 26, 2019 at 14:30
  • I uploaded the two source rasters to make things easier. I think the two reproducible example you created should work, but maybe there’s something we’re missing from the data structure. drive.google.com/open?id=1HW_uVi8Vpm96sJ7sptaoeQSYyQXzNrc0 and drive.google.com/open?id=1Eqt46I3xka01_CYAI6mcnzRXM7Efs7iD
    – Sean Basquill
    Sep 26, 2019 at 16:50
  • ohh sorry i have made the wrong atribution. I have edited the message. Now it should work. if it doesn't I will download your raster when I get home.
    – MarujoRe
    Sep 27, 2019 at 12:19
  • That worked! I am just learning the raster package and building my R skills. I needed this solution for my work and I learned new coding functions and syntax too! Thanks so much.
    – Sean Basquill
    Sep 27, 2019 at 15:43

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.