Geo projections (#52)

Former-commit-id: bc6475b83149db308bf541f1a49b6a8e67729a81

.azure/azure-osx-template.yml

@@ -20,7 +20,7 @@ jobs:
   steps:
   - task: Cache@2
     inputs:
-      key: ${{ parameters.name }} | "${{ parameters.vmImage }}" | $(STACK_YAML) | stack-root | $(Agent.OS) | version4
+      key: ${{ parameters.name }} | "${{ parameters.vmImage }}" | $(STACK_YAML) | stack-root | $(Agent.OS) | version5
       path: $(STACK_ROOT)
       cacheHitVar: CACHE_RESTORED
     displayName: Cache stack root

RESOURCES.md

@@ -2,3 +2,4 @@
 * https://www.solarsystemscope.com/textures/
 * https://cc0textures.com/
 * https://www.publicdomainpictures.net/en/index.php
+* https://www.naturalearthdata.com/

examples/doc_bonneP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth (bonneP (45/180*pi))

examples/doc_bottomleyP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth (bottomleyP (30/180*pi))

examples/doc_collignonP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth collignonP

examples/doc_eckert1P.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth eckert1P

examples/doc_eckert3P.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth eckert3P

examples/doc_eckert5P.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth eckert5P

examples/doc_equirectangularP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth equirectangularP

examples/doc_faheyP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth faheyP

examples/doc_hammerP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth hammerP

examples/doc_lambertP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth lambertP

examples/doc_mercatorP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth mercatorP

examples/doc_mollweideP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth mollweideP

examples/doc_orthoP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth (orthoP 0 0)

examples/doc_sinusoidalP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth sinusoidalP

examples/doc_smallEarth.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.Builtin.Documentation
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ docEnv $ animate $ const $
+  scaleToSize screenWidth screenHeight $
+  embedImage smallEarth

examples/doc_wernerP.hs

@@ -0,0 +1,12 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main (main) where
+
+import           Reanimate
+import           Reanimate.GeoProjection
+import           Reanimate.Builtin.Images
+
+main :: IO ()
+main = reanimate $ animate $
+  const $ scaleToSize screenWidth screenHeight $
+  embedImage $ project smallEarth wernerP

reanimate.cabal

@@ -30,6 +30,7 @@ data-files:           viewer-elm/dist/index.html
                       data/CIE_XYZ.csv
                       data/cone_sensitivity_lms.csv
                       data/*.svg
+                      data/*.jpg
 
 Source-Repository head
     Type:      git
@@ -69,6 +70,7 @@ library
                       Reanimate.Parameters
                       Reanimate.Chiphunk
                       Reanimate.PolyShape
+                      Reanimate.GeoProjection
                       Reanimate.Builtin.Documentation
                       Reanimate.Builtin.Images
   other-modules:      Reanimate.Cache

src/Reanimate/Builtin/CirclePlot.hs

@@ -7,9 +7,9 @@ import           Reanimate.Svg
 import           Reanimate.Constants
 
 circlePlot :: Int -- ^ Pixels in the X-axis.
-           -> (Double -- ^ Angle in radians
-            -> Double -- ^ Radius in percent
-            -> PixelRGBA8) -> Tree
+           -> (Double -> Double -> PixelRGBA8)
+              -- ^ Angle and radius in radians percent respectively.
+           -> Tree
 circlePlot density fn =
     scaleToHeight screenHeight $ flipYAxis $
     embedImage $ generateImage gen density density

src/Reanimate/Builtin/Images.hs

@@ -2,8 +2,10 @@ module Reanimate.Builtin.Images
   ( svgLogo
   , haskellLogo
   , githubIcon
+  , smallEarth
   ) where
 
+import           Codec.Picture
 import qualified Data.ByteString     as B
 import           Graphics.SvgTree    (parseSvgFile)
 import           Paths_reanimate
@@ -19,6 +21,14 @@ embedImage key = do
     Nothing  -> error "Malformed svg"
     Just svg -> return $ embedDocument svg
 
+loadJPG :: FilePath -> Image PixelRGB8
+loadJPG key = unsafePerformIO $ do
+  jpg_file <- getDataFileName key
+  dat <- B.readFile jpg_file
+  case decodeJpeg dat of
+    Left err  -> error err
+    Right img -> return $ convertRGB8 img
+
 -- | <<docs/gifs/doc_svgLogo.gif>>
 svgLogo :: SVG
 svgLogo = unsafePerformIO $ embedImage "data/svg-logo.svg"
@@ -30,3 +40,9 @@ haskellLogo = unsafePerformIO $ embedImage "data/haskell.svg"
 -- | <<docs/gifs/doc_githubIcon.gif>>
 githubIcon :: SVG
 githubIcon = unsafePerformIO $ embedImage "data/github-icon.svg"
+
+-- | 300x150 equirectangular earth
+--
+--   <<docs/gifs/doc_smallEarth.gif>>
+smallEarth :: Image PixelRGB8
+smallEarth = loadJPG "data/small_earth.jpg"

src/Reanimate/GeoProjection.hs

@@ -0,0 +1,671 @@
+{-# LANGUAGE MultiWayIf #-}
+module Reanimate.GeoProjection
+  ( Projection(..)
+  , XYCoord(..)
+  , LonLat(..)
+  , project
+  , interpP
+  , mergeP
+  , isValidP
+  , scaleP
+  , flipYAxisP
+  , moveBottomP
+  , moveTopP
+    -- * Projections
+  , equirectangularP
+  , mercatorP
+  , mollweideP
+  , hammerP
+  , lambertP
+  , bottomleyP
+  , sinusoidalP
+  , wernerP
+  , bonneP
+  , orthoP
+  , cassiniP
+  , augustP
+  , collignonP
+  , eckert1P
+  , eckert3P
+  , eckert5P
+  , faheyP
+  , foucautP
+  , lagrangeP
+  ) where
+
+import           Codec.Picture
+import           Codec.Picture.Types
+import           Control.Monad
+import           Control.Monad.ST
+-- import           Data.List
+import           Data.Maybe
+import           Debug.Trace
+import           Reanimate
+
+-- Constants
+halfPi, sqrtPi, epsilon, tau :: Double
+halfPi = pi/2
+sqrtPi = sqrt pi
+epsilon = 1.0e-12
+tau = pi*2
+
+toRads, cot :: Double -> Double
+toRads dec = dec/180 * pi
+cot = recip . tan
+
+srcPixel :: Image PixelRGB8 -> LonLat -> PixelRGB8
+srcPixel src (LonLat lam phi) =
+    pixelAt src xPx yPx
+  where
+    xPx = round $ ((lam+pi)/tau) * fromIntegral (imageWidth src-1)
+    yPx = round $ (1-((phi+halfPi)/pi)) * fromIntegral (imageHeight src-1)
+
+findValidCoord :: Image PixelRGB8 -> Projection -> XYCoord -> XYCoord
+findValidCoord src p (XYCoord x y) = fromMaybe (XYCoord x y) $ listToMaybe
+    [ XYCoord x' y'
+    | let xi = round $ x * wMax
+          yi = round $ y * hMax
+    , ax <- [xi, xi-1, xi+1]
+    , ay <- [yi, yi-1, yi+1]
+    , (ax >= 0 && ax < w && ay >= 0 && ay < h)
+    , let x' = fromIntegral ax / wMax
+          y' = fromIntegral ay / hMax
+    , validLonLat $ projectionInverse p (XYCoord x' y')
+    ]
+  where
+    w = imageWidth src
+    h = imageHeight src
+    wMax = fromIntegral (w-1)
+    hMax = fromIntegral (h-1)
+
+isInWorld :: Projection -> XYCoord -> Bool
+isInWorld p coord =
+    odd $ length $ isInWorld' p coord
+
+isInWorld' :: Projection -> XYCoord -> [(Double, Double)]
+isInWorld' p (XYCoord x y) =
+    [ (x1, y1)
+    | (XYCoord x1 y1, XYCoord x2 y2) <- world
+    , (y1 > y) /= (y2 > y) -- y is between y1 and y2
+    , x < (x2 - x1) * (y - y1) / (y2 - y1) + x1 -- x is to the left of the line
+    ]
+  where
+    world = worldPolygon p
+
+worldPolygon :: Projection -> [(XYCoord, XYCoord)]
+worldPolygon p =
+    interp (-pi, -halfPi) (-pi, halfPi) ++
+    interp (-pi, halfPi) (pi, halfPi) ++
+    interp (pi, halfPi) (pi, -halfPi) ++
+    interp (pi, -halfPi) (-pi, -halfPi)
+  where
+    apply (lam, phi) = projectionForward p $ LonLat lam phi
+    steps = 100
+    interp (x1, y1) (x2,y2) =
+      [ ( apply (fromToS x1 x2 (n/steps), fromToS y1 y2 (n/steps))
+        , apply (fromToS x1 x2 ((n+1)/steps), fromToS y1 y2 ((n+1)/steps)))
+      | n <- [0..steps-1]]
+
+findNearestPixel :: MutableImage s PixelRGBA8 -> Int -> Int -> Int -> Int -> ST s PixelRGBA8
+findNearestPixel src w h srcX srcY = worker
+    [ (x, y)
+    | n <- [1..]
+    , x <- [srcX-n .. srcY+n]
+    , y <- if x == srcX-n || x == srcX+n then [srcY-n,srcY+n] else [srcY-n .. srcY+n]
+    , x >= 0
+    , y >= 0
+    , x < w
+    , y < h
+    ]
+  where
+    worker [] = undefined
+    worker ((x,y):rest) = do
+      this <- readPixel src x y
+      if this == blank
+        then worker rest
+        else return this
+    blank = PixelRGBA8 0x00 0x00 0x00 0x00
+
+interpP :: Image PixelRGB8 -> Projection -> Projection -> Double -> Image PixelRGBA8
+interpP src p1 p2 t = runST $ do
+    img <- newMutableImage w h
+    let blank = PixelRGBA8 0x00 0x00 0x00 0x00
+    let isBlank pixel = pixel == blank
+    -- forM_ [0..w-1] $ \x ->
+    --   forM_ [0..h-1] $ \y -> do
+    --     let x1 = fromIntegral x / (wMax)
+    --         y1 = 1 - fromIntegral y / (hMax)
+    --     when (isInWorld (mergeP p1 p2 t) (XYCoord x1 y1)) $
+    --       writePixel img x y $ PixelRGBA8 0xFF 0x00 0x00 0xFF
+    let factor = 2
+    forM_ [0..(w*factor)-1] $ \x ->
+      forM_ [0..(h*factor)-1] $ \y -> do
+        let x1' = fromIntegral x / (wMax*fromIntegral factor)
+            y1' = fromIntegral y / (hMax*fromIntegral factor)
+            lonlat = projectionInverse p1 (XYCoord x1' y1')
+            XYCoord x1 y1 = projectionForward p1 lonlat
+            XYCoord x2 y2 = findValidCoord src p2 $ projectionForward p2 lonlat
+            x3 = round $ fromToS x1 x2 t * wMax
+            y3 = round $ (1 - fromToS y1 y2 t) * hMax
+        when (validLonLat lonlat && validXYCoord (XYCoord x2 y2)) $ do
+          when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $ do
+            writePixel img x3 y3 (promotePixel $ srcPixel src lonlat)
+    forM_ [0..(w*factor)-1] $ \x ->
+      forM_ [0..(h*factor)-1] $ \y -> do
+        let x2' = fromIntegral x / (wMax*fromIntegral factor)
+            y2' = fromIntegral y / (hMax*fromIntegral factor)
+            lonlat = projectionInverse p2 (XYCoord x2' y2')
+            XYCoord x2 y2 = projectionForward p2 lonlat
+            XYCoord x1 y1 = findValidCoord src p1 $ projectionForward p1 lonlat
+            -- (x2,y2) = p2 lam phi
+            x3 = round $ fromToS x1 x2 t * wMax
+            y3 = round $ (1 - fromToS y1 y2 t) * hMax
+        when (validLonLat lonlat && validXYCoord (XYCoord x1 y1)) $
+          when (validLonLat lonlat) $
+            when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
+              writePixel img x3 y3 (promotePixel $ srcPixel src lonlat)
+    forM_ [1..w-1] $ \x ->
+      forM_ [0..h-1] $ \y -> do
+        let x1 = fromIntegral x / (wMax)
+            y1 = 1 - fromIntegral y / (hMax)
+        this <- readPixel img x y
+        when (isBlank this) $
+          when (isInWorld (mergeP p1 p2 t) (XYCoord x1 y1)) $
+            writePixel img x y =<< findNearestPixel img w h x y
+    unsafeFreezeImage img
+  where
+    w = imageWidth src
+    h = imageHeight src
+    wMax = fromIntegral (w-1)
+    hMax = fromIntegral (h-1)
+
+eqLonLat :: LonLat -> LonLat -> Bool
+eqLonLat (LonLat x1 y1) (LonLat x2 y2)
+  = eqDouble x1 x2 && eqDouble y1 y2
+{-
+eqCoords :: XYCoord -> XYCoord -> Bool
+eqCoords (XYCoord x1 y1) (XYCoord x2 y2)
+  = eqDouble x1 x2 && eqDouble y1 y2
+-}
+eqDouble :: Double -> Double -> Bool
+eqDouble a b = abs (a-b) < epsilon
+
+data XYCoord = XYCoord Double Double -- 0 to 1
+  deriving (Read,Show,Eq,Ord)
+data LonLat = LonLat Double Double -- -pi to +pi, -halfPi to +halfPi
+  deriving (Read,Show,Eq,Ord)
+data Projection = Projection
+  { projectionForward :: LonLat -> XYCoord
+  , projectionInverse :: XYCoord -> LonLat
+  }
+
+-- FIXME: Verify that 'src' has an aspect ratio of 2:1.
+project :: Image PixelRGB8 -> Projection -> Image PixelRGBA8
+project src (Projection _ pInv) = generateImage fn w h
+  where
+    w = imageWidth src
+    h = imageHeight src
+    fn xPx yPx =
+      let x = (fromIntegral xPx / fromIntegral (w-1))
+          y = 1-(fromIntegral yPx / fromIntegral (h-1))
+          lonlat = pInv (XYCoord x y)
+      in
+      if validLonLat lonlat
+        then promotePixel (srcPixel src lonlat)
+        else PixelRGBA8 0 0 0 0
+
+validLonLat :: LonLat -> Bool
+validLonLat (LonLat lam phi) =
+  lam >= -pi && lam <= pi && phi >= -pi/2 && phi <= pi/2
+
+validXYCoord :: XYCoord -> Bool
+validXYCoord (XYCoord x y) = x >= 0 && x <= 1 && y >= 0 && y <= 1
+
+isValidP :: Projection -> Bool
+isValidP (Projection p pInv) = and
+    [ check x y
+    | x <- [0..w-1::Int]
+    , y <- [0..h-1::Int] ]
+  where
+    w = 100
+    h = 100
+    check xPx yPx =
+      let x = (fromIntegral xPx / fromIntegral (w-1))
+          y = (fromIntegral yPx / fromIntegral (h-1))
+          lonlat = pInv (XYCoord x y)
+          lonlat2 = pInv $ p lonlat
+      in if not (validLonLat lonlat) || eqLonLat lonlat lonlat2
+          then True
+          else trace (show (lonlat, lonlat2)) $ False
+
+moveBottomP :: Double -> Projection -> Projection
+moveBottomP offset (Projection p pInv) = Projection p' pInv'
+  where
+    p' (LonLat lon lat) =
+      case p (LonLat lon lat) of
+        XYCoord x y -> XYCoord x (fromToS offset 1 y)
+    pInv' (XYCoord x y) = pInv (XYCoord x ((y-offset)/(1-offset)))
+
+moveTopP :: Double -> Projection -> Projection
+moveTopP offset = flipYAxisP . moveBottomP offset . flipYAxisP
+
+flipYAxisP :: Projection -> Projection
+flipYAxisP (Projection p pInv) = Projection p' pInv'
+  where
+    p' (LonLat lam phi) =
+      let XYCoord x y = p (LonLat lam (negate phi))
+      in XYCoord x (1-y)
+    pInv' (XYCoord x y) =
+      let LonLat lam phi = pInv (XYCoord x (1-y))
+      in LonLat lam (negate phi)
+
+scaleP :: Double -> Double -> Projection -> Projection
+scaleP xScale yScale (Projection p pInv) = Projection forward inverse
+  where
+    forward lonlat =
+      case p lonlat of
+        XYCoord x y -> XYCoord ((x-0.5)*xScale+0.5) ((y-0.5)*yScale+0.5)
+    inverse (XYCoord x y) =
+      let new = XYCoord ((x-0.5)/xScale+0.5) ((y-0.5)/yScale+0.5)
+      in pInv new
+
+
+mergeP :: Projection -> Projection -> Double -> Projection
+mergeP p1 p2 t = Projection p pInv
+  where
+    p lonlat =
+      let XYCoord x1 y1 = projectionForward p1 lonlat
+          XYCoord x2 y2 = projectionForward p2 lonlat
+      in XYCoord (fromToS x1 x2 t) (fromToS y1 y2 t)
+    pInv coord =
+      let LonLat lon1 lat1 = projectionInverse p1 coord
+          LonLat lon2 lat2 = projectionInverse p2 coord
+      in
+        if | oob lon1 lat1 && oob lon2 lat2 -> LonLat (0/0) (0/0)
+           | otherwise     -> LonLat (fromToS lon1 lon2 t) (fromToS lat1 lat2 t)
+    oob lon lat = lon < (-pi) || lon > pi || lat < (-pi/2) || lat > pi/2
+
+-- | <<docs/gifs/doc_equirectangularP.gif>>
+equirectangularP :: Projection
+equirectangularP = Projection forward inverse
+  where
+    forward (LonLat lam phi) = XYCoord ((lam+pi)/tau) ((phi+pi/2)/pi)
+    inverse (XYCoord x y) = LonLat xPi yPi
+      where
+        xPi = fromToS (-pi) pi x
+        yPi = fromToS (-pi/2) (pi/2) y
+
+-- | <<docs/gifs/doc_mercatorP.gif>>
+mercatorP :: Projection
+mercatorP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+      XYCoord ((lam+pi)/tau)
+        (((log(tan(pi/4+phi/2))) + pi)/tau)
+    inverse (XYCoord x y) = LonLat xPi (atan (sinh yPi))
+      where
+        xPi = fromToS (-pi) pi x
+        yPi = fromToS (-pi) pi y
+
+-- | <<docs/gifs/doc_mollweideP.gif>>
+mollweideP :: Projection
+mollweideP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+        XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
+      where
+        sqrt2 = sqrt 2
+        x = (2*sqrt2)/pi * lam * cos theta
+        y = sqrt2*sin theta
+        theta = find_theta 100
+        find_theta :: Int -> Double
+        find_theta 0 = phi
+        find_theta _ | abs phi == pi/2 = signum phi * pi/2
+        find_theta n =
+          let sub = find_theta (n-1)
+          in sub - (2*sub+sin (2*sub)-pi*sin phi)/(2+2*cos(2*sub))
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-2*sqrt(2)) (2*sqrt(2)) x'
+        y = fromToS (-sqrt 2) (sqrt 2) y'
+        theta = asin (y/sqrt(2))
+        lam = pi*x/(2*sqrt(2)*cos theta)
+        phi = asin ((2*theta+sin(2*theta))/pi)
+
+-- | <<docs/gifs/doc_hammerP.gif>>
+hammerP :: Projection
+hammerP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+        XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
+      where
+        sqrt2 = sqrt 2
+        x = (2*sqrt2*cos phi*sin (lam/2))/(sqrt (1+cos phi*cos (lam/2)))
+        y = (sqrt2*sin phi)/(sqrt (1+cos phi*cos (lam/2)))
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-2*sqrt(2)) (2*sqrt(2)) x'
+        y = fromToS (-sqrt 2) (sqrt 2) y'
+        z = sqrt (1 - (x/4)**2 - (y/2)**2)
+        lam = 2 * atan2 (z*x) (2*(2*z**2-1))
+        phi = asin (z*y)
+
+-- | <<docs/gifs/doc_lambertP.gif>>
+lambertP :: Projection
+lambertP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+      XYCoord ((lam+pi)/tau) ((sin phi+1)/2)
+    inverse (XYCoord x' y') = LonLat x (asin y)
+      where
+        x = fromToS (-pi) (pi) x'
+        y = fromToS (-1) 1 y'
+
+-- | <<docs/gifs/doc_bottomleyP.gif>>
+bottomleyP :: Double -> Projection
+bottomleyP phi_1 = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+        XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
+      where
+        x = (rho * sin e) / sin phi_1
+        y = pi/2 - rho * cos e
+        rho = pi/2 - phi
+        e = lam * sin phi_1 * sin rho / rho
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-pi) (pi) x'
+        y = fromToS (-pi/2) (pi/2) y'
+        x1 = x * sin phi_1
+        y1 = pi/2 - y
+        rho = sqrt (x1*x1 + y1*y1)
+        e = atan2 x1 y1
+        lam = (if rho == 0 then 1 else rho / sin rho) * e/sin phi_1
+        phi = pi/2 - rho
+
+-- | <<docs/gifs/doc_sinusoidalP.gif>>
+sinusoidalP :: Projection
+sinusoidalP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+        XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
+      where
+        x = lam * cos phi
+        y = phi
+    inverse (XYCoord x' y') = LonLat (x/cos y) y
+      where
+        x = fromToS (-pi) (pi) x'
+        y = fromToS (-pi/2) (pi/2) y'
+
+-- | <<docs/gifs/doc_wernerP.gif>>
+wernerP :: Projection
+wernerP = moveTopP 0.23 $ Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+        XYCoord ((x+pi)/tau) ((y+pi/2)/pi)
+      where
+        rho = pi/2 - phi
+        e = lam * sin rho / rho
+        x = rho * sin e
+        y = pi/2 - rho * cos e
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-pi) (pi) x'
+        y = fromToS (-pi/2) (pi/2) y'
+        rho = sqrt (x**2 + (pi/2 - y)**2)
+        e = atan2 x (pi/2 -y)
+        phi = pi/2 - rho
+        lam = e * rho / sin rho
+
+-- FIXME: find right scale and position.
+-- | <<docs/gifs/doc_bonneP.gif>>
+bonneP :: Double -> Projection
+bonneP 0 = sinusoidalP
+bonneP phi_0 = moveTopP (-0.17*factor) $ scaleP 1 (fromToS 1 0.65 factor) $ Projection forward inverse
+  where
+    factor = sin phi_0 / sin (pi/4)
+    forward (LonLat lam phi ) = XYCoord ((x+pi)/tau) ((y+halfPi)/pi)
+      where
+        cotPhi0 = cot phi_0
+        rho = cotPhi0 + phi_0 - phi
+        e = (lam * cos phi) / rho
+        x = rho * sin e
+        y = cotPhi0 - rho * cos e
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-pi) (pi) x'
+        y = fromToS (-pi/2) (pi/2) y'
+        cotPhi0 = cot phi_0
+        rho = sqrt (x*x + (cot phi_0 - y)**2)
+        -- e = atan2 x (cot phi_1 + phi_1 + y)
+        phi = cotPhi0 + phi_0 - rho
+        lam = rho / cos phi * atan2 x (cotPhi0-y)
+
+-- | <<docs/gifs/doc_orthoP.gif>>
+orthoP :: Double -> Double -> Projection
+orthoP lam_0 phi_0 = Projection forward inverse
+  where
+    forward (LonLat lam phi)
+      | (lam+lam_0) < -halfPi || (lam+lam_0) > halfPi ||
+        (phi+phi_0) < -halfPi/2 || (phi+phi_0) > halfPi/2
+        = XYCoord (0/0) (0/0)
+    forward (LonLat lam phi) =
+        XYCoord ((x+(16/9))/(16/9*2)) ((y+1)/2)
+      where
+        x = cos phi * sin (lam - lam_0)
+        y = cos phi_0 * sin phi - sin phi_0 * cos phi * cos (lam - lam_0)
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-16/9) (16/9) x'
+        y = fromToS (-1) (1) y'
+        lam = wrap (-pi) (pi) $
+          lam_0 + atan2 (x * sin c) (rho * cos c * cos phi_0 - y * sin c * sin phi_0)
+        phi = wrap (-pi/2) (pi/2) $
+          asin ((cos c * sin phi_0 + y * sin c * cos phi_0)/rho)
+        rho = sqrt (x**2 + y**2)
+        c = asin rho
+        wrap lower upper v
+          | v > upper = v-upper+lower
+          | v < lower = v+upper-lower
+          | otherwise = v
+
+cassiniP :: Projection
+cassiniP = Projection forward inverse
+  where
+    forward (LonLat lam phi) =
+      XYCoord ((asin (cos phi * sin lam)+halfPi)/pi) ((atan2 (tan phi) (cos lam)+pi)/tau)
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS (-halfPi) halfPi x'
+        y = fromToS (-pi) (pi) y'
+        lam = atan2 (tan x) (cos y)
+        phi = asin (sin y * cos x)
+
+augustP :: Projection
+augustP = scaleP 0.70 0.70 $ Projection forward inverse
+  where
+    xHi = 16/3
+    xLo = -xHi
+    yHi = 8 / 3
+    yLo = -yHi
+    forward (LonLat lam phi) = XYCoord ((xPos-xLo)/(xHi-xLo)) ((yPos-yLo)/(yHi-yLo))
+      where
+        tanPhi = tan (phi/2)
+        k = sqrt (1 - tanPhi * tanPhi)
+        c = 1 + k * cos (lam / 2)
+        x = sin (lam/2) * k / c
+        y = tanPhi / c
+        x2 = x*x
+        y2 = y*y
+        xPos = (4 / 3 * x * (3+x2 - 3*y2))
+        yPos = (4 / 3 * y * (3 + 3*x2 - y2))
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x' * 3 / 8
+        y = fromToS yLo yHi y' * 3 / 8
+        x2 = x*x
+        y2 = y*y
+        s = 1 + x2 + y2
+        sin3Eta = sqrt ((s - sqrt (s*s - 4 * y * y)) / 2)
+        eta = asin (sin3Eta) / 3
+        xi = if sin3Eta /= 0 then acosh (abs (y / sin3Eta)) / 3 else asinh (abs x) / 3
+        cosEta = cos eta
+        coshXi = cosh xi
+        d = coshXi * coshXi - cosEta * cosEta
+        lam = signum x * 2 * atan2 (sinh xi * cosEta) (0.25 - d)
+        phi = signum y * 2 * atan2 (coshXi * sin eta) (0.25 + d)
+
+-- | <<docs/gifs/doc_collignonP.gif>>
+collignonP :: Projection
+collignonP = Projection forward inverse
+  where
+    yHi = sqrtPi
+    yLo = sqrtPi * (1 - sqrt 2)
+    xLo = -pi*(2/sqrtPi)*(sqrt 2)
+    xHi = pi*(2/sqrtPi)*(sqrt 2)
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        alpha = sqrt (1 - sin phi)
+        x = (2 / sqrtPi) * lam * alpha
+        y = sqrtPi * (1 - alpha)
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        l = (y / sqrtPi -1)**2
+        lam = if l > 0 then x * sqrt (pi / l) / 2 else 0
+        phi = asin (1-l)
+
+-- | <<docs/gifs/doc_eckert1P.gif>>
+eckert1P :: Projection
+eckert1P = Projection forward inverse
+  where
+    alpha = sqrt (8 / (3*pi))
+    yLo = -alpha * halfPi
+    yHi = alpha * halfPi
+    xLo = -alpha * pi
+    xHi = alpha * pi
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        x = alpha * lam * (1 - abs phi / pi)
+        y = alpha * phi
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        phi = y / alpha
+        lam = x / (alpha * (1 - abs phi / pi))
+
+-- | <<docs/gifs/doc_eckert3P.gif>>
+eckert3P :: Projection
+eckert3P = Projection forward inverse
+  where
+    k = sqrt (pi * (4 + pi))
+    yLo = negate yHi
+    yHi = 4/k * halfPi
+    xLo = negate xHi
+    xHi = 4/k * pi
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        x = 2 / k * lam * (1 + sqrt (1 - 4*phi*phi/(pi*pi)))
+        y = 4 / k * phi
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        lam = x  * k/2 / (1 + sqrt (1 - y*y* (4+pi)/(4*pi)))
+        phi = y * k/4
+
+-- | <<docs/gifs/doc_eckert5P.gif>>
+eckert5P :: Projection
+eckert5P = Projection forward inverse
+  where
+    k = sqrt (2 + pi)
+    yLo = negate yHi
+    yHi = pi/k
+    xLo = negate xHi
+    xHi = tau/k
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        x = lam * (1 + cos phi) / k
+        y = 2 * phi / k
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        lam = k * x / (1 + cos phi)
+        phi = y * k / 2
+
+-- | <<docs/gifs/doc_faheyP.gif>>
+faheyP :: Projection
+faheyP = Projection forward inverse
+  where
+    faheyK = cos (toRads 35)
+    yLo = negate yHi
+    yHi = 1 + faheyK
+    xLo = negate xHi
+    xHi = pi * faheyK * 16/9
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        t = tan (phi/2)
+        x = lam * faheyK * sqrt (1 - t*t)
+        y = (1 + faheyK) * t
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        t = y / (1 + faheyK)
+        lam = x / (faheyK * sqrt (1 - t*t))
+        phi = 2 * atan2 y (1 + faheyK)
+
+foucautP :: Projection
+foucautP = Projection forward inverse
+  where
+    yLo = negate yHi
+    yHi = sqrtPi * tan (halfPi/2)
+    xLo = negate xHi
+    xHi = tau/sqrtPi
+    forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        k = phi / 2
+        cosk = cos k
+        x = 2 * lam / sqrtPi * cos phi * cosk * cosk
+        y = sqrtPi * tan k
+    inverse (XYCoord x' y') = LonLat lam phi
+      where
+        x = fromToS xLo xHi x'
+        y = fromToS yLo yHi y'
+        k = atan (y / sqrtPi)
+        cosk = cos k
+        phi = 2 * k
+        lam = x * sqrtPi / 2 / (cos phi * cosk * cosk)
+
+lagrangeP :: Projection
+lagrangeP = Projection forward inverse
+  where
+    yLo = negate yHi
+    yHi = 2
+    xLo = negate xHi
+    xHi = 2
+    n = 0.5
+    forward (LonLat lam phi)
+        | abs (abs phi - halfPi) < epsilon = XYCoord 0.5 (if phi < 0 then 0 else 1)
+        | otherwise = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
+      where
+        sinPhi = sin phi
+        v = ((1+sinPhi) / (1 - sinPhi))**(n/2)
+        c = 0.5 * (v + 1/v) + cos (lam*n)
+        x = 2 * sin (lam*n) / c
+        y = (v - 1/v) /c
+    inverse (XYCoord x' y')
+        | abs ((abs y')-1) < epsilon = LonLat 0 (signum y * halfPi)
+        | otherwise = LonLat lam phi
+      where
+        x = fromToS xLo xHi x' / 2
+        y = fromToS yLo yHi y' / 2
+        x2 = x * x
+        y2 = y * y
+        t = 2 * y / (1 + x2 + y2)
+        t' = ((1+t) / (1-t)) ** (1/n)
+        lam = atan2 (2*x) (1-x2-y2) / n
+        phi = asin ((t'-1)/(t'+1))

src/Reanimate/LaTeX.hs

@@ -87,6 +87,7 @@ mkTexScript latexExec latexArgs texHeaders tex = T.unlines $
   [ "% " <> T.pack (unwords (latexExec:latexArgs))
   , "\\documentclass[preview]{standalone}"
   , "\\usepackage{amsmath}"
+  , "\\usepackage{gensymb}"
   ] ++ texHeaders ++
   [ "\\usepackage[english]{babel}"
   , "\\linespread{1}"

src/Reanimate/Raster.hs

@@ -6,6 +6,7 @@ module Reanimate.Raster
   , vectorize
   , vectorize_
   , svgAsPngFile
+  , svgAsPngFile'
   ) where
 
 import           Codec.Picture.Types
@@ -128,8 +129,14 @@ vectorize_ args path = unsafePerformIO $ do
 -- imageAsFile img
 
 svgAsPngFile :: Tree -> FilePath
-svgAsPngFile _ | pNoExternals = "/svgAsPngFile/has/been/disabled"
-svgAsPngFile svg = unsafePerformIO $ cacheFile template $ \pngPath -> do
+svgAsPngFile = svgAsPngFile' width height
+  where
+    width = 2560
+    height = width * 9 `div` 16
+
+svgAsPngFile' :: Int -> Int -> Tree -> FilePath
+svgAsPngFile' _ _ _ | pNoExternals = "/svgAsPngFile/has/been/disabled"
+svgAsPngFile' width height svg = unsafePerformIO $ cacheFile template $ \pngPath -> do
     let svgPath = replaceExtension pngPath "svg"
     -- ffmpeg <- requireExecutable "ffmpeg"
     -- convert <- requireExecutable "convert"
@@ -139,6 +146,4 @@ svgAsPngFile svg = unsafePerformIO $ cacheFile template $ \pngPath -> do
     runCmd inkscape [ svgPath, "--export-png=" ++ pngPath, "--without-gui" ]
   where
     template = show (hash rendered) <.> "png"
-    rendered = renderSvg (Just $ Num width) (Just $ Num height) svg
-    width = 2560
-    height = width * 9 / 16
+    rendered = renderSvg (Just $ Px $ fromIntegral width) (Just $ Px $ fromIntegral height) svg

src/Reanimate/Render.hs

@@ -12,11 +12,11 @@ import           Control.Monad       (forM_)
 import qualified Data.Text           as T
 import qualified Data.Text.IO        as T
 import           Graphics.SvgTree    (Number (..))
+import           Numeric
 import           Reanimate.Animation
 import           Reanimate.Misc
-import           System.FilePath     ((</>))
--- import           System.FilePath    (replaceExtension)
 import           System.Exit
+import           System.FilePath     ((</>))
 import           System.IO
 import           Text.Printf         (printf)
 
@@ -99,14 +99,14 @@ render ani target format width height fps = do
       RenderGif -> withTempFile "png" $ \palette -> do
         runCmd ffmpeg ["-i", template, "-y"
                       ,"-vf", "fps="++show fps++",scale=320:-1:flags=lanczos,palettegen"
-                      ,"-t", show (duration ani)
+                      ,"-t", showFFloat Nothing (duration ani) ""
                       , palette ]
         runCmd ffmpeg ["-i", template, "-y"
                       ,"-i", palette
                       ,"-progress", progress
                       ,"-filter_complex"
                       ,"fps="++show fps++",scale=320:-1:flags=lanczos[x];[x][1:v]paletteuse"
-                      ,"-t", show (duration ani)
+                      ,"-t", showFFloat Nothing (duration ani) ""
                       , target]
       RenderWebm ->
         runCmd ffmpeg ["-r", show fps, "-i", template, "-y"

src/Reanimate/Svg/BoundingBox.hs

@@ -29,12 +29,12 @@ boundingBox t =
 svgHeight :: Tree -> Double
 svgHeight t = h
   where
-    (_w, h, _x, _y) = boundingBox t
+    (_x, _y, _w, h) = boundingBox t
 
 svgWidth :: Tree -> Double
 svgWidth t = w
   where
-    (w, _h, _x, _y) = boundingBox t
+    (_x, _y, w, _h) = boundingBox t
 
 linePoints :: [LineCommand] -> [RPoint]
 linePoints = worker zero
@@ -87,22 +87,22 @@ svgBoundingPoints t = map (Transform.transformPoint m) $
         (Num x, Num y) -> V2 x y
         _ -> error "Reanimate.Svg.svgBoundingPoints: Unrecognized number format."
 
-    circleBoundingPoints circ = 
+    circleBoundingPoints circ =
       let (xnum, ynum) = circ ^. circleCenter
           rnum = circ ^. circleRadius
-      in case mapMaybe unpackNumber [xnum, ynum, rnum] of 
+      in case mapMaybe unpackNumber [xnum, ynum, rnum] of
         [x, y, r] -> [ V2 (x + r * cos angle) (y + r * sin angle) | angle <- [0, pi/10 .. 2 * pi]]
         _  -> []
 
-    ellipseBoundingPoints e = 
+    ellipseBoundingPoints e =
       let (xnum,ynum) = e ^. ellipseCenter
           xrnum = e ^. ellipseXRadius
           yrnum = e ^. ellipseYRadius
-      in case mapMaybe unpackNumber [xnum, ynum, xrnum, yrnum] of 
+      in case mapMaybe unpackNumber [xnum, ynum, xrnum, yrnum] of
         [x,y,xr,yr] -> [V2 (x + xr * cos angle) (y + yr * sin angle) | angle <- [0, pi/10 .. 2 * pi]]
         _ -> []
 
     unpackNumber n =
       case toUserUnit defaultDPI n of
         Num d -> Just d
-        _     -> Nothing
+        _     -> Nothing

src/Reanimate/Svg/Constructors.hs

@@ -42,13 +42,15 @@ module Reanimate.Svg.Constructors
   , flipXAxis
   , flipYAxis
   , aroundCenter
+  , aroundCenterX
+  , aroundCenterY
   , withTransformations
   -- * Other
   , mkColor
   , mkBackground
   , mkBackgroundPixel
   , gridLayout
-  
+
   ) where
 
 import           Codec.Picture                (PixelRGBA8 (..))
@@ -98,6 +100,20 @@ aroundCenter fn t =
     offsetY = -y-h/2
     (x,y,w,h) = boundingBox t
 
+aroundCenterY :: (Tree -> Tree) -> Tree -> Tree
+aroundCenterY fn t =
+    translate 0 (-offsetY) $ fn $ translate 0 offsetY t
+  where
+    offsetY = -y-h/2
+    (_x,y,_w,h) = boundingBox t
+
+aroundCenterX :: (Tree -> Tree) -> Tree -> Tree
+aroundCenterX fn t =
+    translate (-offsetX) 0 $ fn $ translate offsetX 0 t
+  where
+    offsetX = -x-w/2
+    (x,_y,w,_h) = boundingBox t
+
 -- | Scale the image uniformly by given factor along both X and Y axes.
 -- For example @scale 2 image@  makes the image twice as large, while @scale 0.5 image@ makes it half the original size.
 -- Negative values are also allowed, and lead to flipping the image along both X and Y axes.

videos/map-projection/distortion.hs

@@ -0,0 +1,228 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+{-# LANGUAGE ApplicativeDo     #-}
+{-# LANGUAGE OverloadedStrings #-}
+module Main(main) where
+
+import           Codec.Picture
+import           Codec.Picture.Jpg
+import           Codec.Picture.Types
+import           Control.Lens            ((^.))
+import           Control.Monad
+import           Control.Monad.ST
+import           Data.Aeson
+import qualified Data.ByteString         as BS
+import           Data.Foldable
+import           Data.Geospatial         hiding (LonLat)
+import           Data.LinearRing
+import qualified Data.LineString         as Line
+import           Data.Map                (Map)
+import qualified Data.Map                as Map
+import qualified Data.Text               as T
+import           Graphics.SvgTree        (PathCommand (..), Tree (None))
+import           Reanimate
+import           Reanimate.Animation
+import           Reanimate.GeoProjection
+import           Reanimate.Scene
+import           System.IO.Unsafe
+
+
+main :: IO ()
+main = seq equirectangular $ reanimate $ sceneAnimation $ do
+    prevProj <- newVar equirectangularP
+    txtVar <- newVar "Equirectangular"
+    txtS <- newSprite $ renderLabel <$> unVar txtVar
+    txtFade <- spriteVar txtS 1 withGroupOpacity
+    spriteZ txtS 2
+    let pushInterp = pushMerge' (\a b t -> interpP src a b t)
+        pushMerge' fn label proj = do
+          fork $ do
+            tweenVar txtFade 0.2 $ \v -> fromToS v 0 . curveS 2
+            writeVar txtVar label
+            tweenVar txtFade 0.2 $ \v -> fromToS v 1 . curveS 2
+          prev <- readVar prevProj
+          play $ pauseAtEnd waitT $ signalA (curveS 2) $
+            mkAnimation morphT $ \t ->
+              mkGroup $
+              [ scaleToSize screenWidth screenHeight $
+                embedImage $ fn prev proj t
+              , grid $ mergeP prev proj t ]
+          writeVar prevProj proj
+        pushT mkLabel mkProj = do
+          fork $ tweenVar txtVar morphT $ \v t -> if t > 0 then mkLabel t else v
+          play $ pauseAtEnd waitT $ signalA (curveS 2) $
+            mkAnimation morphT $ \t ->
+              mkGroup $
+              [ scaleToSize screenWidth screenHeight $
+                embedImage $ project src $ mkProj t
+              , grid $ mkProj t ]
+          writeVar prevProj (mkProj 1)
+
+    play $ staticFrame (waitT/2) $
+      mkGroup
+      [ scaleToSize screenWidth screenHeight $
+        embedImage $ project src equirectangularP
+      , grid equirectangularP ]
+
+    -- pushInterp "Lambert" lambertP
+    -- 1
+    pushInterp "Web Mercator" mercatorP
+    -- 2
+    pushInterp "Mollweide" mollweideP
+    -- 3
+    pushInterp "Bottomley 30\\degree" (bottomleyP (toRads 30))
+    -- 4
+    pushInterp "Werner" wernerP
+    -- 5
+    pushInterp "Bonne 45\\degree" (bonneP (toRads 45))
+    pushT
+      (\t -> "Bonne " <> T.pack (show $ round $ fromToS 45 0 t) <> "\\degree")
+      (bonneP . toRads . fromToS 45 0)
+    -- 6
+    pushInterp "Eckert I" eckert1P
+    eckert <- newSpriteSVG $ renderLabel "Eckert"
+    spriteZ eckert 2
+    pushInterp "Eckert III" eckert3P
+    pushInterp "Eckert V" eckert5P
+    destroySprite eckert
+    -- 7
+    pushInterp "Fahey" faheyP
+    -- 8
+    pushInterp "August" augustP
+    -- 9
+    pushInterp "Foucaut" foucautP
+    -- 10
+    pushInterp "Lagrange" lagrangeP
+    pushInterp "Equirectangular" equirectangularP
+  where
+    src = equirectangular
+    waitT = 2
+    morphT = 2
+
+renderLabel label =
+  let ref = scale 1.5 $ latex "\\texttt{Tygv123}"
+      glyphs = scale 1.5 $ latex ("\\texttt{" <> label <> "}")
+      svgTxt = mkGroup
+        [ withStrokeColor "black" $ withFillColor "white" $
+          glyphs
+        , withFillColor "white" $
+          glyphs ]
+  in
+    translate (screenWidth*0.01) (screenHeight*0.02) $
+    translate (-screenWidth/2) (-screenHeight/2) $
+    translate 0 (svgHeight ref) svgTxt
+
+equirectangular :: Image PixelRGB8
+equirectangular = unsafePerformIO $ do
+  dat <- BS.readFile "earth.jpg"
+  case decodeJpeg dat of
+    Left err  -> error err
+    Right img -> return $ convertRGB8 img
+
+toRads :: Double -> Double
+toRads dec = dec/180 * pi
+
+
+grid :: Projection -> SVG
+grid p =
+
+  scaleXY
+    (screenWidth)
+    (screenHeight)
+   $
+  translate (-1/2) (-1/2) $
+  withStrokeWidth strokeWidth $
+
+  withFillOpacity 0 $
+  mkGroup
+  [ mkGroup []
+  , withStrokeColorPixel (PixelRGBA8 0x90 0x90 0x90 0x0) $
+    withFillOpacity 0 $ mkGroup
+    [ geometryToSVG p geo
+    | geo <- landBorders
+    ]
+  , withStrokeColorPixel (PixelRGBA8 0x30 0x30 0x30 0x0) $
+    mkGroup $ map mkLinePath (latitudeLines p ++ longitudeLines p)
+  ]
+  where
+    strokeWidth = defaultStrokeWidth * 0.02
+
+worldLine :: Projection -> SVG
+worldLine p =
+  mkLinePath $
+    map apply
+    [ (-pi, -halfPi)
+    , (-pi, halfPi)
+    , (pi, halfPi)
+    , (pi, -halfPi)
+    , (-pi, -halfPi) ]
+  where
+    apply (lam, phi) =
+      let XYCoord x y = projectionForward p $ LonLat lam phi
+      in (x, y)
+
+latitudeLines :: Projection -> [[(Double, Double)]]
+latitudeLines p =
+    [ latitudeLine (fromToS (-pi) pi (n/(latLines*2)))
+    | n <- [0 .. latLines*2]]
+  where
+    latLines = 4
+    segments = 100
+    maxLat = atan (sinh pi)
+    latitudeLine lam =
+      [ (x, y)
+      | n <- [0..segments]
+      , let phi = fromToS (-maxLat) maxLat (n/segments)
+      , let XYCoord x y = projectionForward p $ LonLat lam phi ]
+
+longitudeLines :: Projection -> [[(Double, Double)]]
+longitudeLines p =
+    longitudeLine maxLat :
+    longitudeLine (-maxLat) :
+    [ longitudeLine (fromToS (-halfPi) halfPi (n/(lonLines*2)))
+    | n <- [0 .. lonLines*2]]
+  where
+    lonLines = 4
+    segments = 100
+    maxLat = atan (sinh pi)
+    longitudeLine phi =
+      [ (x, y)
+      | n <- [0..segments]
+      , let lam = fromToS (-pi) pi (n/segments)
+      , let XYCoord x y = projectionForward p $ LonLat lam phi ]
+
+halfPi :: Double
+halfPi = pi/2
+
+landBorders :: [(GeospatialGeometry)]
+landBorders = unsafePerformIO $ do
+  Just geo <- decodeFileStrict "land.geojson"
+  return
+    [ (feature ^. geometry)
+    | feature <- toList $ geo ^. geofeatures
+    , let p = feature ^. properties :: Map String Value
+    ]
+
+geometryToSVG :: Projection -> GeospatialGeometry -> SVG
+geometryToSVG p geometry =
+  case geometry of
+    MultiPolygon mpolygon ->
+      mkGroup $ map (geometryToSVG p . Polygon) $ toList (splitGeoMultiPolygon mpolygon)
+    Polygon poly ->
+      mkGroup
+      [ mkLinePath section
+      | section <- pure
+          [ (x', y')
+          | PointXY x y <- map retrieveXY (fromLinearRing (head (toList (poly^.unGeoPolygon))))
+          , let XYCoord x' y' = projectionForward p $ LonLat (x/180*pi) (y/180*pi)
+          ]
+      ]
+    Line line ->
+      mkLinePath
+      [ (x', y')
+      | PointXY x y <- map retrieveXY (Line.fromLineString (line ^. unGeoLine))
+      , let XYCoord x' y' = projectionForward p $ LonLat (x/180*pi) (y/180*pi)
+      ]
+    MultiLine ml ->
+      mkGroup $ map (geometryToSVG p . Line) $ toList (splitGeoMultiLine ml)
+    _ -> None

videos/map-projection/map-projection.hs

@@ -0,0 +1,6 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+module Main(main) where
+
+main :: IO ()
+main = return ()

videos/map-projection/output.geojson.REMOVED.git-id

@@ -0,0 +1 @@
+088ba82a112f25607c756814500ad4246d7ced4c