{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiWayIf #-}
module Reanimate.GeoProjection
( Projection(..)
, XYCoord(..)
, LonLat(..)
, project
, interpP
, mergeP
, isValidP
, scaleP
, flipYAxisP
, moveBottomP
, moveTopP
, equirectangularP
, mercatorP
, mollweideP
, hammerP
, lambertP
, bottomleyP
, sinusoidalP
, wernerP
, bonneP
, orthoP
, cassiniP
, augustP
, collignonP
, eckert1P
, eckert3P
, eckert5P
, faheyP
, foucautP
, lagrangeP
, drawFeatureCollection
, loadFeatureCollection
, applyProjection
, applyProjection'
, renderGeometry
) where
import Codec.Picture
import Codec.Picture.Types
import Control.Lens ((^.))
import Control.Monad
import Control.Monad.ST
import Data.Aeson
import Data.Foldable
import Data.Geospatial hiding (LonLat)
import Data.LinearRing
import qualified Data.LineString as Line
import Debug.Trace
import Graphics.SvgTree (Tree (None))
import Linear (distance, lerp)
import Linear.V2 hiding (angle)
import Reanimate
import System.IO.Unsafe
halfPi, sqrtPi, sqrt2, epsilon, tau :: Double
halfPi = pi/2
sqrtPi = sqrt pi
sqrt2 = sqrt 2
epsilon = 1.0e-12
tau = pi*2
toRads, cot :: Double -> Double
toRads dec = dec/180 * pi
cot = recip . tan
srcPixel :: Pixel pixel => Image pixel -> LonLat -> pixel
srcPixel src (LonLat lam phi) =
pixelAt src xPx yPx
where
!xPx = round $ ((lam+pi)/tau) * fromIntegral (imageWidth src-1)
!yPx = round $ (1-((phi+halfPi)/pi)) * fromIntegral (imageHeight src-1)
findValidCoord :: Image pixel -> Projection -> XYCoord -> XYCoord
findValidCoord !src !p (XYCoord x y) = foldr const (XYCoord x y)
[ XYCoord x' y'
| let !xi = round $ x * wMax
!yi = round $ y * hMax
, !ax <- [xi, xi-1, xi+1]
, ax >= 0
, ax < w
, let !x' = fromIntegral ax / wMax
, !ay <- [yi, yi-1, yi+1]
, ay >= 0
, ay < h
, let !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)
, x < (x2 - x1) * (y - y1) / (y2 - y1) + x1
]
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 $ take 20
[ (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 [] = pure $ PixelRGBA8 0xFF 0x00 0x00 0xFF
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 PixelRGBA8 -> 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
let factor = 2
let l1 =
loopTo (w*factor) $ \x -> do
loopTo (h*factor) $ \y -> do
let !x1' = fromIntegral x / (wMax*fromIntegral factor)
!y1' = fromIntegral y / (hMax*fromIntegral factor)
!lonlat = projectionInverse p1 $! XYCoord x1' y1'
p = srcPixel src lonlat
when (validLonLat lonlat && pixelOpacity p /= 0) $ do
let 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 (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $
writePixel img x3 y3 p
l2 =
loopTo (w*factor) $ \x ->
loopTo (h*factor) $ \y -> do
let !x2' = fromIntegral x / (wMax*fromIntegral factor)
!y2' = fromIntegral y / (hMax*fromIntegral factor)
!lonlat = projectionInverse p2 (XYCoord x2' y2')
p = srcPixel src lonlat
when (validLonLat lonlat && pixelOpacity p /= 0) $ do
let XYCoord x2 y2 = projectionForward p2 lonlat
XYCoord x1 y1 = findValidCoord src p1 $ projectionForward p1 lonlat
!x3 = round $ fromToS x1 x2 t * wMax
!y3 = round $ (1 - fromToS y1 y2 t) * hMax
when (x3 >= 0 && x3 < w && y3 >= 0 && y3 < h) $ do
writePixel img x3 y3 p
if t < 0.5
then l1 >> l2
else l2 >> l1
when False $
forM_ [0..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
loopTo m fn = go m
where go 0 = return ()
go n = fn (n-1) >> go (n-1)
!w = imageWidth src
!h = imageHeight src
!wMax = fromIntegral (w-1)
!hMax = fromIntegral (h-1)
eqLonLat :: LonLat -> LonLat -> Bool
eqLonLat (LonLat x1 y1) (LonLat x2 y2)
= eqDouble x1 x2 && eqDouble y1 y2
eqDouble :: Double -> Double -> Bool
eqDouble a b = abs (a-b) < epsilon
data XYCoord = XYCoord !Double !Double
deriving (Read,Show,Eq,Ord)
data LonLat = LonLat !Double !Double
deriving (Read,Show,Eq,Ord)
data Projection = Projection
{ projectionForward :: !(LonLat -> XYCoord)
, projectionInverse :: !(XYCoord -> LonLat)
}
project :: Image PixelRGBA8 -> 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 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 not (validLonLat lonlat) || eqLonLat lonlat lonlat2
|| 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
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
mercatorP :: Projection
mercatorP = Projection forward inverse
where
forward (LonLat lam phi) =
XYCoord ((lam+pi)/tau)
(min 1 $ max (0) $ (((log(tan(pi/4+phi/2))) + pi)/tau))
inverse (XYCoord x y) = LonLat xPi (atan (sinh yPi))
where
xPi = fromToS (-pi) pi x
yPi = fromToS (-pi) pi y
mollweideP :: Projection
mollweideP = Projection forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
where
x = (2*sqrt2)/pi * lam * cos theta
y = sqrt2*sin theta
theta = 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*sqrt2) (2*sqrt2) x'
y = fromToS (-sqrt2) sqrt2 y'
theta = asin (y/sqrt2)
lam = pi*x/(2*sqrt2*cos theta)
phi = asin ((2*theta+sin(2*theta))/pi)
hammerP :: Projection
hammerP = Projection forward inverse
where
forward (LonLat lam phi) =
XYCoord ((x+sqrt2*2)/(4*sqrt2)) ((y+sqrt2)/(2*sqrt2))
where
x = (2*sqrt2*cos phi*sin (lam/2))/sqrt (1+cos phi*cos (lam/2))
y = (sqrt2*sin phi)/ sqrt (1+cos phi*cos (lam/2))
inverse (XYCoord x' y') = LonLat lam phi
where
x = fromToS (-2*sqrt2) (2*sqrt2) x'
y = fromToS (-sqrt2) sqrt2 y'
z = sqrt (1 - (x/4)**2 - (y/2)**2)
lam = 2 * atan2 (z*x) (2*(2*z**2-1))
phi = asin (z*y)
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'
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
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'
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 = if rho == 0 then rho else lam * sin rho / rho
x = rho * sin e
y = pi/2 - rho * cos e
inverse (XYCoord x' y') = LonLat lam phi
where
x = fromToS (-pi) pi x'
y = fromToS (-pi/2) (pi/2) y'
rho = sqrt (x**2 + (pi/2 - y)**2)
e = atan2 x (pi/2 -y)
phi = pi/2 - rho
lam = e * rho / sin rho
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)
phi = cotPhi0 + phi_0 - rho
lam = rho / cos phi * atan2 x (cotPhi0-y)
orthoP :: LonLat -> Projection
orthoP (LonLat lam_0 phi_0) = Projection forward inverse
where
forward (LonLat lam phi)
| cosc < 0 =
let ang = atan2 y x
xV = cos ang
yV = sin ang
xPos = 7/32 + ((xV+1)/2 * 9/16)
in
XYCoord xPos ((yV+1)/2)
| otherwise = XYCoord ((x+(16/9))/(16/9*2)) ((y+1)/2)
where
x = cos phi * sin (lam - lam_0)
y = cos phi_0 * sin phi - sin phi_0 * cos phi * cos (lam - lam_0)
cosc = sin phi_0 * sin phi + cos phi_0 * cos phi * cos (lam-lam_0)
inverse (XYCoord x' y') = LonLat lam phi
where
x = fromToS (-16/9) (16/9) x'
y = fromToS (-1) 1 y'
lam = wrap (-pi) pi $
lam_0 + atan2 (x * sin c) (rho * cos c * cos phi_0 - y * sin c * sin phi_0)
phi = wrap (-pi/2) (pi/2) $
asin (cos c * sin phi_0 + (y * sin c * cos phi_0)/rho)
rho = sqrt (x*x + y*y)
c = asin rho
wrap lower upper v
| v > upper = v-upper+lower
| v < lower = v+upper-lower
| otherwise = v
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)
collignonP :: Projection
collignonP = Projection forward inverse
where
yHi = sqrtPi
yLo = sqrtPi * (1 - sqrt2)
xLo = -pi*(2/sqrtPi)*sqrt2
xHi = pi*(2/sqrtPi)*sqrt2
forward (LonLat lam phi) = XYCoord ((x-xLo)/(xHi-xLo)) ((y-yLo)/(yHi-yLo))
where
alpha = sqrt (1 - sin phi)
x = (2 / sqrtPi) * lam * alpha
y = sqrtPi * (1 - alpha)
inverse (XYCoord x' y') = LonLat lam phi
where
x = fromToS xLo xHi x'
y = fromToS yLo yHi y'
l = (y / sqrtPi -1)**2
lam = if l > 0 then x * sqrt (pi / l) / 2 else 0
phi = asin (1-l)
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))
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
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
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))
drawFeatureCollection :: GeoFeatureCollection a -> (a -> SVG -> SVG) -> SVG
drawFeatureCollection geo fn = mkGroup
[ fn (feature ^. properties) $ renderGeometry (feature ^. geometry)
| feature <- toList (geo ^. geofeatures)
]
{-# INLINE loadFeatureCollection #-}
loadFeatureCollection :: FromJSON a => FilePath -> (a -> SVG -> SVG) -> SVG
loadFeatureCollection path = unsafePerformIO $ do
mbGeo <- decodeFileStrict path
case mbGeo of
Nothing -> error $ "loadFeatureCollection: Invalid GeoJSON: " ++ path
Just geo -> return (drawFeatureCollection geo)
renderGeometry :: GeospatialGeometry -> SVG
renderGeometry shape =
case shape of
MultiPolygon mpolygon ->
mkGroup $ map (renderGeometry . Polygon) $ toList (splitGeoMultiPolygon mpolygon)
Polygon poly ->
mkGroup
[ mkLinePath section
| section <- pure
[ (x, y)
| PointXY x y <- map retrieveXY (fromLinearRing (head (toList (poly^.unGeoPolygon))))
]
]
Line line ->
mkLinePath
[ (x, y)
| PointXY x y <- map retrieveXY (Line.fromLineString (line ^. unGeoLine))
]
MultiLine ml ->
mkGroup $ map (renderGeometry . Line) $ toList (splitGeoMultiLine ml)
_ -> None
applyProjection :: Projection -> SVG -> SVG
applyProjection = applyProjection' 1e-2
applyProjection' :: Double -> Projection -> SVG -> SVG
applyProjection' tolerance p = mapSvgLines start
where
start (LineMove x:rest) = LineMove (proj x) : worker x rest
start _ = []
worker a (LineEnd b : rest) =
let (x:xs) = reverse $ drop 1 $ mkChunks a b
in map (\v -> LineBezier [v]) (map proj $ reverse xs) ++ LineEnd (proj x) : start rest
worker a (LineBezier [b] : rest) =
let (x:xs) = reverse $ drop 1 $ mkChunks a b
in map (\v -> LineBezier [v]) (map proj $ reverse xs) ++ LineBezier [proj x] : worker x rest
worker _ (LineBezier ps : rest) =
LineBezier (map proj ps) : worker (last ps) rest
worker _ (LineMove x:rest) = LineMove (proj x) : worker x rest
worker _ [] = []
lowTolerance = tolerance*tolerance
proj (V2 lam phi) =
case projectionForward p $ LonLat lam phi of
XYCoord x y -> V2 x y
mkChunks a b =
let midway = lerp 0.5 a b in
if distance (proj a) (proj b) < tolerance || distance a b < lowTolerance
then [a, b]
else mkChunks a midway ++ drop 1 (mkChunks midway b)