Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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?

like image 241
umbe1987 Avatar asked Feb 08 '23 00:02

umbe1987


2 Answers

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

like image 147
Robert Hijmans Avatar answered Feb 11 '23 16:02

Robert Hijmans


this should work for the 1 values:

A[B == 1]

as well as this, for the 0 values

A[B == 0]
like image 27
Carlos Alberto Avatar answered Feb 11 '23 17:02

Carlos Alberto