Retain coordinate type for empty simplified geometries

This fixes a bug where, in some scenarios, the `Simplify` operation
would drop extended coordinate type information (from XYZ, XYM, or XYZM
back to XY).

The bug happened because some geometry constructors inferred the
geometry type from the input's lowest common denominator coordinate
type. But if the input is empty, then the coordinate type is inferred as
XY. The constructor design is a bit of a sharp edge and is not addressed
here (instead, the bug that the sharp edge made more likely is fixed).

CHANGELOG.md

@@ -7,6 +7,9 @@
 - Fixes a bug where geometry collections with mixed coordinate types were
   erroneously allowed during WKT and WKB parsing.
 
+- Fixes a bug where the `Simplify` method would drop coordinate type to XY in
+  some scenarios where the result is an empty geometry.
+
 ## v0.50.0
 
 2024-05-07

geom/alg_simplify_test.go

@@ -88,6 +88,30 @@ func TestSimplify(t *testing.T) {
 		{"POLYGON Z((2 2 10,2 3 20,3 3 30,3 2 40,2 2 10))", 0.0, "POLYGON Z((2 2 10,2 3 20,3 3 30,3 2 40,2 2 10))"},
 		{"POLYGON M((2 2 10,2 3 20,3 3 30,3 2 40,2 2 10))", 0.0, "POLYGON M((2 2 10,2 3 20,3 3 30,3 2 40,2 2 10))"},
 		{"POLYGON ZM((2 2 10 11,2 3 20 21,3 3 30 31,3 2 40 41,2 2 10 11))", 0.0, "POLYGON ZM((2 2 10 11,2 3 20 21,3 3 30 31,3 2 40 41,2 2 10 11))"},
+
+		// Empty Z, M, ZM
+		{"POINT Z EMPTY", 0.1, "POINT Z EMPTY"},
+		{"POINT M EMPTY", 0.1, "POINT M EMPTY"},
+		{"POINT ZM EMPTY", 0.1, "POINT ZM EMPTY"},
+		{"LINESTRING Z EMPTY", 0.1, "LINESTRING Z EMPTY"},
+		{"LINESTRING M EMPTY", 0.1, "LINESTRING M EMPTY"},
+		{"LINESTRING ZM EMPTY", 0.1, "LINESTRING ZM EMPTY"},
+		{"POLYGON Z EMPTY", 0.1, "POLYGON Z EMPTY"},
+		{"POLYGON M EMPTY", 0.1, "POLYGON M EMPTY"},
+		{"POLYGON ZM EMPTY", 0.1, "POLYGON ZM EMPTY"},
+		{"MULTIPOINT Z EMPTY", 0.1, "MULTIPOINT Z EMPTY"},
+		{"MULTIPOINT M EMPTY", 0.1, "MULTIPOINT M EMPTY"},
+		{"MULTIPOINT ZM EMPTY", 0.1, "MULTIPOINT ZM EMPTY"},
+		{"MULTILINESTRING Z EMPTY", 0.1, "MULTILINESTRING Z EMPTY"},
+		{"MULTILINESTRING M EMPTY", 0.1, "MULTILINESTRING M EMPTY"},
+		{"MULTILINESTRING ZM EMPTY", 0.1, "MULTILINESTRING ZM EMPTY"},
+		{"MULTIPOLYGON Z EMPTY", 0.1, "MULTIPOLYGON Z EMPTY"},
+		{"MULTIPOLYGON M EMPTY", 0.1, "MULTIPOLYGON M EMPTY"},
+		{"MULTIPOLYGON ZM EMPTY", 0.1, "MULTIPOLYGON ZM EMPTY"},
+		{"GEOMETRYCOLLECTION Z EMPTY", 0.1, "GEOMETRYCOLLECTION Z EMPTY"},
+		{"GEOMETRYCOLLECTION M EMPTY", 0.1, "GEOMETRYCOLLECTION M EMPTY"},
+		{"GEOMETRYCOLLECTION ZM EMPTY", 0.1, "GEOMETRYCOLLECTION ZM EMPTY"},
+		{"GEOMETRYCOLLECTION EMPTY", 0.1, "GEOMETRYCOLLECTION EMPTY"},
 	} {
 		t.Run(strconv.Itoa(i), func(t *testing.T) {
 			in := geomFromWKT(t, tc.input)

geom/type_geometry_collection.go

@@ -550,7 +550,7 @@ func (c GeometryCollection) Simplify(threshold float64, nv ...NoValidate) (Geome
 			return GeometryCollection{}, wrapSimplified(err)
 		}
 	}
-	return NewGeometryCollection(geoms), nil
+	return NewGeometryCollection(geoms).ForceCoordinatesType(c.CoordinatesType()), nil
 }
 
 // Densify returns a new GeometryCollection with additional linearly

geom/type_multi_line_string.go

@@ -508,7 +508,7 @@ func (m MultiLineString) Simplify(threshold float64) MultiLineString {
 			lss = append(lss, ls)
 		}
 	}
-	return NewMultiLineString(lss)
+	return NewMultiLineString(lss).ForceCoordinatesType(m.CoordinatesType())
 }
 
 // Densify returns a new MultiLineString with additional linearly interpolated

geom/type_multi_polygon.go

@@ -567,7 +567,7 @@ func (m MultiPolygon) Simplify(threshold float64, nv ...NoValidate) (MultiPolygo
 			return MultiPolygon{}, wrapSimplified(err)
 		}
 	}
-	return simpl, nil
+	return simpl.ForceCoordinatesType(m.CoordinatesType()), nil
 }
 
 // Densify returns a new MultiPolygon with additional linearly interpolated

geom/type_polygon.go

@@ -661,7 +661,7 @@ func (p Polygon) Simplify(threshold float64, nv ...NoValidate) (Polygon, error)
 		return ring.Coordinates().Length() < 4
 	}
 	if hasCollapsed(exterior) {
-		return Polygon{}, nil
+		return Polygon{}.ForceCoordinatesType(p.CoordinatesType()), nil
 	}
 
 	n := p.NumInteriorRings()