I was wondering if there is a bug in terra::freq
function since it returns an error when a polygon is too small and do not overlap enough a pixel (even if it's on top). terra::extract
works fine.
# Make a raster
r <- rast(nrows=10, ncols=10)
set.seed(2)
values(r) <- sample(5, ncell(r), replace=TRUE)
freq(r)
# Make a second raster but change the extent and make it as vector
r2 <- rast(ncols=2, nrows=2)
values(r2) <- 1:ncell(r2)
ext(r2) <- ext(-100, -20, 0, 50) # lots of overlap
p <- as.polygons(r2)
# compare the 2 (they overlap)
plot(r)
plot(p, add=TRUE)
# freq works fine
freq(x = r, zone = p)
# Extract works fine
extract(x = r, y = p, exact=TRUE)
# Make a third raster but change the extent and make it as vector
r3 <- rast(ncols=2, nrows=2)
values(r3) <- 1:ncell(r3)
ext(r3) <- ext(-25, -20, 0, 2) # almost no overlap
p <- as.polygons(r3)
p
plot(r)
plot(p, add=TRUE)
# get the frequency, but see how it fails
freq(x = r, zone = p)
# Error in `$<-.data.frame`(`*tmp*`, "zone", value = 1L) :
# replacement has 1 row, data has 0
# Extract works fine
extract(x = r, y = p, exact=TRUE)
For the situation I have freq
is much faster as I just need the count of pixels per polygon. extract
could work but is much slower and need more steps to get the count.
I don't mind getting an NA for the small polygon, but not an error!
Note that in my situation, I'm using an sf
object that I change to a SpatVector vect(sf::st_as_sf(p))
for the zones. In this case, it doesn't seem to be a resolution issue.
s <- sf::st_as_sf(p)
plot(r)
plot(s, add=TRUE)
freq(x = r, zone = vect(s))
With terra 1.8-56 you do not get this error. And you can also use argument touches=TRUE
.
library(terra)
#terra 1.8.56
r <- rast(nrows=10, ncols=10)
set.seed(2)
values(r) <- sample(5, ncell(r), replace=TRUE)
r3 <- rast(ncols=2, nrows=2)
values(r3) <- 1:ncell(r3)
ext(r3) <- ext(-25, -20, 0, 2) # almost no overlap
p <- as.polygons(r3)
freq(x = r, zone = p)
#NULL
freq(x = r, zone = p, touches=TRUE)
# layer value count zone
#1 1 4 1 1
#2 1 4 1 2
#3 1 4 1 3
#4 1 4 1 4
terra 1.8-56 is currently the development version. You can install it with
install.packages('terra', repos='https://rspatial.r-universe.dev')