ENH: Add fixed precision set operations for GEOS >= 3.9

[brendan-ward]

This PR adds support for the fixed precision variants for set operations introduced in GEOS 3.9.

In ufuncs.c, these were copied and adapted from the YY_Y functions to handle an extra scalar parameter. Since these needed to be version-guarded in C, it seemed reasonable to call into these only if the Python wrapper function determines that the grid_size parameter is non-None and GEOS >= 3.9; otherwise it falls back to the original functions.

The grid_size is scalar for now; this seems like a reasonable simplification because mixed-precision within an array seems a bit awkward.

It looks like there are a few errors due to mixed dimension GeometryCollections not being accepted in the new fixed precision functions. I think this is a known issue in GEOS / JTS, but need to track down the cross-references.

Note: this is built against #257

[brendan-ward]

Based on my understanding of ufunc.reduce(), it doesn't look possible to pass in the grid_size parameter in order to add fixed precision to the set operations that use reduce: difference_all, intersection_all, symmetric_difference_all.

Luckily union_all uses a different approach, so it was easily adaptable to use fixed precision.

May need to postpone adding fixed precision to the other *_all functions to a later PR.

[jorisvandenbossche]

Cool, thanks a lot!

I suppose the easiest is first to get #257 in, so this gets a bit easier to review? (is there actually a reason it built on the other PR, it doesn't seem to depend on it?)

The grid_size is scalar for now; this seems like a reasonable simplification because mixed-precision within an array seems a bit awkward.

+1

May need to postpone adding fixed precision to the other *_all functions to a later PR.

Yes, that seems fine to leave for a follow-up. Having the grid_size option is the plain element-wise ufuncs would already be great!

[jorisvandenbossche]

One user question that comes up: doing a set_precision(geom, grid_size=1) followed by a normal intersection, does that give the same result as doing an intersection(...., grid_size=1) ?

[brendan-ward]

GEOS ticket for mixed dimension overlay errors: https://trac.osgeo.org/geos/ticket/1092

[brendan-ward]

is there actually a reason it built on the other PR, it doesn't seem to depend on it

I was using the other PR to test (manually) if the precision was set on the output geometries after using one of these fixed-precision operations. It is not set (e.g., the precision of the result of union with grid_size=1 is 0, even though results have been rounded). Does this seem like something we should raise with GEOS? The documentation there is not specific on this front.

doing a set_precision(geom, grid_size=1) followed by a normal intersection, does that give the same result as doing an intersection(...., grid_size=1)

Based on simple tests, these are the same. My understanding is that for GEOS >= 3.9, calling set_precision will use the same precision model in GEOS as the fixed-precision overlay functions - so the effect to the inputs of the operation should be the same. It is also my understanding that if a grid_size was set for an input geometry to the set operations, that they would automatically use whichever precision is higher of the inputs, so in this case of both are set to 1 the result should be the same as if the grid_size=1 parameter is used for the set operation.

[brendan-ward]

Benchmarks (*_prec1 denotes grid_size=1 and *_prec2 denotes grid_size=2):

time_difference                                    11.3±0.2ms
time_difference_prec1                                29.1±1ms
time_difference_prec2                              30.1±0.6ms
time_intersection                                  12.6±0.6ms
time_intersection_prec1                              30.4±4ms
time_intersection_prec2                            29.4±0.4ms
time_symmetric_difference                          12.0±0.3ms
time_symmetric_difference_prec1                      30.6±2ms
time_symmetric_difference_prec2                      33.3±5ms
time_union                                         12.1±0.6ms
time_union_all                                       33.6±7ms
time_union_all_prec1                                 76.7±8ms
time_union_all_prec2                                82.0±10ms
time_union_prec1                                   28.0±0.7ms
time_union_prec2                                     31.0±2ms

[brendan-ward]

Benchmark test geometry is randomly generated to look something like this (transparency to show overlap):

[jorisvandenbossche]

Looks good!

Mainly a bunch of minor docstring comments (I made them on the first, but they also apply for the other docstrings), and some comments on the tests.

[jorisvandenbossche]

Maybe also assert that the precision of the result is equal to grid size? (so we actually are sure we properly passed the grid size param, as I suppose this test would pass as well at the moment if the implementation had a mistake there)

(Eg for union_all below you have a few cases with actual input vs expected based on grid_size, which ensures things work as expected)

[brendan-ward]

Unfortunately, GEOS currently returns 0 for the precision of the outputs when grid_size is provided as a parameter here.

I added a test to compare the results to the output of the operation against geometries that have been set to the same grid_size in advance via set_precision.

[jorisvandenbossche]

Unfortunately, GEOS currently returns 0 for the precision of the outputs when grid_size is provided as a parameter here.

Ah, that's interesting to know (it's probably not "unfortunately", but intentionally done that way in GEOS?). So that is then still a difference between doing a set_precision first vs specifying the grid_size keyword. When setting the precision, that gets propagated to the output, while passing the keyword, the grid size is only used "within" the operation, but not set on the output.

[brendan-ward]

Added clarifying note to docstring. Not sure if this is intentional in GEOS.

[jorisvandenbossche]

Just checked on the geos chat channel, and Martin Davies confirmed that it is intentional

[jorisvandenbossche]

Quoting him "The idea is to give more flexibility in using precision for specific operations, without affecting the underlying geometry precision model"

[jorisvandenbossche]

Thanks for the updates!

[jorisvandenbossche]

(needs a merge of master, conflict with #266)

[caspervdw]

Thanks @brendan-ward and @jorisvandenbossche ! Merging right before the 0.9 release