isciences / exactextractr

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

User defined summary function returns summary of coverage fraction #69

Closed amgvandoorn closed 2 years ago

amgvandoorn commented 2 years ago

Hi,

I want to use my own summary function to use with exact_extract. However, this calculates the trait for the coverage_fraction, not for the values of my raster cells. How can I apply it to raster cells please?

My raster has only empty values (NA's) In detail:

  1. exact_extract(raster,poly, "sum")= 0

  2. total=function(values){ sum(values, na.rm=T)} exact_extract(raster, poly, total, summarize_df=T) = 131945.9

Thanks for help !

dbaston commented 2 years ago

Please see the "summary functions" section of the README for some examples: https://github.com/isciences/exactextractr#summary-functions

amgvandoorn commented 2 years ago

Thanks dbaston for your quick reply. I've read it again and if I understand it correctly this example shows how to ignore the coverage fraction, by simply not using it in the function:

Number of distinct raster values within the polygon (coverage fractions are ignored) exact_extract(rast, poly, function(values, coverage_fraction) length(unique(values)))

However, in my situation the following two functions give exactly the same result:

a = exact_extract(rast, poly, function(values, coverage_fraction) sum(values * coverage_fraction, na.rm=TRUE))

b = exact_extract(rast, poly, function(values, coverage_fraction) sum(values , na.rm=TRUE))

But the output that I need is: c = sum(values(crop(raster, poly)), na.rm=T)

and c != b

Do you have a suggestion how to get c by using exact_extract? Thanks!

dbaston commented 2 years ago

c is giving you the sum of values that occupy the bounding rectangle of poly, while a and b are giving you the values that occupy poly. a and b will be identical if the edges of poly follow cell boundaries.

You could get c from exact_extract by constructing a polygon for the bounding rectangle of poly, though it seems easier and more clear to use values(crop()) as you are doing above.