{-# 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 :: Double
halfPi = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2
sqrtPi :: Double
sqrtPi = Double -> Double
forall a. Floating a => a -> a
sqrt Double
forall a. Floating a => a
pi
sqrt2 :: Double
sqrt2 = Double -> Double
forall a. Floating a => a -> a
sqrt Double
2
epsilon :: Double
epsilon = Double
1.0e-12
tau :: Double
tau = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
2

toRads, cot :: Double -> Double
toRads :: Double -> Double
toRads Double
dec = Double
decDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
180 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi
cot :: Double -> Double
cot = Double -> Double
forall a. Fractional a => a -> a
recip (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
tan


srcPixel :: Pixel pixel => Image pixel -> LonLat -> pixel
srcPixel :: Image pixel -> LonLat -> pixel
srcPixel Image pixel
src (LonLat Double
lam Double
phi) =
    -- pixelAt src xPx yPx
    Vector (PixelBaseComponent pixel) -> Int -> pixel
forall a. Pixel a => Vector (PixelBaseComponent a) -> Int -> a
unsafePixelAt (Image pixel -> Vector (PixelBaseComponent pixel)
forall a. Image a -> Vector (PixelBaseComponent a)
imageData Image pixel
src) (Image pixel -> Int -> Int -> Int
forall a. Pixel a => Image a -> Int -> Int -> Int
pixelBaseIndex Image pixel
src Int
xPx Int
yPx)
  where
    !xPx :: Int
xPx = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Image pixel -> Int
forall a. Image a -> Int
imageWidth Image pixel
srcInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    !yPx :: Int
yPx = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-((Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
halfPi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Image pixel -> Int
forall a. Image a -> Int
imageHeight Image pixel
srcInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)


srcPixelFast :: Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast :: Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast Image PixelRGBA8
src Double
wMax Double
hMax (LonLat Double
lam Double
phi) = IO PixelRGBA8 -> ST s PixelRGBA8
forall a s. IO a -> ST s a
unsafeIOToST (IO PixelRGBA8 -> ST s PixelRGBA8)
-> IO PixelRGBA8 -> ST s PixelRGBA8
forall a b. (a -> b) -> a -> b
$
    Vector Word8 -> (Ptr Word8 -> IO PixelRGBA8) -> IO PixelRGBA8
forall a b. Storable a => Vector a -> (Ptr a -> IO b) -> IO b
unsafeWith (Image PixelRGBA8 -> Vector (PixelBaseComponent PixelRGBA8)
forall a. Image a -> Vector (PixelBaseComponent a)
imageData Image PixelRGBA8
src) ((Ptr Word8 -> IO PixelRGBA8) -> IO PixelRGBA8)
-> (Ptr Word8 -> IO PixelRGBA8) -> IO PixelRGBA8
forall a b. (a -> b) -> a -> b
$ \Ptr Word8
ptr -> do
      let ptr' :: Ptr b
ptr' = Ptr Word8 -> Int -> Ptr b
forall a b. Ptr a -> Int -> Ptr b
plusPtr Ptr Word8
ptr Int
idx
      Word8
r <- Ptr Word8 -> IO Word8
forall a. Storable a => Ptr a -> IO a
peek Ptr Word8
forall b. Ptr b
ptr'
      Word8
g <- Ptr Word8 -> IO Word8
forall a. Storable a => Ptr a -> IO a
peek (Ptr Word8 -> IO Word8) -> Ptr Word8 -> IO Word8
forall a b. (a -> b) -> a -> b
$ Ptr Any -> Int -> Ptr Word8
forall a b. Ptr a -> Int -> Ptr b
plusPtr Ptr Any
forall b. Ptr b
ptr' Int
1
      Word8
b <- Ptr Word8 -> IO Word8
forall a. Storable a => Ptr a -> IO a
peek (Ptr Word8 -> IO Word8) -> Ptr Word8 -> IO Word8
forall a b. (a -> b) -> a -> b
$ Ptr Any -> Int -> Ptr Word8
forall a b. Ptr a -> Int -> Ptr b
plusPtr Ptr Any
forall b. Ptr b
ptr' Int
2
      Word8
a <- Ptr Word8 -> IO Word8
forall a. Storable a => Ptr a -> IO a
peek (Ptr Word8 -> IO Word8) -> Ptr Word8 -> IO Word8
forall a b. (a -> b) -> a -> b
$ Ptr Any -> Int -> Ptr Word8
forall a b. Ptr a -> Int -> Ptr b
plusPtr Ptr Any
forall b. Ptr b
ptr' Int
3
      PixelRGBA8 -> IO PixelRGBA8
forall (m :: * -> *) a. Monad m => a -> m a
return (PixelRGBA8 -> IO PixelRGBA8) -> PixelRGBA8 -> IO PixelRGBA8
forall a b. (a -> b) -> a -> b
$ Word8 -> Word8 -> Word8 -> Word8 -> PixelRGBA8
PixelRGBA8 Word8
r Word8
g Word8
b Word8
a
  where
    !idx :: Int
idx = Image PixelRGBA8 -> Int -> Int -> Int
forall a. Pixel a => Image a -> Int -> Int -> Int
pixelBaseIndex Image PixelRGBA8
src Int
xPx Int
yPx
    !xPx :: Int
xPx = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
wMax
    !yPx :: Int
yPx = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-((Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
halfPi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
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 :: Image PixelRGBA8
-> Projection -> Projection -> Double -> Image PixelRGBA8
interpP Image PixelRGBA8
src Projection
p1 Projection
_ Double
0 = Image PixelRGBA8 -> Projection -> Image PixelRGBA8
project Image PixelRGBA8
src Projection
p1
interpP Image PixelRGBA8
src Projection
_ Projection
p2 Double
1 = Image PixelRGBA8 -> Projection -> Image PixelRGBA8
project Image PixelRGBA8
src Projection
p2
interpP !Image PixelRGBA8
src (Projection String
_label1 LonLat -> XYCoord
p1 XYCoord -> LonLat
p1_inv) (Projection String
_label2 LonLat -> XYCoord
p2 XYCoord -> LonLat
p2_inv) !Double
t = (forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8)
-> (forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8
forall a b. (a -> b) -> a -> b
$ do
    !MutableImage s PixelRGBA8
img <- Int -> Int -> ST s (MutableImage (PrimState (ST s)) PixelRGBA8)
forall px (m :: * -> *).
(Pixel px, PrimMonad m) =>
Int -> Int -> m (MutableImage (PrimState m) px)
newMutableImage Int
w Int
h

    let factor :: Int
factor = Int
2
    let l1 :: ST s ()
l1 =
          Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
x ->
            Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
y -> do
              let !x1' :: Double
x1' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
wMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !y1' :: Double
y1' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
hMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !lonlat :: LonLat
lonlat = XYCoord -> LonLat
p1_inv (XYCoord -> LonLat) -> XYCoord -> LonLat
forall a b. (a -> b) -> a -> b
$! Double -> Double -> XYCoord
XYCoord Double
x1' Double
y1'

              Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (LonLat -> Bool
validLonLat LonLat
lonlat) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                PixelRGBA8
p <- Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
forall s.
Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast Image PixelRGBA8
src Double
wMax Double
hMax LonLat
lonlat
                Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (PixelRGBA8 -> PixelBaseComponent PixelRGBA8
forall a. Pixel a => a -> PixelBaseComponent a
pixelOpacity PixelRGBA8
p Word8 -> Word8 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word8
0) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                  let XYCoord Double
x1 Double
y1 = LonLat -> XYCoord
p1 LonLat
lonlat
                      XYCoord Double
x2 Double
y2 = LonLat -> XYCoord
p2 LonLat
lonlat
                      !x3 :: Int
x3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
fromToS Double
x1 Double
x2 Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
wMax
                      !y3 :: Int
y3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double
fromToS Double
y1 Double
y2 Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hMax
                  Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
w Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
h) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$
                    MutableImage (PrimState (ST s)) PixelRGBA8
-> Int -> Int -> PixelRGBA8 -> ST s ()
forall a (m :: * -> *).
(Pixel a, PrimMonad m) =>
MutableImage (PrimState m) a -> Int -> Int -> a -> m ()
writePixel MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img Int
x3 Int
y3 PixelRGBA8
p
        l2 :: ST s ()
l2 =
          Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
x ->
            Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
y -> do
              let !x2' :: Double
x2' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
wMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !y2' :: Double
y2' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
hMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !lonlat :: LonLat
lonlat = XYCoord -> LonLat
p2_inv (Double -> Double -> XYCoord
XYCoord Double
x2' Double
y2')
              Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (LonLat -> Bool
validLonLat LonLat
lonlat) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                PixelRGBA8
p <- Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
forall s.
Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast Image PixelRGBA8
src Double
wMax Double
hMax LonLat
lonlat
                Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (PixelRGBA8 -> PixelBaseComponent PixelRGBA8
forall a. Pixel a => a -> PixelBaseComponent a
pixelOpacity PixelRGBA8
p Word8 -> Word8 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word8
0) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                  let XYCoord Double
x2 Double
y2 = LonLat -> XYCoord
p2 LonLat
lonlat
                      XYCoord Double
x1 Double
y1 = LonLat -> XYCoord
p1 LonLat
lonlat
                      !x3 :: Int
x3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
fromToS Double
x1 Double
x2 Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
wMax
                      !y3 :: Int
y3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double
fromToS Double
y1 Double
y2 Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hMax
                  Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
w Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
h) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$
                    MutableImage (PrimState (ST s)) PixelRGBA8
-> Int -> Int -> PixelRGBA8 -> ST s ()
forall a (m :: * -> *).
(Pixel a, PrimMonad m) =>
MutableImage (PrimState m) a -> Int -> Int -> a -> m ()
writePixel MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img Int
x3 Int
y3 PixelRGBA8
p
    if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5
      then ST s ()
l1 ST s () -> ST s () -> ST s ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> ST s ()
l2
      else ST s ()
l2 ST s () -> ST s () -> ST s ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> ST s ()
l1
    MutableImage (PrimState (ST s)) PixelRGBA8
-> ST s (Image PixelRGBA8)
forall a (m :: * -> *).
(Storable (PixelBaseComponent a), PrimMonad m) =>
MutableImage (PrimState m) a -> m (Image a)
unsafeFreezeImage MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img
  where
    loopTo :: t -> (t -> m a) -> m ()
loopTo t
m t -> m a
fn = t -> m ()
go t
m
      where go :: t -> m ()
go t
0 = () -> m ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
            go t
n = t -> m a
fn (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1) m a -> m () -> m ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> t -> m ()
go (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1)
    !w :: Int
w = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageWidth Image PixelRGBA8
src
    !h :: Int
h = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageHeight Image PixelRGBA8
src
    !wMax :: Double
wMax = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    !hMax :: Double
hMax = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
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 :: Image PixelRGBA8
-> Projection
-> Projection
-> (Double, Double, Double, Double)
-> (Double, Double, Double, Double)
-> Double
-> Image PixelRGBA8
interpBBP !Image PixelRGBA8
src (Projection String
_ LonLat -> XYCoord
p1 XYCoord -> LonLat
p1_inv) (Projection String
_ LonLat -> XYCoord
p2 XYCoord -> LonLat
p2_inv) (Double
fx,Double
fy,Double
fw,Double
fh) (Double
tx, Double
ty, Double
tw, Double
th) !Double
t = (forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8
forall a. (forall s. ST s a) -> a
runST ((forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8)
-> (forall s. ST s (Image PixelRGBA8)) -> Image PixelRGBA8
forall a b. (a -> b) -> a -> b
$ do
    !MutableImage s PixelRGBA8
img <- Int -> Int -> ST s (MutableImage (PrimState (ST s)) PixelRGBA8)
forall px (m :: * -> *).
(Pixel px, PrimMonad m) =>
Int -> Int -> m (MutableImage (PrimState m) px)
newMutableImage Int
w Int
h
    let factor :: Int
factor = Int
2
    let l1 :: ST s ()
l1 =
          Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
x ->
            Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
y -> do
              let !x1' :: Double
x1' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
wMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !y1' :: Double
y1' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
hMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
              Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
x1' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
fx Bool -> Bool -> Bool
&& Double
x1' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
fxDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fw Bool -> Bool -> Bool
&& Double
y1' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
fy Bool -> Bool -> Bool
&& Double
y1' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
fyDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
fh) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                let !lonlat :: LonLat
lonlat = XYCoord -> LonLat
p1_inv (XYCoord -> LonLat) -> XYCoord -> LonLat
forall a b. (a -> b) -> a -> b
$! Double -> Double -> XYCoord
XYCoord Double
x1' Double
y1'

                Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (LonLat -> Bool
validLonLat LonLat
lonlat) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                    PixelRGBA8
p <- Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
forall s.
Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast Image PixelRGBA8
src Double
wMax Double
hMax LonLat
lonlat
                    Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (PixelRGBA8 -> PixelBaseComponent PixelRGBA8
forall a. Pixel a => a -> PixelBaseComponent a
pixelOpacity PixelRGBA8
p Word8 -> Word8 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word8
0) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                      let XYCoord Double
x1 Double
y1 = LonLat -> XYCoord
p1 LonLat
lonlat
                          XYCoord Double
x2 Double
y2 = LonLat -> XYCoord
p2 LonLat
lonlat
                          !x3 :: Int
x3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
fromToS Double
x1 Double
x2 Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
wMax
                          !y3 :: Int
y3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double
fromToS Double
y1 Double
y2 Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hMax
                      Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
w Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
h) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$
                        MutableImage (PrimState (ST s)) PixelRGBA8
-> Int -> Int -> PixelRGBA8 -> ST s ()
forall a (m :: * -> *).
(Pixel a, PrimMonad m) =>
MutableImage (PrimState m) a -> Int -> Int -> a -> m ()
writePixel MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img Int
x3 Int
y3 PixelRGBA8
p
        l2 :: ST s ()
l2 =
          Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
x ->
            Int -> (Int -> ST s ()) -> ST s ()
forall t (m :: * -> *) a.
(Eq t, Num t, Monad m) =>
t -> (t -> m a) -> m ()
loopTo (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
factor) ((Int -> ST s ()) -> ST s ()) -> (Int -> ST s ()) -> ST s ()
forall a b. (a -> b) -> a -> b
$ \Int
y -> do
              let !x2' :: Double
x2' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
wMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
                  !y2' :: Double
y2' = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
hMaxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
factor)
              Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Double
x2' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
tx Bool -> Bool -> Bool
&& Double
x2' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
txDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
tw Bool -> Bool -> Bool
&& Double
y2' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
ty Bool -> Bool -> Bool
&& Double
y2' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
tyDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
th) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                let !lonlat :: LonLat
lonlat = XYCoord -> LonLat
p2_inv (Double -> Double -> XYCoord
XYCoord Double
x2' Double
y2')
                Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (LonLat -> Bool
validLonLat LonLat
lonlat) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                  PixelRGBA8
p <- Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
forall s.
Image PixelRGBA8 -> Double -> Double -> LonLat -> ST s PixelRGBA8
srcPixelFast Image PixelRGBA8
src Double
wMax Double
hMax LonLat
lonlat
                  Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (PixelRGBA8 -> PixelBaseComponent PixelRGBA8
forall a. Pixel a => a -> PixelBaseComponent a
pixelOpacity PixelRGBA8
p Word8 -> Word8 -> Bool
forall a. Eq a => a -> a -> Bool
/= Word8
0) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$ do
                    let XYCoord Double
x2 Double
y2 = LonLat -> XYCoord
p2 LonLat
lonlat
                        XYCoord Double
x1 Double
y1 = LonLat -> XYCoord
p1 LonLat
lonlat
                        !x3 :: Int
x3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double -> Double
fromToS Double
x1 Double
x2 Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
wMax
                        !y3 :: Int
y3 = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double
fromToS Double
y1 Double
y2 Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
hMax
                    Bool -> ST s () -> ST s ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
x3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
w Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 Bool -> Bool -> Bool
&& Int
y3 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
h) (ST s () -> ST s ()) -> ST s () -> ST s ()
forall a b. (a -> b) -> a -> b
$
                      MutableImage (PrimState (ST s)) PixelRGBA8
-> Int -> Int -> PixelRGBA8 -> ST s ()
forall a (m :: * -> *).
(Pixel a, PrimMonad m) =>
MutableImage (PrimState m) a -> Int -> Int -> a -> m ()
writePixel MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img Int
x3 Int
y3 PixelRGBA8
p
    if Double
t Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5
      then ST s ()
l1 ST s () -> ST s () -> ST s ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> ST s ()
l2
      else ST s ()
l2 ST s () -> ST s () -> ST s ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> ST s ()
l1

    MutableImage (PrimState (ST s)) PixelRGBA8
-> ST s (Image PixelRGBA8)
forall a (m :: * -> *).
(Storable (PixelBaseComponent a), PrimMonad m) =>
MutableImage (PrimState m) a -> m (Image a)
unsafeFreezeImage MutableImage s PixelRGBA8
MutableImage (PrimState (ST s)) PixelRGBA8
img
  where
    loopTo :: t -> (t -> m a) -> m ()
loopTo t
m t -> m a
fn = t -> m ()
go t
m
      where go :: t -> m ()
go t
0 = () -> m ()
forall (m :: * -> *) a. Monad m => a -> m a
return ()
            go t
n = t -> m a
fn (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1) m a -> m () -> m ()
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> t -> m ()
go (t
nt -> t -> t
forall a. Num a => a -> a -> a
-t
1)
    !w :: Int
w = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageWidth Image PixelRGBA8
src
    !h :: Int
h = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageHeight Image PixelRGBA8
src
    !wMax :: Double
wMax = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    !hMax :: Double
hMax = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)


eqLonLat :: LonLat -> LonLat -> Bool
eqLonLat :: LonLat -> LonLat -> Bool
eqLonLat (LonLat Double
x1 Double
y1) (LonLat Double
x2 Double
y2)
  = Double -> Double -> Bool
eqDouble Double
x1 Double
x2 Bool -> Bool -> Bool
&& Double -> Double -> Bool
eqDouble Double
y1 Double
y2
{-
eqCoords :: XYCoord -> XYCoord -> Bool
eqCoords (XYCoord x1 y1) (XYCoord x2 y2)
  = eqDouble x1 x2 && eqDouble y1 y2
-}
eqDouble :: Double -> Double -> Bool
eqDouble :: Double -> Double -> Bool
eqDouble Double
a Double
b = Double -> Double
forall a. Num a => a -> a
abs (Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon

-- | XY coordinates on a 2D plane. Valid ranges go from 0 to 1, inclusive.
data XYCoord = XYCoord !Double !Double -- 0 to 1
  deriving (ReadPrec [XYCoord]
ReadPrec XYCoord
Int -> ReadS XYCoord
ReadS [XYCoord]
(Int -> ReadS XYCoord)
-> ReadS [XYCoord]
-> ReadPrec XYCoord
-> ReadPrec [XYCoord]
-> Read XYCoord
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [XYCoord]
$creadListPrec :: ReadPrec [XYCoord]
readPrec :: ReadPrec XYCoord
$creadPrec :: ReadPrec XYCoord
readList :: ReadS [XYCoord]
$creadList :: ReadS [XYCoord]
readsPrec :: Int -> ReadS XYCoord
$creadsPrec :: Int -> ReadS XYCoord
Read,Int -> XYCoord -> ShowS
[XYCoord] -> ShowS
XYCoord -> String
(Int -> XYCoord -> ShowS)
-> (XYCoord -> String) -> ([XYCoord] -> ShowS) -> Show XYCoord
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [XYCoord] -> ShowS
$cshowList :: [XYCoord] -> ShowS
show :: XYCoord -> String
$cshow :: XYCoord -> String
showsPrec :: Int -> XYCoord -> ShowS
$cshowsPrec :: Int -> XYCoord -> ShowS
Show,XYCoord -> XYCoord -> Bool
(XYCoord -> XYCoord -> Bool)
-> (XYCoord -> XYCoord -> Bool) -> Eq XYCoord
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: XYCoord -> XYCoord -> Bool
$c/= :: XYCoord -> XYCoord -> Bool
== :: XYCoord -> XYCoord -> Bool
$c== :: XYCoord -> XYCoord -> Bool
Eq,Eq XYCoord
Eq XYCoord
-> (XYCoord -> XYCoord -> Ordering)
-> (XYCoord -> XYCoord -> Bool)
-> (XYCoord -> XYCoord -> Bool)
-> (XYCoord -> XYCoord -> Bool)
-> (XYCoord -> XYCoord -> Bool)
-> (XYCoord -> XYCoord -> XYCoord)
-> (XYCoord -> XYCoord -> XYCoord)
-> Ord XYCoord
XYCoord -> XYCoord -> Bool
XYCoord -> XYCoord -> Ordering
XYCoord -> XYCoord -> XYCoord
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: XYCoord -> XYCoord -> XYCoord
$cmin :: XYCoord -> XYCoord -> XYCoord
max :: XYCoord -> XYCoord -> XYCoord
$cmax :: XYCoord -> XYCoord -> XYCoord
>= :: XYCoord -> XYCoord -> Bool
$c>= :: XYCoord -> XYCoord -> Bool
> :: XYCoord -> XYCoord -> Bool
$c> :: XYCoord -> XYCoord -> Bool
<= :: XYCoord -> XYCoord -> Bool
$c<= :: XYCoord -> XYCoord -> Bool
< :: XYCoord -> XYCoord -> Bool
$c< :: XYCoord -> XYCoord -> Bool
compare :: XYCoord -> XYCoord -> Ordering
$ccompare :: XYCoord -> XYCoord -> Ordering
$cp1Ord :: Eq XYCoord
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 (ReadPrec [LonLat]
ReadPrec LonLat
Int -> ReadS LonLat
ReadS [LonLat]
(Int -> ReadS LonLat)
-> ReadS [LonLat]
-> ReadPrec LonLat
-> ReadPrec [LonLat]
-> Read LonLat
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [LonLat]
$creadListPrec :: ReadPrec [LonLat]
readPrec :: ReadPrec LonLat
$creadPrec :: ReadPrec LonLat
readList :: ReadS [LonLat]
$creadList :: ReadS [LonLat]
readsPrec :: Int -> ReadS LonLat
$creadsPrec :: Int -> ReadS LonLat
Read,Int -> LonLat -> ShowS
[LonLat] -> ShowS
LonLat -> String
(Int -> LonLat -> ShowS)
-> (LonLat -> String) -> ([LonLat] -> ShowS) -> Show LonLat
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [LonLat] -> ShowS
$cshowList :: [LonLat] -> ShowS
show :: LonLat -> String
$cshow :: LonLat -> String
showsPrec :: Int -> LonLat -> ShowS
$cshowsPrec :: Int -> LonLat -> ShowS
Show,LonLat -> LonLat -> Bool
(LonLat -> LonLat -> Bool)
-> (LonLat -> LonLat -> Bool) -> Eq LonLat
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: LonLat -> LonLat -> Bool
$c/= :: LonLat -> LonLat -> Bool
== :: LonLat -> LonLat -> Bool
$c== :: LonLat -> LonLat -> Bool
Eq,Eq LonLat
Eq LonLat
-> (LonLat -> LonLat -> Ordering)
-> (LonLat -> LonLat -> Bool)
-> (LonLat -> LonLat -> Bool)
-> (LonLat -> LonLat -> Bool)
-> (LonLat -> LonLat -> Bool)
-> (LonLat -> LonLat -> LonLat)
-> (LonLat -> LonLat -> LonLat)
-> Ord LonLat
LonLat -> LonLat -> Bool
LonLat -> LonLat -> Ordering
LonLat -> LonLat -> LonLat
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: LonLat -> LonLat -> LonLat
$cmin :: LonLat -> LonLat -> LonLat
max :: LonLat -> LonLat -> LonLat
$cmax :: LonLat -> LonLat -> LonLat
>= :: LonLat -> LonLat -> Bool
$c>= :: LonLat -> LonLat -> Bool
> :: LonLat -> LonLat -> Bool
$c> :: LonLat -> LonLat -> Bool
<= :: LonLat -> LonLat -> Bool
$c<= :: LonLat -> LonLat -> Bool
< :: LonLat -> LonLat -> Bool
$c< :: LonLat -> LonLat -> Bool
compare :: LonLat -> LonLat -> Ordering
$ccompare :: LonLat -> LonLat -> Ordering
$cp1Ord :: Eq LonLat
Ord)
instance Hashable LonLat where
  hashWithSalt :: Int -> LonLat -> Int
hashWithSalt Int
s (LonLat Double
a Double
b) = Int -> (Double, Double) -> Int
forall a. Hashable a => Int -> a -> Int
hashWithSalt Int
s (Double
a,Double
b)

-- | Projections are named bi-directional mappings between a sphere
--   and a 2D plane.
data Projection = Projection
  { Projection -> String
projectionLabel   :: String -- ^ Name of the projection.
  , Projection -> LonLat -> XYCoord
projectionForward :: !(LonLat -> XYCoord)
    -- ^ Mapping from longitude and latitude on a sphere to
    --   XY coordinates on a 2D plane.
  , Projection -> XYCoord -> LonLat
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 :: Image PixelRGBA8 -> Projection -> Image PixelRGBA8
project Image PixelRGBA8
src (Projection String
_label LonLat -> XYCoord
_ XYCoord -> LonLat
pInv) = (Int -> Int -> PixelRGBA8) -> Int -> Int -> Image PixelRGBA8
forall px. Pixel px => (Int -> Int -> px) -> Int -> Int -> Image px
generateImage Int -> Int -> PixelRGBA8
forall a a. (Integral a, Integral a) => a -> a -> PixelRGBA8
fn Int
w Int
h
  where
    w :: Int
w = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageWidth Image PixelRGBA8
src
    h :: Int
h = Image PixelRGBA8 -> Int
forall a. Image a -> Int
imageHeight Image PixelRGBA8
src
    fn :: a -> a -> PixelRGBA8
fn a
xPx a
yPx =
      let x :: Double
x = (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
xPx Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
          y :: Double
y = Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-(a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
yPx Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
          lonlat :: LonLat
lonlat = XYCoord -> LonLat
pInv (Double -> Double -> XYCoord
XYCoord Double
x Double
y)
      in
      if LonLat -> Bool
validLonLat LonLat
lonlat
        then Image PixelRGBA8 -> LonLat -> PixelRGBA8
forall pixel. Pixel pixel => Image pixel -> LonLat -> pixel
srcPixel Image PixelRGBA8
src LonLat
lonlat
        else Word8 -> Word8 -> Word8 -> Word8 -> PixelRGBA8
PixelRGBA8 Word8
0 Word8
0 Word8
0 Word8
0

validLonLat :: LonLat -> Bool
validLonLat :: LonLat -> Bool
validLonLat (LonLat Double
lam Double
phi) =
  Double
lam Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double
forall a. Floating a => a
pi Bool -> Bool -> Bool
&& Double
lam Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
forall a. Floating a => a
pi Bool -> Bool -> Bool
&& Double
phi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= -Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Bool -> Bool -> Bool
&& Double
phi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2

_validXYCoord :: XYCoord -> Bool
_validXYCoord :: XYCoord -> Bool
_validXYCoord (XYCoord Double
x Double
y) = Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1 Bool -> Bool -> Bool
&& Double
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0 Bool -> Bool -> Bool
&& Double
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1

-- | Returns @True@ iff a projection is consistent and complete.
isValidP :: Projection -> Bool
isValidP :: Projection -> Bool
isValidP (Projection String
_label LonLat -> XYCoord
p XYCoord -> LonLat
pInv) = [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
and
    [ Int -> Int -> Bool
forall a a. (Integral a, Integral a) => a -> a -> Bool
check Int
x Int
y
    | Int
x <- [Int
0..Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1::Int]
    , Int
y <- [Int
0..Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1::Int] ]
  where
    w :: Int
w = Int
100
    h :: Int
h = Int
100
    check :: a -> a -> Bool
check a
xPx a
yPx =
      let x :: Double
x = (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
xPx Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
wInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
          y :: Double
y = (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
yPx Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
hInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
          lonlat :: LonLat
lonlat = XYCoord -> LonLat
pInv (Double -> Double -> XYCoord
XYCoord Double
x Double
y)
          lonlat2 :: LonLat
lonlat2 = XYCoord -> LonLat
pInv (XYCoord -> LonLat) -> XYCoord -> LonLat
forall a b. (a -> b) -> a -> b
$ LonLat -> XYCoord
p LonLat
lonlat
      in Bool -> Bool
not (LonLat -> Bool
validLonLat LonLat
lonlat) Bool -> Bool -> Bool
|| LonLat -> LonLat -> Bool
eqLonLat LonLat
lonlat LonLat
lonlat2
          Bool -> Bool -> Bool
|| String -> Bool -> Bool
forall a. String -> a -> a
trace ((LonLat, LonLat) -> String
forall a. Show a => a -> String
show (LonLat
lonlat, LonLat
lonlat2)) Bool
False

-- | Translate the lower-most point of a projection by an offset.
moveBottomP :: Double -> Projection -> Projection
moveBottomP :: Double -> Projection -> Projection
moveBottomP Double
offset (Projection String
label LonLat -> XYCoord
p XYCoord -> LonLat
pInv) = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
label LonLat -> XYCoord
p' XYCoord -> LonLat
pInv'
  where
    p' :: LonLat -> XYCoord
p' (LonLat Double
lon Double
lat) =
      case LonLat -> XYCoord
p (Double -> Double -> LonLat
LonLat Double
lon Double
lat) of
        XYCoord Double
x Double
y -> Double -> Double -> XYCoord
XYCoord Double
x (Double -> Double -> Double -> Double
fromToS Double
offset Double
1 Double
y)
    pInv' :: XYCoord -> LonLat
pInv' (XYCoord Double
x Double
y) = XYCoord -> LonLat
pInv (Double -> Double -> XYCoord
XYCoord Double
x ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
offset)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
offset)))

-- | Translate the top-most point of a projection by an offset.
moveTopP :: Double -> Projection -> Projection
moveTopP :: Double -> Projection -> Projection
moveTopP Double
offset = Projection -> Projection
flipYAxisP (Projection -> Projection)
-> (Projection -> Projection) -> Projection -> Projection
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Projection -> Projection
moveBottomP Double
offset (Projection -> Projection)
-> (Projection -> Projection) -> Projection -> Projection
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Projection -> Projection
flipYAxisP

-- | Invert the Y axis of projection.
flipYAxisP :: Projection -> Projection
flipYAxisP :: Projection -> Projection
flipYAxisP (Projection String
label LonLat -> XYCoord
p XYCoord -> LonLat
pInv) = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
label LonLat -> XYCoord
p' XYCoord -> LonLat
pInv'
  where
    p' :: LonLat -> XYCoord
p' (LonLat Double
lam Double
phi) =
      let XYCoord Double
x Double
y = LonLat -> XYCoord
p (Double -> Double -> LonLat
LonLat Double
lam (Double -> Double
forall a. Num a => a -> a
negate Double
phi))
      in Double -> Double -> XYCoord
XYCoord Double
x (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y)
    pInv' :: XYCoord -> LonLat
pInv' (XYCoord Double
x Double
y) =
      let LonLat Double
lam Double
phi = XYCoord -> LonLat
pInv (Double -> Double -> XYCoord
XYCoord Double
x (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y))
      in Double -> Double -> LonLat
LonLat Double
lam (Double -> Double
forall a. Num a => a -> a
negate Double
phi)

-- | Scale X and Y axis of projection.
scaleP :: Double -> Double -> Projection -> Projection
scaleP :: Double -> Double -> Projection -> Projection
scaleP Double
xScale Double
yScale (Projection String
label LonLat -> XYCoord
p XYCoord -> LonLat
pInv) = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
label LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward LonLat
lonlat =
      case LonLat -> XYCoord
p LonLat
lonlat of
        XYCoord Double
x Double
y -> Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
0.5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xScaleDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
0.5) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
0.5)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
yScaleDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
0.5)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x Double
y) =
      XYCoord -> LonLat
pInv (XYCoord -> LonLat) -> XYCoord -> LonLat
forall a b. (a -> b) -> a -> b
$ Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
0.5)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
xScaleDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
0.5) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
0.5)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
yScaleDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
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 :: Projection -> Projection -> Double -> Projection
mergeP Projection
p1 Projection
p2 Double
t = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection (Projection -> String
projectionLabel Projection
p1 String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"/" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Projection -> String
projectionLabel Projection
p2) LonLat -> XYCoord
p XYCoord -> LonLat
pInv
  where
    p :: LonLat -> XYCoord
p LonLat
lonlat =
      let XYCoord Double
x1 Double
y1 = Projection -> LonLat -> XYCoord
projectionForward Projection
p1 LonLat
lonlat
          XYCoord Double
x2 Double
y2 = Projection -> LonLat -> XYCoord
projectionForward Projection
p2 LonLat
lonlat
      in Double -> Double -> XYCoord
XYCoord (Double -> Double -> Double -> Double
fromToS Double
x1 Double
x2 Double
t) (Double -> Double -> Double -> Double
fromToS Double
y1 Double
y2 Double
t)
    pInv :: XYCoord -> LonLat
pInv XYCoord
coord =
      let LonLat Double
lon1 Double
lat1 = Projection -> XYCoord -> LonLat
projectionInverse Projection
p1 XYCoord
coord
          LonLat Double
lon2 Double
lat2 = Projection -> XYCoord -> LonLat
projectionInverse Projection
p2 XYCoord
coord
      in
        if Double -> Double -> Bool
forall a a.
(Ord a, Ord a, Floating a, Floating a) =>
a -> a -> Bool
oob Double
lon1 Double
lat1 Bool -> Bool -> Bool
&& Double -> Double -> Bool
forall a a.
(Ord a, Ord a, Floating a, Floating a) =>
a -> a -> Bool
oob Double
lon2 Double
lat2
          then Double -> Double -> LonLat
LonLat (Double
0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0) (Double
0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0)
          else Double -> Double -> LonLat
LonLat (Double -> Double -> Double -> Double
fromToS Double
lon1 Double
lon2 Double
t) (Double -> Double -> Double -> Double
fromToS Double
lat1 Double
lat2 Double
t)
    oob :: a -> a -> Bool
oob a
lon a
lat = a
lon a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< (-a
forall a. Floating a => a
pi) Bool -> Bool -> Bool
|| a
lon a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
forall a. Floating a => a
pi Bool -> Bool -> Bool
|| a
lat a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< (-a
forall a. Floating a => a
pia -> a -> a
forall a. Fractional a => a -> a -> a
/a
2) Bool -> Bool -> Bool
|| a
lat a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
forall a. Floating a => a
pia -> a -> a
forall a. Fractional a => a -> a -> a
/a
2

-- | <<docs/gifs/doc_equirectangularP.gif>>
equirectangularP :: Projection
equirectangularP :: Projection
equirectangularP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"equirectangular" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x Double
y) = Double -> Double -> LonLat
LonLat Double
xPi Double
yPi
      where
        xPi :: Double
xPi = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x
        yPi :: Double
yPi = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double
y

-- | <<docs/gifs/doc_mercatorP.gif>>
mercatorP :: Projection
mercatorP :: Projection
mercatorP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"mercator" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
      Double -> Double -> XYCoord
XYCoord ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau)
        (Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double
forall a. Floating a => a -> a
log(Double -> Double
forall a. Floating a => a -> a
tan(Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
phiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x Double
y) = Double -> Double -> LonLat
LonLat Double
xPi (Double -> Double
forall a. Floating a => a -> a
atan (Double -> Double
forall a. Floating a => a -> a
sinh Double
yPi))
      where
        xPi :: Double
xPi = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x
        yPi :: Double
yPi = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
y

thetas :: V.Vector Double
thetas :: Vector Double
thetas = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
V.fromList ([Double] -> Vector Double) -> [Double] -> Vector Double
forall a b. (a -> b) -> a -> b
$
  (Int -> Double) -> [Int] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double
findThetaFast (Double -> Double) -> (Int -> Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
fromIndex) [Int
0 .. Int
granularity]

granularity :: Int
granularity :: Int
granularity = Int
50000

toIndex :: Double -> Int
toIndex :: Double -> Int
toIndex Double
phi = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
round ((Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
halfPi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
granularity)
fromIndex :: Int -> Double
fromIndex :: Int -> Double
fromIndex Int
x = Double -> Double -> Double -> Double
fromToS (-Double
halfPi) Double
halfPi (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
granularity)
granualize :: Double -> Double
granualize :: Double -> Double
granualize = Int -> Double
fromIndex (Int -> Double) -> (Double -> Int) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Int
toIndex

{-# INLINE mollweideP #-}
-- | <<docs/gifs/doc_mollweideP.gif>>
mollweideP :: Projection
mollweideP :: Projection
mollweideP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"mollweide" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat !Double
lam !Double
phi) =
        Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2))
      where
        x :: Double
x = (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
theta
        y :: Double
y = Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
sin Double
theta
        theta :: Double
theta = Vector Double
thetas Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
V.! Double -> Int
toIndex Double
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 -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2) (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2) Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
1) Double
1 Double
y'
        theta :: Double
theta = Double -> Double
granualize (Double -> Double
forall a. Floating a => a -> a
asin Double
y)
        lam :: Double
lam = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
cos Double
theta)
        phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
asin ((Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
thetaDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double -> Double
forall a. Floating a => a -> a
sin(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
theta))Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)

findThetaFast :: Double -> Double
findThetaFast :: Double -> Double
findThetaFast !Double
phi | Double -> Double
forall a. Num a => a -> a
abs Double
phi Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 = Double -> Double
forall a. Num a => a -> a
signum Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
halfPi
findThetaFast (D# Double#
phi) = Double# -> Double
go Double#
phi
  where
    !(D# Double#
pi#) = Double
forall a. Floating a => a
pi
    go :: Double# -> Double
go Double#
acc =
      let c :: Double#
c = Double# -> Double#
cosDouble# (Double#
acc Double# -> Double# -> Double#
+## Double#
acc)
          s :: Double#
s = Double# -> Double#
sinDouble# (Double#
acc Double# -> Double# -> Double#
+## Double#
acc)
          next :: Double#
next =
            Double#
acc Double# -> Double# -> Double#
-##
            (Double#
acc Double# -> Double# -> Double#
+## Double#
acc Double# -> Double# -> Double#
+## Double#
s Double# -> Double# -> Double#
-## Double#
pi# Double# -> Double# -> Double#
*## Double# -> Double#
sinDouble# Double#
phi)
            Double# -> Double# -> Double#
/## (Double#
2.0## Double# -> Double# -> Double#
+## Double#
c Double# -> Double# -> Double#
+## Double#
c) in
      if Double -> Double
forall a. Num a => a -> a
abs (Double# -> Double
D# (Double#
next Double# -> Double# -> Double#
-## Double#
acc)) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
        then Double# -> Double
D# Double#
next
        else Double# -> Double
go Double#
next


-- | <<docs/gifs/doc_hammerP.gif>>
hammerP :: Projection
hammerP :: Projection
hammerP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"hammer" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
        Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sqrt2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2))
      where
        x :: Double
x = (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
cos Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
sin (Double
lamDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2))Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double -> Double
forall a. Floating a => a -> a
cos Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
cos (Double
lamDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2))
        y :: Double
y = (Double
sqrt2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
sin Double
phi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double -> Double
forall a. Floating a => a -> a
cos Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
cos (Double
lamDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2))
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2) (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2) Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
sqrt2) Double
sqrt2 Double
y'
        z :: Double
z = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4)Double -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
yDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)Double -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)
        lam :: Double
lam = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x) (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
zDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1))
        phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
asin (Double
zDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)

-- | <<docs/gifs/doc_cylindricalEqualAreaP.gif>>
cylindricalEqualAreaP :: Double -> Projection
cylindricalEqualAreaP :: Double -> Projection
cylindricalEqualAreaP Double
phi0 = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"lambert" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    cosPhi0 :: Double
cosPhi0 = Double -> Double
forall a. Floating a => a -> a
cos Double
phi0
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
      Double -> Double -> XYCoord
XYCoord ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double -> Double
forall a. Floating a => a -> a
sin Double
phiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cosPhi0Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cosPhi0)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cosPhi0)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
x (Double -> Double
forall a. Floating a => a -> a
asin Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
asin Double
cosPhi0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
halfPi))
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
1) Double
1 Double
y' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosPhi0

-- | <<docs/gifs/doc_lambertP.gif>>
lambertP :: Projection
lambertP :: Projection
lambertP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"lambert" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
      Double -> Double -> XYCoord
XYCoord ((Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double -> Double
forall a. Floating a => a -> a
sin Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
x (Double -> Double
forall a. Floating a => a -> a
asin Double
y)
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
1) Double
1 Double
y'

-- | <<docs/gifs/doc_bottomleyP.gif>>
bottomleyP :: Double -> Projection
bottomleyP :: Double -> Projection
bottomleyP !Double
phi_1 = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"bottomley" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
        Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)
      where
        x :: Double
x = (Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
e) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sin Double
phi_1
        y :: Double
y = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
e
        rho :: Double
rho = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
phi
        e :: Double
e = Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi_1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
rho Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double
y'
        x1 :: Double
x1 = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi_1
        y1 :: Double
y1 = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y
        rho :: Double
rho = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
x1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y1Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y1)
        e :: Double
e = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
x1 Double
y1
        lam :: Double
lam = (if Double
rho Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 else Double
rho Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sin Double
rho) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
eDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
sin Double
phi_1
        phi :: Double
phi = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho

-- | <<docs/gifs/doc_sinusoidalP.gif>>
sinusoidalP :: Projection
sinusoidalP :: Projection
sinusoidalP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"sinusoidal" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
        Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)
      where
        x :: Double
x = Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi
        y :: Double
y = Double
phi
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat (Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
cos Double
y) Double
y
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double
y'

-- | <<docs/gifs/doc_wernerP.gif>>
wernerP :: Projection
wernerP :: Projection
wernerP = Double -> Projection -> Projection
moveTopP Double
0.23 (Projection -> Projection) -> Projection -> Projection
forall a b. (a -> b) -> a -> b
$ String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"werner" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
        Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)
      where
        rho :: Double
rho = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
phi
        e :: Double
e = if Double
rho Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
rho else Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
rho Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho
        x :: Double
x = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
e
        y :: Double
y = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
e
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double
y'
        rho :: Double
rho = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
xDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)Double -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)
        e :: Double
e = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
x (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y)
        phi :: Double
phi = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho
        lam :: Double
lam = Double
e Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rho Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sin Double
rho

-- FIXME: find right scale and position.
-- | <<docs/gifs/doc_bonneP.gif>>
bonneP :: Double -> Projection
bonneP :: Double -> Projection
bonneP Double
0 = Projection
sinusoidalP
bonneP Double
phi_0 = Double -> Projection -> Projection
moveTopP (-Double
0.17Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
factor) (Projection -> Projection) -> Projection -> Projection
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Projection -> Projection
scaleP Double
1 (Double -> Double -> Double -> Double
fromToS Double
1 Double
0.65 Double
factor) (Projection -> Projection) -> Projection -> Projection
forall a b. (a -> b) -> a -> b
$
    String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"bonne" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    factor :: Double
factor = Double -> Double
forall a. Floating a => a -> a
sin Double
phi_0 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sin (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4)
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi ) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
halfPi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi)
      where
        cotPhi0 :: Double
cotPhi0 = Double -> Double
cot Double
phi_0
        rho :: Double
rho = Double
cotPhi0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
phi
        e :: Double
e = (Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
rho
        x :: Double
x = Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
e
        y :: Double
y = Double
cotPhi0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
e
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double
y'
        cotPhi0 :: Double
cotPhi0 = Double -> Double
cot Double
phi_0
        rho :: Double
rho = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
cot Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y)Double -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2)
        -- e = atan2 x (cot phi_1 + phi_1 + y)
        phi :: Double
phi = Double
cotPhi0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
rho
        lam :: Double
lam = Double
rho Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
x (Double
cotPhi0Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y)

-- | <<docs/gifs/doc_orthoP.gif>>
orthoP :: LonLat -> Projection
orthoP :: LonLat -> Projection
orthoP (LonLat Double
lam_0 Double
phi_0) = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"ortho" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi)
        | Double
cosc Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0  =
            let ang :: Double
ang = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y Double
x
                xV :: Double
xV = Double -> Double
forall a. Floating a => a -> a
cos Double
ang
                yV :: Double
yV = Double -> Double
forall a. Floating a => a -> a
sin Double
ang
                xPos :: Double
xPos = Double
7Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
32 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ((Double
xVDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
9Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
16)
            in -- trace (show (x,y)) $
              Double -> Double -> XYCoord
XYCoord Double
xPos ((Double
yVDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
            --XYCoord (0/0) (0/0)
        | Bool
otherwise = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+(Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
9))Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
9Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
2)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
      where
        x :: Double
x = Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lam_0)
        y :: Double
y = Double -> Double
forall a. Floating a => a -> a
cos Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sin Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lam_0)
        cosc :: Double
cosc = Double -> Double
forall a. Floating a => a -> a
sin Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
lam_0)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
9) (Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
9) Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
1) Double
1 Double
y'
        lam :: Double
lam = Double -> Double -> Double -> Double
forall a. (Ord a, Num a) => a -> a -> a -> a
wrap (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$
          Double
lam_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
c) (Double
rho Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi_0)
        phi :: Double
phi = Double -> Double -> Double -> Double
forall a. (Ord a, Num a) => a -> a -> a -> a
wrap (-Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$
          Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double
forall a. Floating a => a -> a
cos Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
phi_0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi_0)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
rho)
        rho :: Double
rho = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)
        c :: Double
c = Double -> Double
forall a. Floating a => a -> a
asin Double
rho
        wrap :: a -> a -> a -> a
wrap a
lower a
upper a
v
          | a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
upper = a
va -> a -> a
forall a. Num a => a -> a -> a
-a
uppera -> a -> a
forall a. Num a => a -> a -> a
+a
lower
          | a
v a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
lower = a
va -> a -> a
forall a. Num a => a -> a -> a
+a
uppera -> a -> a
forall a. Num a => a -> a -> a
-a
lower
          | Bool
otherwise = a
v

-- | <<docs/gifs/doc_cassiniP.gif>>
cassiniP :: Projection
cassiniP :: Projection
cassiniP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"cassini" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) =
      Double -> Double -> XYCoord
XYCoord ((Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
lam)Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
halfPi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
forall a. Floating a => a
pi) ((Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double -> Double
forall a. Floating a => a -> a
tan Double
phi) (Double -> Double
forall a. Floating a => a -> a
cos Double
lam)Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
tau)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS (-Double
halfPi) Double
halfPi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS (-Double
forall a. Floating a => a
pi) Double
forall a. Floating a => a
pi Double
y'
        lam :: Double
lam = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double -> Double
forall a. Floating a => a -> a
tan Double
x) (Double -> Double
forall a. Floating a => a -> a
cos Double
y)
        phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
asin (Double -> Double
forall a. Floating a => a -> a
sin Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
x)

-- | <<docs/gifs/doc_augustP.gif>>
augustP :: Projection
augustP :: Projection
augustP = Double -> Double -> Projection -> Projection
scaleP Double
0.70 Double
0.70 (Projection -> Projection) -> Projection -> Projection
forall a b. (a -> b) -> a -> b
$ String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"august"  LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    xHi :: Double
xHi = Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
3
    xLo :: Double
xLo = -Double
xHi
    yHi :: Double
yHi = Double
8 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
    yLo :: Double
yLo = -Double
yHi
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xPosDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yPosDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        tanPhi :: Double
tanPhi = Double -> Double
forall a. Floating a => a -> a
tan (Double
phiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
        k :: Double
k = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tanPhi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
tanPhi)
        c :: Double
c = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos (Double
lam Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
        x :: Double
x = Double -> Double
forall a. Floating a => a -> a
sin (Double
lamDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c
        y :: Double
y = Double
tanPhi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c
        x2 :: Double
x2 = Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x
        y2 :: Double
y2 = Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y
        xPos :: Double
xPos = Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y2)
        yPos :: Double
yPos = Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y2)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
3 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
8
        x2 :: Double
x2 = Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x
        y2 :: Double
y2 = Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y
        s :: Double
s = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y2
        sin3Eta :: Double
sin3Eta = Double -> Double
forall a. Floating a => a -> a
sqrt ((Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sqrt (Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2)
        eta :: Double
eta = Double -> Double
forall a. Floating a => a -> a
asin Double
sin3Eta Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
        xi :: Double
xi = if Double
sin3Eta Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 then Double -> Double
forall a. Floating a => a -> a
acosh (Double -> Double
forall a. Num a => a -> a
abs (Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sin3Eta)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3 else Double -> Double
forall a. Floating a => a -> a
asinh (Double -> Double
forall a. Num a => a -> a
abs Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
3
        cosEta :: Double
cosEta = Double -> Double
forall a. Floating a => a -> a
cos Double
eta
        coshXi :: Double
coshXi = Double -> Double
forall a. Floating a => a -> a
cosh Double
xi
        d :: Double
d = Double
coshXi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
coshXi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
cosEta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosEta
        lam :: Double
lam = Double -> Double
forall a. Num a => a -> a
signum Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double -> Double
forall a. Floating a => a -> a
sinh Double
xi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosEta) (Double
0.25 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
d)
        phi :: Double
phi = Double -> Double
forall a. Num a => a -> a
signum Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
coshXi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
eta) (Double
0.25 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
d)

-- | <<docs/gifs/doc_collignonP.gif>>
collignonP :: Projection
collignonP :: Projection
collignonP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"collignon" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    yHi :: Double
yHi = Double
sqrtPi
    yLo :: Double
yLo = Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sqrt2)
    xLo :: Double
xLo = -Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
sqrtPi)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2
    xHi :: Double
xHi = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
sqrtPi)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
sqrt2
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        alpha :: Double
alpha = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
sin Double
phi)
        x :: Double
x = (Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sqrtPi) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
alpha
        y :: Double
y = Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
alpha)
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        l :: Double
l = (Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1)Double -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
2
        lam :: Double
lam = if Double
l Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 then Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
l) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 else Double
0
        phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
asin (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
l)

-- | <<docs/gifs/doc_eckert1P.gif>>
eckert1P :: Projection
eckert1P :: Projection
eckert1P = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"eckert1" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    alpha :: Double
alpha = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
8 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi))
    yLo :: Double
yLo = -Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
halfPi
    yHi :: Double
yHi = Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
halfPi
    xLo :: Double
xLo = -Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi
    xHi :: Double
xHi = Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        x :: Double
x = Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Num a => a -> a
abs Double
phi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
forall a. Floating a => a
pi)
        y :: Double
y = Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
phi
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        phi :: Double
phi = Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
alpha
        lam :: Double
lam = Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Num a => a -> a
abs Double
phi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
forall a. Floating a => a
pi))

-- | <<docs/gifs/doc_eckert3P.gif>>
eckert3P :: Projection
eckert3P :: Projection
eckert3P = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"eckert3" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    k :: Double
k = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi))
    yLo :: Double
yLo = Double -> Double
forall a. Num a => a -> a
negate Double
yHi
    yHi :: Double
yHi = Double
4Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
halfPi
    xLo :: Double
xLo = Double -> Double
forall a. Num a => a -> a
negate Double
xHi
    xHi :: Double
xHi = Double
4Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
forall a. Floating a => a
pi
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        x :: Double
x = Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
phiDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
phiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)))
        y :: Double
y = Double
4 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
phi
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        lam :: Double
lam = Double
x  Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
kDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
forall a. Floating a => a
pi)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi)))
        phi :: Double
phi = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
kDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
4

-- | <<docs/gifs/doc_eckert5P.gif>>
eckert5P :: Projection
eckert5P :: Projection
eckert5P = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"eckert5" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    k :: Double
k = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi)
    yLo :: Double
yLo = Double -> Double
forall a. Num a => a -> a
negate Double
yHi
    yHi :: Double
yHi = Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
k
    xLo :: Double
xLo = Double -> Double
forall a. Num a => a -> a
negate Double
xHi
    xHi :: Double
xHi = Double
tauDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
k
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        x :: Double
x = Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos Double
phi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k
        y :: Double
y = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
phi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
k
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        lam :: Double
lam = Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos Double
phi)
        phi :: Double
phi = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2

{-# INLINE faheyP #-}
-- | <<docs/gifs/doc_faheyP.gif>>
faheyP :: Projection
faheyP :: Projection
faheyP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"fahey" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    faheyK :: Double
faheyK = Double -> Double
forall a. Floating a => a -> a
cos (Double -> Double
toRads Double
35)
    yLo :: Double
yLo = Double -> Double
forall a. Num a => a -> a
negate Double
yHi
    yHi :: Double
yHi = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
faheyK
    xLo :: Double
xLo = Double -> Double
forall a. Num a => a -> a
negate Double
xHi
    xHi :: Double
xHi = Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
faheyK Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
16Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
9
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        t :: Double
t = Double -> Double
forall a. Floating a => a -> a
tan (Double
phiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
        x :: Double
x = Double
lam Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
faheyK Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)
        y :: Double
y = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
faheyK) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        t :: Double
t = Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
faheyK)
        lam :: Double
lam = Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
faheyK Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
tDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t))
        phi :: Double
phi = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 Double
y (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
faheyK)

{-# INLINE foucautP #-}
-- | <<docs/gifs/doc_foucautP.gif>>
foucautP :: Projection
foucautP :: Projection
foucautP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"foucaut" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    yLo :: Double
yLo = Double -> Double
forall a. Num a => a -> a
negate Double
yHi
    yHi :: Double
yHi = Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan (Double
halfPiDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
    xLo :: Double
xLo = Double -> Double
forall a. Num a => a -> a
negate Double
xHi
    xHi :: Double
xHi = Double
tauDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
sqrtPi
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi) = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        k :: Double
k = Double
phi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        cosk :: Double
cosk = Double -> Double
forall a. Floating a => a -> a
cos Double
k
        x :: Double
x = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lam Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosk Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosk
        y :: Double
y = Double
sqrtPi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
tan Double
k
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y') = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x'
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y'
        k :: Double
k = Double -> Double
forall a. Floating a => a -> a
atan (Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
sqrtPi)
        cosk :: Double
cosk = Double -> Double
forall a. Floating a => a -> a
cos Double
k
        phi :: Double
phi = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
k
        lam :: Double
lam = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
sqrtPi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
cos Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosk Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cosk)

{-# INLINE lagrangeP #-}
-- | <<docs/gifs/doc_lagrangeP.gif>>
lagrangeP :: Projection
lagrangeP :: Projection
lagrangeP = String -> (LonLat -> XYCoord) -> (XYCoord -> LonLat) -> Projection
Projection String
"lagrange" LonLat -> XYCoord
forward XYCoord -> LonLat
inverse
  where
    yLo :: Double
yLo = Double -> Double
forall a. Num a => a -> a
negate Double
yHi
    yHi :: Double
yHi = Double
2
    xLo :: Double
xLo = Double -> Double
forall a. Num a => a -> a
negate Double
xHi
    xHi :: Double
xHi = Double
2
    n :: Double
n = Double
0.5
    forward :: LonLat -> XYCoord
forward (LonLat Double
lam Double
phi)
        | Double -> Double
forall a. Num a => a -> a
abs (Double -> Double
forall a. Num a => a -> a
abs Double
phi Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
halfPi) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon = Double -> Double -> XYCoord
XYCoord Double
0.5 (if Double
phi Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 then Double
0 else Double
1)
        | Bool
otherwise = Double -> Double -> XYCoord
XYCoord ((Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
xHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
xLo)) ((Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
yHiDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
yLo))
      where
        sinPhi :: Double
sinPhi = Double -> Double
forall a. Floating a => a -> a
sin Double
phi
        v :: Double
v = ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
sinPhi) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
sinPhi))Double -> Double -> Double
forall a. Floating a => a -> a -> a
**(Double
nDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)
        c :: Double
c = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
v) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
cos (Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
n)
        x :: Double
x = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin (Double
lamDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
n) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
c
        y :: Double
y = (Double
v Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
v) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
c
    inverse :: XYCoord -> LonLat
inverse (XYCoord Double
x' Double
y')
        | Double -> Double
forall a. Num a => a -> a
abs (Double
y'Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
epsilon
          -- = LonLat 0 (signum y * halfPi)
          = Double -> Double -> LonLat
LonLat Double
100 Double
100
        | Bool
otherwise = Double -> Double -> LonLat
LonLat Double
lam Double
phi
      where
        x :: Double
x = Double -> Double -> Double -> Double
fromToS Double
xLo Double
xHi Double
x' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        y :: Double
y = Double -> Double -> Double -> Double
fromToS Double
yLo Double
yHi Double
y' Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        x2 :: Double
x2 = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x
        y2 :: Double
y2 = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y
        t :: Double
t = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y2)
        t' :: Double
t' = ((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
t) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n)
        lam :: Double
lam = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x) (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x2Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
y2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
n
        phi :: Double
phi = Double -> Double
forall a. Floating a => a -> a
asin ((Double
t'Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
t'Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1))















-- | Map for all features and render the geometry.
drawFeatureCollection :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
drawFeatureCollection :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
drawFeatureCollection GeoFeatureCollection a
geo a -> SVG -> SVG
fn = [SVG] -> SVG
mkGroup
  [ a -> SVG -> SVG
fn (GeoFeature a
feature GeoFeature a -> Getting a (GeoFeature a) a -> a
forall s a. s -> Getting a s a -> a
^. Getting a (GeoFeature a) a
forall a1 a2. Lens (GeoFeature a1) (GeoFeature a2) a1 a2
properties) (SVG -> SVG) -> SVG -> SVG
forall a b. (a -> b) -> a -> b
$ GeospatialGeometry -> SVG
renderGeometry (GeoFeature a
feature GeoFeature a
-> Getting GeospatialGeometry (GeoFeature a) GeospatialGeometry
-> GeospatialGeometry
forall s a. s -> Getting a s a -> a
^. Getting GeospatialGeometry (GeoFeature a) GeospatialGeometry
forall a. Lens' (GeoFeature a) GeospatialGeometry
geometry)
  | GeoFeature a
feature <- Seq (GeoFeature a) -> [GeoFeature a]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList (GeoFeatureCollection a
geo GeoFeatureCollection a
-> Getting
     (Seq (GeoFeature a)) (GeoFeatureCollection a) (Seq (GeoFeature a))
-> Seq (GeoFeature a)
forall s a. s -> Getting a s a -> a
^. Getting
  (Seq (GeoFeature a)) (GeoFeatureCollection a) (Seq (GeoFeature a))
forall a1 a2.
Lens
  (GeoFeatureCollection a1)
  (GeoFeatureCollection a2)
  (Seq (GeoFeature a1))
  (Seq (GeoFeature a2))
geofeatures)
  ]

{-# INLINE loadFeatureCollection #-}
-- | Load GeoJSON from a filepath and render the geometry.
loadFeatureCollection :: FromJSON a => FilePath -> (a -> SVG -> SVG) -> SVG
loadFeatureCollection :: String -> (a -> SVG -> SVG) -> SVG
loadFeatureCollection String
path = IO ((a -> SVG -> SVG) -> SVG) -> (a -> SVG -> SVG) -> SVG
forall a. IO a -> a
unsafePerformIO (IO ((a -> SVG -> SVG) -> SVG) -> (a -> SVG -> SVG) -> SVG)
-> IO ((a -> SVG -> SVG) -> SVG) -> (a -> SVG -> SVG) -> SVG
forall a b. (a -> b) -> a -> b
$ do
  Maybe (GeoFeatureCollection a)
mbGeo <- String -> IO (Maybe (GeoFeatureCollection a))
forall a. FromJSON a => String -> IO (Maybe a)
decodeFileStrict String
path
  case Maybe (GeoFeatureCollection a)
mbGeo of
    Maybe (GeoFeatureCollection a)
Nothing  -> String -> IO ((a -> SVG -> SVG) -> SVG)
forall a. HasCallStack => String -> a
error (String -> IO ((a -> SVG -> SVG) -> SVG))
-> String -> IO ((a -> SVG -> SVG) -> SVG)
forall a b. (a -> b) -> a -> b
$ String
"loadFeatureCollection: Invalid GeoJSON: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
path
    Just GeoFeatureCollection a
geo -> ((a -> SVG -> SVG) -> SVG) -> IO ((a -> SVG -> SVG) -> SVG)
forall (m :: * -> *) a. Monad m => a -> m a
return (GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
forall a. GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
drawFeatureCollection GeoFeatureCollection a
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 :: GeospatialGeometry -> SVG
renderGeometry GeospatialGeometry
shape =
  case GeospatialGeometry
shape of
    MultiPolygon GeoMultiPolygon
mpolygon ->
      [SVG] -> SVG
mkGroup ([SVG] -> SVG) -> [SVG] -> SVG
forall a b. (a -> b) -> a -> b
$ (GeoPolygon -> SVG) -> [GeoPolygon] -> [SVG]
forall a b. (a -> b) -> [a] -> [b]
map (GeospatialGeometry -> SVG
renderGeometry (GeospatialGeometry -> SVG)
-> (GeoPolygon -> GeospatialGeometry) -> GeoPolygon -> SVG
forall b c a. (b -> c) -> (a -> b) -> a -> c
. GeoPolygon -> GeospatialGeometry
Polygon) ([GeoPolygon] -> [SVG]) -> [GeoPolygon] -> [SVG]
forall a b. (a -> b) -> a -> b
$ Seq GeoPolygon -> [GeoPolygon]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList (GeoMultiPolygon -> Seq GeoPolygon
splitGeoMultiPolygon GeoMultiPolygon
mpolygon)
    Polygon GeoPolygon
poly ->
      [SVG] -> SVG
mkGroup
      [ [(Double, Double)] -> SVG
mkLinePath [(Double, Double)]
section
      | [(Double, Double)]
section <- [(Double, Double)] -> [[(Double, Double)]]
forall (f :: * -> *) a. Applicative f => a -> f a
pure
          [ (Double
x, Double
y)
          | PointXY Double
x Double
y <- (GeoPositionWithoutCRS -> PointXY)
-> [GeoPositionWithoutCRS] -> [PointXY]
forall a b. (a -> b) -> [a] -> [b]
map GeoPositionWithoutCRS -> PointXY
retrieveXY (LinearRing GeoPositionWithoutCRS -> [GeoPositionWithoutCRS]
forall a. LinearRing a -> [a]
fromLinearRing ([LinearRing GeoPositionWithoutCRS]
-> LinearRing GeoPositionWithoutCRS
forall a. [a] -> a
head (Seq (LinearRing GeoPositionWithoutCRS)
-> [LinearRing GeoPositionWithoutCRS]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList (GeoPolygon
polyGeoPolygon
-> Getting
     (Seq (LinearRing GeoPositionWithoutCRS))
     GeoPolygon
     (Seq (LinearRing GeoPositionWithoutCRS))
-> Seq (LinearRing GeoPositionWithoutCRS)
forall s a. s -> Getting a s a -> a
^.Getting
  (Seq (LinearRing GeoPositionWithoutCRS))
  GeoPolygon
  (Seq (LinearRing GeoPositionWithoutCRS))
Iso' GeoPolygon (Seq (LinearRing GeoPositionWithoutCRS))
unGeoPolygon))))
          ]
      ]
    Line GeoLine
line ->
      [(Double, Double)] -> SVG
mkLinePath
      [ (Double
x, Double
y)
      | PointXY Double
x Double
y <- (GeoPositionWithoutCRS -> PointXY)
-> [GeoPositionWithoutCRS] -> [PointXY]
forall a b. (a -> b) -> [a] -> [b]
map GeoPositionWithoutCRS -> PointXY
retrieveXY (LineString GeoPositionWithoutCRS -> [GeoPositionWithoutCRS]
forall a. LineString a -> [a]
Line.fromLineString (GeoLine
line GeoLine
-> Getting
     (LineString GeoPositionWithoutCRS)
     GeoLine
     (LineString GeoPositionWithoutCRS)
-> LineString GeoPositionWithoutCRS
forall s a. s -> Getting a s a -> a
^. Getting
  (LineString GeoPositionWithoutCRS)
  GeoLine
  (LineString GeoPositionWithoutCRS)
Iso' GeoLine (LineString GeoPositionWithoutCRS)
unGeoLine))
      ]
    MultiLine GeoMultiLine
ml ->
      [SVG] -> SVG
mkGroup ([SVG] -> SVG) -> [SVG] -> SVG
forall a b. (a -> b) -> a -> b
$ (GeoLine -> SVG) -> [GeoLine] -> [SVG]
forall a b. (a -> b) -> [a] -> [b]
map (GeospatialGeometry -> SVG
renderGeometry (GeospatialGeometry -> SVG)
-> (GeoLine -> GeospatialGeometry) -> GeoLine -> SVG
forall b c a. (b -> c) -> (a -> b) -> a -> c
. GeoLine -> GeospatialGeometry
Line) ([GeoLine] -> [SVG]) -> [GeoLine] -> [SVG]
forall a b. (a -> b) -> a -> b
$ Seq GeoLine -> [GeoLine]
forall (t :: * -> *) a. Foldable t => t a -> [a]
toList (GeoMultiLine -> Seq GeoLine
splitGeoMultiLine GeoMultiLine
ml)
    GeospatialGeometry
_ -> SVG
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 :: Projection -> SVG -> SVG
applyProjection = Double -> Projection -> SVG -> SVG
applyProjection' Double
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' :: Double -> Projection -> SVG -> SVG
applyProjection' Double
tolerance Projection
p = ([LineCommand] -> [LineCommand]) -> SVG -> SVG
mapSvgLines [LineCommand] -> [LineCommand]
start
  where
    start :: [LineCommand] -> [LineCommand]
start (LineMove RPoint
x:[LineCommand]
rest) = RPoint -> LineCommand
LineMove (RPoint -> RPoint
proj RPoint
x) LineCommand -> [LineCommand] -> [LineCommand]
forall a. a -> [a] -> [a]
: RPoint -> [LineCommand] -> [LineCommand]
worker RPoint
x [LineCommand]
rest
    start [LineCommand]
_                 = []
    worker :: RPoint -> [LineCommand] -> [LineCommand]
worker RPoint
a (LineEnd RPoint
b : [LineCommand]
rest) =
      let (RPoint
x:[RPoint]
xs) = [RPoint] -> [RPoint]
forall a. [a] -> [a]
reverse ([RPoint] -> [RPoint]) -> [RPoint] -> [RPoint]
forall a b. (a -> b) -> a -> b
$ Int -> [RPoint] -> [RPoint]
forall a. Int -> [a] -> [a]
drop Int
1 ([RPoint] -> [RPoint]) -> [RPoint] -> [RPoint]
forall a b. (a -> b) -> a -> b
$ RPoint -> RPoint -> [RPoint]
mkChunks RPoint
a RPoint
b
      in (RPoint -> LineCommand) -> [RPoint] -> [LineCommand]
forall a b. (a -> b) -> [a] -> [b]
map (\RPoint
v -> [RPoint] -> LineCommand
LineBezier [RPoint -> RPoint
proj RPoint
v]) ([RPoint] -> [RPoint]
forall a. [a] -> [a]
reverse [RPoint]
xs) [LineCommand] -> [LineCommand] -> [LineCommand]
forall a. [a] -> [a] -> [a]
++ RPoint -> LineCommand
LineEnd (RPoint -> RPoint
proj RPoint
x) LineCommand -> [LineCommand] -> [LineCommand]
forall a. a -> [a] -> [a]
: [LineCommand] -> [LineCommand]
start [LineCommand]
rest
    worker RPoint
a (LineBezier [RPoint
b] : [LineCommand]
rest) =
      let (RPoint
x:[RPoint]
xs) = [RPoint] -> [RPoint]
forall a. [a] -> [a]
reverse ([RPoint] -> [RPoint]) -> [RPoint] -> [RPoint]
forall a b. (a -> b) -> a -> b
$ Int -> [RPoint] -> [RPoint]
forall a. Int -> [a] -> [a]
drop Int
1 ([RPoint] -> [RPoint]) -> [RPoint] -> [RPoint]
forall a b. (a -> b) -> a -> b
$ RPoint -> RPoint -> [RPoint]
mkChunks RPoint
a RPoint
b
      in (RPoint -> LineCommand) -> [RPoint] -> [LineCommand]
forall a b. (a -> b) -> [a] -> [b]
map (\RPoint
v -> [RPoint] -> LineCommand
LineBezier [RPoint -> RPoint
proj RPoint
v]) ([RPoint] -> [RPoint]
forall a. [a] -> [a]
reverse [RPoint]
xs) [LineCommand] -> [LineCommand] -> [LineCommand]
forall a. [a] -> [a] -> [a]
++ [RPoint] -> LineCommand
LineBezier [RPoint -> RPoint
proj RPoint
x] LineCommand -> [LineCommand] -> [LineCommand]
forall a. a -> [a] -> [a]
: RPoint -> [LineCommand] -> [LineCommand]
worker RPoint
x [LineCommand]
rest
    worker RPoint
_ (LineBezier [RPoint]
ps : [LineCommand]
rest) =
      [RPoint] -> LineCommand
LineBezier ((RPoint -> RPoint) -> [RPoint] -> [RPoint]
forall a b. (a -> b) -> [a] -> [b]
map RPoint -> RPoint
proj [RPoint]
ps) LineCommand -> [LineCommand] -> [LineCommand]
forall a. a -> [a] -> [a]
: RPoint -> [LineCommand] -> [LineCommand]
worker ([RPoint] -> RPoint
forall a. [a] -> a
last [RPoint]
ps) [LineCommand]
rest
    worker RPoint
_ (LineMove RPoint
x:[LineCommand]
rest) = RPoint -> LineCommand
LineMove (RPoint -> RPoint
proj RPoint
x) LineCommand -> [LineCommand] -> [LineCommand]
forall a. a -> [a] -> [a]
: RPoint -> [LineCommand] -> [LineCommand]
worker RPoint
x [LineCommand]
rest
    worker RPoint
_ [] = []
    lowTolerance :: Double
lowTolerance = Double
toleranceDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
tolerance

    proj :: RPoint -> RPoint
proj (V2 Double
lam Double
phi) =
      case Projection -> LonLat -> XYCoord
projectionForward Projection
p (LonLat -> XYCoord) -> LonLat -> XYCoord
forall a b. (a -> b) -> a -> b
$ Double -> Double -> LonLat
LonLat Double
lam Double
phi of
        XYCoord Double
x Double
y -> Double -> Double -> RPoint
forall a. a -> a -> V2 a
V2 Double
x Double
y
    mkChunks :: RPoint -> RPoint -> [RPoint]
mkChunks RPoint
a RPoint
b =
      let midway :: RPoint
midway = Double -> RPoint -> RPoint -> RPoint
forall (f :: * -> *) a.
(Additive f, Num a) =>
a -> f a -> f a -> f a
lerp Double
0.5 RPoint
a RPoint
b in
      if RPoint -> RPoint -> Double
forall (f :: * -> *) a. (Metric f, Floating a) => f a -> f a -> a
distance (RPoint -> RPoint
proj RPoint
a) (RPoint -> RPoint
proj RPoint
b) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
tolerance Bool -> Bool -> Bool
|| RPoint -> RPoint -> Double
forall (f :: * -> *) a. (Metric f, Floating a) => f a -> f a -> a
distance RPoint
a RPoint
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
lowTolerance
        then [RPoint
a, RPoint
b]
        else RPoint -> RPoint -> [RPoint]
mkChunks RPoint
a RPoint
midway [RPoint] -> [RPoint] -> [RPoint]
forall a. [a] -> [a] -> [a]
++ Int -> [RPoint] -> [RPoint]
forall a. Int -> [a] -> [a]
drop Int
1 (RPoint -> RPoint -> [RPoint]
mkChunks RPoint
midway RPoint
b)