Shapely: Support for precision models other than "float"

[sgillies]

Upcoming, I will require a fixed-precision mode of Shapely. It's something that's been discussed off and on, but now I have an actual need for it at work and am beginning to plan to make it happen.

My initial idea is that Shapely will no longer hide the GEOS context (this struct: https://github.com/libgeos/libgeos/blob/svn-trunk/capi/geos_ts_c.cpp#L137) from users, but will adapt it into a Python context manager (as in Fiona and rasterio). That context manager would be where the precision model would be specified, the usage to be something like this:

import shapely
from shapely.geometry import Point
with shapely.context(precision='fixed', scale=1.0, offset=(0.0, 0.0)):
    # make a point snapped to a fixed-precision grid
    p = Point(0, 0)

Calling shapely.context() will call initGEOS (from libgeos) and instead of creating a default geometry factory (see https://github.com/libgeos/libgeos/blob/svn-trunk/capi/geos_ts_c.cpp#L209) will create a GEOS precision model and a geometry factory based upon the precision model. All geometry operations in that context will use that geometry factory and precision model. When the with block ends, finishGEOS (https://github.com/libgeos/libgeos/blob/svn-trunk/capi/geos_ts_c.cpp#L255) will be called and the precision model and geometry factory will be deleted.

Without a context manager, the default floating point model will be used and older code will continue to work as usual. A context manager will be required to use a fixed model. The specific syntax is yet to be determined and I welcome suggestions.

News from the GEOS project is that exposing the precision model classes via the C API might be a GSoC thing for a student. The good of this is that someone will write that code and it'll be tested and be a driver for a new GEOS release. The downside is waiting for it, and what if it takes all summer? I may develop a Shapely-specific shim using the C++ classes in the meantime and then migrate over to the new GEOS code when it's ready.

[springmeyer]

I remember talking to you about this. Great to see it moving. C++ shim sounds like a good thing to try.

[sgillies]

Here's my post to the geos-devel list: http://lists.osgeo.org/pipermail/geos-devel/2014-February/006796.html.

@springmeyer I'm happy to have a real use for it at last. Pleiades never provided one.

[mabl]

This might also be interesting for my work. As part of my research I develop a CAD system for chip design (nano-optics). There, all geometries are typically exported to file formats which only allow coordinates on a discrete grid. Currently I sometimes have troubles that polygons do touch when rounded to the grid, but not in the float model. As far as I have understood, the issue could be solved with this proposal, correct?

[sgillies]

After more thought, I see a problem with the approach above: it's going to leave fixed-precision geoms twisting in the wind when the context exits. I doubt GEOS is going to be able to migrate data from fixed to float. The following may be more robust:

from shapely.geometry import Point
from shapely.precision import fixed_model

# make a fixed-precision model.
fixed = fixed_model(scale=1.0, offset=(0.0, 0.0))

# make a point snapped to a fixed-precision grid.
p = Point(0, 0, model=fixed)

The object would get a reference to the model to hold for its lifetime. How does that look? I'd like to get requirements into the GEOS wiki ASAP.

cc @Toblerity/owners @Toblerity/pushers @mabl @springmeyer

[sgillies]

Am going with the above.

[mabl]

Any progress on this front? Is there a student working on it?

[sgillies]

@mabl I haven't seen any recent discussion on geos-devel, so I think the GSoC work may have come apart.

[sgillies]

@mabl Also, I've got a new idea for doing the precision model. It's at http://trac.osgeo.org/geos/wiki/GSoC/CAPI_PrecisionModel, but I'll also copy it here.

Requirements for Shapely

Shapely (https://github.com/Toblerity/Shapely) is one of the major Python users of GEOS. I want to let programmers select and use a precision model like this:

from shapely.geometry import Point, factory
from shapely.precision import fixed_model

# make a geometry factory using a fixed-precision model.
fixed = factory(fixed_model(scale=1.0, offset=(0.0, 0.0)))

# make a point snapped to a fixed-precision grid.
p = Point(0, 0, factory=fixed)

# make a float-precision point.
q = Point(0, 0)

The above is a bit of a reversal from what I've previously written about my precision model requirements. A Python API like the one above is going to be very reliable and doesn't leave any geometry objects twisting in the wind when their precision model goes away: they have references to the model which keep it "alive". So, my main requirements for the precision model in the C API are now:

No global state

There shouldn't be a precision model in the global context. Or if there is, it shouldn't preclude precision models outside the global context.

As many precision models as I want

I expect any particular Shapely program to be using only 1-2 at a time, but let's not do anything to arbitrarily limit the number.

[jwass]

Will there also be a way to define a global default precision model? I could see many programs wanting to have all geometry objects created using the same precision model, but not wanting to have factory= everywhere. It would be nice to change it in one place without having to pass it around or have every module grab it from somewhere.

[sgillies]

@jwass globals are the devil, certainly so in a library. No :)

[sgillies]

Well, the Summer of Code project seems to have foundered. I'm going to push this one ahead on the roadmap. Will make a great rallying feature for the next major release.

[vectorien]

+1 for a fixed precision model

[mwtoews]

Here's some ctypes hackery to explore the GEOS precision model capabilities, based on some data from #821:

from shapely import geos, wkt
from ctypes import c_double, c_int, c_void_p

A = wkt.loads('POLYGON ((0 0, 0 55, 16.8 55, 14.8 43, 13 43.3, 5 0, 0 0))')
B = wkt.loads('POLYGON ((4 45, 20 45, 20 35, 4 35, 4 45))')
C = A.intersection(B)
print(C.wkt)
# POLYGON ((15.13333333333333 45, 14.8 43, 13 43.3, 11.46651270207852 35, 4 35, 4 45, 15.13333333333333 45))

# extern GEOSGeometry GEOS_DLL *GEOSGeom_setPrecision_r(
#     GEOSContextHandle_t handle, const GEOSGeometry *g, double gridSize, int flags);
geos._lgeos.GEOSGeom_setPrecision_r.argtypes = (c_void_p, c_void_p, c_double, c_int)
geos._lgeos.GEOSGeom_setPrecision_r.restype = c_void_p

# extern double GEOS_DLL GEOSGeom_getPrecision_r(
#     GEOSContextHandle_t handle, const GEOSGeometry *g)
geos._lgeos.GEOSGeom_getPrecision_r.argtypes = (c_void_p, c_void_p)
geos._lgeos.GEOSGeom_getPrecision_r.restype = c_double

ctx = geos.lgeos.geos_handle

Ap = geos._lgeos.GEOSGeom_setPrecision_r(ctx, A._geom, 0.1, 0)
Bp = geos._lgeos.GEOSGeom_setPrecision_r(ctx, B._geom, 0.1, 0)
Cp = geos.lgeos.GEOSIntersection(Ap, Bp)
print(geos._lgeos.GEOSGeom_getPrecision_r(ctx, Cp))
# 0.1
print(geos.lgeos.GEOSGeomToWKT(Cp))
# POLYGON ((15.1 45.0, 14.8 43.0, 13.0 43.3, 11.5 35.0, 4.0 35.0, 4.0 45.0, 15.1 45.0))

however, I'm not convinced endorsing a fixed precision model will help resolve overlay issues like #821 seeks to get:

Apb = geos.lgeos.GEOSBoundary(Ap)

# using logic in question showing "wrong" output
Cp2 = geos.lgeos.GEOSIntersection(Apb, Cp)
print(geos.lgeos.GEOSGeomToWKT(Cp2))
# MULTILINESTRING ((14.8 43.0, 13.0 43.3), (13.0 43.3, 11.5 35.0))

# "correct" result, avoiding near-parallel lines
Cp3 = geos.lgeos.GEOSIntersection(Apb, Bp)
print(geos.lgeos.GEOSGeomToWKT(Cp3))
# LINESTRING (15.1 45.0, 14.8 43.0, 13.0 43.3, 11.5 35.0)

[mwtoews]

On the flip-side of a good example where a precision model helps, borrowing from the LibGEOS.jl folks:

from shapely.geometry import Polygon

A = Polygon([(75.9, 30.7), (79.9, 30.7), (77.9, 28.7), (75.9, 30.7)])
B = Polygon([(81.0, 26.8), (76.0, 26.8), (81.0, 31.8), (81.0, 26.8)])

A.intersects(B)  # True
A.touches(B)  # False

Ap = geos._lgeos.GEOSGeom_setPrecision_r(ctx, A._geom, 1.0, 0)
Bp = geos._lgeos.GEOSGeom_setPrecision_r(ctx, B._geom, 1.0, 0)

bool(geos.lgeos.GEOSIntersects(Ap, Bp))  # True
bool(geos.lgeos.GEOSTouches(Ap, Bp))  # True