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 Answers
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
-
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! 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
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.
-
-
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
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.small=TRUE
which returns samples for any intersection of the polygon with a cell, even partial. You can addweights=TRUE
to get the area clipped.small=TRUE
misses an intersection gist.github.com/dbaston/3f21c3c399021f3b5aace12e4f09a495