[geos-devel] Using GEOS intersection result

[Maude Giasson]

Hi all,
I'm doing a small application using GEOS.

Given two simple linearRings L1 and L2, my application create a new
linearRing L3 which is similar to L1 but contains, in addition, new vertex
at the points where L1 and L2 cross.

If L1 and L2 cross on a edge of L1, no vertex is added.

If L1 and L2 share a common lineSegment, vertex are added to L3 at the
endPoints of the segment. (if they are not already vertex on L3).

In consequence, every edges of L3 should be within L1.

1- Does a similar function already exists in GEOS (I haven't find any) ?

2- Later in my application, I look if the edges of L3 are within L1 or not
(they should be within L1). It does not work pretty well. The problem seems
to come from GEOS which is not always consistent in it's answers. In fact, I
think the real problem is from my side because I do not know what I can
expect from GEOS answers.
I've made here a little example showing some GEOS results that I don't
understand. Can someone please tell me what is wrong in that example?
Understanding this particular example will help me solving problems in my
application.

In that example, GEOS says that the LineString AB intersects the LineString
XCDY at the point E.
However,  it says that LineString XCDY does not intersect point E.
In consequence, the line segment ED is not within XCDY.

//Coordinate :
geos::Coordinate cA(3.8727033436298370, 34.142325401306152);
geos::Coordinate cB(22.790718823671341, 34.684023857116699);
geos::Coordinate cC(15.851628303527832, 42.153842926025391);
geos::Coordinate cD(4.8595705032348633, 28.923076629638672);
// E : 9.3254120334165922, 34.298458262089483
geos::Coordinate cX(12,55);
geos::Coordinate cY(4,28);
geos::PrecisionModel g_precisionModel;
int g_SRID = 1;
//-- ls AB construction
geos::BasicCoordinateList* bclAB = new geos::BasicCoordinateList();
bclAB->add(cA);
bclAB->add(cB);
geos::LineString* lsAB = new geos::LineString(bclAB, &g_precisionModel,
g_SRID);
//-- ls XCDY construction
geos::BasicCoordinateList* bclXCDY = new geos::BasicCoordinateList();
bclXCDY->add(cX);
bclXCDY->add(cC);
bclXCDY->add(cD);
bclXCDY->add(cY);
geos::LineString* lsXCDY = new geos::LineString(bclXCDY, &g_precisionModel,
g_SRID);

//-- intersection between the 2 lineStrings
geos::Geometry* gABinterXCDY = lsAB->intersection(lsXCDY);

//The intersection between lineStrings AB and XCDY occurs at the single
point of coordinate E.
//E is on AB and on CD. Let's verify that XCDY contains line segment
(lineString) ED.

//-- ls ED construction
geos::BasicCoordinateList* bclED = new geos::BasicCoordinateList();
bclED->add( geos::Coordinate(*(gABinterXCDY->getCoordinate())) );
bclED->add(cD);
geos::LineString* lsED = new geos::LineString(bclED, &g_precisionModel,
g_SRID);

//-- lsXCDY should contain lsED but it does not; why ?
Booleen test3 = lsXCDY->contains(lsED); // why is it false

// gABinterXCDY = lsAB inter lsXYCD
// but lsXCDY does not intersect gABinterXCDY.
// Why does it occur ?
Booleen test1 = lsXCDY->intersects(gABinterXCDY); // why is it false ?


Thanks
Maude Giasson