Search code examples
rdataframerasterr-raster

Extracting pixels from a raster based on specific value of another raster using R

I imported two rasters (raster A and B)in R using the raster function. I would like to extract the pixels of A where B equals 1 into a data frame. I am trying the following, however, all the pixels I obtain have the same values, although they are various in the original dataset.

The two rasters have the same dimensions (ncols, nrows, ncell, resolution, extent, projection).

library(raster)
library(rgdal)

# import inputs
A <- raster('/pat/to/rasterA.tif')
B <- raster('/pat/to/rasterB.tif')

# extract raster values from A over raster B where B == 1
mydata <- data.frame(A[B[B == 1]])

EDIT 1

Might be that when I do A[B[B == 1]], the class of object A and B from RasterLayer becomes numeric, and this creates problems? I discovered this by doing class(A[B[B == 1]]), which gives numeric.

EDIT 2

Ok, this is weird. I tried to do mydata <- data.frame(A[B]) and now the output has the original A only at B == 1 locations. Trying this before it extracted all the pixels from A (as I would expect). I can coinfirm it is right by counting the number of ones in B and the number of elements in mydata, which is the same. It's like if the indexing has skipped all the zeros in B. Can anyone explain this?

Solution

  • Example data:

    library(raster)
r <- raster(nrow=5, ncol=5, xmn=0, xmx=1, ymn=0, ymx=1)
set.seed(1010)
A <- setValues(r, sample(0:5, ncell(r), replace=TRUE))
B <- setValues(r, sample(0:2, ncell(r), replace=TRUE))

    Now you can do:

    s <- stack(A,B)
v <- as.data.frame(s)
v[v[,2] == 1, 1]

    Alternatively:

    A[B==1]

    Or:

    D <- overlay(A, B, fun=function(x,y){ x[y!=0] <- NA; x})
na.omit(values(D))

    Or:

    xy <- rasterToPoints(B, function(x) x == 1)
extract(A, xy[,1:2])

    Or:

    A[B!=1] <- NA
rasterToPoints(A)[, 3]

    etc...

    Now why does this: A[B[B == 1]] not work? Unpack it:

    B[B == 1]
# [1] 1 1 1 1 1 1 1 1 1 1

    The cell values of B where B==1 are, of course, 1. A[B[B == 1]] thus becomes A[c(1,1,1,..)], and returns the value of the first cell many times.

    A[B] is equivalent to A[B!=0] as B is considered a logical statement in this case, and 0 == FALSE and all other values are TRUE