2

I have a question concerning the extraction of spatial points that match a certain value of an additional raster in R. The operation itself should actually be pretty simple, but I'm stuck somehow. (I am using R at a "lower intermediate" level).

I am using two data sets. The first one is a shapefile containing point data.

training1 <- shapefile ("re_27042015projcrop_all.shp")

The second one is a "binary mask" raster object consisting solely of 0 and 1 values.

radar1 <- brick ("binary_01_180115_240415.envi")

How do I extract all points whose coordinates match raster pixels with a value of 1?

I tried some logical expressions like

output <- training1@coords & radar1==1

and various raster functions like crop, mask etc., but none seemed to work.

  • By "points whose coordinates match raster pixels" do you mean the point is somewhere within the raster cell or it is exactly at the centre of the cell? – Spacedman Aug 22 '17 at 17:58
  • I mean somewhere within the raster cell. Basically, the idea of my operation is similar to the clip function in ArcGIS, except that I want to use raster cells of a certain value as the clip feature. – bagadosh Aug 22 '17 at 18:07
  • extract(r, pts) will get the raster values at the points, and then you can subset that for values==1. – Spacedman Aug 22 '17 at 18:51
1

You need to create a logic test to select those points. With an example:

library(raster)
library(sp)

set.seed(123)

# raster creation
r <- raster()

# values 0 and 1
r <- setValues(r, sample(x = 0:1, size = ncell(r), replace = T))

# point creation
pts <- spsample(as(r@extent,"SpatialPolygons"),100,type='random')

r is a binary raster and pts is a point layer with 100 features, let's see raster values of points position:

head(extract(r,pts))
## [1] 0 1 0 0 1 0

We have 0 and 1 values, to select only the points with position over values equal to 0, it's necessary to create a logic vector and select rows with this vector:

ones <- which(extract(r == 1,pts) == 1)

pts2 <- pts[ones,]

Now, let's see values of this new layer:

head(extract(r,pts2))

## [1] 1 1 1 1 1 1

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.