module GIS.Math.Projections (
washingtonDC
, mecca
, littow
, craig
, winkel3
, mercator
, bonne
, albers
, project
, projectMap
) where
import Control.Arrow
import Control.Lens
import GIS.Graphics.Types
import GIS.Math.Utils
import GIS.Types
washingtonDC :: Point
washingtonDC = toRadians (38.9072, -77.0369)
mecca :: Point
mecca = toRadians (21.3891, 39.8579)
littow :: Projection
littow (long, lat) = (sin(long - referenceLong)/cos lat, cos (long - referenceLong) * tan lat)
where referenceLong = radians 0
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
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 (long, lat) = (long - meridian, asinh (tan lat))
where meridian = radians (-98.5795)
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
meridian = radians (-77.0369)
cot = (1/) . tan
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
phi1 = radians 20
phi2 = radians 50
(referenceLong, referenceLat) = referencePoint
project :: Projection -> Polygon -> Polygon
project f = fmap (f . toRadians)
projectMap :: Projection -> Map -> Map
projectMap p = over labelledDistricts (fmap (first (fmap (project p))))