2

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

Improve this question
8
  • 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 '20 at 22:58
  • great idea, thanks a lot! – Robin Kohrs Mar 6 '20 at 7:45
  • @Spacedman won't this test against raster cell centroids only? – dbaston Mar 9 '20 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 '20 at 16:18
  • @Spacedman FWIW, here's a simple example where small=TRUE misses an intersection gist.github.com/dbaston/3f21c3c399021f3b5aace12e4f09a495 – dbaston Mar 10 '20 at 14:36
3

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:

library(sf)
library(raster)

r <-
  raster(
    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
Improve this answer
2
  • 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;) – Robin Kohrs Mar 13 '20 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 '20 at 19:26
3

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 
print(sum(matr[,1][!is.na(matr[,1])]))

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.

Improve this answer
2
  • How did it work better for you? – Simbamangu Mar 14 '20 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 – Robin Kohrs Mar 14 '20 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.