import numpy as np
from ..plot.crosssection import _StructuredCrossSection
from ..plot.vcrosssection import _VertexCrossSection
from ..plot import plotutil

try:
    import matplotlib.pyplot as plt
    import matplotlib.colors
except ImportError:
    plt = None


class PlotCrossSection(object):
    """
    Class to create a cross section of the model.

    Parameters
    ----------
    ax : matplotlib.pyplot axis
        The plot axis.  If not provided it, plt.gca() will be used.
    model : flopy.modflow object
        flopy model object. (Default is None)
    modelgrid : flopy.discretization.Grid object
        can be a StructuredGrid, VertexGrid, or UnstructuredGrid object
    line : dict
        Dictionary with either "row", "column", or "line" key. If key
        is "row" or "column" key value should be the zero-based row or
        column index for cross-section. If key is "line" value should
        be an array of (x, y) tuples with vertices of cross-section.
        Vertices should be in map coordinates consistent with xul,
        yul, and rotation.
    extent : tuple of floats
        (xmin, xmax, ymin, ymax) will be used to specify axes limits.  If None
        then these will be calculated based on grid, coordinates, and rotation.
    geographic_coords : bool
        boolean flag to allow the user to plot cross section lines in
        geographic coordinates. If False (default), cross section is plotted
        as the distance along the cross section line.

    """

    def __init__(self, model=None, modelgrid=None, ax=None,
                 line=None, extent=None, geographic_coords=False):
        if plt is None:
            s = 'Could not import matplotlib.  Must install matplotlib ' + \
                ' in order to use ModelMap method'
            raise ImportError(s)

        if modelgrid is None and model is not None:
            modelgrid = model.modelgrid

        # update this after unstructured grid is finished!
        tmp = modelgrid.grid_type

        if tmp == "structured":
            self.__cls = _StructuredCrossSection(ax=ax, model=model,
                                                 modelgrid=modelgrid,
                                                 line=line, extent=extent,
                                                 geographic_coords=
                                                 geographic_coords)

        elif tmp == "unstructured":
            raise NotImplementedError("Unstructured xc not yet implemented")

        elif tmp == "vertex":
            self.__cls = _VertexCrossSection(ax=ax, model=model,
                                             modelgrid=modelgrid,
                                             line=line, extent=extent,
                                             geographic_coords=
                                             geographic_coords)

        else:
            raise ValueError("Unknown modelgrid type {}".format(tmp))

        self.model = self.__cls.model
        self.mg = self.__cls.mg
        self.ax = self.__cls.ax
        self.direction = self.__cls.direction
        self.pts = self.__cls.pts
        self.xpts = self.__cls.xpts
        self.d = self.__cls.d
        self.ncb = self.__cls.ncb
        self.laycbd = self.__cls.laycbd
        self.active = self.__cls.active
        self.elev = self.__cls.elev
        self.layer0 = self.__cls.layer0
        self.layer1 = self.__cls.layer1
        self.zpts = self.__cls.zpts
        self.xcentergrid = self.__cls.xcentergrid
        self.zcentergrid = self.__cls.zcentergrid
        self.geographic_coords = self.__cls.geographic_coords
        self.extent = self.__cls.extent

    def plot_array(self, a, masked_values=None, head=None, **kwargs):
        """
        Plot a three-dimensional array as a patch collection.

        Parameters
        ----------
        a : numpy.ndarray
            Three-dimensional array to plot.
        masked_values : iterable of floats, ints
            Values to mask.
        head : numpy.ndarray
            Three-dimensional array to set top of patches to the minimum
            of the top of a layer or the head value. Used to create
            patches that conform to water-level elevations.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.collections.PatchCollection

        Returns
        -------
        patches : matplotlib.collections.PatchCollection

        """
        return self.__cls.plot_array(a=a, masked_values=masked_values,
                                     head=head, **kwargs)

    def plot_surface(self, a, masked_values=None, **kwargs):
        """
        Plot a two- or three-dimensional array as line(s).

        Parameters
        ----------
        a : numpy.ndarray
            Two- or three-dimensional array to plot.
        masked_values : iterable of floats, ints
            Values to mask.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.pyplot.plot

        Returns
        -------
        plot : list containing matplotlib.plot objects

        """
        return self.__cls.plot_surface(a=a, masked_values=masked_values,
                                       **kwargs)

    def plot_fill_between(self, a, colors=('blue', 'red'),
                          masked_values=None, head=None, **kwargs):
        """
        Plot a three-dimensional array as lines.

        Parameters
        ----------
        a : numpy.ndarray
            Three-dimensional array to plot.
        colors: list
            matplotlib fill colors, two required
        masked_values : iterable of floats, ints
            Values to mask.
        head : numpy.ndarray
            Three-dimensional array to set top of patches to the minimum
            of the top of a layer or the head value. Used to create
            patches that conform to water-level elevations.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.pyplot.plot

        Returns
        -------
        plot : list containing matplotlib.fillbetween objects

        """
        return self.__cls.plot_fill_between(a=a, colors=colors,
                                            masked_values=masked_values,
                                            head=head, **kwargs)

    def contour_array(self, a, masked_values=None, head=None, **kwargs):
        """
        Contour a three-dimensional array.

        Parameters
        ----------
        a : numpy.ndarray
            Three-dimensional array to plot.
        masked_values : iterable of floats, ints
            Values to mask.
        head : numpy.ndarray
            Three-dimensional array to set top of patches to the minimum
            of the top of a layer or the head value. Used to create
            patches that conform to water-level elevations.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.pyplot.contour

        Returns
        -------
        contour_set : matplotlib.pyplot.contour

        """
        return self.__cls.contour_array(a=a, masked_values=masked_values,
                                        head=head, **kwargs)

    def plot_inactive(self, ibound=None, color_noflow='black', **kwargs):
        """
        Make a plot of inactive cells.  If not specified, then pull ibound
        from the self.ml

        Parameters
        ----------
        ibound : numpy.ndarray
            ibound array to plot.  (Default is ibound in 'BAS6' package.)

        color_noflow : string
            (Default is 'black')

        Returns
        -------
        quadmesh : matplotlib.collections.QuadMesh

        """
        if ibound is None:
            if self.mg.idomain is None:
                raise AssertionError("An idomain array must be provided")
            else:
                ibound = self.mg.idomain

        plotarray = np.zeros(ibound.shape, dtype=np.int)
        idx1 = (ibound == 0)
        plotarray[idx1] = 1
        plotarray = np.ma.masked_equal(plotarray, 0)
        cmap = matplotlib.colors.ListedColormap(['0', color_noflow])
        bounds = [0, 1, 2]
        norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
        patches = self.plot_array(plotarray, cmap=cmap, norm=norm, **kwargs)

        return patches

    def plot_ibound(self, ibound=None, color_noflow='black', color_ch='blue',
                    color_vpt="red", head=None, **kwargs):
        """
        Make a plot of ibound.  If not specified, then pull ibound from the
        self.model

        Parameters
        ----------
        ibound : numpy.ndarray
            ibound array to plot.  (Default is ibound in 'BAS6' package.)
        color_noflow : string
            (Default is 'black')
        color_ch : string
            Color for constant heads (Default is 'blue'.)
        head : numpy.ndarray
            Three-dimensional array to set top of patches to the minimum
            of the top of a layer or the head value. Used to create
            patches that conform to water-level elevations.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.collections.PatchCollection

        Returns
        -------
        patches : matplotlib.collections.PatchCollection

        """
        if ibound is None:
            if self.model is not None:
                if self.model.version == "mf6":
                    color_ch = color_vpt

            if self.mg.idomain is None:
                raise AssertionError("Ibound/Idomain array must be provided")

            ibound = self.mg.idomain

        plotarray = np.zeros(ibound.shape, dtype=np.int)
        idx1 = (ibound == 0)
        idx2 = (ibound < 0)
        plotarray[idx1] = 1
        plotarray[idx2] = 2
        plotarray = np.ma.masked_equal(plotarray, 0)
        cmap = matplotlib.colors.ListedColormap(['none', color_noflow,
                                                 color_ch])
        bounds = [0, 1, 2, 3]
        norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
        # mask active cells
        patches = self.plot_array(plotarray, masked_values=[0], head=head,
                                  cmap=cmap, norm=norm, **kwargs)
        return patches

    def plot_grid(self, **kwargs):
        """
        Plot the grid lines.

        Parameters
        ----------
            kwargs : ax, colors.  The remaining kwargs are passed into the
                the LineCollection constructor.

        Returns
        -------
            lc : matplotlib.collections.LineCollection

        """
        if 'ax' in kwargs:
            ax = kwargs.pop('ax')
        else:
            ax = self.ax

        col = self.get_grid_line_collection(**kwargs)
        if col is not None:
            ax.add_collection(col)
            ax.set_xlim(self.extent[0], self.extent[1])
            ax.set_ylim(self.extent[2], self.extent[3])

        return col

    def plot_bc(self, name=None, package=None, kper=0, color=None,
                head=None, **kwargs):
        """
        Plot boundary conditions locations for a specific boundary
        type from a flopy model

        Parameters
        ----------
        name : string
            Package name string ('WEL', 'GHB', etc.). (Default is None)
        package : flopy.modflow.Modflow package class instance
            flopy package class instance. (Default is None)
        kper : int
            Stress period to plot
        color : string
            matplotlib color string. (Default is None)
        head : numpy.ndarray
            Three-dimensional array (structured grid) or
            Two-dimensional array (vertex grid)
            to set top of patches to the minimum of the top of a\
            layer or the head value. Used to create
            patches that conform to water-level elevations.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.collections.PatchCollection

        Returns
        -------
        patches : matplotlib.collections.PatchCollection

        """
        if 'ftype' in kwargs and name is None:
            name = kwargs.pop('ftype')

        # Find package to plot
        if package is not None:
            p = package
        elif self.model is not None:
            if name is None:
                raise Exception('ftype not specified')
            name = name.upper()
            p = self.model.get_package(name)
        else:
            raise Exception('Cannot find package to plot')

        # trap for mf6 'cellid' vs mf2005 'k', 'i', 'j' convention
        if isinstance(p, list) or p.parent.version == "mf6":
            if not isinstance(p, list):
                p = [p]

            idx = np.array([])
            for pp in p:
                if pp.package_type in ('lak', 'sfr', 'maw', 'uzf'):
                    t = plotutil.advanced_package_bc_helper(pp, self.mg,
                                                            kper)
                else:
                    try:
                        mflist = pp.stress_period_data.array[kper]
                    except Exception as e:
                        raise Exception("Not a list-style boundary package: "
                                        + str(e))
                    if mflist is None:
                        return

                    t = np.array([list(i) for i in mflist['cellid']],
                                 dtype=int).T

                if len(idx) == 0:
                    idx = np.copy(t)
                else:
                    idx = np.append(idx, t, axis=1)

        else:
            # modflow-2005 structured and unstructured grid
            if p.package_type in ('uzf', 'lak'):
                idx = plotutil.advanced_package_bc_helper(p, self.mg, kper)
            else:
                try:
                    mflist = p.stress_period_data[kper]
                except Exception as e:
                    raise Exception("Not a list-style boundary package: "
                                    + str(e))
                if mflist is None:
                    return
                if len(self.mg.shape) == 3:
                    idx = [mflist['k'], mflist['i'], mflist['j']]
                else:
                    idx = mflist['node']

        # Plot the list locations, change this to self.mg.shape
        if len(self.mg.shape) != 3:
            plotarray = np.zeros((self.mg.nlay, self.mg.ncpl), dtype=np.int)
            plotarray[tuple(idx)] = 1
        else:
            plotarray = np.zeros((self.mg.nlay, self.mg.nrow, self.mg.ncol), dtype=np.int)
            plotarray[idx[0], idx[1], idx[2]] = 1

        plotarray = np.ma.masked_equal(plotarray, 0)
        if color is None:
            key = name[:3].upper()
            if key in plotutil.bc_color_dict:
                c = plotutil.bc_color_dict[key]
            else:
                c = plotutil.bc_color_dict['default']
        else:
            c = color
        cmap = matplotlib.colors.ListedColormap(['none', c])
        bounds = [0, 1, 2]
        norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
        patches = self.plot_array(plotarray, masked_values=[0],
                                  head=head, cmap=cmap, norm=norm, **kwargs)

        return patches

    def plot_vector(self, vx, vy, vz, head=None, kstep=1, hstep=1,
                    normalize=False, masked_values=None, **kwargs):
        """
        Plot a vector.

        Parameters
        ----------
        vx : np.ndarray
            x component of the vector to be plotted (non-rotated)
            array shape must be (nlay, nrow, ncol) for a structured grid
            array shape must be (nlay, ncpl) for a unstructured grid
        vy : np.ndarray
            y component of the vector to be plotted (non-rotated)
            array shape must be (nlay, nrow, ncol) for a structured grid
            array shape must be (nlay, ncpl) for a unstructured grid
        vz : np.ndarray
            y component of the vector to be plotted (non-rotated)
            array shape must be (nlay, nrow, ncol) for a structured grid
            array shape must be (nlay, ncpl) for a unstructured grid
        head : numpy.ndarray
            MODFLOW's head array.  If not provided, then the quivers will be
            plotted in the cell center.
        kstep : int
            layer frequency to plot (default is 1)
        hstep : int
            horizontal frequency to plot (default is 1)
        normalize : bool
            boolean flag used to determine if vectors should be normalized
            using the vector magnitude in each cell (default is False)
        masked_values : iterable of floats
            values to mask
        kwargs : matplotlib.pyplot keyword arguments for the
            plt.quiver method

        Returns
        -------
        quiver : matplotlib.pyplot.quiver
            result of the quiver function

        """
        if 'pivot' in kwargs:
            pivot = kwargs.pop('pivot')
        else:
            pivot = 'middle'

        if 'ax' in kwargs:
            ax = kwargs.pop('ax')
        else:
            ax = self.ax

        # this function does not support arbitrary cross-sections, so check it
        arbitrary = False
        if self.mg.grid_type == 'structured':
            if not (self.direction == 'x' or self.direction == 'y'):
                arbitrary = True
        else:
            # check within a tolerance
            pts = self.pts
            xuniform = [True if abs(pts.T[0, 0] - i) < 1
                        else False for i in pts.T[0]]
            yuniform = [True if abs(pts.T[1, 0] - i) < 1
                        else False for i in pts.T[1]]
            if not np.all(xuniform) and not np.all(yuniform):
                arbitrary = True
        if arbitrary:
            err_msg = "plot_specific_discharge() does not " \
                      "support arbitrary cross-sections"
            raise AssertionError(err_msg)

        # get the actual values to plot
        if self.direction == 'x':
            u_tmp = vx
        elif self.direction == 'y':
            u_tmp = -1. * vy
        v_tmp = vz
        if self.mg.grid_type == "structured":
            if isinstance(head, np.ndarray):
                zcentergrid = self.__cls.set_zcentergrid(head)
            else:
                zcentergrid = self.zcentergrid

            if self.geographic_coords:
                xcentergrid = self.__cls.geographic_xcentergrid
            else:
                xcentergrid = self.xcentergrid

            if self.mg.nlay == 1:
                x = []
                z = []
                for k in range(self.mg.nlay):
                    for i in range(xcentergrid.shape[1]):
                        x.append(xcentergrid[k, i])
                        z.append(0.5 * (zcentergrid[k, i] + zcentergrid[k + 1, i]))
                x = np.array(x).reshape((1, xcentergrid.shape[1]))
                z = np.array(z).reshape((1, xcentergrid.shape[1]))
            else:
                x = xcentergrid
                z = zcentergrid

            u = []
            v = []
            xedge, yedge = self.mg.xyedges
            for k in range(self.mg.nlay):
                u.append(plotutil.cell_value_points(self.xpts, xedge,
                                                    yedge, u_tmp[k, :, :]))
                v.append(plotutil.cell_value_points(self.xpts, xedge,
                                                    yedge, v_tmp[k, :, :]))
            u = np.array(u)
            v = np.array(v)
            x = x[::kstep, ::hstep]
            z = z[::kstep, ::hstep]
            u = u[::kstep, ::hstep]
            v = v[::kstep, ::hstep]

            # upts and vpts has a value for the left and right
            # sides of a cell. Sample every other value for quiver
            u = u[:, ::2]
            v = v[:, ::2]

        else:
            # kstep implementation for vertex grid
            projpts = {key: value for key, value in self.__cls.projpts.items()
                       if (key // self.mg.ncpl) % kstep == 0}

            # set x and z centers
            if isinstance(head, np.ndarray):
                # pipe kstep to set_zcentergrid to assure consistent array size
                zcenters = self.__cls.set_zcentergrid(np.ravel(head), kstep=kstep)
            else:
                zcenters = [np.mean(np.array(v).T[1]) for i, v
                            in sorted(projpts.items())]

            u = np.array([u_tmp.ravel()[cell] for cell in sorted(projpts)])

            x = np.array([np.mean(np.array(v).T[0]) for i, v
                          in sorted(projpts.items())])

            z = np.ravel(zcenters)
            v = np.array([v_tmp.ravel()[cell] for cell in sorted(projpts)])

            x = x[::hstep]
            z = z[::hstep]
            u = u[::hstep]
            v = v[::hstep]

        # mask values
        if masked_values is not None:
            for mval in masked_values:
                to_mask = np.logical_or(u==mval, v==mval)
                u[to_mask] = np.nan
                v[to_mask] = np.nan

        # normalize
        if normalize:
            vmag = np.sqrt(u ** 2. + v ** 2.)
            idx = vmag > 0.
            u[idx] /= vmag[idx]
            v[idx] /= vmag[idx]

        # plot with quiver
        quiver = ax.quiver(x, z, u, v, pivot=pivot, **kwargs)

        return quiver

    def plot_specific_discharge(self, spdis, head=None, kstep=1,
                                hstep=1, normalize=False, **kwargs):
        """
        DEPRECATED. Use plot_vector() instead, which should follow after
        postprocessing.get_specific_discharge().

        Use quiver to plot vectors.

        Parameters
        ----------
        spdis : np.recarray
            numpy recarray of specific discharge information. This
            can be grabbed directly from the CBC file if SAVE_SPECIFIC_DISCHARGE
            is used in the MF6 NPF file.
        head : numpy.ndarray
            MODFLOW's head array.  If not provided, then the quivers will be plotted
            in the cell center.
        kstep : int
            layer frequency to plot. (Default is 1.)
        hstep : int
            horizontal frequency to plot. (Default is 1.)
        normalize : bool
            boolean flag used to determine if discharge vectors should
            be normalized using the magnitude of the specific discharge in each
            cell. (default is False)
        kwargs : dictionary
            Keyword arguments passed to plt.quiver()

        Returns
        -------
        quiver : matplotlib.pyplot.quiver
            Vectors

        """
        import warnings
        warnings.warn('plot_specific_discharge() has been deprecated. Use '
                      'plot_vector() instead, which should follow after '
                      'postprocessing.get_specific_discharge()',
                      DeprecationWarning)

        if 'pivot' in kwargs:
            pivot = kwargs.pop('pivot')
        else:
            pivot = 'middle'

        if 'ax' in kwargs:
            ax = kwargs.pop('ax')
        else:
            ax = self.ax

        if isinstance(spdis, list):
            print("Warning: Selecting the final stress period from Specific"
                  " Discharge list")
            spdis = spdis[-1]

        if self.mg.grid_type == "structured":
            ncpl = self.mg.nrow * self.mg.ncol

        else:
            ncpl = self.mg.ncpl

        nlay = self.mg.nlay

        qx = np.zeros((nlay * ncpl))
        qz = np.zeros((nlay * ncpl))
        ib = np.zeros((nlay * ncpl), dtype=bool)

        idx = np.array(spdis['node']) - 1

        # check that vertex grid cross sections are not arbitrary
        # within a tolerance!
        if self.mg.grid_type != 'structured':
            pts = self.pts
            xuniform = [True if abs(pts.T[0, 0] - i) < 1
                        else False for i in pts.T[0]]
            yuniform = [True if abs(pts.T[1, 0] - i) < 1
                        else False for i in pts.T[1]]
            if not np.all(xuniform):
                if not np.all(yuniform):
                    err_msg = "plot_specific_discharge does not " \
                              "support aribtrary cross sections"
                    raise AssertionError(err_msg)

        if self.direction == 'x':
            qx[idx] = spdis['qx']
        elif self.direction == 'y':
            qx[idx] = spdis['qy'] * -1
        else:
            err_msg = 'plot_specific_discharge does not ' \
                      'support arbitrary cross-sections'
            raise AssertionError(err_msg)

        qz[idx] = spdis["qz"]
        ib[idx] = True

        if self.mg.grid_type == "structured":
            qx.shape = (self.mg.nlay, self.mg.nrow, self.mg.ncol)
            qz.shape = (self.mg.nlay, self.mg.nrow, self.mg.ncol)
            ib.shape = (self.mg.nlay, self.mg.nrow, self.mg.ncol)

            if isinstance(head, np.ndarray):
                zcentergrid = self.__cls.set_zcentergrid(head)
            else:
                zcentergrid = self.zcentergrid

            if self.geographic_coords:
                xcentergrid = self.__cls.geographic_xcentergrid
            else:
                xcentergrid = self.xcentergrid

            if nlay == 1:
                x = []
                z = []
                for k in range(nlay):
                    for i in range(xcentergrid.shape[1]):
                        x.append(xcentergrid[k, i])
                        z.append(0.5 * (zcentergrid[k, i] + zcentergrid[k + 1, i]))
                x = np.array(x).reshape((1, xcentergrid.shape[1]))
                z = np.array(z).reshape((1, xcentergrid.shape[1]))
            else:
                x = xcentergrid
                z = zcentergrid

            u = []
            v = []
            ibx = []
            xedge, yedge = self.mg.xyedges
            for k in range(self.mg.nlay):
                u.append(plotutil.cell_value_points(self.xpts, xedge,
                                                    yedge, qx[k, :, :]))
                v.append(plotutil.cell_value_points(self.xpts, xedge,
                                                    yedge, qz[k, :, :]))
                ibx.append(plotutil.cell_value_points(self.xpts, xedge,
                                                      yedge, ib[k, :, :]))
            u = np.array(u)
            v = np.array(v)
            ibx = np.array(ibx)
            x = x[::kstep, ::hstep]
            z = z[::kstep, ::hstep]
            u = u[::kstep, ::hstep]
            v = v[::kstep, ::hstep]
            ib = ibx[::kstep, ::hstep]

            # upts and vpts has a value for the left and right
            # sides of a cell. Sample every other value for quiver
            u = u[:, ::2]
            v = v[:, ::2]
            ib = ib[:, ::2]

        else:
            # kstep implementation for vertex grid
            projpts = {key: value for key, value in self.__cls.projpts.items()
                       if (key // ncpl) % kstep == 0}

            # set x and z centers
            if isinstance(head, np.ndarray):
                # pipe kstep to set_zcentergrid to assure consistent array size
                zcenters = self.__cls.set_zcentergrid(np.ravel(head), kstep=kstep)
            else:
                zcenters = [np.mean(np.array(v).T[1]) for i, v
                            in sorted(projpts.items())]

            u = np.array([qx[cell] for cell in sorted(projpts)])

            x = np.array([np.mean(np.array(v).T[0]) for i, v
                          in sorted(projpts.items())])

            z = np.ravel(zcenters)
            v = np.array([qz[cell] for cell
                          in sorted(projpts)])
            ib = np.array([ib[cell] for cell
                           in sorted(projpts)])

            x = x[::hstep]
            z = z[::hstep]
            u = u[::hstep]
            v = v[::hstep]
            ib = ib[::hstep]

        if normalize:
            vmag = np.sqrt(u ** 2. + v ** 2.)
            idx = vmag > 0.
            u[idx] /= vmag[idx]
            v[idx] /= vmag[idx]

        # mask with an ibound array
        u[~ib] = np.nan
        v[~ib] = np.nan

        quiver = ax.quiver(x, z, u, v, pivot=pivot, **kwargs)

        return quiver

    def plot_discharge(self, frf, fff, flf=None,
                       head=None, kstep=1, hstep=1, normalize=False,
                       **kwargs):
        """
        DEPRECATED. Use plot_vector() instead, which should follow after
        postprocessing.get_specific_discharge().

        Use quiver to plot vectors.

        Parameters
        ----------
        frf : numpy.ndarray
            MODFLOW's 'flow right face'
        fff : numpy.ndarray
            MODFLOW's 'flow front face'
        flf : numpy.ndarray
            MODFLOW's 'flow lower face' (Default is None.)
        head : numpy.ndarray
            MODFLOW's head array.  If not provided, then will assume confined
            conditions in order to calculated saturated thickness.
        kstep : int
            layer frequency to plot. (Default is 1.)
        hstep : int
            horizontal frequency to plot. (Default is 1.)
        normalize : bool
            boolean flag used to determine if discharge vectors should
            be normalized using the magnitude of the specific discharge in each
            cell. (default is False)
        kwargs : dictionary
            Keyword arguments passed to plt.quiver()

        Returns
        -------
        quiver : matplotlib.pyplot.quiver
            Vectors

        """
        import warnings
        warnings.warn('plot_discharge() has been deprecated. Use '
                      'plot_vector() instead, which should follow after '
                      'postprocessing.get_specific_discharge()',
                      DeprecationWarning)

        if self.mg.grid_type != "structured":
            err_msg = "Use plot_specific_discharge for " \
                      "{} grids".format(self.mg.grid_type)
            raise NotImplementedError(err_msg)

        else:
            ib = np.ones((self.mg.nlay, self.mg.nrow, self.mg.ncol))
            if self.mg.idomain is not None:
                ib = self.mg.idomain

            delr = self.mg.delr
            delc = self.mg.delc
            top = self.mg.top
            botm = self.mg.botm
            if not np.all(self.active==1):
                botm = botm[self.active==1]
            nlay = botm.shape[0]
            laytyp = None
            hnoflo = 999.
            hdry = 999.

            if self.model is not None:
                if self.model.laytyp is not None:
                    laytyp = self.model.laytyp

                if self.model.hnoflo is not None:
                    hnoflo = self.model.hnoflo

                if self.model.hdry is not None:
                    hdry = self.model.hdry

            # If no access to head or laytyp, then calculate confined saturated
            # thickness by setting laytyp to zeros
            if head is None or laytyp is None:
                head = np.zeros(botm.shape, np.float32)
                laytyp = np.zeros((nlay), dtype=np.int)
                head[0, :, :] = top
                if nlay > 1:
                    head[1:, :, :] = botm[:-1, :, :]

            sat_thk = plotutil.PlotUtilities. \
                saturated_thickness(head, top, botm,
                                    laytyp, [hnoflo, hdry])

            # Calculate specific discharge
            qx, qy, qz = plotutil.PlotUtilities. \
                centered_specific_discharge(frf, fff, flf,
                                            delr, delc, sat_thk)

            if qz is None:
                qz = np.zeros((qx.shape), dtype=np.float)

            ib = ib.ravel()
            qx = qx.ravel()
            qy = qy.ravel() * -1
            qz = qz.ravel()

            temp = []
            for ix, val in enumerate(ib):
                if val != 0:
                    temp.append((ix + 1, qx[ix], -qy[ix], qz[ix]))

            spdis = np.recarray((len(temp),), dtype=[('node', np.int),
                                                     ("qx", np.float),
                                                     ("qy", np.float),
                                                     ("qz", np.float)])
            for ix, tup in enumerate(temp):
                spdis[ix] = tup

            self.plot_specific_discharge(spdis, head=head, kstep=kstep,
                                         hstep=hstep, normalize=normalize,
                                         **kwargs)

    def get_grid_patch_collection(self, zpts, plotarray, **kwargs):
        """
        Get a PatchCollection of plotarray in unmasked cells

        Parameters
        ----------
        zpts : numpy.ndarray
            array of z elevations that correspond to the x, y, and horizontal
            distance along the cross-section (self.xpts). Constructed using
            plotutil.cell_value_points().
        plotarray : numpy.ndarray
            Three-dimensional array to attach to the Patch Collection.
        **kwargs : dictionary
            keyword arguments passed to matplotlib.collections.PatchCollection

        Returns
        -------
        patches : matplotlib.collections.PatchCollection

        """
        if self.mg.grid_type == "structured":
            return self.__cls.get_grid_patch_collection(zpts=zpts, plotarray=plotarray,
                                                        **kwargs)
        elif self.mg.grid_type == "unstructured":
            raise NotImplementedError()

        else:
            return self.__cls.get_grid_patch_collection(projpts=zpts, plotarray=plotarray,
                                                        **kwargs)

    def get_grid_line_collection(self, **kwargs):
        """
        Get a LineCollection of the grid

        Parameters
        ----------
        **kwargs : dictionary
            keyword arguments passed to matplotlib.collections.LineCollection

        Returns
        -------
        linecollection : matplotlib.collections.LineCollection
        """
        return self.__cls.get_grid_line_collection(**kwargs)


class DeprecatedCrossSection(PlotCrossSection):
    """
    Deprecation handler for the PlotCrossSection class

    Parameters
    ----------
    ax : matplotlib.pyplot.axes object
    model : flopy.modflow.Modflow object
    modelgrid : flopy.discretization.Grid object
    line : dict
        Dictionary with either "row", "column", or "line" key. If key
        is "row" or "column" key value should be the zero-based row or
        column index for cross-section. If key is "line" value should
        be an array of (x, y) tuples with vertices of cross-section.
        Vertices should be in map coordinates consistent with xul,
        yul, and rotation.
    extent : tuple of floats
        (xmin, xmax, ymin, ymax) will be used to specify axes limits.  If None
        then these will be calculated based on grid, coordinates, and rotation.

    """
    def __init__(self, ax=None, model=None, modelgrid=None,
                 line=None, extent=None):
        super(DeprecatedCrossSection, self).__init__(ax=ax, model=model,
                                                     modelgrid=modelgrid,
                                                     line=line,
                                                     extent=extent)
