Repository: shapely/shapely (language: python)

Issue:
shapely crashes on mac with "abort" on very large floats


## Expected behavior and actual behavior.

If passed a very large float, `shapely` should maintain that representation, change to infinity, raise an exception, etc. Instead, I am getting a hard crash

## Steps to reproduce the problem.

Here is a minimalized example:
```python
import shapely
p = shapely.Point(1e+279, 0)
print(p)
```

On my Mac:

```shell
python shapely_thing.py
zsh: abort      poetry run python shapely_thing.py
```
## Operating system

Mac OS X 13.4.1

## Shapely version and provenance

```python
>>> shapely.__version__
'2.0.1'
```

Installed from `PyPI` using `poetry`

## Other notes

1) It appears that the constructor runs fine (commenting out `print(p)` gets rid of the error)
2) If I replace `1e+279` with `np.inf` it also works fine (returns a point with Infinity as one of the coordinates)
3) If I replace `1e+279` with `1e+278` it works fine.
4) it works okay on linux (Ubuntu 20.04.4 LTS), the print statement gives this:

```shell
POINT POINT (1000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0)```


Discussion / hints:
Thanks for the report! I can actually reproduce this on Linux (Ubuntu) as well

> It appears that the constructor runs fine (commenting out print(p) gets rid of the error)

Given this, I assume this is a GEOS issue in the WKT representation. Could you do a report at https://github.com/libgeos/geos/issues?

One way to reproduce this with the GEOS command line:

```
$ geosop -a "POINT (1e+27 0)" -f wkt
POINT (1000000000000000000000000000 0)
$ geosop -a "POINT (1e+279 0)" -f wkt
*** stack smashing detected ***: terminated
Aborted (core dumped)
```


We do have the policy of "covering up" GEOS crashes in shapely, if possible. Any ideas how we could do that here?
If we know the threshold for the coordinate value that triggers it, we could in theory check for that before calling the WKT routine? But then we have to check all the coordinates of a potentially large geometry..

Fix patch:
diff --git a/CHANGES.txt b/CHANGES.txt
index 246db1440..1e913f208 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -4,6 +4,10 @@ Changes
 2.1.0 (unreleased)
 ------------------
 
+Bug fixes:
+
+- Prevent crash when serializing a number > 1e100 to WKT. (#1907)
+
 Improvements:
 
 - Require GEOS >= 3.8, NumPy >= 1.16, and Python >= 3.8 (#1802, #1885)
diff --git a/shapely/_pygeos_api.pxd b/shapely/_pygeos_api.pxd
index 29eb0f9e1..83a1a70a5 100644
--- a/shapely/_pygeos_api.pxd
+++ b/shapely/_pygeos_api.pxd
@@ -27,6 +27,7 @@ cdef extern from "c_api.h":
         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,
diff --git a/src/geos.c b/src/geos.c
index 1189ad5d0..488154e71 100644
--- a/src/geos.c
+++ b/src/geos.c
@@ -304,6 +304,134 @@ 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,
+                                  int dimension) {
+  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))) {
+    return PGERR_GEOS_EXCEPTION;
+  }
+
+  if ((dimension > 2) && GEOSHasZ_r(ctx, geom)) {
+    if (!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
diff --git a/src/geos.h b/src/geos.h
index 1abe9f826..fec584ac7 100644
--- a/src/geos.h
+++ b/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, int dimension);
 #if GEOS_SINCE_3_9_0
 extern char wkt_empty_3d_geometry(GEOSContextHandle_t ctx, GEOSGeometry* geom,
                                   char** wkt);
diff --git a/src/pygeom.c b/src/pygeom.c
index 9d29e9041..2c058c0d0 100644
--- a/src/pygeom.c
+++ b/src/pygeom.c
@@ -70,6 +70,10 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
   char* wkt;
   PyObject* result;
   GEOSGeometry* geom = obj->ptr;
+  char trim = 1;
+  int precision = 3;
+  int dimension = 3;
+  int use_old_3d = 0;
 
   if (geom == NULL) {
     Py_INCREF(Py_None);
@@ -77,6 +81,12 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
   }
 
   GEOS_INIT;
+  if (trim) {
+    errstate = check_to_wkt_trim_compatible(ctx, geom, dimension);
+  }
+  if (errstate != PGERR_SUCCESS) {
+    goto finish;
+  }
 
 #if GEOS_SINCE_3_9_0
   errstate = wkt_empty_3d_geometry(ctx, geom, &wkt);
@@ -101,10 +111,6 @@ static PyObject* GeometryObject_ToWKT(GeometryObject* obj) {
     goto finish;
   }
 
-  char trim = 1;
-  int precision = 3;
-  int dimension = 3;
-  int use_old_3d = 0;
   GEOSWKTWriter_setRoundingPrecision_r(ctx, writer, precision);
   GEOSWKTWriter_setTrim_r(ctx, writer, trim);
   GEOSWKTWriter_setOutputDimension_r(ctx, writer, dimension);
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 2781df774..e0aa6272b 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -2458,9 +2458,9 @@ static void points_func(char** args, const npy_intp* dimensions, const npy_intp*
   GEOS_INIT_THREADS;
 
   char *ip1 = args[0];               
-  npy_intp is1 = steps[0], cs1 = steps[3]; 
-  npy_intp n = dimensions[0], n_c1 = dimensions[1];        
-  npy_intp i;                                        
+  npy_intp is1 = steps[0], cs1 = steps[3];
+  npy_intp n = dimensions[0], n_c1 = dimensions[1];
+  npy_intp i;
   for (i = 0; i < n; i++, ip1 += is1) {
     CHECK_SIGNALS_THREADS(i);
     if (errstate == PGERR_PYSIGNAL) {
@@ -3217,7 +3217,8 @@ static PyUFuncGenericFunction to_wkb_funcs[1] = {&to_wkb_func};
 
 static char to_wkt_dtypes[6] = {NPY_OBJECT, NPY_INT,  NPY_BOOL,
                                 NPY_INT,    NPY_BOOL, NPY_OBJECT};
-static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp* steps, void* data) {
+static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp* steps,
+                        void* data) {
   char *ip1 = args[0], *ip2 = args[1], *ip3 = args[2], *ip4 = args[3], *ip5 = args[4],
        *op1 = args[5];
   npy_intp is1 = steps[0], is2 = steps[1], is3 = steps[2], is4 = steps[3], is5 = steps[4],
@@ -3234,6 +3235,11 @@ static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp*
     return;
   }
 
+  int precision = *(int*)ip2;
+  npy_bool trim = *(npy_bool*)ip3;
+  int dimension = *(int*)ip4;
+  int use_old_3d = *(npy_bool*)ip5;
+
   GEOS_INIT;
 
   /* Create the WKT writer */
@@ -3242,10 +3248,10 @@ static void to_wkt_func(char** args, const npy_intp* dimensions, const npy_intp*
     errstate = PGERR_GEOS_EXCEPTION;
     goto finish;
   }
-  GEOSWKTWriter_setRoundingPrecision_r(ctx, writer, *(int*)ip2);
-  GEOSWKTWriter_setTrim_r(ctx, writer, *(npy_bool*)ip3);
-  GEOSWKTWriter_setOutputDimension_r(ctx, writer, *(int*)ip4);
-  GEOSWKTWriter_setOld3D_r(ctx, writer, *(npy_bool*)ip5);
+  GEOSWKTWriter_setRoundingPrecision_r(ctx, writer, precision);
+  GEOSWKTWriter_setTrim_r(ctx, writer, trim);
+  GEOSWKTWriter_setOutputDimension_r(ctx, writer, dimension);
+  GEOSWKTWriter_setOld3D_r(ctx, writer, use_old_3d);
 
   // Check if the above functions caused a GEOS exception
   if (last_error[0] != 0) {
@@ -3269,6 +3275,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, dimension);
+        if (errstate != PGERR_SUCCESS) {
+          goto finish;
+        }
+      }
 #if GEOS_SINCE_3_9_0
       errstate = wkt_empty_3d_geometry(ctx, in1, &wkt);
       if (errstate != PGERR_SUCCESS) {


Test patch:
diff --git a/shapely/tests/geometry/test_format.py b/shapely/tests/geometry/test_format.py
index 6ae61a34c..689f2e0ba 100644
--- a/shapely/tests/geometry/test_format.py
+++ b/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():
diff --git a/shapely/tests/test_io.py b/shapely/tests/test_io.py
index 6b79d26e6..8508e47c5 100644
--- a/shapely/tests/test_io.py
+++ b/shapely/tests/test_io.py
@@ -7,10 +7,26 @@
 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,
+    point,
+    point_z,
+    polygon_z,
+)
 
 # fmt: off
 POINT11_WKB = b"\x01\x01\x00\x00\x00" + struct.pack("<2d", 1.0, 1.0)
@@ -348,6 +364,68 @@ def test_to_wkt_multipoint_with_point_empty_errors():
         shapely.to_wkt(geom)
 
 
+@pytest.mark.parametrize("geom", [Point(1e100, 0), Point(0, 1e100)])
+def test_to_wkt_large_float_ok(geom):
+    # https://github.com/shapely/shapely/issues/1903
+    shapely.to_wkt(geom)
+    assert "Exception in WKT writer" not in repr(geom)
+
+
+@pytest.mark.parametrize("geom", [Point(1e101, 0), Point(0, 1e101)])
+def test_to_wkt_large_float_err(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)
+
+
+@pytest.mark.parametrize(
+    "geom",
+    [
+        # We implemented our own "GetZMax", so go through all geometry types:
+        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(
+            [(0, 0, 0), (0, 10, 0), (10, 0, 0), (0, 0, 0)],
+            [[(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)]],
+        ),
+        MultiPoint([(0, 0, 0), (0, 0, 1e101)]),
+        MultiLineString(
+            [LineString([(0, 0, 0), (0, 1, 0)]), LineString([(0, 0, 0), (0, 1, 1e101)])]
+        ),
+        MultiPolygon(
+            [polygon_z, Polygon([(0, 0, 0), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)])]
+        ),
+        GeometryCollection([point_z, Point(0, 0, 1e101)]),
+        GeometryCollection([GeometryCollection([Point(0, 0, 1e101)])]),
+        LineString([(0, 0, np.nan), (0, 0, 1e101)]),
+        Polygon([(0, 0, np.nan), (0, 1, 0), (1, 0, 1e101), (0, 0, 0)]),
+        GeometryCollection([Point(0, 0), Point(0, 0, 1e101)]),
+    ],
+)
+def test_to_wkt_large_float_3d_no_crash(geom):
+    # https://github.com/shapely/shapely/issues/1903
+    # just test if there is a crash (detailed behaviour differs per GEOS version)
+    try:
+        shapely.to_wkt(geom)
+    except ValueError as e:
+        assert str(e).startswith("WKT output of coordinates greater than")
+    repr(geom)
+
+
+def test_to_wkt_large_float_skip_z():
+    # https://github.com/shapely/shapely/issues/1903
+    assert shapely.to_wkt(Point(0, 0, 1e101), output_dimension=2) == "POINT (0 0)"
+
+
+def test_to_wkt_large_float_no_trim():
+    # https://github.com/shapely/shapely/issues/1903
+    # don't test the exact number, it is ridiculously large and probably platform dependent
+    assert shapely.to_wkt(Point(1e101, 0), trim=False).startswith("POINT (")
+
+
 def test_repr():
     assert repr(point) == "<POINT (2 3)>"