{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} module Geodetics.TransverseMercator( GridTM (trueOrigin, falseOrigin, gridScale), mkGridTM ) where import Data.Function import Data.Monoid import Geodetics.Ellipsoids import Geodetics.Geodetic import Geodetics.Grid import Numeric.Units.Dimensional.Prelude hiding ((.)) import Prelude () -- | A Transverse Mercator projection gives an approximate mapping of the ellipsoid on to a 2-D grid. It models -- a sheet curved around the ellipsoid so that it touches it at one north-south line (hence making it part of -- a slightly elliptical cylinder). data GridTM e = GridTM { trueOrigin :: Geodetic e, -- ^ A point on the line where the projection touches the ellipsoid (altitude is ignored). falseOrigin :: GridOffset, -- ^ The grid position of the true origin. Used to avoid negative coordinates over -- the area of interest. The altitude gives a vertical offset from the ellipsoid. gridScale :: Dimensionless Double, -- ^ A scaling factor that balances the distortion between the east & west edges and the middle -- of the projection. -- Remaining elements are memoised parameters computed from the ellipsoid underlying the true origin. gridN1, gridN2, gridN3, gridN4 :: Dimensionless Double } deriving (Show) -- | Create a Transverse Mercator grid. mkGridTM :: (Ellipsoid e) => Geodetic e -- ^ True origin. -> GridOffset -- ^ Vector from true origin to false origin. -> Dimensionless Double -- ^ Scale factor. -> GridTM e mkGridTM origin offset sf = GridTM {trueOrigin = origin, falseOrigin = offset, gridScale = sf, gridN1 = _1 + n + (_5/_4) * n^pos2 + (_5/_4) * n^pos3, gridN2 = _3 * n + _3 * n^pos2 + ((21*~one)/_8) * n^pos3, gridN3 = ((15*~one)/_8) * (n^pos2 + n^pos3), gridN4 = ((35*~one)/(24*~one)) * n^pos3 } where f = flattening $ ellipsoid origin n = f / (_2-f) -- Equivalent to (a-b)/(a+b) where b = (1-f)*a -- | Equation C3 from reference [1]. m :: (Ellipsoid e) => GridTM e -> Dimensionless Double -> Length Double m grid lat = bF0 * (gridN1 grid * dLat - gridN2 grid * sin dLat * cos sLat + gridN3 grid * sin (_2 * dLat) * cos (_2 * sLat) - gridN4 grid * sin (_3 * dLat) * cos (_3 * sLat)) where dLat = lat - latitude (trueOrigin grid) sLat = lat + latitude (trueOrigin grid) bF0 = minorRadius (gridEllipsoid grid) * gridScale grid instance (Ellipsoid e) => GridClass (GridTM e) e where fromGrid p = Geodetic (lat' - east' ^ pos2 * tanLat / (_2 * rho * v) -- Term VII + east' ^ pos4 * (tanLat / ((24 *~ one) * rho * v ^ pos3)) * (_5 + _3 * tanLat ^ pos2 + eta2 - _9 * tanLat ^ pos2 * eta2) -- Term VIII - east' * east' ^ pos5 * (tanLat / ((720 *~ one) * rho * v ^ pos5)) * (61 *~ one + (90 *~ one) * tanLat ^ pos2 + (45 *~ one) * tanLat ^ pos4)) -- Term IX (longitude (trueOrigin grid) + east' / (cosLat * v) -- Term X - (east' ^ pos3 / (_6 * cosLat * v ^ pos3)) * (v / rho + _2 * tanLat ^ pos2) -- Term XI + (east' ^ pos5 / ((120 *~ one) * cosLat * v ^ pos5)) * (_5 + (28 *~ one) * tanLat ^ pos2 + (24 *~ one) * tanLat ^ pos4) -- Term XII - (east' ^ pos5 * east' ^ pos2 / ((5040 *~ one) * cosLat * v * v * v ^ pos5)) * ((61 *~ one) + (662 *~ one) * tanLat ^ pos2 + (1320 *~ one) * tanLat ^ pos4 + (720 *~ one) * tanLat * tanLat ^ pos5)) -- Term XIIa (0 *~ meter) (gridEllipsoid grid) where GridPoint east' north' _ _ = falseOrigin grid `applyOffset` p lat' = fst $ head $ dropWhile ((> 0.01 *~ milli meter) . snd) $ tail $ iterate next (latitude $ trueOrigin grid, 1 *~ meter) where next (phi, _) = let delta = north' - m grid phi in (phi + delta / aF0, delta) -- head and tail are safe because iterate returns an infinite list. sinLat = sin lat' cosLat = cos lat' tanLat = tan lat' sinLat2 = sinLat ^ pos2 v = aF0 / sqrt (_1 - e2 * sinLat2) rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one) eta2 = v / rho - _1 aF0 = majorRadius (gridEllipsoid grid) * gridScale grid e2 = eccentricity2 $ gridEllipsoid grid grid = gridBasis p toGrid grid geo = applyOffset (off `mappend` (offsetNegate $ falseOrigin grid)) $ GridPoint _0 _0 _0 grid where v = aF0 / sqrt (_1 - e2 * sinLat2) rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one) eta2 = v / rho - _1 off = GridOffset (dLong * term_IV + dLong ^ pos3 * term_V + dLong ^ pos5 * term_VI) (m grid lat + dLong ^ pos2 * term_II + dLong ^ pos4 * term_III + dLong * dLong ^ pos5 * term_IIIa) (0 *~ meter) -- Terms defined in [1]. term_II = (v/_2) * sinLat * cosLat term_III = (v/(24*~one)) * sinLat * cosLat ^ pos3 * (_5 - tanLat ^ pos2 + _9 * eta2) term_IIIa = (v/(720*~one)) * sinLat * cosLat ^ pos5 * ((61 *~ one) - (58 *~ one) * tanLat ^ pos2 + tanLat ^ pos4) term_IV = v * cosLat term_V = (v/_6) * cosLat ^ pos3 * (v/rho - tanLat ^ pos2) term_VI = (v/(120*~one)) * cosLat ^ pos5 * (_5 - (18*~one) * tanLat ^ pos2 + tanLat ^ pos4 + (14*~one) * eta2 - (58*~one) * tanLat ^ pos2 * eta2) {- -- Trace message for debugging. Uncomment this code for easy access to intermediate values. traceMsg = concat [ "v = ", show v, "\n", "rho = ", show rho, "\n", "eta2 = ", show eta2, "\n", "M = ", show $ m grid lat, "\n", "I = ", show $ m grid lat + deltaNorth (falseOrigin grid), "\n", "II = ", show term_II, "\n", "III = ", show term_III, "\n", "IIIa = ", show term_IIIa, "\n", "IV = ", show term_IV, "\n", "V = ", show term_V, "\n", "VI = ", show term_VI, "\n"] -} -- Common subexpressions lat = latitude geo long = longitude geo dLong = long - longitude (trueOrigin grid) sinLat = sin lat cosLat = cos lat tanLat = tan lat sinLat2 = sinLat ^ pos2 aF0 = (majorRadius $ gridEllipsoid grid) * gridScale grid e2 = eccentricity2 $ gridEllipsoid grid gridEllipsoid = ellipsoid . trueOrigin