Select polygon search tool

[vkbo]

We need to select an appropriate tool that can handle the polygon intersect functionality and which can be run from Python.

[fzeiser]

ne possibility:

• shapely

It is based on the widely deployed GEOS (the engine of PostGIS) and JTS (from which GEOS is ported) libraries. Shapely is not concerned with data formats or coordinate systems, but can be readily integrated with packages that are.

Simple method to get the overlap ratio, e.g.:
poly1.intersection(poly2).area / poly1.area [criticism: might be slower then "only" calculating the overlap ratio]

• Similar functionality: spherical_geometry

• You can use the GDAL/OGR Python bindings for that, GDAL/OGR Cookbook .

Requirements:

• May or may not implement spherical geometry
• Test for more complex polygons, see eg.
• what happens for non-closed linear rings (e.g. due to computational inaccuracy? (linear rings: define cube by 5 points, where p1 and p5 should be the same)

[fzeiser]

Interesting issue here with many suggestions: https://github.com/astropy/regions/issues/46

I see Boost also has some polygon stuff: https://www.boost.org/doc/libs/1_70_0/libs/polygon/doc/index.htm

But again if you look at the spherical polygon code - it's very complex (like in the Google S2 code): https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/strategies/spherical/point_in_poly_winding.hpp

[vkbo]

Shapely seems to be built on a port of JTS, which is the package I was suggesting. So in principle Shapely and JTS should be equivalent.

There's also a package for SolR called Geo3D that works on ellipsoid geometry in addition to Cartesian, but I don't think the polygon areas will be large enough that this is needed. Something to discuss anyway.

[fzeiser]

Used shapely on the kartverket data of #10 to get a feeling (not an exact test) on the timings etc.

• defining a complex polygon: ~ 60 µs

• area of a complex polygon: 5 µs

• intersection of complex polygon: ~ 3 ms

• Some few kommunes actually consist of several polygons (I will therefore call them "true" multipoligons.), e.g. Malvik in Trøndelag. For this specific example, I checked that the sum of area of the intersections of the sub-polygons is the same same as the area of the intersection of the whole kommune with the fylke [which makes sense, right].

• The area of a of a polygon (polygon.area) seems to be calculated when the function is called, and the property is not saved (thus, calling it once again, the it will be recalculated)

Test-file: (ipynb) https://app.zenhub.com/files/376049307/46aa6315-ccc7-4619-ae9f-af36b50de9c4/download

[fzeiser]

There's also a package for SolR called Geo3D that works on ellipsoid geometry in addition to Cartesian, but I don't think the polygon areas will be large enough that this is needed. Something to discuss anyway.

I think we can drop ellipsoid/spherical geometry in the first implementation and potentially open for this later.

[fzeiser]

Update on the timings: As expected, the calculation of the intersection polygon is dependent on the number of vertices. Shapely/GEOS has a simplify method with a given tolerance.

A "slight" simplification (tolerance parameter = 800), will produce following polygons (solid line, vs dashed). The upper map is of Malvik, the lower one is of Trøndelag. The gives a speedup of a factor 10, whilst it is difficult to imagine that this should have a practical impact. (Due to the resolution / scaling, the simplification is not visible by eye on this scale for the Trøndelag picture.)

[vkbo]

Closing this since we have definitely decided to go with shapely. :)