[geos-devel] [GEOS] #717: LINEARRING intersection returns different answers depending on where 'join' is

[GEOS]

#717: LINEARRING intersection returns different answers depending on where 'join'
is
-----------------------+----------------------------------------------------
 Reporter:  messypup   |       Owner:  geos-devel@…                 
     Type:  defect     |      Status:  new                          
 Priority:  minor      |   Milestone:  3.4.3                        
Component:  Default    |     Version:  3.4.2                        
 Severity:  Annoyance  |    Keywords:  GEOSIntersection_r LINEARRING
-----------------------+----------------------------------------------------
 When performing an intersection between a POLYGON and a LINEARRING if the
 first point of the LINEARRING is inside the polygon then what should be
 one LINESTRING returned is split into 2 on that point. Where as if the
 LINEARRING is constructed with the same points but the first point
 ('join') outside the POLYGON then (correctly) only one LINESTRING is
 returned.

 Example:

 #include <stdio.h>
 #include <geos_c.h>

 GEOSContextHandle_t handle;
 GEOSWKTWriter *writer;
 void printGeom(GEOSGeometry*, char*);

 int main(int argc, char **argv) {
     GEOSGeometry *ring1, *ring2, *box,
                  *geom1, *geom2;

     handle = initGEOS_r(NULL, NULL);
     writer = GEOSWKTWriter_create_r(handle);
     GEOSWKTWriter_setRoundingPrecision_r(handle, writer, 1);

     ring1 = GEOSGeomFromWKT_r(handle, "LINEARRING (2 2, 4 2, 4 4, 2 4, 2
 2)");
     ring2 = GEOSGeomFromWKT_r(handle, "LINEARRING (2 4, 2 2, 4 2, 4 4, 2
 4)");
     box = GEOSGeomFromWKT_r(handle, "POLYGON ((3 1, 3 3, 1 3, 1 1, 3
 1))");

     geom1 = GEOSIntersection_r(handle, ring1, box);
     geom2 = GEOSIntersection_r(handle, ring2, box);

     printGeom(ring1, "ring1");
     printGeom(ring2, "ring2");
     printGeom(box, "box");
     printGeom(geom1, "Intersection of ring1 and box");
     printGeom(geom2, "Intersection of ring2 and box");
 }

 void printGeom(GEOSGeometry *geom, char *name) {
     printf("%s is: %s\n", name, GEOSWKTWriter_write_r(handle, writer,
 geom));
 }

 Output:

 ring1 is: LINEARRING (2.0 2.0, 4.0 2.0, 4.0 4.0, 2.0 4.0, 2.0 2.0)
 ring2 is: LINEARRING (2.0 4.0, 2.0 2.0, 4.0 2.0, 4.0 4.0, 2.0 4.0)
 box is: POLYGON ((3.0 1.0, 3.0 3.0, 1.0 3.0, 1.0 1.0, 3.0 1.0))
 Intersection of ring1 and box is: MULTILINESTRING ((2.0 2.0, 3.0 2.0),
 (2.0 3.0, 2.0 2.0))
 Intersection of ring2 and box is: LINESTRING (2.0 3.0, 2.0 2.0, 3.0 2.0)

 I am using geos through the Python package Shapely. I originally submitted
 this as a Shapely issue but the Shapely developer showed the result is
 coming from geos itself. The example above is his.

 (excuse terminology here because I only know the Shapley interface) It is
 also not a trivial matter to work around this by passing the result to
 linemerge because linemerge will not accept a LINESTRING or POINT which
 this intersection could return.

-- 
Ticket URL: <http://trac.osgeo.org/geos/ticket/717>
GEOS <http://trac.osgeo.org/geos>
GEOS (Geometry Engine - Open Source) is a C++ port of the Java Topology Suite (JTS).