<filename>deltametrics/section.py
import abc
import warnings

import numpy as np
from scipy import sparse

import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection

from . import cube
from . import plot
from . import utils


class BaseSectionVariable(np.ndarray):
    """Section variable.

    We subclass the numpy `ndarray`, in order to have `ndarray` methods to
    the object in subclasses (e.g., `add`, slicing).

    This is a really lightweight subclass of the `np.ndarray`, but it allows
    us to add coordinate arrays to the variable, which are always needed to
    display info correctly. We also gain the ability to return the correct
    subclass type following views and slicing of the subclasses.

    .. note::

        Subclasses should implement the ``__init__`` method.

    """
    _spacetime_names = ['full', 'spacetime', 'as spacetime', 'as_spacetime']
    _preserved_names = ['psvd', 'preserved', 'as preserved', 'as_preserved']
    _stratigraphy_names = ['strat', 'strata', 'stratigraphy',
                           'as stratigraphy', 'as_stratigraphy']

    def __new__(cls, _data, _s, _z, _psvd_mask=None, **unused_kwargs):
        # Input array is an already formed ndarray instance
        obj = np.asarray(_data).view(cls)
        if (_psvd_mask is not None):
            _psvd_mask = np.asarray(_psvd_mask)
            if _psvd_mask.shape != obj.shape:
                raise ValueError(
                    'Shape of "_psvd_mask" incompatible with "_data" array.')
        obj._psvd_mask = _psvd_mask
        if (len(_z) != obj.shape[0]) or (len(_s) != obj.shape[1]):
            raise ValueError(
                'Shape of "_s" or "_z" incompatible with "_data" array.')
        obj._s = _s
        obj._z = _z
        obj._S, obj._Z = np.meshgrid(obj._s, obj._z)
        return obj

    def __array_finalize__(self, obj):
        if obj is None:
            return
        self._psvd_mask = getattr(obj, '_psvd_mask', None)
        self._s = getattr(obj, '_s', None)
        self._z = getattr(obj, '_z', None)
        self._S, self._Z = np.meshgrid(self._s, self._z)


class DataSectionVariable(BaseSectionVariable):
    """Variable returned from a DataCube Section.

    Subclasses numpy MaskedArray, so supports arbitrary math.

    """
    _default_data = 'spacetime'

    def __init__(self, _data, _s, _z, _psvd_mask=None, _strat_attr=None):
        """Construct the array from section info.

        Parameters
        ----------
        _data : :obj:`ndarray`
            Slice of underlying data. Generated by slicing the CubeVariable
            with: :code:`cube[var][:, self._y, self._x]`

        _psvd_mask : :obj:`ndarray`
            Mask indicating the *preserved* voxels. Must have same shape as
            `_data`.

        strat_attr : :obj:`dict`
            Dictionary of attributes regarding stratigraphy generated by the
            section, on instantiation. May be a nearly empty dictionary, but
            must always be provided. Only used if
            `Section._knows_stratigraphy`.

        .. note::
            ``__new__`` from the base class is called *before* ``__init__``.
            The ``__new__`` method configures the `_data`, `_s`, `_z`, and
            `_psvd_mask` arguments.
        """
        if not (_strat_attr is None):
            self.strat_attr = _strat_attr
            self._knows_stratigraphy = True
        else:
            self._knows_stratigraphy = False

    @property
    def knows_stratigraphy(self):
        """Whether the data variable knows preservation information."""
        return self._knows_stratigraphy

    def _check_knows_stratigraphy(self):
        """Check whether "knows_stratigraphy".

        Raises
        ------
        AttributeError
            Raises if does not know stratigraphy.
        """
        if not self._knows_stratigraphy:
            raise utils.NoStratigraphyError(obj=self)
        return self._knows_stratigraphy

    def as_preserved(self):
        """Variable with only preserved values.

        Returns
        -------
        ma : :obj:`np.ma.MaskedArray`
            A numpy MaskedArray with non-preserved values masked.
        """
        if self._check_knows_stratigraphy():
            return np.ma.MaskedArray(self, ~self._psvd_mask)

    def as_stratigraphy(self):
        """Variable as preserved stratigraphy.

        .. warning::

            This method returns a sparse array that is not suitable to be
            displayed directly. Use
            :obj:`get_display_arrays(style='stratigraphy')` instead to get
            corresponding x-y coordinates for plotting the array.
        """
        if self._check_knows_stratigraphy():
            # actual data, where preserved
            _psvd_data = self[self.strat_attr['psvd_idx']]
            _sp = sparse.coo_matrix((_psvd_data,
                                     (self.strat_attr['z_sp'],
                                      self.strat_attr['s_sp'])))
            return _sp


class StratigraphySectionVariable(BaseSectionVariable):
    """
    """
    _default_data = 'stratigraphy'

    def __init__(self, _data, _s, _z):
        self._knows_spacetime = False

    @property
    def knows_spacetime(self):
        """Whether the data variable knows preservation information."""
        return self._knows_spacetime

    def _check_knows_spacetime(self):
        """Check whether "knows_spacetime".

        Raises
        ------
        AttributeError
            Raises always when this method is called, because a
            StratigraphySectionVariable will never know spacetime information
            directly.
        """
        raise AttributeError(
            'No "spacetime" or "preserved" information available.')


class BaseSection(abc.ABC):
    """Base section object.

    Defines common attributes and methods of a section object.

    This object should wrap around many of the functions available from
    :obj:`~deltametrics.strat`.

    """

    def __init__(self, section_type, *args, name=None):
        """
        Identify coordinates defining the section.

        Parameters
        ----------
        CubeInstance : :obj:`~deltametrics.cube.Cube` subclass instance, optional
            Connect to this cube. No connection is made if cube is not
            provided.

        Notes
        -----

        If no arguments are passed, an empty section not connected to any cube
        is returned. This cube will will need to be manually connected to have
        any functionality (via the :meth:`connect` method).
        """
        # begin unconnected
        self._s = None
        self._z = None
        self._x = None
        self._y = None
        self._trace = None
        self._shape = None
        self._variables = None
        self.cube = None

        self.section_type = section_type
        self._name = name

        if len(args) > 1:
            raise ValueError('Expected single positional argument to \
                             %s instantiation.'
                             % type(self))

        if len(args) > 0:
            self.connect(args[0])
        else:
            pass

    def connect(self, CubeInstance, name=None):
        """Connect this Section instance to a Cube instance.
        """
        if not issubclass(type(CubeInstance), cube.BaseCube):
            raise TypeError('Expected type is subclass of {_exptype}, '
                            'but received was {_gottype}.'.format(
                                _exptype=type(cube.BaseCube),
                                _gottype=type(CubeInstance)))
        self.cube = CubeInstance
        self._variables = self.cube.variables
        self.name = name  # use the setter to determine the _name
        self._compute_section_coords()
        self._compute_section_attrs()

    @property
    def name(self):
        return self._name

    @name.setter
    def name(self, var):
        if (self._name is None):
            # _name is not yet set
            self._name = var or self.section_type
        else:
            # _name is already set
            if not (var is None):
                warnings.warn(
                    UserWarning("`name` argument supplied to instantiated "
                                "`Section` object. To change the name of "
                                "a Section, you must set the attribute "
                                "directly with `section._name = 'name'`."))
            # do nothing

    @abc.abstractmethod
    def _compute_section_coords(self):
        """Should calculate x-y coordinates of the section.

        Sets the value ``self._x`` and ``self._y`` according to the algorithm
        of each section initialization.

        .. warning::

            When implementing new section types, be sure that ``self._x`` and
            ``self._y`` are *one-dimensional arrays*, or you will get an
            improperly shaped Section array in return.
        """
        ...

    def _compute_section_attrs(self):
        """Compute attrs

        Compute the along-section coordinate array from x-y pts pairs
        definining the section.
        """
        self._s = np.cumsum(np.hstack(
            (0, np.sqrt((self._x[1:] - self._x[:-1])**2
             + (self._y[1:] - self._y[:-1])**2))))
        self._z = self.cube.z
        self._shape = (len(self._z), len(self._s))
        self._trace = np.column_stack((self._x, self._y))

    @property
    def trace(self):
        """Coordinates of the section in the x-y plane.
        """
        return self._trace

    @property
    def s(self):
        """Along-section coordinate."""
        return self._s

    @property
    def z(self):
        """Up-section (vertical) coordinate."""
        return self._z

    @property
    def shape(self):
        """Section shape.

        Simply a `tuple` equivalent to ``(len(z), len(s))``
        """
        return self._shape

    @property
    def variables(self):
        """List of variables.
        """
        return self._variables

    @property
    def strat_attr(self):
        """Stratigraphic attributes data object.

        Raises
        ------
        NoStratigraphyError
            If no stratigraphy information is found for the section.
        """
        if self.cube._knows_stratigraphy:
            return self.cube.strat_attr
        else:
            raise utils.NoStratigraphyError(obj=self, var='strat_attr')

    def __getitem__(self, var):
        """Get a slice of the section.

        Slicing the section instance creates a
        :obj:`~deltametrics.section.SectionVariable` instance from data for
        variable ``var``.

        .. note:: We only support slicing by string.

        Parameters
        ----------
        var : :obj:`str`
            Which variable to slice.

        Returns
        -------
        SectionVariable : :obj:`~deltametrics.section.SectionVariable` instance
            SectionVariable instance for variable ``var``.
        """
        if type(self.cube) is cube.DataCube:
            if self.cube._knows_stratigraphy:
                return DataSectionVariable(
                    _data=self.cube[var].data.values[:, self._y, self._x],
                    _s=self.s, _z=self.z,
                    _psvd_mask=self.cube.strat_attr.psvd_idx[:, self._y, self._x],  # noqa: E501
                    _strat_attr=self.cube.strat_attr(
                        'section', self._y, self._x)
                    )
            else:
                return DataSectionVariable(
                    _data=self.cube[var].data.values[:, self._y, self._x],
                    _s=self.s, _z=self.z
                    )
        elif type(self.cube) is cube.StratigraphyCube:
            return StratigraphySectionVariable(
                _data=self.cube[var].data.values[:, self._y, self._x],
                _s=self.s, _z=self.z
                )
        elif self.cube is None:
            raise AttributeError(
                'No cube connected. Are you sure you ran `.connect()`?')
        else:
            raise TypeError('Unknown Cube type encountered: %s'
                            % type(self.cube))

    def show(self, SectionAttribute, style='shaded', data=None,
             label=False, colorbar=True, colorbar_label=False, ax=None):
        """Show the section.

        Method enumerates convenient routines for visualizing sections of data
        and stratigraphy. Includes support for multiple data `style` and
        mutuple `data` choices as well.

        .. note::

            The colors for `style='lines'` are determined from the left-end
            edge node, and colors for the `style='shaded'` mesh are determined
            from the lower-left-end edge node of the quad.

        Parameters
        ----------

        SectionAttribute : :obj:`str`
            Which attribute to show.

        style : :obj:`str`, optional
            What style to display the section with. Choices are 'mesh' or
            'line'.

        data : :obj:`str`, optional
            Argument passed to
            :obj:`~deltametrics.section.DataSectionVariable.get_display_arrays`
            or
            :obj:`~deltametrics.section.DataSectionVariable.get_display_lines`.
            Supported options are `'spacetime'`, `'preserved'`, and
            `'stratigraphy'`. Default is to display full spacetime plot for
            section generated from a `DataCube`, and stratigraphy for a
            `StratigraphyCube` section.

        label : :obj:`bool`, `str`, optional
            Display a label of the variable name on the plot. Default is
            False, display nothing. If ``label=True``, the label name from the
            :obj:`~deltametrics.plot.VariableSet` is used. Other arguments are
            attempted to coerce to `str`, and the literal is diplayed.

        colorbar : :obj:`bool`, optional
            Whether a colorbar is appended to the axis.

        colorbar_label : :obj:`bool`, `str`, optional
            Display a label of the variable name along the colorbar. Default is
            False, display nothing. If ``label=True``, the label name from the
            :obj:`~deltametrics.plot.VariableSet` is used. Other arguments are
            attempted to coerce to `str`, and the literal is diplayed.

        ax : :obj:`~matplotlib.pyplot.Axes` object, optional
            A `matplotlib` `Axes` object to plot the section. Optional; if not
            provided, a call is made to ``plt.gca()`` to get the current (or
            create a new) `Axes` object.

        Examples
        --------
        *Example 1:* Display the `velocity` spacetime section of a DataCube.

        .. doctest::

            >>> rcm8cube = dm.sample_data.rcm8()
            >>> rcm8cube.register_section(
            ...     'demo', dm.section.StrikeSection(y=5))
            >>> rcm8cube.sections['demo'].show('velocity')

        .. plot:: section/section_demo_spacetime.py

        Note that the last line above is functionally equivalent to
        ``rcm8cube.show_section('demo', 'velocity')``.

        *Example 2:* Display a section, with "quick" stratigraphy, as the
        `depth` attribute, displaying several different section styles.

        .. doctest::

            >>> rcm8cube = dm.sample_data.rcm8()
            >>> rcm8cube.stratigraphy_from('eta')
            >>> rcm8cube.register_section(
            ...     'demo', dm.section.StrikeSection(y=5))

            >>> fig, ax = plt.subplots(4, 1, sharex=True, figsize=(6, 9))
            >>> rcm8cube.sections['demo'].show('depth', data='spacetime',
            ...                                 ax=ax[0], label='spacetime')
            >>> rcm8cube.sections['demo'].show('depth', data='preserved',
            ...                                ax=ax[1], label='preserved')
            >>> rcm8cube.sections['demo'].show('depth', data='stratigraphy',
            ...                                ax=ax[2], label='quick stratigraphy')
            >>> rcm8cube.sections['demo'].show('depth', style='lines', data='stratigraphy',
            ...                                ax=ax[3], label='quick stratigraphy')          # noqa: E501

        .. plot:: section/section_demo_quick_strat.py
        """
        # process arguments and inputs
        if not ax:
            ax = plt.gca()
        _varinfo = self.cube.varset[SectionAttribute] if \
            issubclass(type(self.cube), cube.BaseCube) else \
            plot.VariableSet()[SectionAttribute]
        SectionVariableInstance = self[SectionAttribute]

        # main routines for plot styles
        if style in ['shade', 'shaded']:
            _data, _X, _Y = plot.get_display_arrays(SectionVariableInstance,
                                                    data=data)
            ci = ax.pcolormesh(_X, _Y, _data, cmap=_varinfo.cmap,
                               norm=_varinfo.norm,
                               vmin=_varinfo.vmin, vmax=_varinfo.vmax,
                               rasterized=True, shading='auto')
        elif style in ['line', 'lines']:
            _data, _segments = plot.get_display_lines(SectionVariableInstance,
                                                      data=data)
            lc = LineCollection(_segments, cmap=_varinfo.cmap)
            lc.set_array(_data.flatten())
            lc.set_linewidth(1.25)
            ci = ax.add_collection(lc)
        else:
            raise ValueError('Bad style argument: "%s"' % style)

        # style adjustments
        if colorbar:
            cb = plot.append_colorbar(ci, ax)
            if colorbar_label:
                _colorbar_label = _varinfo.label if (colorbar_label is True) \
                    else str(colorbar_label)  # use custom if passed
                cb.ax.set_ylabel(_colorbar_label, rotation=-90, va="bottom")
        ax.margins(y=0.2)
        if label:
            _label = _varinfo.label if (label is True) else str(
                label)  # use custom if passed
            ax.text(0.99, 0.8, _label, fontsize=10,
                    horizontalalignment='right', verticalalignment='center',
                    transform=ax.transAxes)
        xmin, xmax, ymin, ymax = plot.get_display_limits(
            SectionVariableInstance, data=data)
        ax.set_xlim(xmin, xmax)
        ax.set_ylim(ymin, ymax)

    def show_trace(self, *args, ax=None, **kwargs):
        """Plot section trace (x-y plane path).

        Plot the section trace (:obj:`trace`) onto an x-y planview.

        Parameters
        ----------
        *args
            Passed to `matplotlib` :obj:`~matplotlib.pyplot.plot()`.

        ax : :obj:`~matplotlib.pyplot.Axes` object, optional
            A `matplotlib` `Axes` object to plot the trace. Optional; if not
            provided, a call is made to ``plt.gca()`` to get the current (or
            create a new) `Axes` object.

        **kwargs
            Passed to `matplotlib` :obj:`~matplotlib.pyplot.plot()`.
        """
        if not ax:
            ax = plt.gca()

        _label = kwargs.pop('label', self.name)

        ax.plot(self._x, self._y, label=_label, *args, **kwargs)


class PathSection(BaseSection):
    """Path section object.

    Create a Section along user-specified path. Specify the section location
    as an `(N, 2)` `ndarray` of x-y pairs of coordinates that define the
    verticies of the path. All coordinates along the path will be included in
    the section.

    .. important::

        The vertex coordinates must be specified as cell indices (not
        actual x and y coordinate values). This is a needed patch.

    Parameters
    ----------
    *args : :obj:`DataCube` or `StratigraphyCube`
        The `Cube` object to link for underlying data. This option should be
        ommitted if using the :obj:`register_section` method of a `Cube`.

    path : :obj:`ndarray`
        An `(N, 2)` `ndarray` specifying the x-y pairs of coordinates that
            define the verticies of the path to extract the section from.

    **kwargs
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
        options are `name`.

    Returns
    -------
    section : :obj:`PathSection`
        `PathSection` object with specified parameters. The section is
        automatically connected to the underlying `Cube` data source if the
        :obj:`register_section` method of a `Cube` is used to set up the
        section, or the `Cube` is passed as the first positional argument
        during instantiation.

    Examples
    --------

    To create a `PathSection` that is registered to a `DataCube` at
    specified coordinates:

    .. plot::
        :include-source:

        >>> rcm8cube = dm.sample_data.rcm8()
        >>> rcm8cube.register_section('path', dm.section.PathSection(
        ...     path=np.array([[50, 3], [65, 17], [130, 10]])))
        >>>
        >>> # show the location and the "velocity" variable
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
        >>> rcm8cube.sections['path'].show_trace('r--', ax=ax[0])
        >>> rcm8cube.sections['path'].show('velocity', ax=ax[1])
        >>> plt.show()
    """

    def __init__(self, *args, path, **kwargs):
        """Instantiate.

        Parameters
        ----------
        path : :obj:`ndarray`
            An `(N, 2)` `ndarray` specifying the x-y pairs of coordinates that
            define the verticies of the path to extract the section from.

        .. note::

            :obj:`path` must be supplied as a keyword argument.

        """
        self._input_path = path
        super().__init__('path', *args, **kwargs)

    def _compute_section_coords(self):
        """Calculate coordinates of the strike section.
        """
        # convert the points into segments into lists of cells
        _segs = utils.coordinates_to_segments(self._input_path)
        _cell = utils.segments_to_cells(_segs)

        # determine only unique coordinates along the path
        self._path = np.unique(_cell, axis=0)
        self._vertices = np.unique(self._input_path, axis=0)

        self._x = self._path[:, 0]
        self._y = self._path[:, 1]

    @property
    def path(self):
        """Path of the PathSection.

        Returns same as `trace` property.
        """
        return self.trace


class StrikeSection(BaseSection):
    """Strike section object.

    Section oriented along the delta strike (i.e., perpendicular to an inlet
    channel). Specify the location of the strike section with :obj`y` and
    :obj:`x` keyword parameter options.

    .. important::

        The `y` and `x` parameters must be specified as cell indices (not
        actual x and y coordinate values). This is a needed patch.

    Parameters
    ----------
    *args : :obj:`DataCube` or `StratigraphyCube`
        The `Cube` object to link for underlying data. This option should be
        ommitted if using the :obj:`register_section` method of a `Cube`.

    y : :obj:`int`, optional
        The `y` location of the section. This is the distance to locate the
        section from the domain edge with a channel inlet. Defaults to ``0``
        if no value is given.

    x : :obj:`int`, optional
        The `x` limits for the section. Defaults to the full domain width.
        Specify as a two-element `tuple` or `list` of `int`, giving the lower
        and upper bounds of `x` values to span the section.

    **kwargs
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
        options are `name`.

    Returns
    -------
    section : :obj:`StrikeSection`
        `StrikeSection` object with specified parameters. The section is
        automatically connected to the underlying `Cube` data source if the
        :obj:`register_section` method of a `Cube` is used to set up the
        section, or the `Cube` is passed as the first positional argument
        during instantiation.

    Examples
    --------

    To create a `StrikeSection` that is registered to a `DataCube` at
    specified `y` coordinate ``=10``, and spans the entire model domain:

    .. plot::
        :include-source:

        >>> rcm8cube = dm.sample_data.rcm8()
        >>> rcm8cube.register_section('strike', dm.section.StrikeSection(y=10))
        >>>
        >>> # show the location and the "velocity" variable
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
        >>> rcm8cube.sections['strike'].show_trace('r--', ax=ax[0])
        >>> rcm8cube.sections['strike'].show('velocity', ax=ax[1])
        >>> plt.show()

    Similarly, create a `StrikeSection` that is registered to a
    `StratigraphyCube` at specified `y` coordinate ``=20``, and spans only the
    left side of the model domain:

    .. plot::
        :include-source:

        >>> rcm8cube = dm.sample_data.rcm8()
        >>> sc8cube = dm.cube.StratigraphyCube.from_DataCube(rcm8cube)
        >>> sc8cube.register_section(
        ...     'strike_half', dm.section.StrikeSection(y=20, x=[0, 120]))

        >>> # show the location and the "velocity" variable
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
        >>> sc8cube.sections['strike_half'].show_trace('r--', ax=ax[0])
        >>> sc8cube.sections['strike_half'].show('velocity', ax=ax[1])
        >>> plt.show()
    """

    def __init__(self, *args, y=None, x=None, **kwargs):

        self.y = y  # strike coord scalar
        self._input_xlim = x  # the input x lims
        super().__init__('strike', *args, **kwargs)

    def _compute_section_coords(self):
        """Calculate coordinates of the strike section.
        """
        if self._input_xlim is None:
            _nx = self.cube['eta'].shape[2]
            self._x = np.arange(_nx)
        else:
            self._x = np.arange(self._input_xlim[0], self._input_xlim[1])
            _nx = len(self._x)
        self._y = np.tile(self.y, (_nx))


class DipSection(BaseSection):
    """Dip section object.

    """

    def __init__(self, x=-1):
        raise NotImplementedError
        # choose center point if x=-1


class CircularSection(BaseSection):
    """Circular section object.

    Section drawn as a circular cut, located a along the arc a specified
    radius from specified origin.  Specify the location of the circular section
    with :obj`radius` and :obj:`origin` keyword parameter options. The
    circular section trace is interpolated to the nearest integer model domain
    cells, following the mid-point circle algorithm
    (:obj:`~deltametrics.utils.circle_to_cells`).

    .. important::

        The `radius` and `origin` parameters must be specified as cell indices
        (not actual x and y coordinate values). This is a needed patch.

    .. important::

        The `origin` attempts to detect the land width from bed elevation
        changes, but should use the value of ``L0`` recorded in the netcdf
        file, or defined in the cube.

    Parameters
    ----------
    *args : :obj:`DataCube` or `StratigraphyCube`
        The `Cube` object to link for underlying data. This option should be
        ommitted if using the :obj:`register_section` method of a `Cube`.

    radius : :obj:`float`, `int`, optional
        The `radius` of the section. This is the distance to locate the
        section from the :obj:`origin`. If no value is given, the `radius`
        defaults to half of the minimum model domain edge length if it can be
        determined, otherwise defaults to ``1``.

    origin : :obj:`tuple` or `list` of `int`, optional
        The `origin` of the circular section. This is the center of the
        circle. If no value is given, the origin defaults to the center of the
        x-direction of the model domain, and offsets into the domain a
        distance of ``y == L0``, if these values can be determined. I.e., the
        origin defaults to be centered over the channel inlet. If no value is
        given, and these values cannot be determined, the origin defaults to
        ``(0, 0)``.

    **kwargs
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
        options are `name`.

    Returns
    -------
    section : :obj:`CircularSection`
        `CircularSection` object with specified parameters. The section is
        automatically connected to the underlying `Cube` data source if the
        :obj:`register_section` method of a `Cube` is used to set up the
        section, or the `Cube` is passed as the first positional argument
        during instantiation.

    Examples
    --------

    To create a `CircularSection` that is registered to a `DataCube` with
    radius ``=30``, and using the default `origin` options:

    .. plot::
        :include-source:

        >>> rcm8cube = dm.sample_data.rcm8()
        >>> rcm8cube.register_section(
        ...     'circular', dm.section.CircularSection(radius=30))

        >>> # show the location and the "velocity" variable
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
        >>> rcm8cube.sections['circular'].show_trace('r--', ax=ax[0])
        >>> rcm8cube.sections['circular'].show('velocity', ax=ax[1])
        >>> plt.show()
    """

    def __init__(self, *args, radius=None, origin=None, **kwargs):

        self._input_radius = radius
        self._input_origin = origin
        super().__init__('circular', *args, **kwargs)

    def _compute_section_coords(self):
        if (self._input_radius is None):
            self.radius = int(np.min(self.cube.shape[1:]) / 2)
        else:
            self.radius = self._input_radius

        if (self._input_origin is None):
            if (self.cube.meta is None):
                # try and guess the value (should issue a warning?)
                land_width = np.minimum(utils.guess_land_width_from_land(
                    self.cube['eta'][-1, :, 0]), 5)
            else:
                # extract L0 from the cube
                land_width = self.cube.meta['L0']
            self.origin = (int(self.cube.shape[2] / 2),
                           land_width)
        else:
            self.origin = self._input_origin

        xy = utils.circle_to_cells(self.origin, self.radius)
        self._x = xy[0]
        self._y = xy[1]


class RadialSection(BaseSection):
    """Radial section object.

    Section drawn as a radial cut, located a along the line starting from
    `origin` and proceeding away in direction specified by azimuth. Specify
    the location of the radial section with :obj`azimuth` and :obj:`origin`
    keyword parameter options. The radial section trace is interpolated to the
    nearest integer model domain cells, following the a line-walking algorithm
    (:obj:`~deltametrics.utils.line_to_cells`).

    .. important::

        The `origin` parameter must be specified as cell indices (not actual x
        and y coordinate values). This is a needed patch.

    .. important::

        The `origin` attempts to detect the land width from bed elevation
        changes, but should use the value of ``L0`` recorded in the netcdf
        file, or defined in the cube.

    .. important::

        This Section type will only work for deltas with an inlet along the
        ``y=0`` line. For other delta configurations, specify a radial
        section by defining two end points and instantiating a `Section` with
        the :obj:`PathSection`.

    Parameters
    ----------
    *args : :obj:`DataCube` or `StratigraphyCube`
        The `Cube` object to link for underlying data. This option should be
        ommitted if using the :obj:`register_section` method of a `Cube`.

    azimuth : :obj:`float`, `int`, optional
        The `azimuth` of the section, directed away from the origin. If no
        value is given, the `azimuth` defaults to ``90``.

    origin : :obj:`tuple` or `list` of `int`, optional
        The `origin` of the radial section. This is the "start" of the radial
        line. If no value is given, the origin defaults to the center of the
        x-direction of the model domain, and offsets into the domain a
        distance of ``y == L0``, if these values can be determined. I.e., the
        origin defaults to be centered over the channel inlet. If no value is
        given and these values cannot be determined, the origin defaults to
        ``(0, 0)``.

    length : :obj:`float`, `int`, optional
        The length of the section (note this must be given in pixel length).
        If no value is given, the length defaults to the length required to
        reach a model boundary (if a connection to underlying `Cube` exists).
        Otherwise, length is set to ``1``.

    **kwargs
        Keyword arguments are passed to `BaseSection.__init__()`. Supported
        options are `name`.

    Returns
    -------
    section : :obj:`RadialSection`
        `RadialSection` object with specified parameters. The section is
        automatically connected to the underlying `Cube` data source if the
        :obj:`register_section` method of a `Cube` is used to set up the
        section, or the `Cube` is passed as the first positional argument
        during instantiation.

    Examples
    --------

    To create a `RadialSection` that is registered to a `DataCube` at
    specified `origin` coordinate, and spans the entire model domain:

    .. plot::
        :include-source:

        >>> rcm8cube = dm.sample_data.rcm8()
        >>> rcm8cube.register_section(
        ...     'radial', dm.section.RadialSection(azimuth=45))

        >>> # show the location and the "velocity" variable
        >>> fig, ax = plt.subplots(2, 1, figsize=(8, 4))
        >>> rcm8cube.show_plan('eta', t=-1, ax=ax[0], ticks=True)
        >>> rcm8cube.sections['radial'].show_trace('r--', ax=ax[0])
        >>> rcm8cube.sections['radial'].show('velocity', ax=ax[1])
        >>> plt.show()
    """
    def __init__(self, *args, azimuth=None, origin=None, length=None,
                 **kwargs):
        self._input_azimuth = azimuth
        self._input_origin = origin
        self._input_length = length
        super().__init__('radial', *args, **kwargs)

    def _compute_section_coords(self):

        # determine the azimuth
        if (self._input_azimuth is None):
            self.azimuth = 90
        else:
            self.azimuth = self._input_azimuth

        # determine the origin of the line
        if (self._input_origin is None):
            if (self.cube.meta is None):
                # try and guess the value (should issue a warning?)
                land_width = np.minimum(utils.guess_land_width_from_land(
                    self.cube['eta'][-1, :, 0]), 5)
            else:
                # extract L0 from the cube
                land_width = self.cube.meta['L0']
            self.origin = (int(self.cube.shape[2] / 2),
                           land_width)
        else:
            self.origin = self._input_origin

        # determine the length of the line to travel
        # find the line of the azimuth
        theta = self.azimuth
        m = np.tan(theta * np.pi / 180)
        b = self.origin[1] - m * self.origin[0]
        if (self._input_length is None):
            # if no input
            # find the intersection with an edge
            if self.azimuth <= 90.0 and self.azimuth >= 0:
                dx = (self.cube.W - self.origin[0])
                dy = (np.tan(theta * np.pi / 180) * dx)
                if dy <= self.cube.L:
                    end_y = int(np.minimum(
                        m * (self.cube.W) + b, self.cube.L - 1))
                    end_point = (self.cube.W - 1, end_y)
                else:
                    end_x = int(np.minimum(
                        (self.cube.L - b) / m, self.cube.W - 1))
                    end_point = (end_x, self.cube.L - 1)
            elif self.azimuth > 90 and self.azimuth <= 180:
                dx = (self.origin[0])
                dy = (np.tan(theta * np.pi / 180) * dx)
                if np.abs(dy) <= self.cube.L:
                    end_y = b
                    end_point = (0, end_y)
                else:
                    end_x = int(np.maximum((self.cube.L - b) / m,
                                           0))
                    end_point = (end_x, self.cube.L - 1)
            else:
                raise ValueError('Azimuth must be in range (0, 180).')
        else:
            # if input length
            _len = self._input_length
            # use vector math to determine end point len along azimuth
            #   vector is from (0, b) to (origin)
            vec = np.array([self.origin[0] - 0, self.origin[1] - b])
            vec_norm = vec / np.sqrt(vec**2)
            end_point = (self.origin[0] + _len*vec_norm[0],
                         self.origin[1] + _len*vec_norm[1])

        xy = utils.line_to_cells(self.origin, end_point)
        self._x = xy[0]
        self._y = xy[1]
