6

I have two rasters in R. I want to set values to NA in the first raster where the second raster has values. I think this should be simple, using the raster package, with two RasterLayer objects raster1 and raster2, both the same extent and snapped to each other. They are 29775x29930.

I'm doing:

newraster <- raster1[is.na(raster2)]

But this seems to take an unnecessary amount of memory, and keeps crashing R. My computer has 8GB of memory. Is there a less memory-intensive way to do this?

6

This can easily be solved using the overlay function from the raster package. Objects rst1 and rst2 are replicates of the initial 'volcano' layer, and a random sample of n = 1000 cells in rst2 is set to NA. Afterwards, overlay is applied and the associated function rejects all cells in rst1 that hold a valid value, i.e. different from NA, in rst2.

library(raster)

rst <- raster::raster(volcano)
rst1 <- rst
rst2 <- rst

# artificial gaps
set.seed(123)
id <- sample(1:ncell(rst), 1000)
rst2[id] <- NA

# introduce na in rst1 for all locations that are non-na in rst2
rst1_mod <- overlay(rst1, rst2, fun = function(x, y) {
  x[!is.na(y[])] <- NA
  return(x)
})

plot(rst1_mod)

volcano_incl_na

1

I would recommend raster::calc because according to the documentation:

For large objects calc will compute values chunk by chunk.

raster2<-reclassify(raster2,c(-Inf,Inf,0))
s->stack(raster1,raster2)
rs1<-calc(s, sum)

You might want to look into this vignette for instructions on processing large rasters.

1

The problem is the way you are addressing the raster with newraster <- raster1[is.na(raster2)]. Try is.na(raster2) by itself to see what you get!

Trying to set values of a raster using vector terminology[] doesn't always work the way one expects. Use values() to set the actual data values of the raster:

library(raster)

# Create two small rasters and set them to random values:
r1 <- r2 <- raster(nrows=5, ncols=5)
values(r2) <- values(r1) <- rnorm(length(r1))

# Set the first row of r1 to NA:
r1[1:5] <- NA

Have a look at it:

head(r1)
           1           2          3           4           5
1         NA          NA         NA          NA          NA
2  0.1451180  0.34801689  1.0545334 -1.15284126 -0.04138286
3 -1.1370768  0.05409194 -0.7767229  0.88499661  0.34104942
4 -0.7654063 -1.03248120  1.1414939 -0.07996859 -0.34718092
5 -0.0110303 -1.70203386  0.8650742 -0.69514811 -0.31484591

Now, set the values of r2 that are not NA in r1 to NA:

values(r2)[!is.na(values(r1))] <- NA

Let's check:

head(r2)
           1         2          3        4          5
1 -0.9571857 0.4479621 -0.5601638 1.207951 -0.1459748
2         NA        NA         NA       NA         NA
3         NA        NA         NA       NA         NA
4         NA        NA         NA       NA         NA
5         NA        NA         NA       NA         NA
  • This doesn't work for me, I still go out of memory. – Herman Toothrot Jan 8 '18 at 17:08
0

Use the subs function (also possible with reclassify, calc, etc.) I demonstrate the process with a small raster first:

raster1 <- raster(ncol=99, nrow=99)
# is the raster small enough to use the following functions
canProcessInMemory(raster1)
# TRUE, so the simple method of assigning values should work
setValues(raster1, 100)
raster1[10:50, 10:50] <- NA
newraster <- subs(raster1, data.frame(old=NA, new=TRUE))

The same substitution function should work with a large raster, but you must store it on disk so that it can be processed in chunks. In my experience, any GIS will do this better than R.

raster2 <- raster(ncol=29775, nrow=29930)
canProcessInMemory(raster2)
# for such a large raster
# you will probably get the answer "FALSE" 
# meaning that setValues(raster2, 100) will not work 
# and you should write your raster to a file before working with it
# writing may take awhile
writeRaster(raster2, tempfile())

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.