Hi all,
I'm writing because I have a problem with the GEOSIntersection function.
I'm using this function as follows:
1. You should always provide compilable code, if at all possible, when
asking for help.
Follows the code I'm using.
// Init GEOS
GEOSMessageHandler error_function, notice_function;
initGEOS(notice_function, error_function);
GEOSCoordSeq coordseq = NULL;
GEOSGeom area_1 = NULL, area_2 = NULL, intersection = NULL;
coordseq = (GEOSCoordSeq) GEOSCoordSeq_create(5, 2); //5 points
bi-dimensional
GEOSCoordSeq_setX(coordseq, 0, inputdatat1->Lon_UL); //upper left
GEOSCoordSeq_setY(coordseq, 0, inputdatat1->Lat_UL);
GEOSCoordSeq_setX(coordseq, 1, inputdatat1->Lon_LR); //upper right
GEOSCoordSeq_setY(coordseq, 1, inputdatat1->Lat_UL);
GEOSCoordSeq_setX(coordseq, 2, inputdatat1->Lon_LR); //lower right
GEOSCoordSeq_setY(coordseq, 2, inputdatat1->Lat_LR);
GEOSCoordSeq_setX(coordseq, 3, inputdatat1->Lon_UL); //lower left
GEOSCoordSeq_setY(coordseq, 3, inputdatat1->Lat_LR);
GEOSCoordSeq_setX(coordseq, 4, inputdatat1->Lon_UL); //upper left
GEOSCoordSeq_setY(coordseq, 4, inputdatat1->Lat_UL);
area_1 = GEOSGeom_createLinearRing(coordseq);
if((GEOSisEmpty(area_1) != 0) || (GEOSisValid(area_1) != 1)) {
printf("No valid intersection found.\n");
exit(2); //invalid input parameter
}
GEOSCoordSeq_setX(coordseq, 0, inputdatat2->Lon_UL); //upper left
GEOSCoordSeq_setY(coordseq, 0, inputdatat2->Lat_UL);
GEOSCoordSeq_setX(coordseq, 1, inputdatat2->Lon_LR); //upper right
GEOSCoordSeq_setY(coordseq, 1, inputdatat2->Lat_UL);
GEOSCoordSeq_setX(coordseq, 2, inputdatat2->Lon_LR); //lower right
GEOSCoordSeq_setY(coordseq, 2, inputdatat2->Lat_LR);
GEOSCoordSeq_setX(coordseq, 3, inputdatat2->Lon_UL); //lower left
GEOSCoordSeq_setY(coordseq, 3, inputdatat2->Lat_LR);
GEOSCoordSeq_setX(coordseq, 4, inputdatat2->Lon_UL); //upper left
GEOSCoordSeq_setY(coordseq, 4, inputdatat2->Lat_UL);
area_2 = GEOSGeom_createLinearRing(coordseq);
if((GEOSisEmpty(area_2) != 0) || (GEOSisValid(area_2) != 1)) {
printf("No valid intersection found.\n");
exit(2); //invalid input parameter
}
intersection = GEOSIntersection(area_1, area_2);
if((GEOSisEmpty(intersection) != 0) || (GEOSisValid(intersection) !=
1)) {
printf("No valid intersection found.\n");
exit(2); //invalid input parameter
}
//Getting coords for the vertex
unsigned int num;
double xPoints[4];
double yPoints[4];
double *int_Lon_UL = (double*)malloc(sizeof(double)*1);
double *int_Lat_UL = (double*)malloc(sizeof(double)*1);
double *int_Lon_LR = (double*)malloc(sizeof(double)*1);
double *int_Lat_LR = (double*)malloc(sizeof(double)*1);
GEOSGeom geom;
num = GEOSGetNumGeometries(intersection);
printf("Geometries: %d\n",num);
GEOSCoordSeq_destroy(coordseq);
coordseq = (GEOSCoordSeq) GEOSCoordSeq_create(2, 2); //2 points
bi-dimensional
for(i=0; i < num; i++) {
geom = (GEOSGeom) GEOSGetGeometryN(intersection, i);
coordseq = (GEOSCoordSeq) GEOSGeom_getCoordSeq(geom);
GEOSCoordSeq_getX(coordseq, 0, &xPoints[i]);
GEOSCoordSeq_getY(coordseq, 0, &yPoints[i]);
}
//Filling the values in the right place
*int_Lon_UL = xPoints[0];
*int_Lat_UL = yPoints[0];
*int_Lon_LR = xPoints[2];
*int_Lat_LR = yPoints[2];
// Finalizzo GEOS
finishGEOS();
printf("First Area:\t%.2lf %.2lf,%.2lf %.2lf,%.2lf %.2lf,%.2lf
%.2lf\n",inputdatat1->Lat_UL, inputdatat1->Lon_UL, inputdatat1->Lat_UL,
inputdatat1->Lon_LR, inputdatat1->Lat_LR, inputdatat1->Lon_LR,
inputdatat1->Lat_LR, inputdatat1->Lon_UL);
printf("Second Area:\t%.2lf %.2lf,%.2lf %.2lf,%.2lf %.2lf,%.2lf
%.2lf\n",inputdatat2->Lat_UL, inputdatat2->Lon_UL, inputdatat2->Lat_UL,
inputdatat2->Lon_LR, inputdatat2->Lat_LR, inputdatat2->Lon_LR,
inputdatat2->Lat_LR, inputdatat2->Lon_UL);
printf("intersection:\t%.2lf %.2lf,%.2lf %.2lf,%.2lf %.2lf,%.2lf
%.2lf\n",*int_Lat_UL, *int_Lon_UL, *int_Lat_UL, *int_Lon_LR,
*int_Lat_LR, *int_Lon_LR, *int_Lat_LR, *int_Lon_UL);
And the output is:
First Area: 42.46 -131.80,42.46 -112.91,21.96 -112.91,21.96 -131.80
Second Area: 43.22 -125.52,43.22 -106.47,22.71 -106.47,22.71 -125.52
intersection: 43.22 -125.52,43.22 -106.47,22.71 -106.47,22.71 -125.52
As you can see the intersection's vertex are the same of the second area.
2. I think GEOSGeom_createLinearRing does not copy the input argument
--
Giacomo
_______________________________________________
geos-devel mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/geos-devel
