# (C) British Crown Copyright 2011 - 2019, Met Office
#
# This file is part of cartopy.
#
# cartopy is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# cartopy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with cartopy.  If not, see <https://www.gnu.org/licenses/>.

from __future__ import (absolute_import, division, print_function)

import numpy as np
import pytest
import shapely.geometry as sgeom
import shapely.wkt


import cartopy.crs as ccrs


class TestBoundary(object):
    def test_no_polygon_boundary_reversal(self):
        # Check that polygons preserve their clockwise or counter-clockwise
        # ordering when they are attached to the boundary.
        # Failure to do so will result in invalid polygons (their boundaries
        # cross-over).
        polygon = sgeom.Polygon([(-10, 30), (10, 60), (10, 50)])
        projection = ccrs.Robinson(170.5)
        multi_polygon = projection.project_geometry(polygon)
        for polygon in multi_polygon:
            assert polygon.is_valid

    def test_polygon_boundary_attachment(self):
        # Check the polygon is attached to the boundary even when no
        # intermediate point for one of the crossing segments would normally
        # exist.
        polygon = sgeom.Polygon([(-10, 30), (10, 60), (10, 50)])
        projection = ccrs.Robinson(170.6)
        # This will raise an exception if the polygon/boundary intersection
        # fails.
        projection.project_geometry(polygon)

    def test_out_of_bounds(self):
        # Check that a polygon that is completely out of the map boundary
        # doesn't produce an empty result.
        projection = ccrs.TransverseMercator(central_longitude=0, approx=True)

        polys = [
            # All valid
            ([(86, -1), (86, 1), (88, 1), (88, -1)], 1),
            # One out of backwards projection range
            ([(86, -1), (86, 1), (130, 1), (88, -1)], 1),
            # An out of backwards projection range segment
            ([(86, -1), (86, 1), (130, 1), (130, -1)], 1),
            # All out of backwards projection range
            ([(120, -1), (120, 1), (130, 1), (130, -1)], 0),
        ]

        # Try all four combinations of valid/NaN vs valid/NaN.
        for coords, expected_polys in polys:
            polygon = sgeom.Polygon(coords)
            multi_polygon = projection.project_geometry(polygon)
            assert len(multi_polygon) == expected_polys


class TestMisc(object):
    def test_misc(self):
        projection = ccrs.TransverseMercator(central_longitude=-90,
                                             approx=False)
        polygon = sgeom.Polygon([(-10, 30), (10, 60), (10, 50)])
        projection.project_geometry(polygon)

    def test_small(self):
        projection = ccrs.Mercator()
        polygon = sgeom.Polygon([
            (-179.7933201090486079, -16.0208822567412312),
            (-180.0000000000000000, -16.0671326636424396),
            (-179.9173693847652942, -16.5017831356493616),
        ])
        multi_polygon = projection.project_geometry(polygon)
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].exterior.coords) == 4

    def test_former_infloop_case(self):
        # test a polygon which used to get stuck in an infinite loop
        # see https://github.com/SciTools/cartopy/issues/60
        coords = [(260.625, 68.90383337092122), (360.0, 79.8556091996901),
                  (360.0, 77.76848175458498), (0.0, 88.79068047337279),
                  (210.0, 90.0), (135.0, 88.79068047337279),
                  (260.625, 68.90383337092122)]
        geom = sgeom.Polygon(coords)

        target_projection = ccrs.PlateCarree()
        source_crs = ccrs.Geodetic()

        multi_polygon = target_projection.project_geometry(geom, source_crs)
        # check the result is non-empty
        assert not multi_polygon.is_empty

    def test_project_previous_infinite_loop(self):
        mstring1 = shapely.wkt.loads(
            'MULTILINESTRING ('
            '(-179.9999990464349651 -80.2000000000000171, '
            '-179.5000000001111005 -80.2000000000000171, '
            '-179.5000000001111005 -79.9000000000000199, '
            '-179.9999995232739138 -79.9499999523163041, '
            '-179.8000000001110550 -80.0000000000000000, '
            '-179.8000000001110550 -80.0999999999999943, '
            '-179.9999999047436177 -80.0999999999999943), '
            '(179.9999995231628702 -79.9499999523163041, '
            '179.5000000000000000 -79.9000000000000199, '
            '179.5000000000000000 -80.0000000000000000, '
            '179.9999995231628702 -80.0499999523162842, '
            '179.5000000000000000 -80.0999999999999943, '
            '179.5000000000000000 -80.2000000000000171, '
            '179.9999990463256836 -80.2000000000000171))')
        mstring2 = shapely.wkt.loads(
            'MULTILINESTRING ('
            '(179.9999996185302678 -79.9999999904632659, '
            '179.5999999999999943 -79.9899999999999949, '
            '179.5999999999999943 -79.9399999999999977, '
            '179.9999996185302678 -79.9599999809265114), '
            '(-179.9999999047436177 -79.9600000000000080, '
            '-179.9000000001110777 -79.9600000000000080, '
            '-179.9000000001110777 -80.0000000000000000, '
            '-179.9999999047436177 -80.0000000000000000))')
        multi_line_strings = [mstring1, mstring2]

        src = ccrs.PlateCarree()
        src._attach_lines_to_boundary(multi_line_strings, True)

    @pytest.mark.parametrize('proj',
                             [ccrs.InterruptedGoodeHomolosine, ccrs.Mollweide])
    def test_infinite_loop_bounds(self, proj):
        # test a polygon which used to get stuck in an infinite loop but is now
        # erroneously clipped.
        # see https://github.com/SciTools/cartopy/issues/1131

        # These names are for IGH; effectively the same for Mollweide.
        bottom = [0., 70.]
        right = [0., 90.]
        top = [-180., 90.]
        left = [-180., 70.]
        verts = np.array([
            bottom,
            right,
            top,
            left,
            bottom,
        ])
        bad_path = sgeom.Polygon(verts)

        target = proj()
        source = ccrs.PlateCarree()

        projected = target.project_geometry(bad_path, source)

        # When transforming segments was broken, the resulting path did not
        # close, and either filled most of the domain, or a smaller portion
        # than it should. Check that the bounds match the individual points at
        # the expected edges.
        projected_left = target.transform_point(left[0], left[1], source)
        assert projected.bounds[0] == pytest.approx(projected_left[0],
                                                    rel=target.threshold)
        projected_bottom = target.transform_point(bottom[0], bottom[1], source)
        assert projected.bounds[1] == pytest.approx(projected_bottom[1],
                                                    rel=target.threshold)
        projected_right = target.transform_point(right[0], right[1], source)
        assert projected.bounds[2] == pytest.approx(projected_right[0],
                                                    rel=target.threshold,
                                                    abs=1e-8)
        projected_top = target.transform_point(top[0], top[1], source)
        assert projected.bounds[3] == pytest.approx(projected_top[1],
                                                    rel=target.threshold)

    def test_3pt_poly(self):
        projection = ccrs.OSGB(approx=True)
        polygon = sgeom.Polygon([(-1000, -1000),
                                 (-1000, 200000),
                                 (200000, -1000)])
        multi_polygon = projection.project_geometry(polygon,
                                                    ccrs.OSGB(approx=True))
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].exterior.coords) == 4

    def test_self_intersecting_1(self):
        # Geometry comes from a matplotlib contourf (see #537)
        wkt = ('POLYGON ((366.22000122 -9.71489298, '
               '366.73212393 -9.679999349999999, '
               '366.77412634 -8.767753000000001, '
               '366.17762962 -9.679999349999999, '
               '366.22000122 -9.71489298), '
               '(366.22000122 -9.692636309999999, '
               '366.32998657 -9.603356099999999, '
               '366.74765799 -9.019999500000001, '
               '366.5094086 -9.63175386, '
               '366.22000122 -9.692636309999999))')
        geom = shapely.wkt.loads(wkt)
        source, target = ccrs.RotatedPole(198.0, 39.25), ccrs.EuroPP()
        projected = target.project_geometry(geom, source)
        # Before handling self intersecting interiors, the area would be
        # approximately 13262233761329.
        area = projected.area
        assert 2.2e9 < area < 2.3e9, \
            'Got area {}, expecting ~2.2e9'.format(area)

    def test_self_intersecting_2(self):
        # Geometry comes from a matplotlib contourf (see #509)
        wkt = ('POLYGON ((343 20, 345 23, 342 25, 343 22, '
               '340 25, 341 25, 340 25, 343 20), (343 21, '
               '343 22, 344 23, 343 21))')
        geom = shapely.wkt.loads(wkt)
        source = target = ccrs.RotatedPole(193.0, 41.0)
        projected = target.project_geometry(geom, source)
        # Before handling self intersecting interiors, the area would be
        # approximately 64808.
        assert 7.9 < projected.area < 8.1

    def test_tiny_point_between_boundary_points(self):
        # Geometry comes from #259.
        target = ccrs.Orthographic(0, -75)
        source = ccrs.PlateCarree()
        wkt = 'POLYGON ((132 -40, 133 -6, 125.3 1, 115 -6, 132 -40))'
        geom = shapely.wkt.loads(wkt)

        target = ccrs.Orthographic(central_latitude=90., central_longitude=0)
        source = ccrs.PlateCarree()
        projected = target.project_geometry(geom, source)
        area = projected.area
        # Before fixing, this geometry used to fill the whole disk. Approx
        # 1.2e14.
        assert 81330 < area < 81340, \
            'Got area {}, expecting ~81336'.format(area)

    def test_same_points_on_boundary_1(self):
        source = ccrs.PlateCarree()
        target = ccrs.PlateCarree(central_longitude=180)

        geom = sgeom.Polygon([(-20, -20), (20, -20), (20, 20), (-20, 20)],
                             [[(-10, 0), (0, 20), (10, 0), (0, -20)]])
        projected = target.project_geometry(geom, source)

        assert abs(1200 - projected.area) < 1e-5

    def test_same_points_on_boundary_2(self):
        source = ccrs.PlateCarree()
        target = ccrs.PlateCarree(central_longitude=180)

        geom = sgeom.Polygon([(-20, -20), (20, -20), (20, 20), (-20, 20)],
                             [[(0, 0), (-10, 10), (0, 20), (10, 10)],
                              [(0, 0), (10, -10), (0, -20), (-10, -10)]])
        projected = target.project_geometry(geom, source)

        assert abs(1200 - projected.area) < 1e-5


class TestQuality(object):
    def setup_class(self):
        projection = ccrs.RotatedPole(pole_longitude=177.5,
                                      pole_latitude=37.5)
        polygon = sgeom.Polygon([
            (177.5, -57.38460319),
            (180.0, -57.445077),
            (175.0, -57.19913331),
        ])
        self.multi_polygon = projection.project_geometry(polygon)
        # from cartopy.tests.mpl import show
        # show(projection, self.multi_polygon)

    def test_split(self):
        # Start simple ... there should be two projected polygons.
        assert len(self.multi_polygon) == 2

    def test_repeats(self):
        # Make sure we don't have repeated points at the boundary, because
        # they mess up the linear extrapolation to the boundary.

        # Make sure there aren't any repeated points.
        xy = np.array(self.multi_polygon[0].exterior.coords)
        same = (xy[1:] == xy[:-1]).all(axis=1)
        assert not any(same), 'Repeated points in projected geometry.'

    def test_symmetry(self):
        # Make sure the number of points added on the way towards the
        # boundary is similar to the number of points added on the way away
        # from the boundary.

        # Identify all the contiguous sets of non-boundary points.
        xy = np.array(self.multi_polygon[0].exterior.coords)
        boundary = np.logical_or(xy[:, 1] == 90, xy[:, 1] == -90)
        regions = (boundary[1:] != boundary[:-1]).cumsum()
        regions = np.insert(regions, 0, 0)

        # For each region, check if the number of increasing steps is roughly
        # equal to the number of decreasing steps.
        for i in range(int(boundary[0]), regions.max(), 2):
            indices = np.where(regions == i)
            x = xy[indices, 0]
            delta = np.diff(x)
            num_incr = np.count_nonzero(delta > 0)
            num_decr = np.count_nonzero(delta < 0)
            assert abs(num_incr - num_decr) < 3, 'Too much asymmetry.'


class PolygonTests(object):
    def _assert_bounds(self, bounds, x1, y1, x2, y2, delta=1):
        assert abs(bounds[0] - x1) < delta
        assert abs(bounds[1] - y1) < delta
        assert abs(bounds[2] - x2) < delta
        assert abs(bounds[3] - y2) < delta


class TestWrap(PolygonTests):
    # Test that Plate Carree projection "does the right thing"(tm) with
    # source data tha extends outside the [-180, 180] range.
    def test_plate_carree_no_wrap(self):
        proj = ccrs.PlateCarree()
        poly = sgeom.box(0, 0, 10, 10)
        multi_polygon = proj.project_geometry(poly, proj)
        # Check the structure
        assert len(multi_polygon) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, 0, 0, 10, 10)

    def test_plate_carree_partial_wrap(self):
        proj = ccrs.PlateCarree()
        poly = sgeom.box(170, 0, 190, 10)
        multi_polygon = proj.project_geometry(poly, proj)
        # Check the structure
        assert len(multi_polygon) == 2
        # Check the rough shape
        poly1, poly2 = multi_polygon
        # The order of these polygons is not guaranteed, so figure out
        # which is appropriate
        if 170.0 not in poly1.bounds:
            poly1, poly2 = poly2, poly1
        self._assert_bounds(poly1.bounds, 170, 0, 180, 10)
        self._assert_bounds(poly2.bounds, -180, 0, -170, 10)

    def test_plate_carree_wrap(self):
        proj = ccrs.PlateCarree()
        poly = sgeom.box(200, 0, 220, 10)
        multi_polygon = proj.project_geometry(poly, proj)
        # Check the structure
        assert len(multi_polygon) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, -160, 0, -140, 10)


def ring(minx, miny, maxx, maxy, ccw):
    box = sgeom.box(minx, miny, maxx, maxy, ccw)
    return np.array(box.exterior.coords)


class TestHoles(PolygonTests):
    def test_simple(self):
        proj = ccrs.PlateCarree()
        poly = sgeom.Polygon(ring(-40, -40, 40, 40, True),
                             [ring(-20, -20, 20, 20, False)])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].interiors) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, -40, -47, 40, 47)
        self._assert_bounds(polygon.interiors[0].bounds, -20, -21, 20, 21)

    def test_wrapped_poly_simple_hole(self):
        proj = ccrs.PlateCarree(-150)
        poly = sgeom.Polygon(ring(-40, -40, 40, 40, True),
                             [ring(-20, -20, 20, 20, False)])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 2

        poly1, poly2 = multi_polygon
        # The order of these polygons is not guaranteed, so figure out
        # which is appropriate
        if not len(poly1.interiors) == 1:
            poly1, poly2 = poly2, poly1

        assert len(poly1.interiors) == 1
        assert len(poly2.interiors) == 0
        # Check the rough shape
        self._assert_bounds(poly1.bounds, 110, -47, 180, 47)
        self._assert_bounds(poly1.interiors[0].bounds, 130, -21, 170, 21)
        self._assert_bounds(poly2.bounds, -180, -43, -170, 43)

    def test_wrapped_poly_wrapped_hole(self):
        proj = ccrs.PlateCarree(-180)
        poly = sgeom.Polygon(ring(-40, -40, 40, 40, True),
                             [ring(-20, -20, 20, 20, False)])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 2
        assert len(multi_polygon[0].interiors) == 0
        assert len(multi_polygon[1].interiors) == 0
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, 140, -47, 180, 47)
        polygon = multi_polygon[1]
        self._assert_bounds(polygon.bounds, -180, -47, -140, 47)

    def test_inverted_poly_simple_hole(self):
        proj = ccrs.NorthPolarStereo()
        poly = sgeom.Polygon([(0, 0), (-90, 0), (-180, 0), (-270, 0)],
                             [[(0, -30), (90, -30), (180, -30), (270, -30)]])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].interiors) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, -2.4e7, -2.4e7, 2.4e7, 2.4e7, 1e6)
        self._assert_bounds(polygon.interiors[0].bounds,
                            - 1.2e7, -1.2e7, 1.2e7, 1.2e7, 1e6)

    def test_inverted_poly_clipped_hole(self):
        proj = ccrs.NorthPolarStereo()
        poly = sgeom.Polygon([(0, 0), (-90, 0), (-180, 0), (-270, 0)],
                             [[(-135, -60), (-45, -60),
                               (45, -60), (135, -60)]])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].interiors) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, -5.0e7, -5.0e7, 5.0e7, 5.0e7, 1e6)
        self._assert_bounds(polygon.interiors[0].bounds,
                            - 1.2e7, -1.2e7, 1.2e7, 1.2e7, 1e6)
        assert abs(polygon.area - 7.30e15) < 1e13

    def test_inverted_poly_removed_hole(self):
        proj = ccrs.NorthPolarStereo(globe=ccrs.Globe(ellipse='WGS84'))
        poly = sgeom.Polygon([(0, 0), (-90, 0), (-180, 0), (-270, 0)],
                             [[(-135, -75), (-45, -75),
                               (45, -75), (135, -75)]])
        multi_polygon = proj.project_geometry(poly)
        # Check the structure
        assert len(multi_polygon) == 1
        assert len(multi_polygon[0].interiors) == 1
        # Check the rough shape
        polygon = multi_polygon[0]
        self._assert_bounds(polygon.bounds, -5.0e7, -5.0e7, 5.0e7, 5.0e7, 1e6)
        self._assert_bounds(polygon.interiors[0].bounds,
                            - 1.2e7, -1.2e7, 1.2e7, 1.2e7, 1e6)
        assert abs(polygon.area - 7.34e15) < 1e13

    def test_multiple_interiors(self):
        exterior = ring(0, 0, 12, 12, True)
        interiors = [ring(1, 1, 2, 2, False), ring(1, 8, 2, 9, False)]

        poly = sgeom.Polygon(exterior, interiors)

        target = ccrs.PlateCarree()
        source = ccrs.Geodetic()

        assert len(list(target.project_geometry(poly, source))) == 1