Merge 79f2222 into 42392bb

shapely/tests/geometry/test_format.py

@@ -18,13 +18,10 @@ def test_format_invalid():
             format(pt, format_spec)
 
 
-def test_format_point():
-    # example coordinate data
+def get_tst_format_point_params():
     xy1 = (0.12345678901234567, 1.2345678901234567e10)
     xy2 = (-169.910918, -18.997564)
     xyz3 = (630084, 4833438, 76)
-
-    # list of tuples to test; see structure at top of the for-loop
     test_list = [
         (".0f", xy1, "POINT (0 12345678901)", True),
         (".1f", xy1, "POINT (0.1 12345678901.2)", True),
@@ -53,26 +50,32 @@ def test_format_point():
         ("F", xyz3, "POINT Z ({:.16f} {:.16f} {:.16f})".format(*xyz3), False),
         ("g", xyz3, "POINT Z (630084 4833438 76)", False),
     ]
-    for format_spec, coords, expt_wkt, same_python_float in test_list:
-        pt = Point(*coords)
-        # basic checks
-        assert f"{pt}" == pt.wkt
-        assert format(pt, "") == pt.wkt
-        assert format(pt, "x") == pt.wkb_hex.lower()
-        assert format(pt, "X") == pt.wkb_hex
-        # check formatted WKT to expected
-        assert format(pt, format_spec) == expt_wkt, format_spec
-        # check Python's format consistency
-        text_coords = expt_wkt[expt_wkt.index("(") + 1 : expt_wkt.index(")")]
-        is_same = []
-        for coord, expt_coord in zip(coords, text_coords.split()):
-            py_fmt_float = format(float(coord), format_spec)
-            if same_python_float:
-                assert py_fmt_float == expt_coord, format_spec
-            else:
-                is_same.append(py_fmt_float == expt_coord)
-        if not same_python_float:
-            assert not all(is_same), f"{format_spec!r} with {expt_wkt}"
+    return test_list
+
+
+@pytest.mark.parametrize(
+    "format_spec, coords, expt_wkt, same_python_float", get_tst_format_point_params()
+)
+def test_format_point(format_spec, coords, expt_wkt, same_python_float):
+    pt = Point(*coords)
+    # basic checks
+    assert f"{pt}" == pt.wkt
+    assert format(pt, "") == pt.wkt
+    assert format(pt, "x") == pt.wkb_hex.lower()
+    assert format(pt, "X") == pt.wkb_hex
+    # check formatted WKT to expected
+    assert format(pt, format_spec) == expt_wkt, format_spec
+    # check Python's format consistency
+    text_coords = expt_wkt[expt_wkt.index("(") + 1 : expt_wkt.index(")")]
+    is_same = []
+    for coord, expt_coord in zip(coords, text_coords.split()):
+        py_fmt_float = format(float(coord), format_spec)
+        if same_python_float:
+            assert py_fmt_float == expt_coord, format_spec
+        else:
+            is_same.append(py_fmt_float == expt_coord)
+    if not same_python_float:
+        assert not all(is_same), f"{format_spec!r} with {expt_wkt}"
 
 
 def test_format_polygon():

shapely/tests/test_io.py

@@ -7,10 +7,28 @@
 import pytest
 
 import shapely
-from shapely import GeometryCollection, LineString, Point, Polygon
+from shapely import (
+    GeometryCollection,
+    LinearRing,
+    LineString,
+    MultiLineString,
+    MultiPoint,
+    MultiPolygon,
+    Point,
+    Polygon,
+)
 from shapely.errors import UnsupportedGEOSVersionError
 from shapely.testing import assert_geometries_equal
-from shapely.tests.common import all_types, empty_point, empty_point_z, point, point_z
+from shapely.tests.common import (
+    all_types,
+    empty_point,
+    empty_point_z,
+    line_string,
+    linear_ring,
+    point,
+    point_z,
+    polygon,
+)
 
 # fmt: off
 POINT11_WKB = b"\x01\x01\x00\x00\x00" + struct.pack("<2d", 1.0, 1.0)
@@ -348,6 +366,43 @@ def test_to_wkt_multipoint_with_point_empty_errors():
         shapely.to_wkt(geom)
 
 
+def test_to_wkt_large_float_ok():
+    # https://github.com/shapely/shapely/issues/1903
+    shapely.to_wkt(Point([1e100, 0.0]))
+
+
+def test_to_wkt_large_float_not_ok():
+    # https://github.com/shapely/shapely/issues/1903
+    with pytest.raises(ValueError, match="WKT output of coordinates greater than.*"):
+        shapely.to_wkt(Point([1e101, 0.0]))
+
+
+@pytest.mark.parametrize(
+    "geom",
+    [
+        Point(1e101, 0, 0),
+        Point(0, 1e101, 0),
+        Point(0, 0, 1e101),
+        LineString([(0, 0, 0), (0, 0, 1e101)]),
+        LinearRing([(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)]),
+        Polygon([(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)]),
+        Polygon(linear_ring, [[(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)]]),
+        MultiPoint([(0, 0, 0), (0, 0, 1e101)]),
+        MultiLineString([line_string, LineString([(0, 0, 0), (0, 0, 1e101)])]),
+        MultiPolygon(
+            [polygon, Polygon([(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)])]
+        ),
+        GeometryCollection([point, Point(0, 0, 1e101)]),
+        GeometryCollection([GeometryCollection([Point(0, 0, 1e101)])]),
+    ],
+)
+def test_to_wkt_large_float(geom):
+    # https://github.com/shapely/shapely/issues/1903
+    with pytest.raises(ValueError, match="WKT output of coordinates greater than.*"):
+        shapely.to_wkt(geom)
+    assert "Exception in WKT writer" in repr(geom)
+
+
 def test_repr():
     assert repr(point) == "<POINT (2 3)>"
 

src/geos.c

@@ -304,6 +304,128 @@ char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom) {
   }
 }
 
+/* Checks whether the geometry contains a coordinate greater than 1E+100
+ *
+ * See also:
+ *
+ * https://github.com/shapely/shapely/issues/1903
+ */
+char get_zmax(GEOSContextHandle_t, const GEOSGeometry*, double*);
+char get_zmax_simple(GEOSContextHandle_t, const GEOSGeometry*, double*);
+char get_zmax_polygon(GEOSContextHandle_t, const GEOSGeometry*, double*);
+char get_zmax_collection(GEOSContextHandle_t, const GEOSGeometry*, double*);
+
+char get_zmax(GEOSContextHandle_t ctx, const GEOSGeometry* geom, double* zmax) {
+  int type = GEOSGeomTypeId_r(ctx, geom);
+  if ((type == 0) || (type == 1) || (type == 2)) {
+    return get_zmax_simple(ctx, geom, zmax);
+  } else if (type == 3) {
+    return get_zmax_polygon(ctx, geom, zmax);
+  } else if ((type >= 4) && (type <= 7)) {
+    return get_zmax_collection(ctx, geom, zmax);
+  } else {
+    return 0;
+  }
+}
+
+char get_zmax_simple(GEOSContextHandle_t ctx, const GEOSGeometry* geom, double* zmax) {
+  const GEOSCoordSequence* seq;
+  unsigned int n, i;
+  double coord;
+
+  seq = GEOSGeom_getCoordSeq_r(ctx, geom);
+  if (seq == NULL) {
+    return 0;
+  }
+  if (GEOSCoordSeq_getSize_r(ctx, seq, &n) == 0) {
+    return 0;
+  }
+
+  for (i = 0; i < n; i++) {
+    if (!GEOSCoordSeq_getZ_r(ctx, seq, i, &coord)) {
+      return 0;
+    }
+    if (npy_isfinite(coord) && (coord > *zmax)) {
+      *zmax = coord;
+    }
+  }
+  return 1;
+}
+
+char get_zmax_polygon(GEOSContextHandle_t ctx, const GEOSGeometry* geom, double* zmax) {
+  const GEOSGeometry* ring;
+  int n, i;
+
+  ring = GEOSGetExteriorRing_r(ctx, geom);
+  if (ring == NULL) {
+    return 0;
+  }
+  if (!get_zmax_simple(ctx, ring, zmax)) {
+    return 0;
+  }
+  n = GEOSGetNumInteriorRings_r(ctx, geom);
+  if (n == -1) {
+    return 0;
+  }
+
+  for (i = 0; i < n; i++) {
+    ring = GEOSGetInteriorRingN_r(ctx, geom, i);
+    if (ring == NULL) {
+      return 0;
+    }
+    if (!get_zmax_simple(ctx, ring, zmax)) {
+      return 0;
+    }
+  }
+  return 1;
+}
+
+char get_zmax_collection(GEOSContextHandle_t ctx, const GEOSGeometry* geom,
+                         double* zmax) {
+  const GEOSGeometry* elem;
+  int n, i;
+
+  n = GEOSGetNumGeometries_r(ctx, geom);
+  if (n == -1) {
+    return 0;
+  }
+
+  for (i = 0; i < n; i++) {
+    elem = GEOSGetGeometryN_r(ctx, geom, i);
+    if (elem == NULL) {
+      return 0;
+    }
+    if (!get_zmax(ctx, elem, zmax)) {
+      return 0;
+    }
+  }
+  return 1;
+}
+
+char check_to_wkt_trim_compatible(GEOSContextHandle_t ctx, const GEOSGeometry* geom) {
+  char is_empty;
+  double xmax = 0.0;
+  double ymax = 0.0;
+  double zmax = 0.0;
+
+  if (GEOSisEmpty_r(ctx, geom)) {
+    return PGERR_SUCCESS;
+  }
+
+  // use max coordinates to check if any coordinate is too large
+  if (!(GEOSGeom_getXMax_r(ctx, geom, &xmax) && GEOSGeom_getYMax_r(ctx, geom, &ymax) &&
+        get_zmax(ctx, geom, &zmax))) {
+    return PGERR_GEOS_EXCEPTION;
+  }
+
+  if ((npy_isfinite(xmax) && (xmax > 1E100)) || (npy_isfinite(ymax) && (ymax > 1E100)) ||
+      (npy_isfinite(zmax) && (zmax > 1E100))) {
+    return PGERR_COORD_OUT_OF_BOUNDS;
+  }
+
+  return PGERR_SUCCESS;
+}
+
 #if GEOS_SINCE_3_9_0
 
 /* Checks whether the geometry is a 3D empty geometry and, if so, create the WKT string

src/geos.h

@@ -52,6 +52,7 @@ enum ShapelyErrorCode {
   PGERR_NO_MALLOC,
   PGERR_GEOMETRY_TYPE,
   PGERR_MULTIPOINT_WITH_POINT_EMPTY,
+  PGERR_COORD_OUT_OF_BOUNDS,
   PGERR_EMPTY_GEOMETRY,
   PGERR_GEOJSON_EMPTY_POINT,
   PGERR_LINEARRING_NCOORDS,
@@ -90,6 +91,11 @@ enum ShapelyErrorCode {
                       "WKT output of multipoints with an empty point is unsupported on " \
                       "this version of GEOS.");                                          \
       break;                                                                             \
+    case PGERR_COORD_OUT_OF_BOUNDS:                                                      \
+      PyErr_SetString(PyExc_ValueError,                                                  \
+                      "WKT output of coordinates greater than 1E+100 is unsupported on " \
+                      "this version of GEOS.");                                          \
+      break;                                                                             \
     case PGERR_EMPTY_GEOMETRY:                                                           \
       PyErr_SetString(PyExc_ValueError, "One of the Geometry inputs is empty.");         \
       break;                                                                             \
@@ -167,6 +173,7 @@ extern char has_point_empty(GEOSContextHandle_t ctx, GEOSGeometry* geom);
 extern GEOSGeometry* point_empty_to_nan_all_geoms(GEOSContextHandle_t ctx,
                                                   GEOSGeometry* geom);
 extern char check_to_wkt_compatible(GEOSContextHandle_t ctx, GEOSGeometry* geom);
+extern char check_to_wkt_trim_compatible(GEOSContextHandle_t ctx, const GEOSGeometry* geom);
 #if GEOS_SINCE_3_9_0
 extern char wkt_empty_3d_geometry(GEOSContextHandle_t ctx, GEOSGeometry* geom,
                                   char** wkt);

src/pygeom.c

@@ -77,6 +77,10 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
   }
 
   GEOS_INIT;
+  errstate = check_to_wkt_trim_compatible(ctx, geom);
+  if (errstate != PGERR_SUCCESS) {
+    goto finish;
+  }
 
 #if GEOS_SINCE_3_9_0
   errstate = wkt_empty_3d_geometry(ctx, geom, &wkt);

src/ufuncs.c

@@ -3234,6 +3234,8 @@ static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp*
     return;
   }
 
+  npy_bool trim = *(npy_bool*)ip3;
+
   GEOS_INIT;
 
   /* Create the WKT writer */
@@ -3243,7 +3245,7 @@ static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp*
     goto finish;
   }
   GEOSWKTWriter_setRoundingPrecision_r(ctx, writer, *(int*)ip2);
-  GEOSWKTWriter_setTrim_r(ctx, writer, *(npy_bool*)ip3);
+  GEOSWKTWriter_setTrim_r(ctx, writer, trim);
   GEOSWKTWriter_setOutputDimension_r(ctx, writer, *(int*)ip4);
   GEOSWKTWriter_setOld3D_r(ctx, writer, *(npy_bool*)ip5);
 
@@ -3269,6 +3271,12 @@ static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp*
       Py_INCREF(Py_None);
       *out = Py_None;
     } else {
+      if (trim) {
+        errstate = check_to_wkt_trim_compatible(ctx, in1);
+        if (errstate != PGERR_SUCCESS) {
+          goto finish;
+        }
+      }
 #if GEOS_SINCE_3_9_0
       errstate = wkt_empty_3d_geometry(ctx, in1, &wkt);
       if (errstate != PGERR_SUCCESS) {