rspatialterra

terra freq fails when a polygon is too small (no overlaping pixel) between SpatRaster and SpatVector


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!

EDIT:

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))

Solution

  • 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')