# Copyright Crown and Cartopy Contributors
#
# This file is part of Cartopy and is released under the BSD 3-clause license.
# See LICENSE in the root of the repository for full licensing details.

import copy
from io import BytesIO
import os
from pathlib import Path
import pickle
import warnings

import numpy as np
from numpy.testing import assert_almost_equal, assert_array_equal
from numpy.testing import assert_array_almost_equal as assert_arr_almost_eq
import pyproj
import pytest
import shapely.geometry as sgeom

import cartopy.crs as ccrs


class TestCRS:
    def test_hash(self):
        stereo = ccrs.Stereographic(90)
        north = ccrs.NorthPolarStereo()
        assert stereo == north
        assert not stereo != north
        assert hash(stereo) == hash(north)

        assert ccrs.Geodetic() == ccrs.Geodetic()

    @pytest.mark.parametrize('approx', [True, False])
    def test_osni(self, approx):
        osni = ccrs.OSNI(approx=approx)
        ll = ccrs.Geodetic()

        # results obtained by nearby.org.uk.
        lon, lat = np.array([-5.54159863617957, 54.5622169298669],
                            dtype=np.double)
        east, north = np.array([359000, 371000], dtype=np.double)

        assert_arr_almost_eq(osni.transform_point(lon, lat, ll),
                             np.array([east, north]),
                             -1)
        assert_arr_almost_eq(ll.transform_point(east, north, osni),
                             np.array([lon, lat]),
                             3)

    def _check_osgb(self, osgb):
        precision = 1

        if os.environ.get('PROJ_NETWORK') != 'ON':
            grid_name = 'uk_os_OSTN15_NTv2_OSGBtoETRS.tif'
            available = (
                Path(pyproj.datadir.get_data_dir(), grid_name).exists() or
                Path(pyproj.datadir.get_user_data_dir(), grid_name).exists()
            )
            if not available:
                import warnings
                warnings.warn(f'{grid_name} is unavailable; '
                              'testing OSGB at reduced precision')
                precision = -1

        ll = ccrs.Geodetic()

        # results obtained by streetmap.co.uk.
        lon, lat = np.array([-3.478831, 50.462023], dtype=np.double)
        east, north = np.array([295132.1, 63512.6], dtype=np.double)

        # note the handling of precision here...
        assert_almost_equal(osgb.transform_point(lon, lat, ll), [east, north],
                            decimal=precision)
        assert_almost_equal(ll.transform_point(east, north, osgb), [lon, lat],
                            decimal=2)

        r_lon, r_lat = ll.transform_point(east, north, osgb)
        r_inverted = np.array(osgb.transform_point(r_lon, r_lat, ll))
        assert_arr_almost_eq(r_inverted, [east, north], 3)

        r_east, r_north = osgb.transform_point(lon, lat, ll)
        r_inverted = np.array(ll.transform_point(r_east, r_north, osgb))
        assert_arr_almost_eq(r_inverted, [lon, lat])

    @pytest.mark.parametrize('approx', [True, False])
    def test_osgb(self, approx):
        self._check_osgb(ccrs.OSGB(approx=approx))

    def test_epsg(self):
        uk = ccrs.epsg(27700)
        assert uk.epsg_code == 27700
        expected_x = (-104009.357, 688806.007)
        expected_y = (-8908.37, 1256558.45)
        expected_threshold = 7928.15
        if pyproj.__proj_version__ >= '9.2.0':
            expected_x = (-104728.764, 688806.007)
            expected_y = (-8908.36, 1256616.32)
            expected_threshold = 7935.34
        assert_almost_equal(uk.x_limits,
                            expected_x, decimal=3)
        assert_almost_equal(uk.y_limits, expected_y, decimal=2)
        assert_almost_equal(uk.threshold, expected_threshold, decimal=2)
        self._check_osgb(uk)

    def test_epsg_compound_crs(self):
        projection = ccrs.epsg(5973)
        assert projection.epsg_code == 5973

    def test_europp(self):
        europp = ccrs.EuroPP()
        proj4_init = europp.proj4_init
        # Transverse Mercator, UTM zone 32,
        assert '+proj=utm' in proj4_init
        assert '+zone=32' in proj4_init
        # International 1924 ellipsoid.
        assert '+ellps=intl' in proj4_init

    def test_transform_points_nD(self):
        rlons = np.array([[350., 352., 354.], [350., 352., 354.]])
        rlats = np.array([[-5., -0., 1.], [-4., -1., 0.]])

        src_proj = ccrs.RotatedGeodetic(pole_longitude=178.0,
                                        pole_latitude=38.0)
        target_proj = ccrs.Geodetic()
        res = target_proj.transform_points(x=rlons, y=rlats,
                                           src_crs=src_proj)
        unrotated_lon = res[..., 0]
        unrotated_lat = res[..., 1]

        # Solutions derived by proj direct.
        solx = np.array([[-16.42176094, -14.85892262, -11.90627520],
                         [-16.71055023, -14.58434624, -11.68799988]])
        soly = np.array([[46.00724251, 51.29188893, 52.59101488],
                         [46.98728486, 50.30706042, 51.60004528]])
        assert_arr_almost_eq(unrotated_lon, solx)
        assert_arr_almost_eq(unrotated_lat, soly)

    def test_transform_points_1D(self):
        rlons = np.array([350., 352., 354., 356.])
        rlats = np.array([-5., -0., 5., 10.])

        src_proj = ccrs.RotatedGeodetic(pole_longitude=178.0,
                                        pole_latitude=38.0)
        target_proj = ccrs.Geodetic()
        res = target_proj.transform_points(x=rlons, y=rlats,
                                           src_crs=src_proj)
        unrotated_lon = res[..., 0]
        unrotated_lat = res[..., 1]

        # Solutions derived by proj direct.
        solx = np.array([-16.42176094, -14.85892262,
                         -12.88946157, -10.35078336])
        soly = np.array([46.00724251, 51.29188893,
                         56.55031485, 61.77015703])

        assert_arr_almost_eq(unrotated_lon, solx)
        assert_arr_almost_eq(unrotated_lat, soly)

    def test_transform_points_xyz(self):
        # Test geodetic transforms when using z value
        rx = np.array([2574.32516e3])
        ry = np.array([837.562e3])
        rz = np.array([5761.325e3])

        src_proj = ccrs.Geocentric()
        target_proj = ccrs.Geodetic()

        res = target_proj.transform_points(x=rx, y=ry, z=rz,
                                           src_crs=src_proj)

        glat = res[..., 0]
        glon = res[..., 1]
        galt = res[..., 2]

        # Solution generated by pyproj
        solx = np.array([18.0224043189])
        soly = np.array([64.9796515089])
        solz = np.array([5048.03893734])

        assert_arr_almost_eq(glat, solx)
        assert_arr_almost_eq(glon, soly)
        assert_arr_almost_eq(galt, solz)

    def test_transform_points_180(self):
        # Test that values less than -180 and more than 180
        # get mapped to the -180, 180 interval
        x = np.array([-190, 190])
        y = np.array([0, 0])

        proj = ccrs.PlateCarree()

        res = proj.transform_points(x=x, y=y, src_crs=proj)
        assert_array_equal(res[..., :2], [[170, 0], [-170, 0]])

    def test_globe(self):
        # Ensure the globe affects output.
        rugby_globe = ccrs.Globe(semimajor_axis=9000000,
                                 semiminor_axis=9000000,
                                 ellipse=None)
        footy_globe = ccrs.Globe(semimajor_axis=1000000,
                                 semiminor_axis=1000000,
                                 ellipse=None)

        rugby_moll = ccrs.Mollweide(globe=rugby_globe)
        footy_moll = ccrs.Mollweide(globe=footy_globe)

        rugby_pt = rugby_moll.transform_point(
            10, 10, rugby_moll.as_geodetic(),
        )
        footy_pt = footy_moll.transform_point(
            10, 10, footy_moll.as_geodetic(),
        )

        assert_arr_almost_eq(rugby_pt, (1400915, 1741319), decimal=0)
        assert_arr_almost_eq(footy_pt, (155657, 193479), decimal=0)

    def test_project_point(self):
        point = sgeom.Point([0, 45])
        multi_point = sgeom.MultiPoint([point, sgeom.Point([180, 45])])

        pc = ccrs.PlateCarree()
        pc_rotated = ccrs.PlateCarree(central_longitude=180)

        result = pc_rotated.project_geometry(point, pc)
        assert_arr_almost_eq(result.xy, [[-180.], [45.]])

        result = pc_rotated.project_geometry(multi_point, pc)
        assert isinstance(result, sgeom.MultiPoint)
        assert len(result.geoms) == 2
        assert_arr_almost_eq(result.geoms[0].xy, [[-180.], [45.]])
        assert_arr_almost_eq(result.geoms[1].xy, [[0], [45.]])

    def test_utm(self):
        utm30n = ccrs.UTM(30)
        ll = ccrs.Geodetic()
        lon, lat = np.array([-3.0, 51.5], dtype=np.double)
        east, north = np.array([500000, 5705429.2], dtype=np.double)
        assert_arr_almost_eq(utm30n.transform_point(lon, lat, ll),
                             [east, north],
                             decimal=1)
        assert_arr_almost_eq(ll.transform_point(east, north, utm30n),
                             [lon, lat],
                             decimal=1)
        utm38s = ccrs.UTM(38, southern_hemisphere=True)
        lon, lat = np.array([47.5, -18.92], dtype=np.double)
        east, north = np.array([763316.7, 7906160.8], dtype=np.double)
        assert_arr_almost_eq(utm38s.transform_point(lon, lat, ll),
                             [east, north],
                             decimal=1)
        assert_arr_almost_eq(ll.transform_point(east, north, utm38s),
                             [lon, lat],
                             decimal=1)


@pytest.fixture(params=[
    [ccrs.PlateCarree, {}],
    [ccrs.PlateCarree, dict(central_longitude=1.23)],
    [ccrs.NorthPolarStereo, dict(central_longitude=42.5,
                                 globe=ccrs.Globe(ellipse="helmert"))],
    [ccrs.CRS, dict(proj4_params="3088")],
    [ccrs.epsg, dict(code="3088")]
])
def proj_to_copy(request):
    cls, kwargs = request.param
    return cls(**kwargs)


def test_pickle(proj_to_copy):
    # check that we can pickle a simple CRS
    fh = BytesIO()
    pickle.dump(proj_to_copy, fh)
    fh.seek(0)
    pickled_prj = pickle.load(fh)
    assert proj_to_copy == pickled_prj


def test_deepcopy(proj_to_copy):
    prj_cp = copy.deepcopy(proj_to_copy)
    assert proj_to_copy.proj4_params == prj_cp.proj4_params
    assert proj_to_copy == prj_cp


def test_PlateCarree_shortcut():
    central_lons = [[0, 0], [0, 180], [0, 10], [10, 0], [-180, 180], [
        180, -180]]

    target = [([[-180, -180], [-180, 180]], 0),
              ([[-180, 0], [0, 180]], 180),
              ([[-180, -170], [-170, 180]], 10),
              ([[-180, 170], [170, 180]], -10),
              ([[-180, 180], [180, 180]], 360),
              ([[-180, -180], [-180, 180]], -360),
              ]

    assert len(target) == len(central_lons)

    for expected, (s_lon0, t_lon0) in zip(target, central_lons):
        expected_bboxes, expected_offset = expected

        src = ccrs.PlateCarree(central_longitude=s_lon0)
        target = ccrs.PlateCarree(central_longitude=t_lon0)

        bbox, offset = src._bbox_and_offset(target)

        assert offset == expected_offset
        assert bbox == expected_bboxes


def test_transform_points_empty():
    """Test CRS.transform_points with empty array."""
    crs = ccrs.Stereographic()
    result = crs.transform_points(ccrs.PlateCarree(),
                                  np.array([]), np.array([]))
    assert_array_equal(result, np.array([], dtype=np.float64).reshape(0, 3))


def test_transform_points_outside_domain():
    """Test CRS.transform_points with out of domain arrays."""
    # Length-1 arrays error out with a bad status code, while
    # greater than 1 arrays put infinity into the return array
    # where the bad values occur
    crs = ccrs.Orthographic()
    result = crs.transform_points(ccrs.PlateCarree(),
                                  np.array([-120]), np.array([80]))
    assert np.all(np.isnan(result))
    result = crs.transform_points(ccrs.PlateCarree(),
                                  np.array([-120]), np.array([80]),
                                  trap=True)
    assert np.all(np.isnan(result))
    # A length-2 array of the same transform produces "inf" rather
    # than nan due to PROJ never returning nan itself.
    result = crs.transform_points(ccrs.PlateCarree(),
                                  np.array([-120, -120]), np.array([80, 80]))
    assert np.all(~np.isfinite(result[..., :2]))

    # Test singular transform to make sure it is producing all nan's
    # the same as the transform_points call with a length-1 array
    result = crs.transform_point(-120, 80, ccrs.PlateCarree())
    assert np.all(np.isnan(result))


def test_projection__from_string():
    crs = ccrs.Projection("NAD83 / Pennsylvania South")
    assert crs.as_geocentric().datum.name == "North American Datum 1983"
    assert_almost_equal(
        crs.bounds,
        [361633.1351868, 859794.6690229, 45575.5693199, 209415.9845754],
    )


def test_crs__from_pyproj_crs():
    assert ccrs.CRS(pyproj.CRS("EPSG:4326")) == "EPSG:4326"


def test_transform_point_no_warning():
    # Make sure we aren't warning on single-point numpy arrays
    # see https://github.com/SciTools/cartopy/pull/2194
    p = ccrs.PlateCarree()
    p2 = ccrs.Mercator()
    with warnings.catch_warnings():
        warnings.simplefilter("error")
        p2.transform_point(1, 2, p)


def test_geographic_bounds_no_area_of_use():
    # Should default to global bounds if no area of use
    wkt = ('GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",'
           'SPHEROID["WGS_1984",6378137.0,298.257223563]],'
           'PRIMEM["Greenwich",0.0],UNIT["Degree",0.017453292519943295]]')
    p = pyproj.CRS.from_wkt(wkt)
    assert p.is_geographic
    assert p.area_of_use is None
    p_cartopy = ccrs.Projection(p)
    assert p_cartopy.bounds == (-180, 180, -90, 90)


def test_geographic_bounds_with_area_of_use():
    # Should default to area of use bounds if available
    p = pyproj.CRS.from_epsg(4267)
    assert p.is_geographic
    assert p.area_of_use is not None
    p_cartopy = ccrs.Projection(p)
    x0, x1, y0, y1 = p_cartopy.bounds
    assert_array_equal((x1, y0, x0, y1), p.area_of_use.bounds)