content: {-# LANGUAGE BangPatterns #-}
  {-# LANGUAGE MagicHash    #-}
  {-# LANGUAGE PatternSynonyms   #-}
  {-|
  Module      : Reanimate.GeoProjection
  Copyright   : Written by David Himmelstrup
  License     : Unlicense
  Maintainer  : lemmih@gmail.com
  Stability   : experimental
  Portability : POSIX
  
  This module provides functions for mapping the surface of
  a sphere on to a 2D plane. It also has convenience functions
  for loading GeoJSON data.
  
  -}
  module Reanimate.GeoProjection
    ( Projection(..)
    , XYCoord(..)
    , LonLat(..)
    , project
    , interpP
    , interpBBP
    , mergeP
    , isValidP
    , scaleP
    , flipYAxisP
    , moveBottomP
    , moveTopP
      -- * Projections
    , equirectangularP
    , mercatorP
    , mollweideP
    , hammerP
    , cylindricalEqualAreaP
    , lambertP
    , bottomleyP
    , sinusoidalP
    , wernerP
    , bonneP
    , orthoP
    , cassiniP
    , augustP
    , collignonP
    , eckert1P
    , eckert3P
    , eckert5P
    , faheyP
    , foucautP
    , lagrangeP
      -- * GeoJSON helpers
    , drawFeatureCollection -- :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
    , loadFeatureCollection -- :: FromJSON a => FilePath -> (a -> SVG -> SVG) -> SVG
    , applyProjection -- :: Projection -> SVG -> SVG
    , applyProjection' -- :: Double -> Projection -> SVG -> SVG
    , renderGeometry
    ) where
  
  import           Codec.Picture
  import           Codec.Picture.Types
  import           Control.Lens            ((^.))
  import           Control.Monad
  import           Control.Monad.ST
  import           Control.Monad.ST.Unsafe
  import           Data.Aeson
  import           Data.Foldable
  import           Data.Geospatial         hiding (LonLat)
  import           Data.Hashable
  import           Data.LinearRing
  import qualified Data.LineString         as Line
  import           Data.Vector.Storable    (unsafeWith)
  import qualified Data.Vector.Unboxed     as V
  import           Debug.Trace
  import           Foreign
  import           GHC.Exts                (Double (..), cosDouble#, sinDouble#,
                                            (*##), (+##), (-##), (/##))
  import           Graphics.SvgTree        (pattern None)
  import           Linear                  (distance, lerp)
  import           Linear.V2               hiding (angle)
  import           Reanimate
  import           System.IO.Unsafe
  
  {-# INLINE halfPi #-}
  {-# INLINE sqrtPi #-}
  {-# INLINE sqrt2 #-}
  {-# INLINE epsilon #-}
  {-# INLINE tau #-}
  -- Constants
  halfPi, sqrtPi, sqrt2, epsilon, tau :: Double
  halfPi = pi/2
  sqrtPi = sqrt pi
  sqrt2 = sqrt 2
  epsilon = 1.0e-12
  tau = pi*2
  
  toRads, cot :: Double -> Double
  toRads dec = dec/180 * pi
  cot = recip . tan
  
  
  srcPixel :: Pixel pixel => Image pixel -> LonLat -> pixel
  srcPixel src (LonLat lam phi) =
      -- pixelAt src xPx yPx
      unsafePixelAt (imageData src) (pixelBaseIndex src xPx yPx)
    where
      !xPx = round $ ((lam+pi)/tau) * fromIntegral (imageWidth src-1)
      !yPx = round $ (1-((phi+halfPi)/pi)) * fromIntegral (imageHeight src-1)
  
  
  srcPixelFast :: Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
  srcPixelFast src wMax hMax (LonLat lam phi) = unsafeIOToST $
      unsafeWith (imageData src) $ \ptr -> do
        let ptr' = plusPtr ptr idx
        r <- peek ptr'
        g <- peek $ plusPtr ptr' 1
        b <- peek $ plusPtr ptr' 2
        a <- peek $ plusPtr ptr' 3
        return $ PixelRGBA8 r g b a
    where
      !idx = pixelBaseIndex src xPx yPx
      !xPx = round $ ((lam+pi)/tau) * wMax
      !yPx = round $ (1-((phi+halfPi)/pi)) * hMax
  
  -- | Interpolate between two projections and apply the result to an image in
  --   equirectangular format. The source image must have an aspect ratio of 2:1.
  interpP :: Image PixelRGBA8 -> Projection -> Projection -> Double -> Image PixelRGBA8
  interpP src p1 _ 0 = project src p1
  interpP src _ p2 1 = project src p2
  interpP !src (Projection _label1 p1 p1_inv) (Projection _label2 p2 p2_inv) !t = runST $ do
      !img <- newMutableImage w h
  
      let factor = 2
      let l1 =
            loopTo (w*factor) $ \x ->
              loopTo (h*factor) $ \y -> do
                let !x1' = fromIntegral x / (wMax*fromIntegral factor)
                    !y1' = fromIntegral y / (hMax*fromIntegral factor)
                    !lonlat = p1_inv $! XYCoord x1' y1'
  
                when (validLonLat lonlat) $ do
                  p <- srcPixelFast src wMax hMax lonlat
                  when (pixelOpacity p /= 0) $ do
                    let XYCoord x1 y1 = p1 lonlat
                        XYCoord x2 y2 = p2 lonlat
                        !x3 = round $ fromToS x1 x2 t * wMax
                        !y3 = round $ (1 - fromToS y1 y2 t) * hMax
                    when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
                      writePixel img x3 y3 p
          l2 =
            loopTo (w*factor) $ \x ->
              loopTo (h*factor) $ \y -> do
                let !x2' = fromIntegral x / (wMax*fromIntegral factor)
                    !y2' = fromIntegral y / (hMax*fromIntegral factor)
                    !lonlat = p2_inv (XYCoord x2' y2')
                when (validLonLat lonlat) $ do
                  p <- srcPixelFast src wMax hMax lonlat
                  when (pixelOpacity p /= 0) $ do
                    let XYCoord x2 y2 = p2 lonlat
                        XYCoord x1 y1 = p1 lonlat
                        !x3 = round $ fromToS x1 x2 t * wMax
                        !y3 = round $ (1 - fromToS y1 y2 t) * hMax
                    when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
                      writePixel img x3 y3 p
      if t < 0.5
        then l1 >> l2
        else l2 >> l1
      unsafeFreezeImage img
    where
      loopTo m fn = go m
        where go 0 = return ()
              go n = fn (n-1) >> go (n-1)
      !w = imageWidth src
      !h = imageHeight src
      !wMax = fromIntegral (w-1)
      !hMax = fromIntegral (h-1)
  
  -- | Interpolate between two projections and apply the result to an image in
  --   equirectangular format. The source image must have an aspect ratio of 2:1.
  --   Only the areas inside of the two bounding boxes (applying to the source and
  --   target projection, respectively) are mapped. Pixels outside of these bounding-boxes
  --   are undefined.
  interpBBP :: Image PixelRGBA8 -> Projection -> Projection ->
              (Double,Double,Double,Double) -> (Double,Double,Double,Double) -> Double -> Image PixelRGBA8
  interpBBP !src (Projection _ p1 p1_inv) (Projection _ p2 p2_inv) (fx,fy,fw,fh) (tx, ty, tw, th) !t = runST $ do
      !img <- newMutableImage w h
      let factor = 2
      let l1 =
            loopTo (w*factor) $ \x ->
              loopTo (h*factor) $ \y -> do
                let !x1' = fromIntegral x / (wMax*fromIntegral factor)
                    !y1' = fromIntegral y / (hMax*fromIntegral factor)
                when (x1' >= fx && x1' <= fx+fw && y1' >= fy && y1' <= fy+fh) $ do
                  let !lonlat = p1_inv $! XYCoord x1' y1'
  
                  when (validLonLat lonlat) $ do
                      p <- srcPixelFast src wMax hMax lonlat
                      when (pixelOpacity p /= 0) $ do
                        let XYCoord x1 y1 = p1 lonlat
                            XYCoord x2 y2 = p2 lonlat
                            !x3 = round $ fromToS x1 x2 t * wMax
                            !y3 = round $ (1 - fromToS y1 y2 t) * hMax
                        when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
                          writePixel img x3 y3 p
          l2 =
            loopTo (w*factor) $ \x ->
              loopTo (h*factor) $ \y -> do
                let !x2' = fromIntegral x / (wMax*fromIntegral factor)
                    !y2' = fromIntegral y / (hMax*fromIntegral factor)
                when (x2' >= tx && x2' <= tx+tw && y2' >= ty && y2' <= ty+th) $ do
                  let !lonlat = p2_inv (XYCoord x2' y2')
                  when (validLonLat lonlat) $ do
                    p <- srcPixelFast src wMax hMax lonlat
                    when (pixelOpacity p /= 0) $ do
                      let XYCoord x2 y2 = p2 lonlat
                          XYCoord x1 y1 = p1 lonlat
                          !x3 = round $ fromToS x1 x2 t * wMax
                          !y3 = round $ (1 - fromToS y1 y2 t) * hMax
                      when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
                        writePixel img x3 y3 p
      if t < 0.5
        then l1 >> l2
        else l2 >> l1
  
      unsafeFreezeImage img
    where
      loopTo m fn = go m
        where go 0 = return ()
              go n = fn (n-1) >> go (n-1)
      !w = imageWidth src
      !h = imageHeight src
      !wMax = fromIntegral (w-1)
      !hMax = fromIntegral (h-1)
  
  
  eqLonLat :: LonLat -> LonLat -> Bool
  eqLonLat (LonLat x1 y1) (LonLat x2 y2)
    = eqDouble x1 x2 && eqDouble y1 y2
  {-
  eqCoords :: XYCoord -> XYCoord -> Bool
  eqCoords (XYCoord x1 y1) (XYCoord x2 y2)
    = eqDouble x1 x2 && eqDouble y1 y2
  -}
  eqDouble :: Double -> Double -> Bool
  eqDouble a b = abs (a-b) < epsilon
  
  -- | XY coordinates on a 2D plane. Valid ranges go from 0 to 1, inclusive.
  data XYCoord = XYCoord !Double !Double -- 0 to 1
    deriving (Read,Show,Eq,Ord)
  
  -- | Longitude and latitude. Valid range for longitude is -pi to +pi.
  --   Valid range for latitude is -pi/2 to +pi/2.
  data LonLat = LonLat !Double !Double -- -pi to +pi, -halfPi to +halfPi
    deriving (Read,Show,Eq,Ord)
  instance Hashable LonLat where
    hashWithSalt s (LonLat a b) = hashWithSalt s (a,b)
  
  -- | Projections are named bi-directional mappings between a sphere
  --   and a 2D plane.
  data Projection = Projection
    { projectionLabel   :: String -- ^ Name of the projection.
    , projectionForward :: !(LonLat -> XYCoord)
      -- ^ Mapping from longitude and latitude on a sphere to
      --   XY coordinates on a 2D plane.
    , projectionInverse :: !(XYCoord -> LonLat)
      -- ^ Mapping from XY coordinates on a 2D plane to longitude
      --   and latitude on a sphere.
    }
  
  -- FIXME: Verify that 'src' has an aspect ratio of 2:1.
  -- | Apply on an image in equirectangular format. The source image
  --   therfore must have an aspect ratio of 2:1.
  project :: Image PixelRGBA8 -> Projection -> Image PixelRGBA8
  project src (Projection _label _ pInv) = generateImage fn w h
    where
      w = imageWidth src
      h = imageHeight src
      fn xPx yPx =
        let x = (fromIntegral xPx / fromIntegral (w-1))
            y = 1-(fromIntegral yPx / fromIntegral (h-1))
            lonlat = pInv (XYCoord x y)
        in
        if validLonLat lonlat
          then srcPixel src lonlat
          else PixelRGBA8 0 0 0 0
  
  validLonLat :: LonLat -> Bool
  validLonLat (LonLat lam phi) =
    lam >= -pi && lam <= pi && phi >= -pi/2 && phi <= pi/2
  
  _validXYCoord :: XYCoord -> Bool
  _validXYCoord (XYCoord x y) = x >= 0 && x <= 1 && y >= 0 && y <= 1
  
  -- | Returns @True@ iff a projection is consistent and complete.
  isValidP :: Projection -> Bool
  isValidP (Projection _label p pInv) = and
      [ check x y
      | x <- [0..w-1::Int]
      , y <- [0..h-1::Int] ]
    where
      w = 100
      h = 100
      check xPx yPx =
        let x = (fromIntegral xPx / fromIntegral (w-1))
            y = (fromIntegral yPx / fromIntegral (h-1))
            lonlat = pInv (XYCoord x y)
            lonlat2 = pInv $ p lonlat
        in not (validLonLat lonlat) || eqLonLat lonlat lonlat2
            || trace (show (lonlat, lonlat2)) False
  
  -- | Translate the lower-most point of a projection by an offset.
  moveBottomP :: Double -> Projection -> Projection
  moveBottomP offset (Projection label p pInv) = Projection label p' pInv'
    where
      p' (LonLat lon lat) =
        case p (LonLat lon lat) of
          XYCoord x y -> XYCoord x (fromToS offset 1 y)
      pInv' (XYCoord x y) = pInv (XYCoord x ((y-offset)/(1-offset)))
  
  -- | Translate the top-most point of a projection by an offset.
  moveTopP :: Double -> Projection -> Projection
  moveTopP offset = flipYAxisP . moveBottomP offset . flipYAxisP
  
  -- | Invert the Y axis of projection.
  flipYAxisP :: Projection -> Projection
  flipYAxisP (Projection label p pInv) = Projection label p' pInv'
    where
      p' (LonLat lam phi) =
        let XYCoord x y = p (LonLat lam (negate phi))
        in XYCoord x (1-y)
      pInv' (XYCoord x y) =
        let LonLat lam phi = pInv (XYCoord x (1-y))
        in LonLat lam (negate phi)
  
  -- | Scale X and Y axis of projection.
  scaleP :: Double -> Double -> Projection -> Projection
  scaleP xScale yScale (Projection label p pInv) = Projection label forward inverse
    where
      forward lonlat =
        case p lonlat of
          XYCoord x y -> XYCoord ((x-0.5)*xScale+0.5) ((y-0.5)*yScale+0.5)
      inverse (XYCoord x y) =
        pInv $ XYCoord ((x-0.5)/xScale+0.5) ((y-0.5)/yScale+0.5)
  
  -- | Attempt to smoothly interpolate two projections. The result may not be continuous
  --   and 'interpP' may give prettier results.
  mergeP :: Projection -> Projection -> Double -> Projection
  mergeP p1 p2 t = Projection (projectionLabel p1 ++ "/" ++ projectionLabel p2) p pInv
    where
      p lonlat =
        let XYCoord x1 y1 = projectionForward p1 lonlat
            XYCoord x2 y2 = projectionForward p2 lonlat
        in XYCoord (fromToS x1 x2 t) (fromToS y1 y2 t)
      pInv coord =
        let LonLat lon1 lat1 = projectionInverse p1 coord
            LonLat lon2 lat2 = projectionInverse p2 coord
        in
          if oob lon1 lat1 && oob lon2 lat2
            then LonLat (0/0) (0/0)
            else LonLat (fromToS lon1 lon2 t) (fromToS lat1 lat2 t)
      oob lon lat = lon < (-pi) || lon > pi || lat < (-pi/2) || lat > pi/2
  
  -- | <<docs/gifs/doc_equirectangularP.gif>>
  equirectangularP :: Projection
  equirectangularP = Projection "equirectangular" forward inverse
    where
      forward (LonLat lam phi) = XYCoord ((lam+pi)/tau) ((phi+pi/2)/pi)
      inverse (XYCoord x y) = LonLat xPi yPi
        where
          xPi = fromToS (-pi) pi x
          yPi = fromToS (-pi/2) (pi/2) y
  
  -- | <<docs/gifs/doc_mercatorP.gif>>
  mercatorP :: Projection
  mercatorP = Projection "mercator" forward inverse
    where
      forward (LonLat lam phi) =
        XYCoord ((lam+pi)/tau)
          (min 1 $ max 0 $ (log(tan(pi/4+phi/2)) + pi)/tau)
      inverse (XYCoord x y) = LonLat xPi (atan (sinh yPi))
        where
          xPi = fromToS (-pi) pi x
          yPi = fromToS (-pi) pi y
  
  thetas :: V.Vector Double
  thetas = V.fromList $
    map (findThetaFast . fromIndex) [0 .. granularity]
  
  granularity :: Int
  granularity = 50000
  
  toIndex :: Double -> Int
  toIndex phi = round ((phi+halfPi)/pi * fromIntegral granularity)
  fromIndex :: Int -> Double
  fromIndex x = fromToS (-halfPi) halfPi (fromIntegral x / fromIntegral granularity)
  granualize :: Double -> Double
  granualize = fromIndex . toIndex
  
  {-# INLINE mollweideP #-}
  -- | <<docs/gifs/doc_mollweideP.gif>>
  mollweideP :: Projection
  mollweideP = Projection "mollweide" forward inverse
    where
      forward (LonLat !lam !phi) =
          XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
        where
          x = (2*sqrt2)/pi * lam * cos theta
          y = sqrt2*sin theta
          theta = thetas V.! toIndex phi
          -- find_theta :: Int -> Double
          -- find_theta 0 = phi
          -- find_theta _ | abs phi == pi/2 = signum phi * pi/2
          -- find_theta n =
          --   let !sub = find_theta (n-1)
          --   in sub - (2*sub+sin (2*sub)-pi*sin phi)/(2+2*cos(2*sub))
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-2*sqrt2) (2*sqrt2) x'
          y = fromToS (-1) 1 y'
          theta = granualize (asin y)
          lam = pi*x/(2*sqrt2*cos theta)
          phi = asin ((2*theta+sin(2*theta))/pi)
  
  findThetaFast :: Double -> Double
  findThetaFast !phi | abs phi == pi/2 = signum phi * halfPi
  findThetaFast (D# phi) = go phi
    where
      !(D# pi#) = pi
      go acc =
        let c = cosDouble# (acc +## acc)
            s = sinDouble# (acc +## acc)
            next =
              acc -##
              (acc +## acc +## s -## pi# *## sinDouble# phi)
              /## (2.0## +## c +## c) in
        if abs (D# (next -## acc)) < epsilon
          then D# next
          else go next
  
  
  -- | <<docs/gifs/doc_hammerP.gif>>
  hammerP :: Projection
  hammerP = Projection "hammer" forward inverse
    where
      forward (LonLat lam phi) =
          XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
        where
          x = (2*sqrt2*cos phi*sin (lam/2))/sqrt (1+cos phi*cos (lam/2))
          y = (sqrt2*sin phi)/ sqrt (1+cos phi*cos (lam/2))
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-2*sqrt2) (2*sqrt2) x'
          y = fromToS (-sqrt2) sqrt2 y'
          z = sqrt (1 - (x/4)**2 - (y/2)**2)
          lam = 2 * atan2 (z*x) (2*(2*z**2-1))
          phi = asin (z*y)
  
  -- | <<docs/gifs/doc_cylindricalEqualAreaP.gif>>
  cylindricalEqualAreaP :: Double -> Projection
  cylindricalEqualAreaP phi0 = Projection "lambert" forward inverse
    where
      cosPhi0 = cos phi0
      forward (LonLat lam phi) =
        XYCoord ((lam+pi)/tau) ((sin phi/cosPhi0+1/cosPhi0)/2/cosPhi0)
      inverse (XYCoord x' y') = LonLat x (asin y / (asin cosPhi0 / halfPi))
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-1) 1 y' * cosPhi0
  
  -- | <<docs/gifs/doc_lambertP.gif>>
  lambertP :: Projection
  lambertP = Projection "lambert" forward inverse
    where
      forward (LonLat lam phi) =
        XYCoord ((lam+pi)/tau) ((sin phi+1)/2)
      inverse (XYCoord x' y') = LonLat x (asin y)
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-1) 1 y'
  
  -- | <<docs/gifs/doc_bottomleyP.gif>>
  bottomleyP :: Double -> Projection
  bottomleyP !phi_1 = Projection "bottomley" forward inverse
    where
      forward (LonLat lam phi) =
          XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
        where
          x = (rho * sin e) / sin phi_1
          y = pi/2 - rho * cos e
          rho = pi/2 - phi
          e = lam * sin phi_1 * sin rho / rho
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-pi/2) (pi/2) y'
          x1 = x * sin phi_1
          y1 = pi/2 - y
          rho = sqrt (x1*x1 + y1*y1)
          e = atan2 x1 y1
          lam = (if rho == 0 then 1 else rho / sin rho) * e/sin phi_1
          phi = pi/2 - rho
  
  -- | <<docs/gifs/doc_sinusoidalP.gif>>
  sinusoidalP :: Projection
  sinusoidalP = Projection "sinusoidal" forward inverse
    where
      forward (LonLat lam phi) =
          XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
        where
          x = lam * cos phi
          y = phi
      inverse (XYCoord x' y') = LonLat (x/cos y) y
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-pi/2) (pi/2) y'
  
  -- | <<docs/gifs/doc_wernerP.gif>>
  wernerP :: Projection
  wernerP = moveTopP 0.23 $ Projection "werner" forward inverse
    where
      forward (LonLat lam phi) =
          XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
        where
          rho = pi/2 - phi
          e = if rho == 0 then rho else lam * sin rho / rho
          x = rho * sin e
          y = pi/2 - rho * cos e
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-pi/2) (pi/2) y'
          rho = sqrt (x**2 + (pi/2 - y)**2)
          e = atan2 x (pi/2 -y)
          phi = pi/2 - rho
          lam = e * rho / sin rho
  
  -- FIXME: find right scale and position.
  -- | <<docs/gifs/doc_bonneP.gif>>
  bonneP :: Double -> Projection
  bonneP 0 = sinusoidalP
  bonneP phi_0 = moveTopP (-0.17*factor) $ scaleP 1 (fromToS 1 0.65 factor) $
      Projection "bonne" forward inverse
    where
      factor = sin phi_0 / sin (pi/4)
      forward (LonLat lam phi ) = XYCoord ((x+pi)/tau) ((y+halfPi)/pi)
        where
          cotPhi0 = cot phi_0
          rho = cotPhi0 + phi_0 - phi
          e = (lam * cos phi) / rho
          x = rho * sin e
          y = cotPhi0 - rho * cos e
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-pi) pi x'
          y = fromToS (-pi/2) (pi/2) y'
          cotPhi0 = cot phi_0
          rho = sqrt (x*x + (cot phi_0 - y)**2)
          -- e = atan2 x (cot phi_1 + phi_1 + y)
          phi = cotPhi0 + phi_0 - rho
          lam = rho / cos phi * atan2 x (cotPhi0-y)
  
  -- | <<docs/gifs/doc_orthoP.gif>>
  orthoP :: LonLat -> Projection
  orthoP (LonLat lam_0 phi_0) = Projection "ortho" forward inverse
    where
      forward (LonLat lam phi)
          | cosc < 0  =
              let ang = atan2 y x
                  xV = cos ang
                  yV = sin ang
                  xPos = 7/32 + ((xV+1)/2 * 9/16)
              in -- trace (show (x,y)) $
                XYCoord xPos ((yV+1)/2)
              --XYCoord (0/0) (0/0)
          | otherwise = XYCoord ((x+(16/9))/(16/9*2)) ((y+1)/2)
        where
          x = cos phi * sin (lam - lam_0)
          y = cos phi_0 * sin phi - sin phi_0 * cos phi * cos (lam - lam_0)
          cosc = sin phi_0 * sin phi + cos phi_0 * cos phi * cos (lam-lam_0)
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-16/9) (16/9) x'
          y = fromToS (-1) 1 y'
          lam = wrap (-pi) pi $
            lam_0 + atan2 (x * sin c) (rho * cos c * cos phi_0 - y * sin c * sin phi_0)
          phi = wrap (-pi/2) (pi/2) $
            asin (cos c * sin phi_0 + (y * sin c * cos phi_0)/rho)
          rho = sqrt (x*x + y*y)
          c = asin rho
          wrap lower upper v
            | v > upper = v-upper+lower
            | v < lower = v+upper-lower
            | otherwise = v
  
  -- | <<docs/gifs/doc_cassiniP.gif>>
  cassiniP :: Projection
  cassiniP = Projection "cassini" forward inverse
    where
      forward (LonLat lam phi) =
        XYCoord ((asin (cos phi * sin lam)+halfPi)/pi) ((atan2 (tan phi) (cos lam)+pi)/tau)
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS (-halfPi) halfPi x'
          y = fromToS (-pi) pi y'
          lam = atan2 (tan x) (cos y)
          phi = asin (sin y * cos x)
  
  -- | <<docs/gifs/doc_augustP.gif>>
  augustP :: Projection
  augustP = scaleP 0.70 0.70 $ Projection "august"  forward inverse
    where
      xHi = 16/3
      xLo = -xHi
      yHi = 8 / 3
      yLo = -yHi
      forward (LonLat lam phi) = XYCoord ((xPos-xLo)/(xHi-xLo)) ((yPos-yLo)/(yHi-yLo))
        where
          tanPhi = tan (phi/2)
          k = sqrt (1 - tanPhi * tanPhi)
          c = 1 + k * cos (lam / 2)
          x = sin (lam/2) * k / c
          y = tanPhi / c
          x2 = x*x
          y2 = y*y
          xPos = 4 / 3 * x * (3+x2 - 3*y2)
          yPos = 4 / 3 * y * (3 + 3*x2 - y2)
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x' * 3 / 8
          y = fromToS yLo yHi y' * 3 / 8
          x2 = x*x
          y2 = y*y
          s = 1 + x2 + y2
          sin3Eta = sqrt ((s - sqrt (s*s - 4 * y * y)) / 2)
          eta = asin sin3Eta / 3
          xi = if sin3Eta /= 0 then acosh (abs (y / sin3Eta)) / 3 else asinh (abs x) / 3
          cosEta = cos eta
          coshXi = cosh xi
          d = coshXi * coshXi - cosEta * cosEta
          lam = signum x * 2 * atan2 (sinh xi * cosEta) (0.25 - d)
          phi = signum y * 2 * atan2 (coshXi * sin eta) (0.25 + d)
  
  -- | <<docs/gifs/doc_collignonP.gif>>
  collignonP :: Projection
  collignonP = Projection "collignon" forward inverse
    where
      yHi = sqrtPi
      yLo = sqrtPi * (1 - sqrt2)
      xLo = -pi*(2/sqrtPi)*sqrt2
      xHi = pi*(2/sqrtPi)*sqrt2
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          alpha = sqrt (1 - sin phi)
          x = (2 / sqrtPi) * lam * alpha
          y = sqrtPi * (1 - alpha)
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          l = (y / sqrtPi -1)**2
          lam = if l > 0 then x * sqrt (pi / l) / 2 else 0
          phi = asin (1-l)
  
  -- | <<docs/gifs/doc_eckert1P.gif>>
  eckert1P :: Projection
  eckert1P = Projection "eckert1" forward inverse
    where
      alpha = sqrt (8 / (3*pi))
      yLo = -alpha * halfPi
      yHi = alpha * halfPi
      xLo = -alpha * pi
      xHi = alpha * pi
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          x = alpha * lam * (1 - abs phi / pi)
          y = alpha * phi
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          phi = y / alpha
          lam = x / (alpha * (1 - abs phi / pi))
  
  -- | <<docs/gifs/doc_eckert3P.gif>>
  eckert3P :: Projection
  eckert3P = Projection "eckert3" forward inverse
    where
      k = sqrt (pi * (4 + pi))
      yLo = negate yHi
      yHi = 4/k * halfPi
      xLo = negate xHi
      xHi = 4/k * pi
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          x = 2 / k * lam * (1 + sqrt (1 - 4*phi*phi/(pi*pi)))
          y = 4 / k * phi
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          lam = x  * k/2 / (1 + sqrt (1 - y*y* (4+pi)/(4*pi)))
          phi = y * k/4
  
  -- | <<docs/gifs/doc_eckert5P.gif>>
  eckert5P :: Projection
  eckert5P = Projection "eckert5" forward inverse
    where
      k = sqrt (2 + pi)
      yLo = negate yHi
      yHi = pi/k
      xLo = negate xHi
      xHi = tau/k
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          x = lam * (1 + cos phi) / k
          y = 2 * phi / k
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          lam = k * x / (1 + cos phi)
          phi = y * k / 2
  
  {-# INLINE faheyP #-}
  -- | <<docs/gifs/doc_faheyP.gif>>
  faheyP :: Projection
  faheyP = Projection "fahey" forward inverse
    where
      faheyK = cos (toRads 35)
      yLo = negate yHi
      yHi = 1 + faheyK
      xLo = negate xHi
      xHi = pi * faheyK * 16/9
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          t = tan (phi/2)
          x = lam * faheyK * sqrt (1 - t*t)
          y = (1 + faheyK) * t
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          t = y / (1 + faheyK)
          lam = x / (faheyK * sqrt (1 - t*t))
          phi = 2 * atan2 y (1 + faheyK)
  
  {-# INLINE foucautP #-}
  -- | <<docs/gifs/doc_foucautP.gif>>
  foucautP :: Projection
  foucautP = Projection "foucaut" forward inverse
    where
      yLo = negate yHi
      yHi = sqrtPi * tan (halfPi/2)
      xLo = negate xHi
      xHi = tau/sqrtPi
      forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          k = phi / 2
          cosk = cos k
          x = 2 * lam / sqrtPi * cos phi * cosk * cosk
          y = sqrtPi * tan k
      inverse (XYCoord x' y') = LonLat lam phi
        where
          x = fromToS xLo xHi x'
          y = fromToS yLo yHi y'
          k = atan (y / sqrtPi)
          cosk = cos k
          phi = 2 * k
          lam = x * sqrtPi / 2 / (cos phi * cosk * cosk)
  
  {-# INLINE lagrangeP #-}
  -- | <<docs/gifs/doc_lagrangeP.gif>>
  lagrangeP :: Projection
  lagrangeP = Projection "lagrange" forward inverse
    where
      yLo = negate yHi
      yHi = 2
      xLo = negate xHi
      xHi = 2
      n = 0.5
      forward (LonLat lam phi)
          | abs (abs phi - halfPi) < epsilon = XYCoord 0.5 (if phi < 0 then 0 else 1)
          | otherwise = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
        where
          sinPhi = sin phi
          v = ((1+sinPhi) / (1 - sinPhi))**(n/2)
          c = 0.5 * (v + 1/v) + cos (lam*n)
          x = 2 * sin (lam*n) / c
          y = (v - 1/v) /c
      inverse (XYCoord x' y')
          | abs (y'-1) < epsilon
            -- = LonLat 0 (signum y * halfPi)
            = LonLat 100 100
          | otherwise = LonLat lam phi
        where
          x = fromToS xLo xHi x' / 2
          y = fromToS yLo yHi y' / 2
          x2 = x * x
          y2 = y * y
          t = 2 * y / (1 + x2 + y2)
          t' = ((1+t) / (1-t)) ** (1/n)
          lam = atan2 (2*x) (1-x2-y2) / n
          phi = asin ((t'-1)/(t'+1))
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  -- | Map for all features and render the geometry.
  drawFeatureCollection :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
  drawFeatureCollection geo fn = mkGroup
    [ fn (feature ^. properties) $ renderGeometry (feature ^. geometry)
    | feature <- toList (geo ^. geofeatures)
    ]
  
  {-# INLINE loadFeatureCollection #-}
  -- | Load GeoJSON from a filepath and render the geometry.
  loadFeatureCollection :: FromJSON a => FilePath -> (a -> SVG -> SVG) -> SVG
  loadFeatureCollection path = unsafePerformIO $ do
    mbGeo <- decodeFileStrict path
    case mbGeo of
      Nothing  -> error $ "loadFeatureCollection: Invalid GeoJSON: " ++ path
      Just geo -> return (drawFeatureCollection geo)
  
  -- drawFeatureCollection :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
  -- loadFeatureColection :: FromJSON a => FilePath -> (a -> SVG -> SVG) -> SVG
  -- modifyPoints :: ((Double,Double) -> (Double, Double)) -> SVG -> SVG
  -- pointsToRadians :: SVG -> SVG
  -- applyProjection :: Projection -> SVG -> SVG
  
  -- | Render GeoJSON geometry as SVG.
  renderGeometry :: GeospatialGeometry -> SVG
  renderGeometry shape =
    case shape of
      MultiPolygon mpolygon ->
        mkGroup $ map (renderGeometry . Polygon) $ toList (splitGeoMultiPolygon mpolygon)
      Polygon poly ->
        mkGroup
        [ mkLinePath section
        | section <- pure
            [ (x, y)
            | PointXY x y <- map retrieveXY (fromLinearRing (head (toList (poly^.unGeoPolygon))))
            ]
        ]
      Line line ->
        mkLinePath
        [ (x, y)
        | PointXY x y <- map retrieveXY (Line.fromLineString (line ^. unGeoLine))
        ]
      MultiLine ml ->
        mkGroup $ map (renderGeometry . Line) $ toList (splitGeoMultiLine ml)
      _ -> None
  
  
  -- | Apply a projection to an SVG image. This is a lossy transformation
  --   but the default tolerance is low enough that inaccuracies should not
  --   be visible.
  applyProjection :: Projection -> SVG -> SVG
  applyProjection = applyProjection' 1e-2
  
  -- | Apply a projection to an SVG image with a specified tolerance.
  --   Projections may turn straight lines into disjointed curves and
  --   the tolerance argument determined the accuracy of this transformation.
  applyProjection' :: Double -> Projection -> SVG -> SVG
  applyProjection' tolerance p = mapSvgLines start
    where
      start (LineMove x:rest) = LineMove (proj x) : worker x rest
      start _                 = []
      worker a (LineEnd b : rest) =
        let (x:xs) = reverse $ drop 1 $ mkChunks a b
        in map (\v -> LineBezier [proj v]) (reverse xs) ++ LineEnd (proj x) : start rest
      worker a (LineBezier [b] : rest) =
        let (x:xs) = reverse $ drop 1 $ mkChunks a b
        in map (\v -> LineBezier [proj v]) (reverse xs) ++ LineBezier [proj x] : worker x rest
      worker _ (LineBezier ps : rest) =
        LineBezier (map proj ps) : worker (last ps) rest
      worker _ (LineMove x:rest) = LineMove (proj x) : worker x rest
      worker _ [] = []
      lowTolerance = tolerance*tolerance
  
      proj (V2 lam phi) =
        case projectionForward p $ LonLat lam phi of
          XYCoord x y -> V2 x y
      mkChunks a b =
        let midway = lerp 0.5 a b in
        if distance (proj a) (proj b) < tolerance || distance a b < lowTolerance
          then [a, b]
          else mkChunks a midway ++ drop 1 (mkChunks midway b)
  