[postgis-users] ST_Intersection/GEOS: should it return POLYGON EMPTY or GEOMETRYCOLLECTION EMPTY ?

Jakub Wartak Jakub.Wartak at tomtom.com
Fri Apr 30 00:17:20 PDT 2021


Hello postgis community,

The https://postgis.net/docs/ST_Intersection.html states that "If the geometries do not share any space (are disjoint), then an empty geometry >> collection << is returned." however we are observing that the following SQL reproducer:

select st_astext(st_intersection(st_geomfromtext('POLYGON ((3.8553818 51.0083094, 3.8560454 51.0083533, 3.8565504 51.008428, 3.8568357 51.0084649, 3.8572921 51.0084757, 3.8575745 51.0084765, 3.8584016 51.0084513, 3.8585614 51.0084387, 3.8587326 51.0084305, 3.8589394 51.0084251, 3.8589694 51.008417, 3.858975 51.0084089, 3.8589736 51.0083936, 3.8589466 51.008381, 3.8588924 51.0083738, 3.8587939 51.0083711, 3.8577014 51.0084073, 3.8576829 51.0083235, 3.8575702 51.0084054, 3.8572607 51.0084063, 3.8569797 51.0083902, 3.8566417 51.0083587, 3.8566161 51.0083019, 3.8566161 51.0082632, 3.8565832 51.0082461, 3.8565576 51.0082461, 3.8565376 51.0082749, 3.8565361 51.0083307, 3.8562823 51.0083092, 3.856046 51.008303, 3.8555091 51.0082619, 3.8555905 51.0080904, 3.8555001 51.0080694, 3.8553818 51.0083094))'), 'POLYGON ((3.8232422 50.9985352, 3.8232422 51.0095215, 3.8452149 51.0095215, 3.8452149 50.9985352, 3.8232422 50.9985352))'))

returns "POLYGON EMPTY" in case of using newer GEOS library, however in earlier GEOS releases "GEOMETRYCOLLECTION EMPTY" was returned instead. Is this bug in the code or the PostGIS documentation is misleading or maybe it's documented change in behavior of GEOS 3.8+? (I was searching in Trac, but I couldn't find anything specific). It's not big deal as we can workaround that by using ST_IsEmpty() I think, but still thanks for any information which one is proper one.

It is PostGIS difference, however here I've GEOS reproducer:

## geos39-3.9.0-1.rhel7.x86_64
[root at box ~]# c++ -I/usr/geos39/include -L/usr/geos39/lib64 -lgeos-3.9.0 -lgeos_c t.cpp -o t [root at box ~]# ./t
Intersection: POLYGON EMPTY

## geos38-3.8.1-2.rhel7.x86_64
[root at box ~]# c++ -I/usr/geos38/include -L/usr/geos38/lib64 -lgeos -lgeos_c t.cpp -o t [root at box ~]# ./t
Intersection: POLYGON EMPTY

## geos37-3.7.2-2.rhel7.x86_64
[root at box ~]# c++ -I/usr/geos37/include -L/usr/geos37/lib64 -lgeos -lgeos_c t.cpp -o t [root at box ~]# ./t
Intersection: GEOMETRYCOLLECTION EMPTY

[root at box ~]# cat t.cpp
#include <geos_c.h>
#include <string>
#include <limits>
#include <cstdarg>
#include <cstdio>

using namespace std;

static void printGEOSNotice( const char *fmt, ... ) {
        va_list ap;
        char buffer[1024];
        va_start(ap, fmt);
        vsnprintf(buffer, sizeof buffer, fmt, ap);
        printf("%s\n", buffer);
        va_end(ap);
        exit(1);
}

int main() {
        initGEOS(&printGEOSNotice, &printGEOSNotice);
        GEOSWKTReader *reader = GEOSWKTReader_create();
        GEOSGeometry *geomA = GEOSWKTReader_read(reader, "POLYGON ((3.8553818 51.0083094, 3.8560454 51.0083533, 3.8565504 51.008428, 3.8568357 51.0084649, 3.8572921 51.0084757, 3.8575745 51.0084765, 3.8584016 51.0084513, 3.8585614 51.0084387, 3.8587326 51.0084305, 3.8589394 51.0084251, 3.8589694 51.008417, 3.858975 51.0084089, 3.8589736 51.0083936, 3.8589466 51.008381, 3.8588924 51.0083738, 3.8587939 51.0083711, 3.8577014 51.0084073, 3.8576829 51.0083235, 3.8575702 51.0084054, 3.8572607 51.0084063, 3.8569797 51.0083902, 3.8566417 51.0083587, 3.8566161 51.0083019, 3.8566161 51.0082632, 3.8565832 51.0082461, 3.8565576 51.0082461, 3.8565376 51.0082749, 3.8565361 51.0083307, 3.8562823 51.0083092, 3.856046 51.008303, 3.8555091 51.0082619, 3.8555905 51.0080904, 3.8555001 51.0080694, 3.8553818 51.0083094))");
        GEOSGeometry *geomB = GEOSWKTReader_read(reader, "POLYGON ((3.8232422 50.9985352, 3.8232422 51.0095215, 3.8452149 51.0095215, 3.8452149 50.9985352, 3.8232422 50.9985352))");
        GEOSGeometry *geomResult = GEOSIntersection(geomA, geomB);
        char *ptr = GEOSGeomToWKT(geomResult);
        printf("Intersection: %s\n", ptr);
        return 0;
}
[root at box ~]#

-Jakub Wartak.


More information about the postgis-users mailing list