Projections.hs

-- | Module containing several useful projections for generating maps.
module GIS.Math.Projections ( -- * Reference points
                              washingtonDC
                            , mecca
                            -- * Projections
                            , littow
                            , craig
                            , winkel3
                            , mercator
                            , bonne
                            , albers
                            -- * Helper functions
                            , project
                            , projectMap
                            ) where

import           Control.Arrow
import           Control.Lens
import           GIS.Graphics.Types
import           GIS.Math.Utils
import           GIS.Types

-- | For use as a reference point
washingtonDC :: Point
washingtonDC = toRadians (38.9072, -77.0369)

-- | For use as a reference point
mecca :: Point
mecca = toRadians (21.3891, 39.8579)

-- | Littow retroazimuthal + conformal projection
littow :: Projection
littow (long, lat) = (sin(long - referenceLong)/cos lat, cos (long - referenceLong) * tan lat)
    where referenceLong = radians 0

-- | Craig retroazimuthal projection
craig :: Point -> Projection
craig referencePoint (long, lat) = (long - referenceLong, y)
    where (referenceLong, referenceLat) = referencePoint
          expr = sin lat * cos (long - referenceLong) - tan referenceLat * cos lat
          y | long - referenceLong == 0 = expr
            | otherwise = (long - referenceLong) / sin (long - referenceLong) * expr

-- | Winkel Tripel projection
winkel3 :: Projection
winkel3 (long, lat) = ((lambda * cos phi1 + (2 * cos lat * sin (lambda/2)/sinc alpha))/2, (lat + sin lat/sinc alpha)/2)
    where lambda = long - lambda0
          phi1 = acos $ 2 / pi
          alpha = acos $ cos lat * cos (lambda/2)
          lambda0 = radians (-77.0369)

-- | Mercator projection.
mercator :: Projection
mercator (long, lat) = (long - meridian, asinh (tan lat))
    where meridian = radians (-98.5795)

-- | Bonne projection with standard parallel at 45 N and central meridian
-- centered at Washington DC.
bonne :: Projection
bonne (long, lat) = (rho * sin e, cot phi1 - rho * cos e)
    where rho = cot phi1 + phi1 - lat
          e = (long - meridian) * cos lat / rho
          phi1 = radians 45 -- standard parallel @ 45 N
          meridian = radians (-77.0369) -- central meridian @ dc
          cot = (1/) . tan

-- | Albers projection for a given reference point.
--
-- > ablers washingtonDC
albers :: Point -> Projection
albers referencePoint (long, lat) = (rho * sin theta, rho' - rho * cos theta)
    where n = (sin phi1 + sin phi2) / 2
          theta = n * (long - referenceLong)
          c = cos phi1^(2 :: Int) + 2 * n * sin phi1
          rho = sqrt (c - 2 * n * sin lat) / n
          rho' = sqrt (c - 2 * n * sin referenceLat) / n
          -- standard parallels @ 20, 50 degrees
          phi1 = radians 20
          phi2 = radians 50
          (referenceLong, referenceLat) = referencePoint

-- | Project given a `Polygon`.
project :: Projection -> Polygon -> Polygon
project f = fmap (f . toRadians)

-- | Apply a given projection a `Map`.
projectMap :: Projection -> Map -> Map
projectMap p = over labelledDistricts (fmap (first (fmap (project p))))