Repository: shapely/shapely (language: python)

Issue:
BUG: creating LineString from array with invalid dimensions uses wrong coordinates
This changed in the last pygeos release (and thus also happens on shapely main):

```python
In [1]: import pygeos

# and using GEOS 3.9.1
In [2]: pygeos.__version__
Out[2]: '0.10.2'

In [3]: pygeos.linestrings(np.array([(1, 2, 3, 4), (4, 5, 6, 7)]))
---------------------------------------------------------------------------
GEOSException                             Traceback (most recent call last)
<ipython-input-3-9a89aa951713> in <module>
----> 1 pygeos.linestrings(np.array([(1, 2, 3, 4), (4, 5, 6, 7)]))

~/miniconda3/envs/geo/lib/python3.9/site-packages/pygeos/decorators.py in wrapped(*args, **kwargs)
     78             for arr in array_args:
     79                 arr.flags.writeable = False
---> 80             return func(*args, **kwargs)
     81         finally:
     82             for arr, old_flag in zip(array_args, old_flags):

~/miniconda3/envs/geo/lib/python3.9/site-packages/pygeos/creation.py in linestrings(coords, y, z, indices, **kwargs)
    101     coords = _xyz_to_coords(coords, y, z)
    102     if indices is None:
--> 103         return lib.linestrings(coords, **kwargs)
    104     else:
    105         return simple_geometries_1d(coords, indices, GeometryType.LINESTRING)

GEOSException: IllegalArgumentException: Unknown ordinate index 3
```

vs

```python
In [1]: import pygeos

# and uses GEOS 3.10.1
In [2]: pygeos.__version__
Out[2]: '0.12.0'

In [3]: pygeos.linestrings(np.array([(1, 2, 3, 4), (4, 5, 6, 7)]))
Out[3]: <pygeos.Geometry LINESTRING (1 2, 3 4)>
```

I suppose this is related to https://github.com/pygeos/pygeos/pull/436


Discussion / hints:
For (n, 3) it's still working correctly:

```
In [2]: shapely.linestrings(np.array([(1, 2, 3), (4, 5, 6)]))
Out[2]: <shapely.geometry.LineString LINESTRING Z (1 2 3, 4 5 6)>
```


Fix patch:
diff --git a/src/geos.c b/src/geos.c
index a3c3da16e..1577f2767 100644
--- a/src/geos.c
+++ b/src/geos.c
@@ -882,6 +882,12 @@ GEOSGeometry* PyGEOSForce3D(GEOSContextHandle_t ctx, GEOSGeometry* geom, double
   return force_dims(ctx, geom, 3, z);
 }
 
+/* Create a GEOSCoordSequence from an array
+ *
+ * Note: this function assumes that the dimension of the buffer is already
+ * checked before calling this function, so the buffer and the dims argument
+ * is only 2D or 3D.
+ */
 GEOSCoordSequence* coordseq_from_buffer(GEOSContextHandle_t ctx, const double* buf,
                                         unsigned int size, unsigned int dims, char ring_closure,
                                         npy_intp cs1, npy_intp cs2) {
diff --git a/src/ufuncs.c b/src/ufuncs.c
index 4f8fe13bf..79500b10a 100644
--- a/src/ufuncs.c
+++ b/src/ufuncs.c
@@ -2191,6 +2191,14 @@ static void linestrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
   GEOSCoordSequence* coord_seq = NULL;
   GEOSGeometry** geom_arr;
 
+  // check the ordinate dimension before calling coordseq_from_buffer
+  if (dimensions[2] < 2 || dimensions[2] > 3) {
+    PyErr_Format(PyExc_ValueError,
+                 "The ordinate (last) dimension should be 2 or 3, got %ld",
+                 dimensions[2]);
+    return;
+  }
+
   // allocate a temporary array to store output GEOSGeometry objects
   geom_arr = malloc(sizeof(void*) * dimensions[0]);
   CHECK_ALLOC(geom_arr);
@@ -2233,6 +2241,14 @@ static void linearrings_func(char** args, npy_intp* dimensions, npy_intp* steps,
   char ring_closure = 0;
   double first_coord, last_coord;
 
+  // check the ordinate dimension before calling coordseq_from_buffer
+  if (dimensions[2] < 2 || dimensions[2] > 3) {
+    PyErr_Format(PyExc_ValueError,
+                 "The ordinate (last) dimension should be 2 or 3, got %ld",
+                 dimensions[2]);
+    return;
+  }
+
   // allocate a temporary array to store output GEOSGeometry objects
   geom_arr = malloc(sizeof(void*) * dimensions[0]);
   CHECK_ALLOC(geom_arr);


Test patch:
diff --git a/shapely/tests/geometry/test_linestring.py b/shapely/tests/geometry/test_linestring.py
index 4b5346a30..0e6ccf637 100644
--- a/shapely/tests/geometry/test_linestring.py
+++ b/shapely/tests/geometry/test_linestring.py
@@ -102,8 +102,6 @@ def test_numpy_object_array():
     assert ar[0] == geom
 
 
-# TODO(shapely-2.0)
-@pytest.mark.xfail(strict=False, reason="Not yet implemented for Shapely 2.0")
 def test_from_invalid_dim():
     # TODO(shapely-2.0) better error message?
     # pytest.raises(ValueError, match="at least 2 coordinate tuples|at least 2 coordinates"):
@@ -111,21 +109,22 @@ def test_from_invalid_dim():
         LineString([(1, 2)])
 
     with pytest.raises(
-        ValueError,
-        match="Inconsistent coordinate dimensionality|Input operand 0 does not have enough dimensions",
+        ValueError, match="Input operand 0 does not have enough dimensions"
     ):
         LineString([(1, 2, 3), (4, 5)])
 
     with pytest.raises(
-        ValueError,
-        match="Inconsistent coordinate dimensionality|Input operand 0 does not have enough dimensions",
+        ValueError, match="Input operand 0 does not have enough dimensions"
     ):
         LineString([(1, 2), (3, 4, 5)])
 
-    # TODO better error, right now raises AssertionError
-    with pytest.raises(Exception):
+    msg = r"The ordinate \(last\) dimension should be 2 or 3, got {}"
+    with pytest.raises(ValueError, match=msg.format(4)):
         LineString([(1, 2, 3, 4), (4, 5, 6, 7)])
 
+    with pytest.raises(ValueError, match=msg.format(1)):
+        LineString([(1,), (4,)])
+
 
 def test_from_single_coordinate():
     """Test for issue #486"""
diff --git a/shapely/tests/test_creation.py b/shapely/tests/test_creation.py
index 541cb4a30..fe10a16eb 100644
--- a/shapely/tests/test_creation.py
+++ b/shapely/tests/test_creation.py
@@ -88,6 +88,23 @@ def test_linestrings_from_xyz():
     assert_geometries_equal(actual, shapely.Geometry("LINESTRING Z (0 2 0, 1 3 0)"))
 
 
+@pytest.mark.parametrize("dim", [2, 3])
+def test_linestrings_buffer(dim):
+    coords = np.random.randn(10, 3, dim)
+    coords1 = np.asarray(coords, order="C")
+    result1 = shapely.linestrings(coords1)
+
+    coords2 = np.asarray(coords1, order="F")
+    result2 = shapely.linestrings(coords2)
+    assert_geometries_equal(result1, result2)
+
+    # creating (.., 8, 8*3) strided array so it uses copyFromArrays
+    coords3 = np.asarray(np.swapaxes(np.swapaxes(coords, 0, 2), 1, 0), order="F")
+    coords3 = np.swapaxes(np.swapaxes(coords3, 0, 2), 1, 2)
+    result3 = shapely.linestrings(coords3)
+    assert_geometries_equal(result1, result3)
+
+
 def test_linestrings_invalid_shape_scalar():
     with pytest.raises(ValueError):
         shapely.linestrings((1, 1))
@@ -106,6 +123,28 @@ def test_linestrings_invalid_shape(shape):
         shapely.linestrings(np.ones(shape))
 
 
+def test_linestrings_invalid_ndim():
+    msg = r"The ordinate \(last\) dimension should be 2 or 3, got {}"
+
+    coords = np.ones((10, 2, 4), order="C")
+    with pytest.raises(ValueError, match=msg.format(4)):
+        shapely.linestrings(coords)
+
+    coords = np.ones((10, 2, 4), order="F")
+    with pytest.raises(ValueError, match=msg.format(4)):
+        shapely.linestrings(coords)
+
+    coords = np.swapaxes(np.swapaxes(np.ones((10, 2, 4)), 0, 2), 1, 0)
+    coords = np.swapaxes(np.swapaxes(np.asarray(coords, order="F"), 0, 2), 1, 2)
+    with pytest.raises(ValueError, match=msg.format(4)):
+        shapely.linestrings(coords)
+
+    # too few ordinates
+    coords = np.ones((10, 2, 1))
+    with pytest.raises(ValueError, match=msg.format(1)):
+        shapely.linestrings(coords)
+
+
 def test_linearrings():
     actual = shapely.linearrings(box_tpl(0, 0, 1, 1))
     assert_geometries_equal(
@@ -157,6 +196,23 @@ def test_linearrings_invalid_shape(shape):
         shapely.linearrings(coords)
 
 
+def test_linearrings_invalid_ndim():
+    msg = r"The ordinate \(last\) dimension should be 2 or 3, got {}"
+
+    coords1 = np.random.randn(10, 3, 4)
+    with pytest.raises(ValueError, match=msg.format(4)):
+        shapely.linearrings(coords1)
+
+    coords2 = np.hstack((coords1, coords1[:, [0], :]))
+    with pytest.raises(ValueError, match=msg.format(4)):
+        shapely.linearrings(coords2)
+
+    # too few ordinates
+    coords3 = np.random.randn(10, 3, 1)
+    with pytest.raises(ValueError, match=msg.format(1)):
+        shapely.linestrings(coords3)
+
+
 def test_linearrings_all_nan():
     coords = np.full((4, 2), np.nan)
     with pytest.raises(shapely.GEOSException):