gridintersect.py

import numpy as np

try:
    import matplotlib.pyplot as plt
except (ImportError, RuntimeError):
    plt = None

from .geometry import transform
from .geospatial_utils import GeoSpatialUtil

try:
    from shapely.geometry import (
        MultiPoint,
        Point,
        Polygon,
        box,
        GeometryCollection,
        MultiPolygon,
    )
    from shapely.strtree import STRtree
    from shapely.affinity import translate, rotate
    from shapely.prepared import prep

    shply = True
except:
    shply = False

import contextlib
import warnings
from distutils.version import LooseVersion

try:
    import shapely

    SHAPELY_GE_20 = str(shapely.__version__) >= LooseVersion("2.0")
except:
    shapely = None
    SHAPELY_GE_20 = False

try:
    from shapely.errors import ShapelyDeprecationWarning as shapely_warning
except:
    shapely_warning = None

if shapely_warning is not None and not SHAPELY_GE_20:

    @contextlib.contextmanager
    def ignore_shapely_warnings_for_object_array():
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore", category=shapely_warning)
            yield


else:

    @contextlib.contextmanager
    def ignore_shapely_warnings_for_object_array():
        yield


def parse_shapely_ix_result(collection, ix_result, shptyps=None):
    """Recursive function for parsing shapely intersection results. Returns a
    list of shapely shapes matching shptyp.

    Parameters
    ----------
    collection : list
        state variable for storing result, generally
        an empty list
    ix_result : shapely.geometry type
        any shapely intersection result
    shptyp : str, list of str, or None, optional
        if None (default), return all types of shapes.
        if str, return shapes of that type, if list of str,
        return all types in list

    Returns
    -------
    collection : list
        list containing shapely geometries of type shptyp
    """
    # convert shptyps to list if needed
    if isinstance(shptyps, str):
        shptyps = [shptyps]
    elif shptyps is None:
        shptyps = [None]

    # if empty
    if ix_result.is_empty:
        return collection
    # base case: geom_type is partial or exact match to shptyp
    elif ix_result.geom_type in shptyps:
        collection.append(ix_result)
        return collection
    # recursion for collections
    elif hasattr(ix_result, "geoms"):
        for ishp in ix_result.geoms:
            parse_shapely_ix_result(collection, ishp, shptyps=shptyps)
    # if collecting all types
    elif shptyps[0] is None:
        return collection.append(ix_result)
    return collection


class GridIntersect:
    """Class for intersecting shapely shapes (Point, Linestring, Polygon, or
    their Multi variants) with MODFLOW grids. Contains optimized search
    routines for structured grids.

    Notes
    -----
     - The STR-tree query is based on the bounding box of the shape or
       collection, if the bounding box of the shape covers nearly the entire
       grid, the query won't be able to limit the search space much, resulting
       in slower performance. Therefore, it can sometimes be faster to
       intersect each individual shape in a collection than it is to intersect
       with the whole collection at once.
     - Building the STR-tree can take a while for large grids. Once built the
       intersect routines (for individual shapes) should be pretty fast. It
       is possible to perform intersects without building the STR-tree by
       setting `rtree=False`.
     - The optimized routines for structured grids will often outperform
       the shapely routines because of the reduced overhead of building and
       parsing the STR-tree. However, for polygons the STR-tree implementation
       is often faster than the optimized structured routines, especially
       for larger grids.
    """

    def __init__(self, mfgrid, method=None, rtree=True):
        """Intersect shapes (Point, Linestring, Polygon) with a modflow grid.

        Parameters
        ----------
        mfgrid : flopy modflowgrid
            MODFLOW grid as implemented in flopy
        method : str, optional
            default is None, which determines intersection method based on
            the grid type. Options are either 'vertex' which uses shapely
            interesection operations or 'structured' which uses optimized
            methods that only work for structured grids
        rtree : bool, optional
            whether to build an STR-Tree, default is True. If False no
            STR-tree is built (which saves some time), but intersects will
            loop through all model gridcells (which is generally slower).
            Only read when `method='vertex'`.
        """
        if not shply:
            msg = (
                "Shapely is needed for grid intersect operations! "
                "Please install shapely if you need to use grid intersect "
                "functionality."
            )
            raise ModuleNotFoundError(msg)

        self.mfgrid = mfgrid
        if method is None:
            # determine method from grid_type
            self.method = self.mfgrid.grid_type
        else:
            # set method
            self.method = method
        self.rtree = rtree

        if self.method == "vertex":
            # set method to get gridshapes depending on grid type
            self._set_method_get_gridshapes()

            # build STR-tree if specified
            if self.rtree:
                self.strtree = STRtree(self._get_gridshapes())

        elif self.method == "structured" and mfgrid.grid_type == "structured":
            pass

        else:
            raise ValueError(
                "Method '{0}' not recognized or "
                "not supported "
                "for grid_type '{1}'!".format(
                    self.method, self.mfgrid.grid_type
                )
            )

    def intersect(self, shp, **kwargs):
        """
        Method to intersect a shape with a model grid

        Parameters
        ----------
        shp : shapely.geometry, geojson object, shapefile.Shape,
              or flopy geomerty object
        sort_by_cellid : bool
            Sort results by cellid
        keepzerolengths : bool
            boolean method to keep zero length intersections for
            linestring intersection

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        gu = GeoSpatialUtil(shp)
        shp = gu.shapely
        sort_by_cellid = kwargs.pop("sort_by_cellid", True)
        keepzerolengths = kwargs.pop("keepzerolengths", False)

        if gu.shapetype in ("Point", "MultiPoint"):
            if (
                self.method == "structured"
                and self.mfgrid.grid_type == "structured"
            ):
                rec = self._intersect_point_structured(shp)
            else:
                rec = self._intersect_point_shapely(shp, sort_by_cellid)
        elif gu.shapetype in ("LineString", "MultiLineString"):
            if (
                self.method == "structured"
                and self.mfgrid.grid_type == "structured"
            ):
                rec = self._intersect_linestring_structured(
                    shp, keepzerolengths
                )
            else:
                rec = self._intersect_linestring_shapely(
                    shp, keepzerolengths, sort_by_cellid
                )
        elif gu.shapetype in ("Polygon", "MultiPolygon"):
            if (
                self.method == "structured"
                and self.mfgrid.grid_type == "structured"
            ):
                rec = self._intersect_polygon_structured(shp)
            else:
                rec = self._intersect_polygon_shapely(shp, sort_by_cellid)
        else:
            err = "Shapetype {} is not supported".format(gu.shapetype)
            raise TypeError(err)

        return rec

    def _set_method_get_gridshapes(self):
        """internal method, set self._get_gridshapes to the certain method for
        obtaining gridcells."""
        # Set method for obtaining grid shapes
        if self.mfgrid.grid_type == "structured":
            self._get_gridshapes = self._rect_grid_to_shape_generator
        elif self.mfgrid.grid_type == "vertex":
            self._get_gridshapes = self._vtx_grid_to_shape_generator
        elif self.mfgrid.grid_type == "unstructured":
            raise NotImplementedError()

    def _rect_grid_to_shape_generator(self):
        """internal method, generator yielding shapely polygons for structured
        grid cells.

        Returns
        -------
        generator :
            generator of shapely Polygons
        """
        for i in range(self.mfgrid.nrow):
            for j in range(self.mfgrid.ncol):
                xy = self.mfgrid.get_cell_vertices(i, j)
                p = Polygon(xy)
                p.name = (i, j)
                yield p

    def _usg_grid_to_shape_generator(self):
        """internal method, convert unstructred grid to list of shapely
        polygons.

        Returns
        -------
        list
            list of shapely Polygons
        """
        raise NotImplementedError()

    def _vtx_grid_to_shape_generator(self):
        """internal method, generator yielding shapely polygons for vertex
        grids.

        Returns
        -------
        generator :
            generator of shapely Polygons
        """
        # for cell2d rec-arrays
        if isinstance(self.mfgrid._cell2d, np.recarray):
            for icell in self.mfgrid._cell2d.icell2d:
                points = []
                icverts = [
                    "icvert_{}".format(i)
                    for i in range(self.mfgrid._cell2d["ncvert"][icell])
                ]
                for iv in self.mfgrid._cell2d[icverts][icell]:
                    points.append(
                        (
                            self.mfgrid._vertices.xv[iv],
                            self.mfgrid._vertices.yv[iv],
                        )
                    )
                # close the polygon, if necessary
                if points[0] != points[-1]:
                    points.append(points[0])
                p = Polygon(points)
                p.name = icell
                yield p
        # for cell2d lists
        elif isinstance(self.mfgrid._cell2d, list):
            for icell in range(len(self.mfgrid._cell2d)):
                points = []
                for iv in self.mfgrid._cell2d[icell][4:]:
                    points.append(
                        (
                            self.mfgrid._vertices[iv][1],
                            self.mfgrid._vertices[iv][2],
                        )
                    )
                # close the polygon, if necessary
                if points[0] != points[-1]:
                    points.append(points[0])
                p = Polygon(points)
                p.name = icell
                yield p

    def _rect_grid_to_shape_list(self):
        """internal method, list of shapely polygons for structured grid cells.

        Returns
        -------
        list :
            list of shapely Polygons
        """
        return list(self._rect_grid_to_shape_generator())

    def _usg_grid_to_shape_list(self):
        """internal method, convert unstructred grid to list of shapely
        polygons.

        Returns
        -------
        list
            list of shapely Polygons
        """
        raise NotImplementedError()

    def _vtx_grid_to_shape_list(self):
        """internal method, list of shapely polygons for vertex grids.

        Returns
        -------
        list :
            list of shapely Polygons
        """
        return list(self._vtx_grid_to_shape_generator())

    def query_grid(self, shp):
        """Perform spatial query on grid with shapely geometry. If no spatial
        query is possible returns all grid cells.

        Parameters
        ----------
        shp : shapely.geometry
            shapely geometry

        Returns
        -------
        list or generator expression
            list or generator containing grid cells in query result
        """
        if self.rtree:
            result = self.strtree.query(shp)
        else:
            # no spatial query
            result = self._get_gridshapes()
        return result

    @staticmethod
    def filter_query_result(qresult, shp):
        """Filter query result to obtain grid cells that intersect with shape.
        Used to (further) reduce query result to cells that definitely
        intersect with shape.

        Parameters
        ----------
        qresult : iterable
            query result, iterable of polygons
        shp : shapely.geometry
            shapely geometry that is prepared and used to filter
            query result

        Returns
        -------
        qfiltered
            filter or generator containing polygons that intersect with shape
        """
        # prepare shape for efficient batch intersection check
        prepshp = prep(shp)
        # get only gridcells that intersect
        qfiltered = filter(prepshp.intersects, qresult)
        return qfiltered

    @staticmethod
    def sort_gridshapes(shape_iter):
        """Sort query result by node id.

        Parameters
        ----------
        shape_iter : iterable
            list or iterable of gridcells

        Returns
        -------
        list
            sorted list of gridcells
        """
        if not isinstance(shape_iter, list):
            shapelist = list(shape_iter)
        else:
            shapelist = shape_iter

        def sort_key(o):
            return o.name

        shapelist.sort(key=sort_key)
        return shapelist

    def _intersect_point_shapely(self, shp, sort_by_cellid=True):
        """intersect grid with Point or MultiPoint.

        Parameters
        ----------
        shp : Point or MultiPoint
            shapely Point or MultiPoint to intersect with grid. Note,
            it is generally faster to loop over a MultiPoint and intersect
            per point than to intersect a MultiPoint directly.
        sort_by_cellid : bool, optional
            flag whether to sort cells by id, used to ensure node
            with lowest id is returned, by default True

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        # query grid
        qresult = self.query_grid(shp)
        # prepare shape for efficient batch intersection check
        prepshp = prep(shp)
        # get only gridcells that intersect
        qfiltered = filter(prepshp.intersects, qresult)

        # sort cells to ensure lowest cell ids are returned
        if sort_by_cellid:
            qfiltered = self.sort_gridshapes(qfiltered)

        isectshp = []
        cellids = []
        vertices = []
        parsed_points = []  # for keeping track of points

        # loop over cells returned by filtered spatial query
        for r in qfiltered:
            name = r.name
            # do intersection
            intersect = shp.intersection(r)
            # parse result per Point
            collection = parse_shapely_ix_result(
                [], intersect, shptyps=["Point"]
            )
            # loop over intersection result and store information
            cell_verts = []
            cell_shps = []
            for c in collection:
                verts = c.__geo_interface__["coordinates"]
                # avoid returning multiple cells for points on boundaries
                if verts in parsed_points:
                    continue
                parsed_points.append(verts)
                cell_shps.append(c)  # collect only new points
                cell_verts.append(verts)
            # if any new ix found
            if len(cell_shps) > 0:
                # combine new points in MultiPoint
                isectshp.append(
                    MultiPoint(cell_shps)
                    if len(cell_shps) > 1
                    else cell_shps[0]
                )
                vertices.append(tuple(cell_verts))
                cellids.append(name)

        rec = np.recarray(
            len(isectshp),
            names=["cellids", "vertices", "ixshapes"],
            formats=["O", "O", "O"],
        )
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = isectshp
        rec.vertices = vertices
        rec.cellids = cellids

        return rec

    def _intersect_linestring_shapely(
        self, shp, keepzerolengths=False, sort_by_cellid=True
    ):
        """intersect with LineString or MultiLineString.

        Parameters
        ----------
        shp : shapely.geometry.LineString or MultiLineString
            LineString to intersect with the grid
        keepzerolengths : bool, optional
            keep linestrings with length zero, default is False
        sort_by_cellid : bool, optional
            flag whether to sort cells by id, used to ensure node
            with lowest id is returned, by default True

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        # query grid
        qresult = self.query_grid(shp)
        # filter result further if possible (only strtree and filter methods)
        qfiltered = self.filter_query_result(qresult, shp)
        # sort cells to ensure lowest cell ids are returned
        if sort_by_cellid:
            qfiltered = self.sort_gridshapes(qfiltered)

        # initialize empty lists for storing results
        isectshp = []
        cellids = []
        vertices = []
        lengths = []

        # loop over cells returned by filtered spatial query
        for r in qfiltered:
            name = r.name
            # do intersection
            intersect = shp.intersection(r)
            # parse result
            collection = parse_shapely_ix_result(
                [], intersect, shptyps=["LineString", "MultiLineString"]
            )
            # loop over intersection result and store information
            for c in collection:
                verts = c.__geo_interface__["coordinates"]
                # test if linestring was already processed (if on boundary)
                if verts in vertices:
                    continue
                # if keep zero don't check length
                if not keepzerolengths:
                    if c.length == 0.0:
                        continue
                isectshp.append(c)
                lengths.append(c.length)
                vertices.append(verts)
                cellids.append(name)

        rec = np.recarray(
            len(isectshp),
            names=["cellids", "vertices", "lengths", "ixshapes"],
            formats=["O", "O", "f8", "O"],
        )
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = isectshp
        rec.vertices = vertices
        rec.lengths = lengths
        rec.cellids = cellids

        return rec

    def _intersect_polygon_shapely(self, shp, sort_by_cellid=True):
        """intersect with Polygon or MultiPolygon.

        Parameters
        ----------
        shp : shapely.geometry.Polygon or MultiPolygon
            shape to intersect with the grid
        sort_by_cellid : bool, optional
            flag whether to sort cells by id, used to ensure node
            with lowest id is returned, by default True

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        # query grid
        qresult = self.query_grid(shp)
        # filter result further if possible (only strtree and filter methods)
        qfiltered = self.filter_query_result(qresult, shp)
        # sort cells to ensure lowest cell ids are returned
        if sort_by_cellid:
            qfiltered = self.sort_gridshapes(qfiltered)

        isectshp = []
        cellids = []
        vertices = []
        areas = []

        # loop over cells returned by filtered spatial query
        for r in qfiltered:
            name = r.name
            # do intersection
            intersect = shp.intersection(r)
            # parse result
            collection = parse_shapely_ix_result(
                [], intersect, shptyps=["Polygon", "MultiPolygon"]
            )
            if len(collection) > 1:
                collection = [MultiPolygon(collection)]
            # loop over intersection result and store information
            for c in collection:
                # don't store intersections with 0 area
                if c.area == 0.0:
                    continue
                verts = c.__geo_interface__["coordinates"]
                isectshp.append(c)
                areas.append(c.area)
                vertices.append(verts)
                cellids.append(name)

        rec = np.recarray(
            len(isectshp),
            names=["cellids", "vertices", "areas", "ixshapes"],
            formats=["O", "O", "f8", "O"],
        )
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = isectshp
        rec.vertices = vertices
        rec.areas = areas
        rec.cellids = cellids

        return rec

    def intersects(self, shp):
        """Return cellIDs for shapes that intersect with shape.

        Parameters
        ----------
        shp : shapely.geometry, geojson geometry, shapefile.shape,
              or flopy geometry object
            shape to intersect with the grid

        Returns
        -------
        rec : numpy.recarray
            a record array containing cell IDs of the gridcells
            the shape intersects with
        """
        # query grid
        shp = GeoSpatialUtil(shp).shapely

        qresult = self.query_grid(shp)
        # filter result further if possible (only strtree and filter methods)
        qfiltered = self.filter_query_result(qresult, shp)
        # get cellids
        cids = [cell.name for cell in qfiltered]
        # build rec-array
        rec = np.recarray(len(cids), names=["cellids"], formats=["O"])
        rec.cellids = cids
        return rec

    def _intersect_point_structured(self, shp):
        """intersection method for intersecting points with structured grids.

        Parameters
        ----------
        shp : shapely.geometry.Point or MultiPoint
            point shape to intersect with grid

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        nodelist = []

        Xe, Ye = self.mfgrid.xyedges

        if isinstance(shp, Point):
            shp = [shp]
        elif isinstance(shp, MultiPoint):
            shp = list(shp.geoms)
        else:
            raise ValueError("expected Point or MultiPoint")

        ixshapes = []
        for p in shp:
            # if grid is rotated or offset transform point to local coords
            if (
                self.mfgrid.angrot != 0.0
                or self.mfgrid.xoffset != 0.0
                or self.mfgrid.yoffset != 0.0
            ):
                rx, ry = transform(
                    p.x,
                    p.y,
                    self.mfgrid.xoffset,
                    self.mfgrid.yoffset,
                    self.mfgrid.angrot_radians,
                    inverse=True,
                )
            else:
                rx = p.x
                ry = p.y

            # two dimensional point
            jpos = ModflowGridIndices.find_position_in_array(Xe, rx)
            ipos = ModflowGridIndices.find_position_in_array(Ye, ry)

            if jpos is not None and ipos is not None:
                # three dimensional point
                if p._ndim == 3:
                    # find k
                    kpos = ModflowGridIndices.find_position_in_array(
                        self.mfgrid.botm[:, ipos, jpos], p.z
                    )
                    if kpos is not None:
                        nodelist.append((kpos, ipos, jpos))
                        ixshapes.append(p)
                else:
                    nodelist.append((ipos, jpos))
                    ixshapes.append(p)

        # remove duplicates
        tempnodes = []
        tempshapes = []
        for node, ixs in zip(nodelist, ixshapes):
            if node not in tempnodes:
                tempnodes.append(node)
                tempshapes.append(ixs)
            else:
                # TODO: not sure if this is correct
                tempshapes[-1] = MultiPoint([tempshapes[-1], ixs])

        ixshapes = tempshapes
        nodelist = tempnodes

        rec = np.recarray(
            len(nodelist), names=["cellids", "ixshapes"], formats=["O", "O"]
        )
        rec.cellids = nodelist
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = ixshapes
        return rec

    def _intersect_linestring_structured(self, shp, keepzerolengths=False):
        """method for intersecting linestrings with structured grids.

        Parameters
        ----------
        shp : shapely.geometry.Linestring or MultiLineString
            linestring to intersect with grid
        keepzerolengths : bool, optional
            if True keep intersection results with length=0, in
            other words, grid cells the linestring does not cross
            but does touch, by default False

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """
        # get local extent of grid
        if (
            self.mfgrid.angrot != 0.0
            or self.mfgrid.xoffset != 0.0
            or self.mfgrid.yoffset != 0.0
        ):
            xmin = np.min(self.mfgrid.xyedges[0])
            xmax = np.max(self.mfgrid.xyedges[0])
            ymin = np.min(self.mfgrid.xyedges[1])
            ymax = np.max(self.mfgrid.xyedges[1])
        else:
            xmin, xmax, ymin, ymax = self.mfgrid.extent
        pl = box(xmin, ymin, xmax, ymax)

        # rotate and translate linestring to local coords
        if self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0:
            shp = translate(
                shp, xoff=-self.mfgrid.xoffset, yoff=-self.mfgrid.yoffset
            )
        if self.mfgrid.angrot != 0.0:
            shp = rotate(shp, -self.mfgrid.angrot, origin=(0.0, 0.0))

        # clip line to mfgrid bbox
        lineclip = shp.intersection(pl)

        if lineclip.length == 0.0:  # linestring does not intersect modelgrid
            return np.recarray(
                0,
                names=["cellids", "vertices", "lengths", "ixshapes"],
                formats=["O", "O", "f8", "O"],
            )
        if lineclip.geom_type == "MultiLineString":  # there are multiple lines
            nodelist, lengths, vertices = [], [], []
            ixshapes = []
            for ls in lineclip.geoms:
                n, l, v, ixs = self._get_nodes_intersecting_linestring(ls)
                nodelist += n
                lengths += l
                # if necessary, transform coordinates back to real
                # world coordinates
                if (
                    self.mfgrid.angrot != 0.0
                    or self.mfgrid.xoffset != 0.0
                    or self.mfgrid.yoffset != 0.0
                ):
                    v_realworld = []
                    for pt in v:
                        pt = np.array(pt)
                        rx, ry = transform(
                            pt[:, 0],
                            pt[:, 1],
                            self.mfgrid.xoffset,
                            self.mfgrid.yoffset,
                            self.mfgrid.angrot_radians,
                            inverse=False,
                        )
                        v_realworld.append(list(zip(rx, ry)))
                    ixs_realworld = []
                    for ix in ixs:
                        ix_realworld = rotate(
                            ix, self.mfgrid.angrot, origin=(0.0, 0.0)
                        )
                        ix_realworld = translate(
                            ix_realworld,
                            self.mfgrid.xoffset,
                            self.mfgrid.yoffset,
                        )
                        ixs_realworld.append(ix_realworld)
                else:
                    v_realworld = v
                    ixs_realworld = ixs
                vertices += v_realworld
                ixshapes += ixs_realworld
        else:  # linestring is fully within grid
            (
                nodelist,
                lengths,
                vertices,
                ixshapes,
            ) = self._get_nodes_intersecting_linestring(lineclip)
            # if necessary, transform coordinates back to real
            # world coordinates
            if (
                self.mfgrid.angrot != 0.0
                or self.mfgrid.xoffset != 0.0
                or self.mfgrid.yoffset != 0.0
            ):
                v_realworld = []
                for pt in vertices:
                    pt = np.array(pt)
                    rx, ry = transform(
                        pt[:, 0],
                        pt[:, 1],
                        self.mfgrid.xoffset,
                        self.mfgrid.yoffset,
                        self.mfgrid.angrot_radians,
                        inverse=False,
                    )
                    v_realworld.append(list(zip(rx, ry)))
                vertices = v_realworld

                ix_shapes_realworld = []
                for ixs in ixshapes:
                    ixs = rotate(ixs, self.mfgrid.angrot, origin=(0.0, 0.0))
                    ixs = translate(
                        ixs, self.mfgrid.xoffset, self.mfgrid.yoffset
                    )
                    ix_shapes_realworld.append(ixs)
                ixshapes = ix_shapes_realworld

        # bundle linestrings in same cell
        tempnodes = []
        templengths = []
        tempverts = []
        tempshapes = []
        unique_nodes = list(set(nodelist))
        if len(unique_nodes) < len(nodelist):
            for inode in unique_nodes:
                templengths.append(
                    sum([l for l, i in zip(lengths, nodelist) if i == inode])
                )
                tempverts.append(
                    [v for v, i in zip(vertices, nodelist) if i == inode]
                )
                tempshapes.append(
                    [ix for ix, i in zip(ixshapes, nodelist) if i == inode]
                )

            nodelist = unique_nodes
            lengths = templengths
            vertices = tempverts
            ixshapes = tempshapes

        # eliminate any nodes that have a zero length
        if not keepzerolengths:
            tempnodes = []
            templengths = []
            tempverts = []
            tempshapes = []
            for i, _ in enumerate(nodelist):
                if lengths[i] > 0:
                    tempnodes.append(nodelist[i])
                    templengths.append(lengths[i])
                    tempverts.append(vertices[i])
                    tempshapes.append(ixshapes[i])
            nodelist = tempnodes
            lengths = templengths
            vertices = tempverts
            ixshapes = tempshapes

        rec = np.recarray(
            len(nodelist),
            names=["cellids", "vertices", "lengths", "ixshapes"],
            formats=["O", "O", "f8", "O"],
        )
        rec.vertices = vertices
        rec.lengths = lengths
        rec.cellids = nodelist
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = ixshapes

        return rec

    def _get_nodes_intersecting_linestring(self, linestring):
        """helper function, intersect the linestring with the a structured grid
        and return a list of node indices and the length of the line in that
        node.

        Parameters
        ----------
        linestring: shapely.geometry.LineString or MultiLineString
            shape to intersect with the grid

        Returns
        -------
        nodelist, lengths, vertices: lists
            lists containing node ids, lengths of intersects and the
            start and end points of the intersects
        """
        nodelist = []
        lengths = []
        vertices = []
        ixshapes = []

        # start at the beginning of the line
        x, y = linestring.xy

        # linestring already in local coords but
        # because intersect_point does transform again
        # we transform back to real world here if necessary
        if (
            self.mfgrid.angrot != 0.0
            or self.mfgrid.xoffset != 0.0
            or self.mfgrid.yoffset != 0.0
        ):
            x0, y0 = transform(
                [x[0]],
                [y[0]],
                self.mfgrid.xoffset,
                self.mfgrid.yoffset,
                self.mfgrid.angrot_radians,
                inverse=False,
            )
        else:
            x0 = [x[0]]
            y0 = [y[0]]

        (i, j) = self.intersect(Point(x0[0], y0[0])).cellids[0]
        Xe, Ye = self.mfgrid.xyedges
        xmin = Xe[j]
        xmax = Xe[j + 1]
        ymax = Ye[i]
        ymin = Ye[i + 1]
        pl = box(xmin, ymin, xmax, ymax)
        intersect = linestring.intersection(pl)
        # if linestring starts in cell, exits, and re-enters
        # a MultiLineString is returned.
        ixshapes.append(intersect)
        length = intersect.length
        lengths.append(length)
        if hasattr(intersect, "geoms"):
            x, y = [], []
            for igeom in intersect.geoms:
                x.append(igeom.xy[0])
                y.append(igeom.xy[1])
            x = np.concatenate(x)
            y = np.concatenate(y)
        else:
            x = intersect.xy[0]
            y = intersect.xy[1]
        verts = [(ixy[0], ixy[1]) for ixy in zip(x, y)]
        vertices.append(verts)
        nodelist.append((i, j))

        n = 0
        while True:
            (i, j) = nodelist[n]
            (
                node,
                length,
                verts,
                ixshape,
            ) = self._check_adjacent_cells_intersecting_line(
                linestring, (i, j), nodelist
            )

            for inode, ilength, ivert, ix in zip(node, length, verts, ixshape):
                if inode is not None:
                    if ivert not in vertices:
                        nodelist.append(inode)
                        lengths.append(ilength)
                        vertices.append(ivert)
                        ixshapes.append(ix)

            if n == len(nodelist) - 1:
                break
            n += 1

        return nodelist, lengths, vertices, ixshapes

    def _check_adjacent_cells_intersecting_line(
        self, linestring, i_j, nodelist
    ):
        """helper method that follows a line through a structured grid.

        Parameters
        ----------
        linestring : shapely.geometry.LineString
            shape to intersect with the grid
        i_j : tuple
            tuple containing (nrow, ncol)
        nodelist : list of tuples
            list of node ids that have already been added
            as intersections

        Returns
        -------
        node, length, verts: lists
            lists containing nodes, lengths and vertices of
            intersections with adjacent cells relative to the
            current cell (i, j)
        """
        i, j = i_j

        Xe, Ye = self.mfgrid.xyedges

        node = []
        length = []
        verts = []
        ixshape = []

        # check to left
        if j > 0:
            ii = i
            jj = j - 1
            if (ii, jj) not in nodelist:
                xmin = Xe[jj]
                xmax = Xe[jj + 1]
                ymax = Ye[ii]
                ymin = Ye[ii + 1]
                pl = box(xmin, ymin, xmax, ymax)
                if linestring.intersects(pl):
                    intersect = linestring.intersection(pl)
                    ixshape.append(intersect)
                    length.append(intersect.length)
                    if hasattr(intersect, "geoms"):
                        x, y = [], []
                        for igeom in intersect.geoms:
                            x.append(igeom.xy[0])
                            y.append(igeom.xy[1])
                        x = np.concatenate(x)
                        y = np.concatenate(y)
                    else:
                        x = intersect.xy[0]
                        y = intersect.xy[1]
                    verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)])
                    node.append((ii, jj))

        # check to right
        if j < self.mfgrid.ncol - 1:
            ii = i
            jj = j + 1
            if (ii, jj) not in nodelist:
                xmin = Xe[jj]
                xmax = Xe[jj + 1]
                ymax = Ye[ii]
                ymin = Ye[ii + 1]
                pl = box(xmin, ymin, xmax, ymax)
                if linestring.intersects(pl):
                    intersect = linestring.intersection(pl)
                    ixshape.append(intersect)
                    length.append(intersect.length)
                    if hasattr(intersect, "geoms"):
                        x, y = [], []
                        for igeom in intersect.geoms:
                            x.append(igeom.xy[0])
                            y.append(igeom.xy[1])
                        x = np.concatenate(x)
                        y = np.concatenate(y)
                    else:
                        x = intersect.xy[0]
                        y = intersect.xy[1]
                    verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)])
                    node.append((ii, jj))

        # check to back
        if i > 0:
            ii = i - 1
            jj = j
            if (ii, jj) not in nodelist:
                xmin = Xe[jj]
                xmax = Xe[jj + 1]
                ymax = Ye[ii]
                ymin = Ye[ii + 1]
                pl = box(xmin, ymin, xmax, ymax)
                if linestring.intersects(pl):
                    intersect = linestring.intersection(pl)
                    ixshape.append(intersect)
                    length.append(intersect.length)
                    if hasattr(intersect, "geoms"):
                        x, y = [], []
                        for igeom in intersect.geoms:
                            x.append(igeom.xy[0])
                            y.append(igeom.xy[1])
                        x = np.concatenate(x)
                        y = np.concatenate(y)
                    else:
                        x = intersect.xy[0]
                        y = intersect.xy[1]
                    verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)])
                    node.append((ii, jj))

        # check to front
        if i < self.mfgrid.nrow - 1:
            ii = i + 1
            jj = j
            if (ii, jj) not in nodelist:
                xmin = Xe[jj]
                xmax = Xe[jj + 1]
                ymax = Ye[ii]
                ymin = Ye[ii + 1]
                pl = box(xmin, ymin, xmax, ymax)
                if linestring.intersects(pl):
                    intersect = linestring.intersection(pl)
                    ixshape.append(intersect)
                    length.append(intersect.length)
                    if hasattr(intersect, "geoms"):
                        x, y = [], []
                        for igeom in intersect.geoms:
                            x.append(igeom.xy[0])
                            y.append(igeom.xy[1])
                        x = np.concatenate(x)
                        y = np.concatenate(y)
                    else:
                        x = intersect.xy[0]
                        y = intersect.xy[1]
                    verts.append([(ixy[0], ixy[1]) for ixy in zip(x, y)])
                    node.append((ii, jj))

        return node, length, verts, ixshape

    def _intersect_rectangle_structured(self, rectangle):
        """intersect a rectangle with a structured grid to retrieve node ids of
        intersecting grid cells.

        Note: only works in local coordinates (i.e. non-rotated grid
        with origin at (0, 0))

        Parameters
        ----------
        rectangle : list of tuples
            list of lower-left coordinate and upper-right
            coordinate: [(xmin, ymin), (xmax, ymax)]

        Returns
        -------
        nodelist: list of tuples
            list of tuples containing node ids with which
            the rectangle intersects
        """

        nodelist = []

        # return if rectangle does not contain any cells
        if (
            self.mfgrid.angrot != 0.0
            or self.mfgrid.xoffset != 0.0
            or self.mfgrid.yoffset != 0.0
        ):
            minx = np.min(self.mfgrid.xyedges[0])
            maxx = np.max(self.mfgrid.xyedges[0])
            miny = np.min(self.mfgrid.xyedges[1])
            maxy = np.max(self.mfgrid.xyedges[1])
            local_extent = [minx, maxx, miny, maxy]
        else:
            local_extent = self.mfgrid.extent

        xmin, xmax, ymin, ymax = local_extent
        bgrid = box(xmin, ymin, xmax, ymax)
        (rxmin, rymin), (rxmax, rymax) = rectangle
        b = box(rxmin, rymin, rxmax, rymax)

        if not b.intersects(bgrid):
            # return with nodelist as an empty list
            return []

        Xe, Ye = self.mfgrid.xyedges

        jmin = ModflowGridIndices.find_position_in_array(Xe, xmin)
        if jmin is None:
            if xmin <= Xe[0]:
                jmin = 0
            elif xmin >= Xe[-1]:
                jmin = self.mfgrid.ncol - 1

        jmax = ModflowGridIndices.find_position_in_array(Xe, xmax)
        if jmax is None:
            if xmax <= Xe[0]:
                jmax = 0
            elif xmax >= Xe[-1]:
                jmax = self.mfgrid.ncol - 1

        imin = ModflowGridIndices.find_position_in_array(Ye, ymax)
        if imin is None:
            if ymax >= Ye[0]:
                imin = 0
            elif ymax <= Ye[-1]:
                imin = self.mfgrid.nrow - 1

        imax = ModflowGridIndices.find_position_in_array(Ye, ymin)
        if imax is None:
            if ymin >= Ye[0]:
                imax = 0
            elif ymin <= Ye[-1]:
                imax = self.mfgrid.nrow - 1

        for i in range(imin, imax + 1):
            for j in range(jmin, jmax + 1):
                nodelist.append((i, j))

        return nodelist

    def _intersect_polygon_structured(self, shp):
        """intersect polygon with a structured grid. Uses bounding box of the
        Polygon to limit search space.

        Notes
        -----
        If performance is slow, try setting the method to 'vertex'
        in the GridIntersect object. For polygons this is often
        faster.

        Parameters
        ----------
        shp : shapely.geometry.Polygon
            polygon to intersect with the grid

        Returns
        -------
        numpy.recarray
            a record array containing information about the intersection
        """

        # initialize the result lists
        nodelist = []
        areas = []
        vertices = []
        ixshapes = []

        # transform polygon to local grid coordinates
        if self.mfgrid.xoffset != 0.0 or self.mfgrid.yoffset != 0.0:
            shp = translate(
                shp, xoff=-self.mfgrid.xoffset, yoff=-self.mfgrid.yoffset
            )
        if self.mfgrid.angrot != 0.0:
            shp = rotate(shp, -self.mfgrid.angrot, origin=(0.0, 0.0))

        # use the bounds of the polygon to restrict the cell search
        minx, miny, maxx, maxy = shp.bounds
        rectangle = ((minx, miny), (maxx, maxy))
        nodes = self._intersect_rectangle_structured(rectangle)

        for (i, j) in nodes:
            if (
                self.mfgrid.angrot != 0.0
                or self.mfgrid.xoffset != 0.0
                or self.mfgrid.yoffset != 0.0
            ):
                cell_coords = [
                    (self.mfgrid.xyedges[0][j], self.mfgrid.xyedges[1][i]),
                    (self.mfgrid.xyedges[0][j + 1], self.mfgrid.xyedges[1][i]),
                    (
                        self.mfgrid.xyedges[0][j + 1],
                        self.mfgrid.xyedges[1][i + 1],
                    ),
                    (self.mfgrid.xyedges[0][j], self.mfgrid.xyedges[1][i + 1]),
                ]
            else:
                cell_coords = self.mfgrid.get_cell_vertices(i, j)
            node_polygon = Polygon(cell_coords)
            if shp.intersects(node_polygon):
                intersect = shp.intersection(node_polygon)
                if intersect.area > 0.0:
                    nodelist.append((i, j))
                    areas.append(intersect.area)

                    # if necessary, transform coordinates back to real
                    # world coordinates
                    if (
                        self.mfgrid.angrot != 0.0
                        or self.mfgrid.xoffset != 0.0
                        or self.mfgrid.yoffset != 0.0
                    ):
                        v_realworld = []
                        if intersect.geom_type.startswith("Multi"):
                            for ipoly in intersect:
                                v_realworld += (
                                    self._transform_geo_interface_polygon(
                                        ipoly
                                    )
                                )
                        else:
                            v_realworld += (
                                self._transform_geo_interface_polygon(
                                    intersect
                                )
                            )
                        intersect_realworld = rotate(
                            intersect, self.mfgrid.angrot, origin=(0.0, 0.0)
                        )
                        intersect_realworld = translate(
                            intersect_realworld,
                            self.mfgrid.xoffset,
                            self.mfgrid.yoffset,
                        )
                    else:
                        v_realworld = intersect.__geo_interface__[
                            "coordinates"
                        ]
                        intersect_realworld = intersect
                    ixshapes.append(intersect_realworld)
                    vertices.append(v_realworld)

        rec = np.recarray(
            len(nodelist),
            names=["cellids", "vertices", "areas", "ixshapes"],
            formats=["O", "O", "f8", "O"],
        )
        rec.vertices = vertices
        rec.areas = areas
        rec.cellids = nodelist
        with ignore_shapely_warnings_for_object_array():
            rec.ixshapes = ixshapes

        return rec

    def _transform_geo_interface_polygon(self, polygon):
        """Internal method, helper function to transform geometry
        __geo_interface__.

        Used for translating intersection result coordinates back into
        real-world coordinates.

        Parameters
        ----------
        polygon : shapely.geometry.Polygon
            polygon to transform coordinates for

        Returns
        -------
        geom_list : list
            list containing transformed coordinates in same structure as
            the original __geo_interface__.
        """

        if polygon.geom_type.startswith("Multi"):
            raise TypeError("Does not support Multi geometries!")

        geom_list = []
        for coords in polygon.__geo_interface__["coordinates"]:
            geoms = []
            try:
                # test depth of list/tuple
                _ = coords[0][0][0]
                if len(coords) == 2:
                    shell, holes = coords
                else:
                    raise ValueError("Cannot parse __geo_interface__")
            except TypeError:
                shell = coords
                holes = None
            except Exception as e:
                raise e
            # transform shell coordinates
            shell_pts = []
            for pt in shell:
                rx, ry = transform(
                    [pt[0]],
                    [pt[1]],
                    self.mfgrid.xoffset,
                    self.mfgrid.yoffset,
                    self.mfgrid.angrot_radians,
                    inverse=False,
                )
                shell_pts.append((rx, ry))
            geoms.append(shell_pts)
            # transform holes coordinates if necessary
            if holes:
                holes_pts = []
                for pt in holes:
                    rx, ry = transform(
                        [pt[0]],
                        [pt[1]],
                        self.mfgrid.xoffset,
                        self.mfgrid.yoffset,
                        self.mfgrid.angrot_radians,
                        inverse=False,
                    )
            # append (shells, holes) to transformed coordinates list
            geom_list.append(tuple(geoms))
        return geom_list

    @staticmethod
    def plot_polygon(rec, ax=None, **kwargs):
        """method to plot the polygon intersection results from the resulting
        numpy.recarray.

        Note: only works when recarray has 'intersects' column!

        Parameters
        ----------
        rec : numpy.recarray
            record array containing intersection results
            (the resulting shapes)
        ax : matplotlib.pyplot.axes, optional
            axes to plot onto, if not provided, creates a new figure
        **kwargs:
            passed to the plot function

        Returns
        -------
        ax: matplotlib.pyplot.axes
            returns the axes handle
        """
        try:
            from descartes import PolygonPatch
        except ImportError:
            msg = "descartes package needed for plotting polygons"
            if plt is None:
                msg = (
                    "matplotlib and descartes packages needed for "
                    "plotting polygons"
                )
            raise ImportError(msg)

        if plt is None:
            msg = "matplotlib package needed for plotting polygons"
            raise ImportError(msg)

        if ax is None:
            _, ax = plt.subplots()

        for i, ishp in enumerate(rec.ixshapes):
            if "facecolor" in kwargs:
                fc = kwargs.pop("facecolor")
            else:
                fc = "C{}".format(i % 10)
            ppi = PolygonPatch(ishp, facecolor=fc, **kwargs)
            ax.add_patch(ppi)

        return ax

    @staticmethod
    def plot_linestring(rec, ax=None, **kwargs):
        """method to plot the linestring intersection results from the
        resulting numpy.recarray.

        Note: only works when recarray has 'intersects' column!

        Parameters
        ----------
        rec : numpy.recarray
            record array containing intersection results
            (the resulting shapes)
        ax : matplotlib.pyplot.axes, optional
            axes to plot onto, if not provided, creates a new figure
        **kwargs:
            passed to the plot function

        Returns
        -------
        ax: matplotlib.pyplot.axes
            returns the axes handle
        """
        if plt is None:
            msg = "matplotlib package needed for plotting polygons"
            raise ImportError(msg)

        if ax is None:
            _, ax = plt.subplots()

        for i, ishp in enumerate(rec.ixshapes):
            if "c" in kwargs:
                c = kwargs.pop("c")
            elif "color" in kwargs:
                c = kwargs.pop("color")
            else:
                c = "C{}".format(i % 10)
            if ishp.type == "MultiLineString":
                for part in ishp:
                    ax.plot(part.xy[0], part.xy[1], ls="-", c=c, **kwargs)
            else:
                ax.plot(ishp.xy[0], ishp.xy[1], ls="-", c=c, **kwargs)

        return ax

    @staticmethod
    def plot_point(rec, ax=None, **kwargs):
        """method to plot the point intersection results from the resulting
        numpy.recarray.

        Note: only works when recarray has 'intersects' column!

        Parameters
        ----------
        rec : numpy.recarray
            record array containing intersection results
        ax : matplotlib.pyplot.axes, optional
            axes to plot onto, if not provided, creates a new figure
        **kwargs:
            passed to the scatter function

        Returns
        -------
        ax: matplotlib.pyplot.axes
            returns the axes handle
        """
        if plt is None:
            msg = "matplotlib package needed for plotting polygons"
            raise ImportError(msg)

        if ax is None:
            _, ax = plt.subplots()

        x, y = [], []
        geo_coll = GeometryCollection(list(rec.ixshapes))
        collection = parse_shapely_ix_result([], geo_coll, ["Point"])
        for c in collection:
            x.append(c.x)
            y.append(c.y)
        ax.scatter(x, y, **kwargs)

        return ax


class ModflowGridIndices:
    """Collection of methods that can be used to find cell indices for a
    structured, but irregularly spaced MODFLOW grid."""

    @staticmethod
    def find_position_in_array(arr, x):
        """If arr has x positions for the left edge of a cell, then return the
        cell index containing x.

        Parameters
        ----------
        arr : A one dimensional array (such as Xe) that contains
            coordinates for the left cell edge.

        x : float
            The x position to find in arr.
        """
        jpos = None

        if x == arr[-1]:
            return len(arr) - 2

        if x < min(arr[0], arr[-1]):
            return None

        if x > max(arr[0], arr[-1]):
            return None

        # go through each position
        for j in range(len(arr) - 1):
            xl = arr[j]
            xr = arr[j + 1]
            frac = (x - xl) / (xr - xl)
            if 0.0 <= frac <= 1.0:
                # if min(xl, xr) <= x < max(xl, xr):
                jpos = j
                return jpos

        return jpos

    @staticmethod
    def kij_from_nodenumber(nodenumber, nlay, nrow, ncol):
        """Convert the modflow node number to a zero-based layer, row and
        column format.  Return (k0, i0, j0).

        Parameters
        ----------
        nodenumber: int
            The cell nodenumber, ranging from 1 to number of
            nodes.
        nlay: int
            The number of layers.
        nrow: int
            The number of rows.
        ncol: int
            The number of columns.
        """
        if nodenumber > nlay * nrow * ncol:
            raise Exception("Error in function kij_from_nodenumber...")
        n = nodenumber - 1
        k = int(n / nrow / ncol)
        i = int((n - k * nrow * ncol) / ncol)
        j = n - k * nrow * ncol - i * ncol
        return (k, i, j)

    @staticmethod
    def nodenumber_from_kij(k, i, j, nrow, ncol):
        """Calculate the nodenumber using the zero-based layer, row, and column
        values.  The first node has a value of 1.

        Parameters
        ----------
        k : int
            The model layer number as a zero-based value.
        i : int
            The model row number as a zero-based value.
        j : int
            The model column number as a zero-based value.
        nrow : int
            The number of model rows.
        ncol : int
            The number of model columns.
        """
        return k * nrow * ncol + i * ncol + j + 1

    @staticmethod
    def nn0_from_kij(k, i, j, nrow, ncol):
        """Calculate the zero-based nodenumber using the zero-based layer, row,
        and column values.  The first node has a value of 0.

        Parameters
        ----------
        k : int
            The model layer number as a zero-based value.
        i : int
            The model row number as a zero-based value.
        j : int
            The model column number as a zero-based value.
        nrow : int
            The number of model rows.
        ncol : int
            The number of model columns.
        """
        return k * nrow * ncol + i * ncol + j

    @staticmethod
    def kij_from_nn0(n, nlay, nrow, ncol):
        """Convert the node number to a zero-based layer, row and column
        format.  Return (k0, i0, j0).

        Parameters
        ----------
        nodenumber : int
            The cell nodenumber, ranging from 0 to number of
            nodes - 1.
        nlay : int
            The number of layers.
        nrow : int
            The number of rows.
        ncol : int
            The number of columns.
        """
        if n > nlay * nrow * ncol:
            raise Exception("Error in function kij_from_nodenumber...")
        k = int(n / nrow / ncol)
        i = int((n - k * nrow * ncol) / ncol)
        j = n - k * nrow * ncol - i * ncol
        return (k, i, j)