[geos-devel] GEOSIntersection returns the input

[Giacomo Piva]

>> 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