[geos-devel] GEOSIntersection returns the input
Giacomo Piva
piva at meeo.it
Mon Aug 24 07:56:36 EDT 2009
Hi all,
I'm writing because I have a problem with the GEOSIntersection function.
I'm using this function as follows:
GEOSCoordSeq coordseq = NULL, coordseq_I = NULL;
GEOSGeom area_1 = NULL, area_2 = NULL, intersection = NULL;
coordseq = (GEOSCoordSeq) GEOSCoordSeq_create(5, 2); //5 coords to
close the ring
GEOSCoordSeq_setX(coordseq, 0, ...);
GEOSCoordSeq_setY(coordseq, 0, ...);
[...] //Fill the Seq with the coords value for the first ring
area_1 = GEOSGeom_createLinearRing(coordseq);
//Re-fill the coords for the second ring ...
GEOSCoordSeq_setX(coordseq, 0, ...);
GEOSCoordSeq_setY(coordseq, 0, ...);
[...]
area_2 = GEOSGeom_createLinearRing(coordseq);
intersection = GEOSIntersection(area_1, area_2);
//checks on valid and non empty geometry returns no error...
Now I get the coordinates of the intersection in this way:
GEOSGeom geom;
num = GEOSGetNumGeometries(intersection); //num is 4, because there is
4 lines
for(i=0; i < num; i++) {
geom = (GEOSGeom) GEOSGetGeometryN(intersection, i); //Here I
get i-th line...
coordseq_I = (GEOSCoordSeq) GEOSCoordSeq_create(2, 2); //each
line has 2 point of 2 dimensions
coordseq_I = (GEOSCoordSeq) GEOSGeom_getCoordSeq(geom);
GEOSCoordSeq_getX(coordseq_I, 0, alpha); //the first point is
in 0-th position
GEOSCoordSeq_getY(coordseq_I, 0, beta);
printf("[%d] %.2lf %.2lf\n",i,*beta, *alpha);
GEOSCoordSeq_getX(coordseq_I, 1, alpha); //the second in the first
GEOSCoordSeq_getY(coordseq_I, 1, beta);
printf("[%d] %.2lf %.2lf\n\n",i,*beta, *alpha);
}
Now, if I print the coordinates of each point of the 4 geometries
(lines) composing the intersection ring, I get the same value of the
"area_2" ring.
Second Area:
[A] 43.22 -125.52
[B] 43.22 -106.47
[C] 22.71 -106.47
[D] 22.71 -125.52
Geometries: 4
[A] 43.22 -125.52
[B] 43.22 -106.47
[B] 43.22 -106.47
[C] 22.71 -106.47
[C] 22.71 -106.47
[D] 22.71 -125.52
[D] 22.71 -125.52
[A] 43.22 -125.52
Really I can't figure out what is wrong.
Someone can help me?
Thanks.
--
Giacomo Piva
More information about the geos-devel
mailing list