Empty geometries not all created the same way

[mwtoews]

Currently there are two ways to create an empty geometry (e.g. Polygon):

from shapely import wkt
from shapely.geometry import mapping, Polygon

# Method 1: using an empty constructor
g1 = Polygon()
print('{}, {}, {}'.format(mapping(g1), g1.geom_type, g1.wkt))
# {'type': 'Polygon', 'coordinates': ()}, GeometryCollection, GEOMETRYCOLLECTION EMPTY

# Method 2: loading empty WKT
g2 = wkt.loads('POLYGON EMPTY')
print('{}, {}, {}'.format(mapping(g2), g2.geom_type, g2.wkt))
# {'type': 'Polygon', 'coordinates': ()}, Polygon, POLYGON EMPTY

(all other geometry types behave the same way).

In both cases, the Python objects know they are Polygons based on the Python object type. However the GEOS objects (referenced by ._geom) are either GeometryCollection or Polygon, based on how the geometries were constructed. This seems inconsistent to me, as I'd expect the two methods to be identical.

[snorfalorpagus]

I have come across this issue before, see: #338 #352

@snorfalorpagus
I'm leaning towards leaving things the way they are, because of the issue with there being no representation of POINT EMPTY in WKB, i.e. to keep consistency with PostGIS.

Although it's possible to have Point().wkt return POINT EMPTY it breaks Point().wkb.

[mwtoews]

Thanks @snorfalorpagus, this empty discussion did seem familiar. I think PR #352 should be reopened to make shapely more like modern PostGIS.

It was correct that older PostGIS (1.x?) treated all empty geometries as GEOMETRYCOLLECTION EMPTY, however modern versions have many more empty types.

From PostGIS 2.4 there are 60 distinct empty geometry type combinations, and that's distinct with respect to WKT and WKB:

SELECT ST_AsText(geom), geom::text
FROM (
  SELECT i, (postgis_typmod_type(i) || ' EMPTY')::geometry AS geom
  FROM generate_series(4, 63) AS i) f
ORDER BY i;

          st_astext          |                                    geom
-----------------------------+----------------------------------------------------------------------------
 POINT EMPTY                 | 0101000000000000000000F87F000000000000F87F
 POINT M EMPTY               | 0101000040000000000000F87F000000000000F87F000000000000F87F
 POINT Z EMPTY               | 0101000080000000000000F87F000000000000F87F000000000000F87F
 POINT ZM EMPTY              | 01010000C0000000000000F87F000000000000F87F000000000000F87F000000000000F87F
 LINESTRING EMPTY            | 010200000000000000
 LINESTRING M EMPTY          | 010200004000000000
 LINESTRING Z EMPTY          | 010200008000000000
 LINESTRING ZM EMPTY         | 01020000C000000000
 POLYGON EMPTY               | 010300000000000000
 POLYGON M EMPTY             | 010300004000000000
 POLYGON Z EMPTY             | 010300008000000000
 POLYGON ZM EMPTY            | 01030000C000000000
 MULTIPOINT EMPTY            | 010400000000000000
 MULTIPOINT M EMPTY          | 010400004000000000
 MULTIPOINT Z EMPTY          | 010400008000000000
 MULTIPOINT ZM EMPTY         | 01040000C000000000
 MULTILINESTRING EMPTY       | 010500000000000000
 MULTILINESTRING M EMPTY     | 010500004000000000
 MULTILINESTRING Z EMPTY     | 010500008000000000
 MULTILINESTRING ZM EMPTY    | 01050000C000000000
 MULTIPOLYGON EMPTY          | 010600000000000000
 MULTIPOLYGON M EMPTY        | 010600004000000000
 MULTIPOLYGON Z EMPTY        | 010600008000000000
 MULTIPOLYGON ZM EMPTY       | 01060000C000000000
 GEOMETRYCOLLECTION EMPTY    | 010700000000000000
 GEOMETRYCOLLECTION M EMPTY  | 010700004000000000
 GEOMETRYCOLLECTION Z EMPTY  | 010700008000000000
 GEOMETRYCOLLECTION ZM EMPTY | 01070000C000000000
 CIRCULARSTRING EMPTY        | 010800000000000000
 CIRCULARSTRING M EMPTY      | 010800004000000000
 CIRCULARSTRING Z EMPTY      | 010800008000000000
 CIRCULARSTRING ZM EMPTY     | 01080000C000000000
 COMPOUNDCURVE EMPTY         | 010900000000000000
 COMPOUNDCURVE M EMPTY       | 010900004000000000
 COMPOUNDCURVE Z EMPTY       | 010900008000000000
 COMPOUNDCURVE ZM EMPTY      | 01090000C000000000
 CURVEPOLYGON EMPTY          | 010A00000000000000
 CURVEPOLYGON M EMPTY        | 010A00004000000000
 CURVEPOLYGON Z EMPTY        | 010A00008000000000
 CURVEPOLYGON ZM EMPTY       | 010A0000C000000000
 MULTICURVE EMPTY            | 010B00000000000000
 MULTICURVE M EMPTY          | 010B00004000000000
 MULTICURVE Z EMPTY          | 010B00008000000000
 MULTICURVE ZM EMPTY         | 010B0000C000000000
 MULTISURFACE EMPTY          | 010C00000000000000
 MULTISURFACE M EMPTY        | 010C00004000000000
 MULTISURFACE Z EMPTY        | 010C00008000000000
 MULTISURFACE ZM EMPTY       | 010C0000C000000000
 POLYHEDRALSURFACE EMPTY     | 010F00000000000000
 POLYHEDRALSURFACE M EMPTY   | 010F00004000000000
 POLYHEDRALSURFACE Z EMPTY   | 010F00008000000000
 POLYHEDRALSURFACE ZM EMPTY  | 010F0000C000000000
 TRIANGLE EMPTY              | 011100000000000000
 TRIANGLE M EMPTY            | 011100004000000000
 TRIANGLE Z EMPTY            | 011100008000000000
 TRIANGLE ZM EMPTY           | 01110000C000000000
 TIN EMPTY                   | 011000000000000000
 TIN M EMPTY                 | 011000004000000000
 TIN Z EMPTY                 | 011000008000000000
 TIN ZM EMPTY                | 01100000C000000000
(60 rows)

(I'm not suggesting to implement all of these, just showing the scene report from modern PostGIS).

Where I'm going with this is, suppose you have a geopandas GeoSeries of geometries, which can be a sequence of geometries of any type. But if you wanted to check if they were all (e.g.) Polygon via the geom_type property, you'd get a mixture of Polygon and/or GeometryCollection, depending on how they were constructed. It would be much easier for the user to see a consistent type, regardless of the emptiness.

Point is a special case, and probably a separate discussion. GEOS does support reading a native POINT EMPTY, but it can't seem to form the WKB, as it spits out a log message and returns a zero length string:

from shapely import wkb
p = wkb.loads('0101000000000000000000F87F000000000000F87F', hex=True)
print(p.wkt)  # POINT EMPTY
print(repr(p.wkb_hex))  # ''
# logger writes: IllegalArgumentException: Empty Points cannot be represented in WKB

in this case, Shapely should just return the same WKB that PostGIS returns from the above table dump.

[snorfalorpagus]

This explanation of why representing POINT EMPTY in WKB is awkward from Paul Ramsey here: https://trac.osgeo.org/postgis/ticket/3031

Every type except point as a WKB representation that goes

endian[byte] + typenumber[integer] + numpoints[integer]

those are the patterns you're seeing for the other types. endian, type, and then a zero, indicating the object has zero elements or zero points
trouble is, the wkb for point is

endian[byte] + typenumber[integer] + x[double] + y[double]

there's nowhere to put in the information that the point is actually empty

It looks like the workaround was to represent an empty point as a point with NaN for both coordinates.

From the geopackage specification: http://www.geopackage.org/spec/

Well-Known Binary as defined in OGC 06-103r4 [9] does not provide a standardized encoding for an empty point set (i.e., 'Point Empty' in Well-Known Text). In GeoPackages these points SHALL be encoded as a Point where each coordinate value is set to an IEEE-754 quiet NaN value. GeoPackages SHALL use big endian 0x7ff8000000000000 or little endian 0x000000000000f87f as the binary encoding of the NaN values.

Implementation of POINT EMPTY as POINT (NaN NaN) in PostGIS: https://trac.osgeo.org/postgis/ticket/3181

I see no reason not to follow suit and return this new representation, instead of raising the existing exception. Do we need some special handling in the WKB writer for handling of SRID?

[jorisvandenbossche]

Similar example with Points:

In [9]: import shapely.geometry

In [10]: import shapely.wkt 

In [11]: empty1 = shapely.wkt.loads("POINT EMPTY") 

In [12]: type(empty1)   
Out[12]: shapely.geometry.point.Point

In [13]: empty1.geom_type 
Out[13]: 'Point'

In [14]: print(empty1)   
POINT EMPTY

In [17]: empty2 = shapely.geometry.Point() 

In [18]: type(empty2)  
Out[18]: shapely.geometry.point.Point

In [19]: empty2.geom_type 
Out[19]: 'GeometryCollection'

In [20]: print(empty2) 
GEOMETRYCOLLECTION EMPTY

where Point() creates an inconsistent Point object as its type is actually a GeometryCollection

[jorisvandenbossche]

I explored this a bit more, see https://gist.github.com/jorisvandenbossche/266aedbd7c3e07bc16b817f3f6c9f5d8

I also tried to look if GEOS actually supports more empty geometry types like PostGIS does (thinks like POLYGON Z EMPTY), but it seems that the extra dimension is lost when creating empty geometries (although this might also be a bug in the "has_z" implementation, might need to look in the GEOS code to be more certain about GEOS' support of this).

[karolzlot]

Another discrepancy:

geo_dict = {
    "ogc_fid": [1, 2, 3, 4, 5, 6],
    "abcdef": ["11b", "22b", "33b", "44b", "55b", "66b"],
    "geometry": [
        wkt.loads(f"POINT Z (30 10 5)"),
        wkt.loads(f"POLYGON Z EMPTY"),
        wkt.loads(f"LINESTRING Z EMPTY"),
        wkt.loads(f"POINT Z EMPTY"),
        wkt.loads(f"POLYGON Z EMPTY"),
        wkt.loads(f"LINESTRING Z EMPTY"),
    ],
}

gdf=gpd.GeoDataFrame(geo_dict,crs="EPSG:4326")

print(gdf)

gdf.to_postgis(
    table_name,
    engine,
    schema=schema,
)

Above works ok, but when I change any geometry above to MULTI (e.g. "MULTIPOINT Z EMPTY" instead of "POINT Z EMPTY") then I will get an error during to_postgis:

Column has Z dimension but geometry does not

[karolzlot]

Also, empty geometries have has_z always as False:

empty_geometry = wkt.loads(f"POINT Z EMPTY")
print(empty_geometry.has_z)
empty_geometry = wkt.loads(f"LINESTRING Z EMPTY")
print(empty_geometry.has_z)
empty_geometry = wkt.loads(f"POLYGON Z EMPTY")
print(empty_geometry.has_z)

empty_geometry = wkt.loads(f"MULTIPOLYGON Z EMPTY")
print(empty_geometry.has_z)
empty_geometry = wkt.loads(f"MULTIPOINT Z EMPTY")
print(empty_geometry.has_z)
empty_geometry = wkt.loads(f"MULTILINESTRING Z EMPTY")
print(empty_geometry.has_z)

Results in:

False
False
False
False
False
False

[jorisvandenbossche]

Also, empty geometries have has_z always as False:

Yes, that's a known limitation. GEOS doesn't really properly support empty 3D geometries. See also #1345 (comment)

In the documentation on the main branch we have a note about this in has_z:

shapely/shapely/predicates.py

Lines 39 to 43 in fc55006

	def has_z(geometry, **kwargs):
	"""Returns True if a geometry has a Z coordinate.
	Note that this function returns False if the (first) Z coordinate equals NaN or
	if the geometry is empty.

[jorisvandenbossche]

Small update on one of the mentioned issues here: the WKB serialization for POINT EMPTY now works and uses the format of PONIT (NaN NaN) (as was already done before by eg PostGIS as well). GEOS itself now does this since version 3.9.0, and with shapely 2.0 we also ensure consistent behaviour for this across all GEOS versions.

[theroggy]

Everything reported in this item has been resolved, except for Z values in multi-geometries (#742 (comment)), and that's an upstream issue in GEOS, tracked by this issue: libgeos/geos#888

To avoid confusion with the many old pieces of information in this issue, I'm closing this and opening a new one with the one problem that is left: #2390