isciences / exactextractr

R package for fast and accurate raster zonal statistics
https://isciences.gitlab.io/exactextractr/
272 stars 26 forks source link

Error when polygon does not overlap raster if include_area = TRUE #68

Closed see24 closed 2 years ago

see24 commented 2 years ago

There is an error if the polygon does not overlap with the raster when the include_area = TRUE but if include_area = FALSE it returns a 0 row data frame. It would be great if the behaviour was consistent.

library(raster)
library(exactextractr)
library(sf)

rast <- raster(matrix(0, nrow = 100, ncol = 100))

nonoverlap_poly <- st_sf(st_sfc(st_polygon(list(matrix(c(0, 0, 1, 0, 1,
                                                         -0.25, 0, -0.25, 0, 0),
                                                       ncol = 2, byrow = TRUE)))))

# returns 0 row dataframe without include_area

exact_extract(rast, nonoverlap_poly)
# [[1]]
# [1] value             coverage_fraction
# <0 rows> (or 0-length row.names)

# but error with include_area = TRUE

exact_extract(rast, nonoverlap_poly, include_area = TRUE)

# Error in CPP_exact_extract(x, weights, wkb, default_value, default_weight,  : 
#                              logical subsetting requires vectors of identical size
dbaston commented 2 years ago

Thanks for the report!

MatthieuStigler commented 2 years ago

Somehow related to this issue, I get sometimes the error: Error: [crop] extents do not overlap, even using the version just updated (0.8).

When does this error arise? Put differently, how come this error does not arise int he case above?

Thanks!

dbaston commented 2 years ago

Thanks for the feedback, @MatthieuStigler . The current development version (0.8) tries to provide reasonable performance for many different types of inputs. A downside is that it has many different code paths to do so. I suspect you get the error with a multi-layer input raster that is stored on disk that intersects none of the input geometries. Can you please let me know if 3eae9a0d0b4e98b9edeb0017e1e849ebfe75f89f fixes the issue for you?

Looks like this fix relies on the most recent version of terra. I've updated the DESCRIPTION to match.

MatthieuStigler commented 2 years ago

Thanks for your rapid response @dbaston !

So the problem was still there using the version with your fix, but simply disappeared once I updated to the latest version of terra (from 1.4.22 to 1.5-17). Sorry, I should have checked that before!

Just for the record, this happened on a single-layer raster, but a quite huge one. Weirdly enough, it did happen when I was querying exact_extract with multiple geometries, but not with individual ones. Here get_1() is just a wrapper for exact_extract(), where the sole argument is indexing my shapefile.

get_1(4674:(4674+1))

 Error: [crop] extents do not overlap

> get_1(4675)
# A tibble: 0 × 5
# … with 5 variables: MUN_CODE <dbl>, prp_typ <chr>, value <dbl>, n_pixels <int>, sum_area <dbl>
> get_1(4674)
# A tibble: 0 × 5
# … with 5 variables: MUN_CODE <dbl>, prp_typ <chr>, value <dbl>, n_pixels <int>, sum_area <dbl>
dbaston commented 2 years ago

Thanks for checking in, and please let me know if you run into any other issues. It's really helpful to have testing of pre-CRAN versions!