
# -*- coding: utf-8 -*-

u'''Universal Polar Stereographic (UPS) classes L{Ups} and L{UPSError}
and functions L{parseUPS5}, L{toUps8} and L{upsZoneBand5}.

A pure Python implementation, partially transcribed from C++ class U{PolarStereographic
<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1PolarStereographic.html>}
by I{Charles Karney}.

The U{UPS<https://WikiPedia.org/wiki/Universal_polar_stereographic_coordinate_system>}
system is used in conjuction with U{UTM
<https://WikiPedia.org/wiki/Universal_Transverse_Mercator_coordinate_system>}
to location on the polar regions of the earth.  UPS covers areas south of 79.5°S
and north of 83.5°N (slightly overlapping the UTM range from 80°S to 84°N).

@newfield example: Example, Examples
'''

from pygeodesy.datum import Datums, _TOL
from pygeodesy.dms import clipDMS, degDMS, parseDMS2, _parseUTMUPS, RangeError
from pygeodesy.fmath import EPS, hypot, hypot1
from pygeodesy.lazily import _ALL_LAZY
from pygeodesy.named import EasNor2Tuple, UtmUps5Tuple, UtmUps8Tuple, \
                            UtmUpsLatLon5Tuple, _xnamed
from pygeodesy.utily import degrees90, degrees180, property_RO, sincos2d
from pygeodesy.utmupsBase import _LLEB, _hemi, _to4lldn, _to3zBhp, _to3zll, \
                                 _UPS_LAT_MAX, _UPS_LAT_MIN, _UPS_ZONE, \
                                 _UPS_ZONE_STR, UtmUpsBase

from math import atan, atan2, radians, sqrt, tan

# all public contants, classes and functions
__all__ = _ALL_LAZY.ups
__version__ = '20.01.22'

_Bands   = 'A', 'B', 'Y', 'Z'    #: (INTERNAL) Polar bands.
_Falsing = 2000e3  #: (INTERNAL) False easting and northing (C{meter}).
_K0      = 0.994   #: (INTERNAL) Central UPS scale factor.
_K1      = 1.0     #: (INTERNAL) Rescale point scale factor.


class UPSError(ValueError):
    '''Universal Polar Stereographic (UPS) parse or other L{Ups} issue.
    '''
    pass


def _Band(a, b):
    # determine the polar band letter
    return _Bands[(0 if a < 0 else 2) + (0 if b < 0 else 1)]


def _scale(E, rho, tau):
    # compute the point scale factor, ala Karney
    t = hypot1(tau)
    return (rho / E.a) * t * sqrt(E.e12 + E.e2 / t**2)


class Ups(UtmUpsBase):
    '''Universal Polar Stereographic (UPS) coordinate.
    '''
    _band        = ''    #: (INTERNAL) Polar band ('A', 'B', 'Y' or 'Z').
    _latlon_args = True  #: (INTERNAL) unfalse from _latlon (C{bool}).
    _pole        = ''    #: (INTERNAL) UPS projection top/center ('N' or 'S').
    _scale       = None  #: (INTERNAL) Point scale factor (C{scalar}).
    _scale0      = _K0   #: (INTERNAL) Central scale factor (C{scalar}).
    _utm         = None  #: (INTERNAL) toUtm cache (L{Utm}).

    def __init__(self, zone, pole, easting, northing, band='',  # PYCHOK expected
                                   datum=Datums.WGS84, falsed=True,
                                   convergence=None, scale=None, name=''):
        '''New L{Ups} UPS coordinate.

           @param zone: UPS zone (C{int}, zero) or zone with/-out Band
                        letter (C{str}, '00', '00A', '00B', '00Y' or '00Z').
           @param pole: Top/center of (stereographic) projection
                        (C{str}, C{'N[orth]'} or C{'S[outh]'}).
           @param easting: Easting, see B{C{falsed}} (C{meter}).
           @param northing: Northing, see B{C{falsed}} (C{meter}).
           @keyword band: Optional, polar Band (C{str}, 'A'|'B'|'Y'|'Z').
           @keyword datum: Optional, this coordinate's datum (L{Datum}).
           @keyword falsed: Both B{C{easting}} and B{C{northing}} are
                            falsed (C{bool}).
           @keyword convergence: Optionally, save gamma meridian
                                 convergence (C{degrees}).
           @keyword scale: Optionally, save computed k scale (C{scalar}).
           @keyword name: Optional name (C{str}).

           @raise UPSError: Invalid B{C{zone}}, B{C{pole}} or B{C{band}}.
        '''
        if name:
            self.name = name

        try:
            z, B, p = _to3zBhp(zone, band=band, hemipole=pole)
            if z != _UPS_ZONE or (B and B not in _Bands):
                raise ValueError
        except ValueError:
            raise UPSError('%s, %s or %s invalid: %r' %
                           ('zone', 'pole', 'band', (zone, pole, band)))

        self._band        = B
        self._convergence = convergence
        self._easting     = float(easting)
        self._falsed      = falsed
        self._northing    = float(northing)
        self._pole        = p
        self._datum       = datum
        self._scale       = scale

    def __eq__(self, other):
        return isinstance(other, Ups) and other.zone     == self.zone \
                                      and other.pole     == self.pole \
                                      and other.easting  == self.easting \
                                      and other.northing == self.northing \
                                      and other.band     == self.band \
                                      and other.datum    == self.datum

    @property_RO
    def band(self):
        '''Get the polar band letter ('A', 'B', 'Y' or 'Z').
        '''
        if not self._band:
            self.toLatLon(unfalse=True)
        return self._band

    @property_RO
    def falsed2(self):
        '''Get the easting and northing falsing (L{EasNor2Tuple}C{(easting, northing)}).
        '''
        f = _Falsing if self.falsed else 0
        return EasNor2Tuple(f, f)

    @property_RO
    def hemisphere(self):
        '''Get the hemisphere (C{str}, 'N'|'S').
        '''
        if not self._hemisphere:
            self.toLatLon(unfalse=True)
        return self._hemisphere

    def parseUPS(self, strUPS):
        '''Parse a string to a UPS coordinate.

           @return: The coordinate (L{Ups}).

           @see: Function L{parseUPS5} in this module L{ups}.
        '''
        return parseUPS5(strUPS, datum=self.datum, Ups=self.classof)

    @property_RO
    def pole(self):
        '''Get the top/center of (stereographic) projection (C{'N'|'S'} or C{""}).
        '''
        return self._pole

    def rescale0(self, lat, scale0=_K0):
        '''Set the central scale factor for this UPS projection.

           @param lat: Northern latitude (C{degrees}).
           @param scale0: UPS k0 scale at B{C{lat}} latitude (C{scalar}).

           @raise RangeError: If B{C{lat}} outside the valid range
                              and L{rangerrors} set to C{True}.

           @raise ValueError: Invalid B{C{scale}}.
        '''
        try:
            s0 = float(scale0)
            if not 0 < s0:  # <= 1.003 or 1.0016?
                raise ValueError
        except (TypeError, ValueError):
            raise ValueError('%s invalid: %r' % ('scale', scale0))

        lat = clipDMS(lat, 90)  # clip and force N
        u = toUps8(abs(lat), 0, datum=self.datum, Ups=_UpsK1)
        k = s0 / u.scale
        if self.scale0 != k:
            self._band = ''  # force re-compute
            self._latlon = self._epsg = self._mgrs = self._utm = None
            self._scale0 = k

    def toLatLon(self, LatLon=None, unfalse=True):
        '''Convert this UPS coordinate to an (ellipsoidal) geodetic point.

           @keyword LatLon: Optional, ellipsoidal (sub-)class to return
                            the point (C{LatLon}) or C{None}.
           @keyword unfalse: Unfalse B{C{easting}} and B{C{northing}}
                             if falsed (C{bool}).

           @return: This UPS coordinate as (B{C{LatLon}}) or A
                    L{LatLonDatum5Tuple}C{(lat, lon, datum,
                    convergence, scale)} if B{C{LatLon}} is C{None}.

           @raise TypeError: If B{C{LatLon}} is not ellipsoidal.

           @raise UPSError: Invalid meridional radius or H-value.
        '''
        if self._latlon and self._latlon_args == unfalse:
            return self._latlon5(LatLon)

        E = self.datum.ellipsoid  # XXX vs LatLon.datum.ellipsoid

        x, y = self.to2en(falsed=not unfalse)

        r = hypot(x, y)
        t = (r / (2 * self.scale0 * E.a / E.es_c)) if r > 0 else EPS**2
        t = E.es_tauf((1 / t - t) * 0.5)
        if self._pole == 'N':
            a, b, c = atan(t), atan2(x, -y), 1
        else:
            a, b, c = -atan(t), atan2(x, y), -1

        a, b = degrees90(a), degrees180(b)
        if not self._band:
            self._band = _Band(a, b)
        if not self._hemisphere:
            self._hemisphere = _hemi(a)

        ll = _LLEB(a, b, datum=self._datum, name=self.name)
        ll._convergence = b * c  # gamma
        ll._scale = _scale(E, r, t) if r > 0 else self.scale0

        self._latlon_to(ll, unfalse)
        return self._latlon5(LatLon)

    def _latlon_to(self, ll, unfalse):
        '''(INTERNAL) See C{.toLatLon}, C{toUps8}.
        '''
        self._latlon, self._latlon_args = ll, unfalse

    def toMgrs(self):
        '''Convert this UPS coordinate to an MGRS grid reference.

           @return: The MGRS grid reference (L{Mgrs}).

           @see: Methods L{Ups.toUtm} and L{Utm.toMgrs}.
        '''
        if self._mgrs is None:
            self._mgrs = self.toUtm(None).toMgrs()  # via .toUtm
        return self._mgrs

    def toStr(self, prec=0, sep=' ', B=False, cs=False):  # PYCHOK expected
        '''Return a string representation of this UPS coordinate.

           Note that UPS coordinates are rounded, not truncated
           (unlike MGRS grid references).

           @keyword prec: Optional number of decimals, unstripped (C{int}).
           @keyword sep: Optional separator to join (C{str}).
           @keyword B: Optionally, include and polar band letter (C{bool}).
           @keyword cs: Optionally, include gamma meridian convergence
                        and point scale factor (C{bool}).

           @return: This UPS as a string with C{00[Band] pole, easting,
                    northing, [convergence, scale]} as C{"00[B] N|S
                    meter meter"} plus C{" DMS float"} if B{C{cs}} is C{True},
                    where C{[Band]} is present and C{'A'|'B'|'Y'|'Z'} only
                    if B{C{B}} is C{True} and convergence C{DMS} is in
                    I{either} degrees, minutes I{or} seconds (C{str}).

           @note: Zone zero (C{"00"}) for UPS follows Karney's U{zone UPS
                  <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1UTMUPS.html>}.
        '''
        return self._toStr4_6(self.pole, B, cs, prec, sep)  # PYCHOK pole

    def toStr2(self, prec=0, fmt='[%s]', sep=', ', B=False, cs=False):  # PYCHOK expected
        '''Return a string representation of this UPS coordinate.

           Note that UPS coordinates are rounded, not truncated
           (unlike MGRS grid references).

           @keyword prec: Optional number of decimals, unstripped (C{int}).
           @keyword fmt: Optional, enclosing backets format (C{str}).
           @keyword sep: Optional separator between name:value pairs (C{str}).
           @keyword B: Optionally, include polar band letter (C{bool}).
           @keyword cs: Optionally, include gamma meridian convergence
                        and point scale factor (C{bool}).

           @return: This UPS as a string with C{00[Band] pole, easting,
                    northing, [convergence, scale]} as C{"[Z:00[Band],
                    P:N|S, E:meter, N:meter]"} plus C{", C:DMS, S:float"}
                    if B{C{cs}} is C{True}, where C{[Band]} is present and
                    C{'A'|'B'|'Y'|'Z'} only if B{C{B}} is C{True} and
                    convergence C{DMS} is in I{either} degrees, minutes
                    I{or} seconds (C{str}).

           @note: Pseudo zone zero (C{"00"}) for UPS follows Karney's U{zone UPS
                  <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1UTMUPS.html>}.
        '''
        return self._toStr2(prec=prec, fmt=fmt, sep=sep, B=B, cs=cs)

    def toUps(self, pole='', **unused):
        '''Duplicate this UPS coordinate.

           @keyword pole: Optional top/center of the UPS projection,
                          (C{str}, 'N[orth]'|'S[outh]').

           @return: A copt of this UPS coordinate (L{Ups}).

           @raise UPSError: Invalid B{C{pole}} or attempt to transfer
                            the projection top/center.
        '''
        if self.pole == pole or not pole:  # PYCHOK pole
            return self.copy()
        raise UPSError('%s transfer invalid: %r to %r' %
                       ('pole', self.pole, pole))  # PYCHOK pole

    def toUtm(self, zone, falsed=True, **unused):
        '''Convert this UPS coordinate to a UTM coordinate.

           @param zone: The UTM zone (C{int}).
           @keyword falsed: False both easting and northing (C{bool}).

           @return: The UTM coordinate (L{Utm}).
        '''
        u = self._utm
        if u is None or u.zone != zone or falsed != u.falsed:
            from pygeodesy.utm import toUtm8, Utm  # PYCHOK recursive import
            ll = self.toLatLon(LatLon=None, unfalse=True)
            self._utm = toUtm8(ll, Utm=Utm, falsed=falsed, name=self.name, zone=zone)
        return self._utm

    @property_RO
    def zone(self):
        '''Get the polar pseudo zone (C{0}), like Karney's U{zone UPS<https://
           GeographicLib.SourceForge.io/html/classGeographicLib_1_1UTMUPS.html>}.
        '''
        return _UPS_ZONE


class _UpsK1(Ups):
    '''(INTERNAL) For method L{Ups.rescale}.
    '''
    _scale0 = _K1


def parseUPS5(strUPS, datum=Datums.WGS84, Ups=Ups, falsed=True, name=''):
    '''Parse a string representing a UPS coordinate, consisting of
       C{"[zone][band] pole easting northing"} where B{C{zone}} is
       pseudo zone C{"00"|"0"|""} and C{band} is C{'A'|'B'|'Y'|'Z'|''}.

       @param strUPS: A UPS coordinate (C{str}).
       @keyword datum: Optional datum to use (L{Datum}).
       @keyword Ups: Optional (sub-)class to return the UPS
                     coordinate (L{Ups}) or C{None}.
       @keyword falsed: Both B{C{easting}} and B{C{northing}}
                        are falsed (C{bool}).
       @keyword name: Optional B{C{Ups}} name (C{str}).

       @return: The UPS coordinate (B{C{Ups}}) or a
                L{UtmUps5Tuple}C{(zone, hemipole, easting, northing,
                band)} if B{C{Ups}} is C{None}.  The C{hemipole} is
                the C{'N'|'S'} pole, the UPS projection top/center.

       @raise UPSError: Invalid B{C{strUPS}}.
    '''
    try:
        u = strUPS.lstrip()
        if not u.startswith(_UPS_ZONE_STR):
            raise ValueError

        z, p, e, n, B = _parseUTMUPS(u)
        if z != _UPS_ZONE or (B and B not in _Bands):
            raise ValueError
    except (AttributeError, TypeError, ValueError):
        raise UPSError('%s invalid: %r' % ('strUPS', strUPS))

    if Ups is None:
        r = UtmUps5Tuple(z, p, e, n, B)
    else:
        r = Ups(z, p, e, n, band=B, falsed=falsed, datum=datum)
    return _xnamed(r, name)


def toUps8(latlon, lon=None, datum=None, Ups=Ups, pole='',
                             falsed=True, strict=True, name=''):
    '''Convert a lat-/longitude point to a UPS coordinate.

       @param latlon: Latitude (C{degrees}) or an (ellipsoidal)
                      geodetic C{LatLon} point.
       @keyword lon: Optional longitude (C{degrees}) or C{None}
                     if B{C{latlon}} is a C{LatLon}.
       @keyword datum: Optional datum for this UPS coordinate,
                       overriding B{C{latlon}}'s datum (C{Datum}).
       @keyword Ups: Optional (sub-)class to return the UPS
                     coordinate (L{Ups}) or C{None}.
       @keyword pole: Optional top/center of (stereographic) projection
                      (C{str}, C{'N[orth]'} or C{'S[outh]'}).
       @keyword falsed: False both easting and northing (C{bool}).
       @keyword strict: Restrict B{C{lat}} to UPS ranges (C{bool}).
       @keyword name: Optional B{C{Ups}} name (C{str}).

       @return: The UPS coordinate (B{C{Ups}}) or a
                L{UtmUps8Tuple}C{(zone, hemipole, easting, northing,
                band, datum, convergence, scale)} if B{C{Ups}} is
                C{None}.  The C{hemipole} is the C{'N'|'S'} pole,
                the UPS projection top/center.

       @raise RangeError: If B{C{strict}} and B{C{lat}} outside the
                          valid UPS bands or if B{C{lat}} or B{C{lon}}
                          outside the valid range and L{rangerrors}
                          set to C{True}.

       @raise TypeError: If B{C{latlon}} is not ellipsoidal.

       @raise ValueError: If B{C{lon}} value is missing or if B{C{latlon}}
                          is invalid.

       @see: Karney's C++ class U{UPS
             <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1UPS.html>}.
    '''
    lat, lon, d, name = _to4lldn(latlon, lon, datum, name)
    z, B, p, lat, lon = upsZoneBand5(lat, lon, strict=strict)  # PYCHOK UtmUpsLatLon5Tuple

    E = d.ellipsoid

    p = str(pole or p)[:1].upper()
    N = p == 'N'  # is north

    a = lat if N else -lat
    A = abs(a - 90) < _TOL  # at pole

    t = tan(radians(a))
    T = E.es_taupf(t)
    r = hypot1(T) + abs(T)
    if T >= 0:
        r = 0 if A else 1 / r

    k0 = getattr(Ups, '_scale0', _K0)  # Ups is class or None
    r *= 2 * k0 * E.a / E.es_c

    k = k0 if A else _scale(E, r, t)
    c = lon  # [-180, 180) from .upsZoneBand5
    x, y = sincos2d(c)
    x *= r
    y *= r
    if N:
        y = -y
    else:
        c = -c

    if falsed:
        x += _Falsing
        y += _Falsing

    if Ups is None:
        r = UtmUps8Tuple(z, p, x, y, B, d, c, k)
    else:
        if z != _UPS_ZONE and not strict:
            z = _UPS_ZONE  # ignore UTM zone
        r = Ups(z, p, x, y, band=B, datum=d, falsed=falsed,
                                    convergence=c, scale=k)
        r._hemisphere = _hemi(lat)
        if isinstance(latlon, _LLEB) and d is latlon.datum:
            r._latlon_to(latlon, falsed)  # XXX weakref(latlon)?
    return _xnamed(r, name)


def upsZoneBand5(lat, lon, strict=True):
    '''Return the UTM/UPS zone number, (polar) Band letter, pole and
       clipped lat- and longitude for a given location.

       @param lat: Latitude in degrees (C{scalar} or C{str}).
       @param lon: Longitude in degrees (C{scalar} or C{str}).
       @keyword strict: Restrict B{C{lat}} to UPS ranges (C{bool}).

       @return: A L{UtmUpsLatLon5Tuple}C{(zone, band, hemipole,
                lat, lon)} where C{hemipole} is the C{'N'|'S'} pole,
                the UPS projection top/center.

       @raise RangeError: If B{C{strict}} and B{C{lat}} in the UTM
                          and not the UPS range or if B{C{lat}} or
                          B{C{lon}} outside the valid range and
                          L{rangerrors} set to C{True}.

       @raise ValueError: Invalid B{C{lat}} or B{C{lon}}.
    '''
    z, lat, lon = _to3zll(*parseDMS2(lat, lon))
    if lat < _UPS_LAT_MIN:  # includes 30' overlap
        z, B, p = _UPS_ZONE, _Band(lat, lon), 'S'

    elif lat > _UPS_LAT_MAX:  # includes 30' overlap
        z, B, p = _UPS_ZONE, _Band(lat, lon), 'N'

    elif strict:
        x = '%s [%s, %s]' % ('range', _UPS_LAT_MIN, _UPS_LAT_MAX)
        raise RangeError('%s inside UTM %s: %s' % ('lat', x, degDMS(lat)))

    else:
        B, p = '', _hemi(lat)
    return UtmUpsLatLon5Tuple(z, B, p, lat, lon)

# **) MIT License
#
# Copyright (C) 2016-2020 -- mrJean1 at Gmail -- All Rights Reserved.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included
# in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
# OTHER DEALINGS IN THE SOFTWARE.
