From 18a36795d334baa965bacab3b2aa547268548b33 Mon Sep 17 00:00:00 2001
From: w-bonelli <wpb3hw@virginia.edu>
Date: Fri, 2 Sep 2022 15:55:01 -0400
Subject: [PATCH] refactor(CrossSectionPlot, GeoSpatialUtil): (#1521)

* accept flopy & shapely LineString for cross section line via GeoSpatialUtil
* check when creating new GeoSpatialUtil that linestrings have 2+ coordinates
---
 autotest/test_geospatial_util.py |  7 +++
 autotest/test_plot.py            | 78 ++++++++++++++++++++++++++++++++
 flopy/plot/crosssection.py       |  9 +++-
 flopy/utils/geospatial_utils.py  |  3 ++
 4 files changed, 96 insertions(+), 1 deletion(-)

diff --git a/autotest/test_geospatial_util.py b/autotest/test_geospatial_util.py
index 309b31b25..037667b05 100644
--- a/autotest/test_geospatial_util.py
+++ b/autotest/test_geospatial_util.py
@@ -543,3 +543,10 @@ def test_mixed_collection(
                     raise AssertionError(
                         "GeoSpatialCollection conversion error"
                     )
+
+
+@requires_pkg("shapely", "geojson")
+def test_create_linestring_with_single_point_fails():
+    point = [0, 0]
+    with pytest.raises(ValueError):
+        GeoSpatialUtil(point, shapetype="linestring")
diff --git a/autotest/test_plot.py b/autotest/test_plot.py
index d600eebe6..2385dac21 100644
--- a/autotest/test_plot.py
+++ b/autotest/test_plot.py
@@ -28,6 +28,7 @@
 from flopy.modpath import Modpath6, Modpath6Bas
 from flopy.plot import PlotCrossSection, PlotMapView
 from flopy.utils import CellBudgetFile, HeadFile, PathlineFile
+from flopy.utils.geometry import LineString
 
 
 def test_map_view():
@@ -523,3 +524,80 @@ def test_cross_section_with_quasi3d_layers(quasi3d_model):
     cs.plot_bc("wel")
     cs.plot_vector(frf, fff, flf, head=head)
     plt.savefig(os.path.join(str(quasi3d_model.model_ws), "plt02.png"))
+    plt.close()
+
+
+def structured_square_grid(side: int = 10, thick: int = 10):
+    """
+    Creates a basic 1-layer structured grid with the given thickness and number of cells per side
+    Parameters
+    ----------
+    side : The number of cells per side
+    thick : The thickness of the grid's single layer
+    Returns
+    -------
+    A single-layer StructuredGrid of the given size and thickness
+    """
+
+    from flopy.discretization.structuredgrid import StructuredGrid
+
+    delr = np.ones(side)
+    delc = np.ones(side)
+    top = np.ones((side, side)) * thick
+    botm = np.ones((side, side)) * (top - thick).reshape(1, side, side)
+    return StructuredGrid(delr=delr, delc=delc, top=top, botm=botm)
+
+
+@pytest.mark.parametrize(
+    "line",
+    [
+        (0, 0),
+        [0, 0],
+    ],
+)
+def test_cross_section_line_as_point_fails(line):
+    grid = structured_square_grid(side=10)
+    with pytest.raises(ValueError):
+        flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line})
+
+
+@requires_pkg("shapely")
+@pytest.mark.parametrize(
+    "line",
+    [
+        # diagonal line
+        [(0, 0), (10, 10)],  # test list of tuples
+        ([0, 0], [10, 10]),  # and tuple of lists
+        # horizontal line
+        ([0, 5.5], [10, 5.5]),
+        [(0, 5.5), (10, 5.5)],
+        # vertical line
+        [(0, 5.5), (10, 5.5)],
+        ([0, 5.5], [10, 5.5]),
+        # line with multiple segments
+        [(0, 0), (4, 6), (10, 10)],
+        ([0, 0], [4, 6], [10, 10]),
+    ],
+)
+def test_cross_section_line_representations(line):
+    import shapely
+
+    grid = structured_square_grid(side=10)
+
+    fls = LineString(line)
+    sls = shapely.geometry.LineString(line)
+
+    # use raw, flopy.utils.geometry and shapely.geometry representations
+    lxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": line})
+    fxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": fls})
+    sxc = flopy.plot.PlotCrossSection(modelgrid=grid, line={"line": sls})
+
+    # make sure parsed points are identical for all line representations
+    assert np.allclose(lxc.pts, fxc.pts) and np.allclose(lxc.pts, sxc.pts)
+    assert (
+        set(lxc.xypts.keys()) == set(fxc.xypts.keys()) == set(sxc.xypts.keys())
+    )
+    for k in lxc.xypts.keys():
+        assert np.allclose(lxc.xypts[k], fxc.xypts[k]) and np.allclose(
+            lxc.xypts[k], sxc.xypts[k]
+        )
diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py
index c0752ad5f..8558f8b35 100644
--- a/flopy/plot/crosssection.py
+++ b/flopy/plot/crosssection.py
@@ -7,6 +7,7 @@
 from matplotlib.patches import Polygon
 
 from ..utils import geometry
+from ..utils.geospatial_utils import GeoSpatialUtil
 from . import plotutil
 
 warnings.simplefilter("always", PendingDeprecationWarning)
@@ -135,7 +136,13 @@ def __init__(
                     (xcenter[int(line[onkey])], yedge[-1] - eps),
                 ]
         else:
-            verts = line[onkey]
+            ln = line[onkey]
+            gu = GeoSpatialUtil(ln, shapetype="linestring")
+            assert (
+                gu.shapetype == "LineString"
+            )  # unnecessary if GSU guarantees shapetype is same as requested
+
+            verts = gu.points
             xp = []
             yp = []
             for [v1, v2] in verts:
diff --git a/flopy/utils/geospatial_utils.py b/flopy/utils/geospatial_utils.py
index 0998f37e0..1ad14f630 100644
--- a/flopy/utils/geospatial_utils.py
+++ b/flopy/utils/geospatial_utils.py
@@ -75,6 +75,9 @@ def __init__(self, obj, shapetype=None):
                 )
                 raise AssertionError(err)
 
+            if shapetype == "linestring" and len(np.array(obj).shape) == 1:
+                raise ValueError(f"linestring must have 2 or more coordinates")
+
             self.__geo_interface = {
                 "type": shape_types[shapetype],
                 "coordinates": list(obj),