precision model or setPrecision

[jaakkor2]

This example is problematic with float precision

using LibGEOS
p1 = LibGEOS.Polygon([[[75.9, 30.7], [79.9, 30.7], [77.9, 28.7], [75.9, 30.7]]]);
p2 = LibGEOS.Polygon([[[81.0, 26.8], [76.0, 26.8], [81.0, 31.8], [81.0, 26.8]]]);

julia> intersects(p1,p2) # should be true
true
julia> touches(p1,p2) # should be true
false
julia> is12 = intersection(p1,p2);
julia> area(is12)
3.552713678800501e-15
julia> writegeom(is12)
"POLYGON ((79.9000000000000057 30.6999999999999993, 77.9000000000000057 28.6999999999999993, 77.9000000000000057 28.7000000000000028, 79.9000000000000057 30.6999999999999993))"
julia> LibGEOS.GEOSGeom_getPrecision(p1.ptr)
-1.0
julia> LibGEOS.GEOSGeom_getPrecision(p2.ptr)
-1.0

using Plots
plot([p1,p2])

I cannot figure out usage of setPrecision that does not seem to do anything, this is what I tried

LibGEOS.GEOSGeom_setPrecision(p2.ptr,0.1,0);

There seems to be a precision model in GEOS
https://github.com/libgeos/geos/blob/master/src/geom/PrecisionModel.cpp
but is this exposed in LibGEOS.jl ?

[visr]

Hmm I'm not sure how this is supposed to work. Is there an example somewhere using the C API that we can learn from? Unfortunately https://geos.osgeo.org/doxygen/geos__c_8h_source.html is not very helpful here.

With your example:

using LibGEOS
p1 = LibGEOS.Polygon([[[75.9, 30.7], [79.9, 30.7], [77.9, 28.7], [75.9, 30.7]]])
p2 = LibGEOS.Polygon([[[81.0, 26.8], [76.0, 26.8], [81.0, 31.8], [81.0, 26.8]]])
LibGEOS.GEOSGeom_setPrecision(p2.ptr, 0.1, 0)  # returns C_NULL, so something went wrong

I tried to use this method (is the other one deprecated already?), and this did give a valid pointer, but did not change the results:

p1_grid = Polygon(LibGEOS.GEOSGeom_setPrecision_r(LibGEOS._context.ptr, p1.ptr, 0.1, 0))
p2_grid = Polygon(LibGEOS.GEOSGeom_setPrecision_r(LibGEOS._context.ptr, p2.ptr, 0.1, 0))
intersects(p1_grid, p2_grid)
touches(p1_grid, p2_grid)

Looking at https://shapely.readthedocs.io/en/stable/manual.html#object.touches it says:

Returns True if the objects have at least one point in common and their interiors do not intersect with any part of the other.

So is false not the correct answer here? What does shapely return here? That should be the same.

[jaakkor2]

Compare with a bit more fortunate geometry

using LibGEOS
p1m = LibGEOS.Polygon([[[76.0, 30.0], [80.0, 30.0], [78.0, 28.0], [76.0, 30.0]]]);
p2m = LibGEOS.Polygon([[[81.0, 26.0], [76.0, 26.0], [81.0, 31.0], [81.0, 26.0]]]);
intersects(p1m, p2m) # true
touches(p1m, p2m) # **true**
is12m = intersection(p1m, p2m)
area(is12m) # 0.0

Here the intersection is a LineString(Ptr{Nothing} @0x000000002f450be0), what I would expect, and touches returns true.

In the grid-example, intersection(p1_grid, p2_grid) gives LineString(Ptr{Nothing} @0x000000002f452a60) like what I would expect. Only the touches is still giving false. I think I can live with LibGEOS.GEOSGeom_setPrecision_r ! So thank you!

To compare with shapely

from shapely.geometry import Polygon
p1 = Polygon([(75.9, 30.7), (79.9, 30.7), (77.9, 28.7), (75.9, 30.7)])
p2 = Polygon([(81.0, 26.8), (76.0, 26.8), (81.0, 31.8), (81.0, 26.8)])
p1.intersects(p2) # True
p1.touches(p2) # False
is12 = p1.intersection(p2)
is12.area

So shapely is saying p1.touches(p2) is False, and the intersection is a <shapely.geometry.polygon.Polygon object at 0x000002BE8D7A0668>. Just like the LibGEOS.jl-version.

There is discussion about GEOS precision models for shapely shapely/shapely#110 but maybe best not to hold one's breath.

[visr]

Ok then at least we are not doing something wrong in this package.

Of course it would be nice to create an easier interface for this functionality, right now it's still hard to use with passing the pointers directly. Also we should probably make the _r functions the default and pass the default context as a keyword argument.

Btw it might be nice to know there is currently also work going on to set up this library: https://github.com/pygeos/pygeos