module GIS.Math.Projections where
import Control.Lens
import Control.Lens.Tuple
import GIS.Types
import GIS.Math.Utils
import GIS.Graphics.Types
washingtonDC :: Point
washingtonDC = over _2 radians $ over _1 radians (38.9072, 77.0369)
mecca :: Point
mecca = over _2 radians $ over _1 radians (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
y = if long referenceLong == 0 then expr else (long referenceLong)/(sin(long referenceLong)) * expr
expr = ((sin lat) * (cos (long referenceLong)) (tan referenceLat) * (cos lat))
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 + 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 (over _1 (fmap (project p))))