isciences / exactextractr

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

Error in weighted.mean #22

Closed ssaxe-usgs closed 4 years ago

ssaxe-usgs commented 4 years ago

When using weighted.mean as the summary function with multiple raster layers in: exact_extract(ras, shp, fun = weighted.mean)

where:

I receive the error:

Error in weighted.mean.default(vals, cov_fracs, ...) : 'x' and 'w' must have the same length

Looking into the methods for getMethod("exact_extract" , signature = c( x = "Raster" , y = "sfc_POLYGON")) shows that the error is occurring because weighted.mean(x = vals, w = cov_fracs) requires that input w have the same "length" as x. Thus, when object vals is created in the function through raster::getValuesBlock() for a "RasterStack", weighted.mean fails because it can't handle multiple columns. This can be fixed on the users end by defining a custom version of weighted.mean() using apply():

apply.weighted.mean <- function(x, y){apply(X = x, MARGIN = 2, FUN = weighted.mean, w = y)}

v <- exact_extract(x = ras, y = shp, fun = apply.weighted.mean)

dbaston commented 4 years ago

I guess this is the problem with having one method that does a lot of different things. The behavior for a RasterStack is to provide the summary function with the values of all layers in the stack at simultaneously. From the README:

If exact_extract is called with a RasterStack instead of a RasterLayer, the R summary function will be called with a data frame of raster values and a vector of coverage fractions as arguments. Each column in the data frame represents values from one layer in the stack, and the columns are named using the names of the layers in the stack.

It looks like you were able to define a function that does what you need, but do you think something should be changed in the package code/documentation?

Side note: I'd expect you'd get better performance for your use case with the equivalent:

v <- exact_extract(x = ras, y = shp, fun = 'mean')

ssaxe-usgs commented 4 years ago

Ah, I think my head was a little muddled from a long day, I completely missed the paragraph you quoted. It might be helpful to future users if there was an example summarizing a RasterStack in the README. Otherwise, fantastic package. The processing speeds are incredible.