Use GEOSGridIntersectionFractions_r() directly in st_interpolate_aw

[edzer]

st_interpolate_aw() is for area-weighted interpolation, estimating attribute values from geometry set A to an arbitrary geometry set B, by using averaging (intensive variables) or summing/redistributing (extensive variables). The naive implementation computes all A-B intersections, and uses their areas for weighting. Since GEOS 3.14, GEOSGridIntersectionFractions_r (migrated from @dbaston's exactextracter) computes such weights much more efficiently in case A or B is (regularly) gridded; and tbh that is probably a very common, if not dominant, use case. The current PR implements this both for x (source geoms) being a regularly gridded stars object, as well as to (target geoms).

The st_interpolated_aw.stars() has been added to sf in the current CRAN release for that reason, and needs to be removed from stars (which happened in r-spatial/stars@50ad7e8).

[edzer]

@dbaston when doing so, I discovered that when calling GEOSGridIntersectionFractions_r repeatedly (once for every geometry), I had to set the float output buffer values all to zero before every call; cells not intersecting seem to keep their old value. it would be good if this went into the docs, or got fixed.

[edzer]

@dbaston found out the hard way

[dbaston]

cells not intersecting seem to keep their old value

It's actually cells outside the geometry's bounding box that keep their old value. This behavior is wrong, and a fix is in libgeos/geos#1532

Still, you may want to adjust your usage to crop the raster down to each bbox before the call. Otherwise, a lot of effort is spent copying zeros around.