[GH-3329] Preserve dimensions of empty Points, LineStrings, and Polygons

python/sedona/spark/utils/geometry_serde_general.py

@@ -36,6 +36,11 @@
 from shapely.wkb import dumps as wkb_dumps
 from shapely.wkt import loads as wkt_loads
 
+try:
+    from shapely import geos_version
+except ImportError:
+    from shapely.geos import geos_version
+
 CoordType = Union[
     Tuple[float, float], Tuple[float, float, float], Tuple[float, float, float, float]
 ]
@@ -85,6 +90,19 @@ def type_of(geom) -> int:
         else:
             raise ValueError(f"Invalid coordinate dimension: {geom._ndim}")
 
+    @staticmethod
+    def type_of_empty(geom) -> int:
+        # GEOS < 3.9 cannot write empty Points as WKB. Keep the fallback's
+        # existing XY encoding on those versions without calling the writer.
+        if isinstance(geom, Point) and geos_version < (3, 9, 0):
+            return CoordinateType.XY
+        # Shapely 1.x reports _ndim == 2 even for explicit XYZ empty geometries.
+        # Their WKB still records Z, so read its type flag instead.
+        wkb = wkb_dumps(geom)
+        byte_order = "<I" if wkb[0] else ">I"
+        geometry_type = struct.unpack_from(byte_order, wkb, 1)[0]
+        return CoordinateType.XYZ if geometry_type & 0x80000000 else CoordinateType.XY
+
     @staticmethod
     def bytes_per_coord(coord_type: int) -> int:
         return CoordinateType.BYTES_PER_COORDINATE[coord_type - 1]
@@ -163,6 +181,12 @@ def read_coordinate(self) -> CoordType:
         self.coords_offset += self.bytes_per_coord
         return coord
 
+    def read_empty(self, geometry_type: str) -> BaseGeometry:
+        # Constructors such as Point() create empty GeometryCollections in
+        # Shapely 1.x. WKT preserves both the primitive type and its dimension.
+        dimension = " Z" if self.coord_type == CoordinateType.XYZ else ""
+        return wkt_loads(f"{geometry_type}{dimension} EMPTY")
+
     def read_int(self) -> int:
         value = struct.unpack_from("i", self.buffer, self.ints_offset)[0]
         if value > len(self.buffer):
@@ -327,15 +351,14 @@ def serialize_point(geom: Point) -> bytes:
         coords = coords[0]
         return struct.pack(pack_format, preamble_byte, 0, 0, 0, 1, *coords)
     else:
-        return struct.pack("BBBBi", 18, 0, 0, 0, 0)
+        return generate_header_bytes(
+            GeometryTypeID.POINT, CoordinateType.type_of_empty(geom), 0
+        )
 
 
 def deserialize_point(geom_buffer: GeometryBuffer) -> Point:
     if geom_buffer.num_coords == 0:
-        # Here we don't call Point() directly since it would create an empty GeometryCollection
-        # in shapely 1.x. You'll find similar code for creating empty geometries in other
-        # deserialization functions.
-        return wkt_loads("POINT EMPTY")
+        return geom_buffer.read_empty("POINT")
     coord = geom_buffer.read_coordinate()
     return Point(coord)
 
@@ -385,12 +408,14 @@ def serialize_linestring(geom: LineString) -> bytes:
         )
         return header + array.array("d", [x for c in coords for x in c]).tobytes()
     else:
-        return generate_header_bytes(GeometryTypeID.LINESTRING, 1, 0)
+        return generate_header_bytes(
+            GeometryTypeID.LINESTRING, CoordinateType.type_of_empty(geom), 0
+        )
 
 
 def deserialize_linestring(geom_buffer: GeometryBuffer) -> LineString:
     if geom_buffer.num_coords == 0:
-        return wkt_loads("LINESTRING EMPTY")
+        return geom_buffer.read_empty("LINESTRING")
     coords = geom_buffer.read_coordinates(geom_buffer.num_coords)
     return LineString(coords)
 
@@ -437,7 +462,9 @@ def serialize_polygon(geom: Polygon) -> bytes:
     num_rings = struct.unpack_from(int_format, wkb_string, 5)[0]
 
     if num_rings == 0:
-        return generate_header_bytes(GeometryTypeID.POLYGON, CoordinateType.XY, 0)
+        return generate_header_bytes(
+            GeometryTypeID.POLYGON, CoordinateType.type_of_empty(geom), 0
+        )
 
     coord_bytes = b""
     ring_lengths = []
@@ -464,7 +491,7 @@ def serialize_polygon(geom: Polygon) -> bytes:
 
 def deserialize_polygon(geom_buffer: GeometryBuffer) -> Polygon:
     if geom_buffer.num_coords == 0:
-        return wkt_loads("POLYGON EMPTY")
+        return geom_buffer.read_empty("POLYGON")
     return geom_buffer.read_polygon()
 
 

python/src/geom_buf.c

@@ -147,7 +147,10 @@ static SedonaErrorCode copy_coord_seq_to_buffer(
 static SedonaErrorCode copy_buffer_to_coord_seq(
     GEOSContextHandle_t handle, double *buf, int num_coords, int has_z,
     int has_m, GEOSCoordSequence **p_coord_seq) {
-  if (dyn_GEOSCoordSeq_copyFromBuffer_r != NULL) {
+  /* Older GEOS versions infer XYZ for an empty buffer. Use the explicit
+   * dimension constructor below for empty XY/XYZ sequences. M layouts need
+   * copyFromBuffer, which preserves their dimensions on GEOS >= 3.12. */
+  if (dyn_GEOSCoordSeq_copyFromBuffer_r != NULL && (num_coords > 0 || has_m)) {
     /* fast path for libgeos >= 3.10.0 */
     GEOSCoordSequence *coord_seq = dyn_GEOSCoordSeq_copyFromBuffer_r(
         handle, buf, num_coords, has_z, has_m);

python/src/geomserde.c

@@ -65,9 +65,7 @@ static SedonaErrorCode sedona_deserialize_point(GEOSContextHandle_t handle,
                                                 CoordinateSequenceInfo *cs_info,
                                                 GEOSGeometry **p_geom) {
   GEOSGeometry *geom = NULL;
-  if (cs_info->num_coords == 0) {
-    geom = dyn_GEOSGeom_createEmptyPoint_r(handle);
-  } else if (cs_info->dims == 2) {
+  if (cs_info->num_coords > 0 && cs_info->dims == 2) {
     /* fast path for 2D points */
     double x = *geom_buf->buf_coord++;
     double y = *geom_buf->buf_coord++;
@@ -122,15 +120,9 @@ static SedonaErrorCode sedona_serialize_linestring(
 static SedonaErrorCode sedona_deserialize_linestring(
     GEOSContextHandle_t handle, int srid, GeomBuffer *geom_buf,
     CoordinateSequenceInfo *cs_info, GEOSGeometry **p_geom) {
-  if (cs_info->num_coords == 0) {
-    GEOSGeometry *geom = dyn_GEOSGeom_createEmptyLineString_r(handle);
-    if (geom == NULL) {
-      return SEDONA_GEOS_ERROR;
-    }
-    *p_geom = geom;
-    return SEDONA_SUCCESS;
-  }
-
+  /* Preserve the stored dimensions for empty LineStrings too. The default
+   * GEOS empty constructor can add or drop dimensions depending on the version.
+   */
   GEOSCoordSequence *coord_seq = NULL;
   SedonaErrorCode err =
       geom_buf_read_coords(geom_buf, handle, cs_info, &coord_seq);
@@ -185,8 +177,22 @@ static SedonaErrorCode sedona_deserialize_polygon(
     GEOSContextHandle_t handle, int srid, GeomBuffer *geom_buf,
     CoordinateSequenceInfo *cs_info, GEOSGeometry **p_geom) {
   if (cs_info->num_coords == 0) {
-    GEOSGeometry *geom = dyn_GEOSGeom_createEmptyPolygon_r(handle);
+    /* An explicit empty shell preserves the stored Z/M layout. The default
+     * empty polygon constructor always creates an XY polygon. */
+    GEOSCoordSequence *coord_seq = NULL;
+    SedonaErrorCode err =
+        geom_buf_read_coords(geom_buf, handle, cs_info, &coord_seq);
+    if (err != SEDONA_SUCCESS) {
+      return err;
+    }
+    GEOSGeometry *shell = dyn_GEOSGeom_createLinearRing_r(handle, coord_seq);
+    if (shell == NULL) {
+      dyn_GEOSCoordSeq_destroy_r(handle, coord_seq);
+      return SEDONA_GEOS_ERROR;
+    }
+    GEOSGeometry *geom = dyn_GEOSGeom_createPolygon_r(handle, shell, NULL, 0);
     if (geom == NULL) {
+      dyn_GEOSGeom_destroy_r(handle, shell);
       return SEDONA_GEOS_ERROR;
     }
     *p_geom = geom;

python/tests/geopandas/test_match_geopandas_series.py

@@ -153,6 +153,7 @@ def setup_method(self):
             ),
             # Collection predicates must retain collection semantics with one member.
             GeometryCollection([LineString([(0, 0), (1, 1), (0, 0)])]),
+            GeometryCollection([Point(), LineString(), Polygon()]),
         ]
 
         self.geoms = [
@@ -168,7 +169,8 @@ def setup_method(self):
 
         self.pairs = [
             (self.points, self.multipolygons),
-            (self.geomcollection, self.polygons),
+            # Keep equal-length inputs so the align=False cases still run.
+            (self.geomcollection[: len(self.polygons)], self.polygons),
             (self.linestrings, self.multipoints),
             (self.linearrings, self.multilinestrings),
         ]
@@ -569,14 +571,6 @@ def test_to_arrow(self):
         import pyarrow as pa
 
         for geom in self.geoms:
-            # LINEARRING EMPTY and LineString EMPTY
-            # result in 01EA03000000000000 instead of 010200000000000000.
-            # Sedona returns the right result, so this bug is likely in pyarrow or geoarrow
-            # Below we set the modify the failing case as a workaround to pass the test
-            # Occurs in python 3.9, but fixed by python 3.10.
-            if geom[0] in [LineString(), LinearRing()]:
-                geom[0] = LineString([(0, 0), (1, 1)])
-
             sgpd_result = pa.array(GeoSeries(geom).to_arrow())
             gpd_result = pa.array(gpd.GeoSeries(geom).to_arrow())
             assert sgpd_result == gpd_result

python/tests/utils/test_geometry_serde.py

@@ -34,6 +34,27 @@
 
 
 class TestGeometrySerde(TestBase):
+    @pytest.mark.parametrize(
+        "wkt",
+        [
+            "POINT EMPTY",
+            "POINT Z EMPTY",
+            "LINESTRING EMPTY",
+            "LINESTRING Z EMPTY",
+            "POLYGON EMPTY",
+            "POLYGON Z EMPTY",
+            "GEOMETRYCOLLECTION (POINT EMPTY, LINESTRING EMPTY, POLYGON EMPTY)",
+            "GEOMETRYCOLLECTION Z (POINT Z EMPTY, LINESTRING Z EMPTY, POLYGON Z EMPTY)",
+        ],
+    )
+    def test_spark_empty_geometry_dimensions(self, wkt):
+        geometry = wkt_loads(wkt)
+        actual = self.spark.createDataFrame(
+            [(geometry,)], StructType().add("geom", GeometryType())
+        ).first()[0]
+
+        assert actual.wkb == geometry.wkb
+
     @pytest.mark.parametrize(
         "geom",
         [