Opened 8 years ago

Closed 4 years ago

#717 closed defect (wontfix)

LINEARRING intersection returns different answers depending on where 'join' is

Reported by: messypup Owned by: geos-devel@…
Priority: minor Milestone: 3.6.3
Component: Default Version: 3.4.2
Severity: Annoyance Keywords: GEOSIntersection_r LINEARRING
Cc: messypup

Description

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.

Change History (5)

comment:1 by messypup, 8 years ago

Cc: messypup added

comment:2 by strk, 6 years ago

Milestone: 3.4.33.6.1

Ticket retargeted after milestone deleted

comment:3 by strk, 5 years ago

Milestone: 3.6.13.6.2

Ticket retargeted after milestone closed

comment:4 by strk, 5 years ago

Milestone: 3.6.23.6.3

Ticket retargeted after milestone closed

comment:5 by robe, 4 years ago

Resolution: wontfix
Status: newclosed

this is too old. Anyone feel free to reopen if you plan to fix.

Note: See TracTickets for help on using tickets.