Search code examples
rgispolygonrasterterra

terra::extract gives NaN when polygon is smaller than raster cell

In the following example, the extract function tells correctly us that the mean value of r within the polygon x2 is 5.14. However, for polygons like x1 what are smaller than the raster, extract returns a value of "NaN"

r <- rast(nrows = 10, ncol = 10, nlyrs = 1, vals = sample(1:10, 100, replace = TRUE), names = "temp")

x1 <- rbind(c(-145,-10), c(-145,-5), c(-140, -5), c(-140,-10))
x2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55))
z <- rbind(cbind(object=1, part=1, x1, hole=0),
           cbind(object=3, part=1, x2, hole=0))
colnames(z)[3:4] <- c('x', 'y')
p <- vect(z, "polygons")

plot(r)
plot(p, add = T)

test <- terra::extract(r, p, fun = mean, cell = TRUE)

test
  ID     temp
1  1      NaN
2  2 5.142857

How can I get the value of r at x1enter image description here?

Solution

  • You can use exact=TRUE

    Example data

    library(terra)
r <- rast(nrows = 10, ncols = 10, nlyrs = 1, vals =1:100, names = "temp")
x1 <- rbind(c(-145,-10), c(-145,-5), c(-140, -5), c(-140,-10))
x2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55))
z <- rbind(cbind(object=1, part=1, x1), cbind(object=2, part=1, x2))
colnames(z)[3:4] <- c('x', 'y')
p <- vect(z, "polygons")

    The default

    extract(r, p, fun = mean)
#  ID temp
#1  1  NaN
#2  2   53

    With touches=TRUE you get all cells that are touched

    extract(r, p, fun = mean, touches=TRUE)
#  ID     temp
#1  1 51.50000
#2  2 52.62069

    Or you can do

    e <- extract(r, p, exact=TRUE)
head(e)
#  ID temp    fraction
#1  1   51 0.007539715
#2  1   52 0.030169815
#3  2   19 0.104219078
#4  2   28 0.282198174
#5  2   29 0.883159178
#6  2   30 0.043386000

x <- by(e[,2:3], e[,1], function(x) weighted.mean(x[,1], x[,2]))
as.vector(x)
# [1] 51.80006 52.21312

    (or use dplyr or data.table if you are familiar with that syntax)

    With the development version (1.3.11), available from install.packages('terra', repos='https://rspatial.r-universe.dev'), you get:

    extract(r, p, fun=mean)
#  ID temp
#1  1 51.5
#2  2 53.0

    And you can do

    extract(r, p, fun=mean, exact=TRUE)
#     ID     temp
#[1,]  1 51.80006
#[2,]  2 52.21312