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

[Martin Davis]

This is an intentional change in the semantics of the overlay operations.
An empty result is now returned as an empty geometry of appropriate type.
This is to make it easier to store and use the result in further
computation.  So yes, the correct technique to check for this is to use
ST_IsEmpty.

Thanks for the alert about the manual - I will fix that.

On Fri, Apr 30, 2021 at 12:17 AM Jakub Wartak <Jakub.Wartak at tomtom.com>
wrote:

> 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.
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/postgis-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20210430/40c09dbc/attachment.html>