I have a rasterstack and a vector-geometry (sf-object). The raster's values are either ǸA or some integer. I would like to check if the raster (in the best case just as a binary response) has any intersection with the shape in the points that are non-NA in the raster. And I don't want to use the bounding box of the raster. I thought about polygonzing the raster where the values are not NA. But this seems a lot of effort and space. I'm not asking for any specific code, but maybe an approach that could be used in R

  • 2
    Can you do extract(raster, polygons) and if you get all NA values then its not intersecting any data values in your raster. That will test against the polygon even clipping a tiny part of a raster cell.
    – Spacedman
    Mar 5, 2020 at 22:58
  • great idea, thanks a lot! Mar 6, 2020 at 7:45
  • @Spacedman won't this test against raster cell centroids only?
    – dbaston
    Mar 9, 2020 at 15:38
  • 1
    The default for extract(raster, polygons) has small=TRUE which returns samples for any intersection of the polygon with a cell, even partial. You can add weights=TRUE to get the area clipped.
    – Spacedman
    Mar 9, 2020 at 16:18
  • @Spacedman FWIW, here's a simple example where small=TRUE misses an intersection gist.github.com/dbaston/3f21c3c399021f3b5aace12e4f09a495
    – dbaston
    Mar 10, 2020 at 14:36

2 Answers 2


You can use extract from the raster package to get a list of overlapping raster and polygon values, and then can check the extracted values - if there are any NAs then the sum of the values will be NA:

Reproducible example, with a raster half integer 1s and half NAs:


r <-
    nrows = 50, ncols = 50,
    xmn = 0, xmx = 5, ymn = 0, ymx = 5,
    crs = CRS("+init=epsg:4326"),
    vals = c(rep(1, 1250), rep(NA, 1250))

# Polygon without overlap of NA area:
p1 <-
  data.frame(id = 1, wkt = 'POLYGON((2 3, 2.5 3, 2.5 3.5, 2 3.5, 2 3))')
p1 <- st_as_sf(p1, wkt = 'wkt', crs = 4326)

sum(unlist(extract(r, p1))) # Result [1] 25

# Polygon totally overlapping with NA area:
p2 <-
  data.frame(id = 1, wkt = 'POLYGON((0 0, 0 0.5, 0.5 0.5, 0.5 0, 0 0))')
p2 <- st_as_sf(p2, wkt = 'wkt', crs = 4326)

sum(unlist(extract(r, p2))) # Result [1] NA
  • sorry a little late, but do you know the velox-package? It's amazing for exactly this task of extracting raster-values based on some vector-data. It's lightning fast;) Mar 13, 2020 at 9:11
  • I wasn't aware of it - if it provides a better solution than this answer, go ahead and answer the question yourself with a good code example!
    – Simbamangu
    Mar 13, 2020 at 19:26

I decided to go the following way. There is a package called velox which does extraction of raster values based on polygons really fast. For some reason it worked better for me using sp than sf-objects. So here's more or less what I did:

rast.norm = raster(rasts[i]) # load it first as RasterLayer
rast.trim = trim(rast.norm) # trim that raster to delete outer NAs
shape.crp = st_crop(shape, rast.trim) 
shape.sp = sf::as_Spatial(shape.crp)
rast.vx = velox(rast.trim) # load runout raster as velox raster
ex.vx = rast.vx$extract(shape.sp) 
matr = do.call(rbind, ex.vx) # rowbind all matrices 

If the sum is 0 or the crop-step doesn't work it means there is no intersection. I think maybe isn't the fastest way, but for me it worked.

  • How did it work better for you?
    – Simbamangu
    Mar 14, 2020 at 8:23
  • It was faster I guess. But actually that's just my own specific solution to the problem at hand. I don't know if it's "better" than any other solution Mar 14, 2020 at 11:22

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.