isciences / exactextractr

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

Cannot use weights="area" with weighted_frac #107

Closed goergen95 closed 2 months ago

goergen95 commented 2 months ago

Hi,

thank you for this amazing work!!

I am trying to obtain a weighted_frac summary with weights = "area". The docs state:

https://github.com/isciences/exactextractr/blob/efa5d50d0301df1477a883ac6f6c504c1cf5d8a1/R/exact_extract.R#L172-L175

and:

https://github.com/isciences/exactextractr/blob/efa5d50d0301df1477a883ac6f6c504c1cf5d8a1/R/exact_extract.R#L134-L136

So, given the repex below, I assume that we currently cannot obtain area weighted summaries for weighted_frac? If that is actually the case, would it be possible to add that feature?

library(exactextractr)
library(sf)
#> Linking to GEOS 3.12.2, GDAL 3.9.1, PROJ 9.4.1; sf_use_s2() is TRUE
library(terra)
#> terra 1.7.78

f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
s <- scales::pretty_breaks(n = 5)
breaks <- s(values(r))
r2 <- classify(r, breaks)
v <- vect(ext(r2))
v <- st_as_sf(v)
st_crs(v) <- st_crs(r)

exact_extract(r2, v, fun =  "frac")
#>       frac_0    frac_1    frac_2    frac_3     frac_4
#> 1 0.02083333 0.2858073 0.4292535 0.2419705 0.02213542
exact_extract(r2, v, fun =  "weighted_frac", weights = "area")
#>   weighted_frac_0 weighted_frac_1 weighted_frac_2 weighted_frac_3
#> 1      0.02089308       0.2863966       0.4296582       0.2410372
#>   weighted_frac_4
#> 1      0.02201497

Created on 2024-08-10 with reprex v2.1.0

dbaston commented 2 months ago

This appears to be working as expected. There isn't much variation in cell area over this polygon, so taking it into account doesn't change the results by much.

goergen95 commented 2 months ago

You are right, I now see that this indeed works as documented. I guess what I am asking is if there is a way to calculate the area per category without using an UDF? Initially, I thought combining weighted_frac with weights = "area" would do this. Sorry for the noise! Please feel free to close this issue anytime.

dbaston commented 2 months ago

Maybe simpler than using a UDF would be to convert the fraction to an area after the fact, like

library(dplyr)
library(stringr)

v$area <- st_area(v)
exact_extract(r2, v, fun='frac', append_cols='area') |>
  mutate(across(starts_with('frac'),
                ~ .x * area)) |>
  rename_with(str_replace, 'frac', 'area', .cols=starts_with('frac')) |>
  select(-area)

#           area_0           area_1           area_2           area_3          area_4
# 1 98683299 [m^2] 1353811417 [m^2] 2033287011 [m^2] 1146165331 [m^2] 104850999 [m^2]
goergen95 commented 2 months ago

Yes, this does work, however, requires additional code. Over at mapme.biodiversity we frequently are interested in calculating the area per category and I was hoping this was possible in a more generalized way. Would love to see this implemented here (e.g. something like frac_coverage or similar), but I also understand if you consider this out of scope.

dbaston commented 2 months ago

I've added an upstream issue for this.

goergen95 commented 2 months ago

Awesome, thanks!