# pylint: disable=invalid-name
"""
Model domains
=============

Regular latitude/longitude
--------------------------

For regular latitude/longitude grids a simple bounding box test is all
that is required to determine if a point lies inside the domain.

Irregular boundaries
--------------------

Regional ocean models with irregular boundaries can perform an additional
point in polygon check to determine if a point is inside the domain.

Global models
-------------

Ocean models of global extent typically have a southern boundary, since
Antarctica is a land mass covering the South Pole. A North/South extent
check may be sufficient to determine whether a point belongs to the domain or
not.
"""
from itertools import izip
import numpy as np


def inside(grid_longitudes,
           grid_latitudes,
           observed_longitudes,
           observed_latitudes,
           kind="regular"):
    """Detect points inside a model domain

    Apply geometry techniques to determine if an observation is inside a model
    domain. Strategies included can handle regular lon/lat, global latitude
    bands and polygonal grids.

    ======= ===============================================
    Key     Description
    ======= ===============================================
    regular model grids bounded by a longitude/latitude box
    band    model grid bounded north/south
    polygon general case where model boundary treated as
            polygon
    ======= ===============================================

    :param grid_longitudes: 2D array shaped (X, Y)
    :param grid_latitudes: 2D array shaped (X, Y)
    :param observed_longitudes: 1D array shaped (N,)
    :param observed_latitudes: 1D array shaped (N,)
    :param kind: string indicating which kind of model grid is being used
    :returns: boolean array where True indicates observation is inside domain
    """
    if kind.lower() in ["regular", "box", "lonlat"]:
        return Box.from2d(grid_longitudes,
                          grid_latitudes).inside(observed_longitudes,
                                                 observed_latitudes)
    elif kind.lower() in ["band", "latitude_band"]:
        return LatitudeBand.from_latitudes(grid_latitudes).inside(observed_latitudes)
    elif kind.lower() in ["polygon", "irregular", "rotated"]:
        return Polygon(grid_longitudes,
                       grid_latitudes).inside(observed_longitudes,
                                              observed_latitudes)


class LatitudeBand(object):
    """Latitude band domain

    Positions within two latitude parallels are considered inside the domain.

    :param minimum: southernmost latitude
    :param maximum: northernmost latitude
    """
    def __init__(self, minimum, maximum):
        self.minimum = minimum
        self.maximum = maximum

    @classmethod
    def from_latitudes(cls, latitudes):
        """Construct latitude band from latitude array"""
        return cls(np.ma.min(latitudes),
                   np.ma.max(latitudes))

    def inside(self, latitudes):
        """Determine observations inside grid"""
        return ((latitudes >= self.minimum) &
                (latitudes <= self.maximum))


class Polygon(object):
    """General ocean grid definition

    Arbitrary shaped non-intersecting polygon domain. Uses
    a combination of bounding box and point in polygon
    algorithms to decide if a point is inside the domain.

    :param longitudes: 2D array shaped (X, Y)
    :param latitudes: 2D array shaped (X, Y)
    """
    def __init__(self, longitudes, latitudes):
        self.bounding_box = Box.from2d(longitudes,
                                       latitudes)
        self.point_in_polygon = PointInPolygon.from2d(longitudes,
                                                      latitudes)

    def inside(self, longitudes, latitudes):
        """Check observations are contained within domain

        :param longitudes: x-coordinate to test
        :param latitudes: y-coordinate to test
        :returns: logical indicating coordinates contained in polygon
        """
        longitudes = np.asarray(longitudes, dtype="d")
        latitudes = np.asarray(latitudes, dtype="d")

        if longitudes.ndim == 0:
            return self.point_in_polygon.inside(longitudes,
                                                latitudes)
        # Optimal vector algorithm
        result = np.zeros_like(longitudes, dtype=np.bool)
        in_box = self.bounding_box.inside(longitudes, latitudes)
        result[in_box] = self.point_in_polygon.inside(longitudes[in_box],
                                                      latitudes[in_box])
        return result


class Box(object):
    """Bounding box surrounding collection of vertices"""
    def __init__(self, xmin, xmax, ymin, ymax):
        self.xmin = xmin
        self.xmax = xmax
        self.ymin = ymin
        self.ymax = ymax

    @classmethod
    def from2d(cls, longitudes, latitudes):
        """Construct from 2D coordinates"""
        return cls(np.min(longitudes),
                   np.max(longitudes),
                   np.min(latitudes),
                   np.max(latitudes))

    def inside(self, x, y):
        """Check point lies inside box

        :param x: x-coordinate to test
        :param y: y-coordinate to test
        :returns: logical indicating coordinates contained in box
        """
        return ((x <= self.xmax) &
                (x >= self.xmin) &
                (y <= self.ymax) &
                (y >= self.ymin))


def boundary(values):
    """Extract boundary from grid

    A boundary of a grid shaped N x M consists of 2N + 2M - 4 points.
    2 rows, 2 columns minus 4 corners that have been double counted.

    :param values: 2D array shaped (N, M)
    :returns: array shaped (B, 2) where B is the number of points on the