{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} module Geodetics.Grid ( -- ** Grid types GridClass (..), GridPoint (..), GridOffset (..), -- ** Grid operations polarOffset, offsetScale, offsetNegate, applyOffset, offsetDistance, offsetDistanceSq, offsetBearing, gridOffset, -- ** Unsafe conversion unsafeGridCoerce, -- ** Utility functions for grid references fromGridDigits, toGridDigits ) where import Data.Char import Data.Function import Data.Monoid (Monoid) import Data.Semigroup (Semigroup, (<>)) import Geodetics.Altitude import Geodetics.Geodetic import Numeric.Units.Dimensional.Prelude hiding ((.)) import qualified Prelude as P -- | A Grid is a two-dimensional projection of the ellipsoid onto a plane. Any given type of grid can -- usually be instantiated with parameters such as a tangential point or line, and these parameters -- will include the terrestrial reference frame ("Ellipsoid" in this library) used as a foundation. -- Hence conversion from a geodetic to a grid point requires the \"basis\" for the grid in question, -- and grid points carry that basis with them because without it there is no defined relationship -- between the grid points and terrestrial positions. class GridClass r e | r->e where fromGrid :: GridPoint r -> Geodetic e toGrid :: r -> Geodetic e -> GridPoint r gridEllipsoid :: r -> e -- | A point on the specified grid. data GridPoint r = GridPoint { eastings, northings, altGP :: Length Double, gridBasis :: r } deriving (Show) instance Eq (GridPoint r) where p1 == p2 = eastings p1 == eastings p2 && northings p1 == northings p2 && altGP p1 == altGP p2 instance HasAltitude (GridPoint g) where altitude = altGP setAltitude h gp = gp{altGP = h} -- | A vector relative to a point on a grid. -- Operations that use offsets will only give -- meaningful results if all the points come from the same grid. -- -- The monoid instance is the sum of offsets. data GridOffset = GridOffset { deltaEast, deltaNorth, deltaAltitude :: Length Double } deriving (Eq, Show) instance Semigroup GridOffset where g1 <> g2 = GridOffset (deltaEast g1 + deltaEast g2) (deltaNorth g1 + deltaNorth g2) (deltaAltitude g1 + deltaAltitude g2) instance Monoid GridOffset where mempty = GridOffset _0 _0 _0 mappend = (<>) -- | An offset defined by a distance and a bearing to the right of North. -- -- There is no elevation parameter because we are using a plane to approximate an ellipsoid, -- so elevation would not provide a useful result. If you want to work with elevations -- then "rayPath" will give meaningful results. polarOffset :: Length Double -> Angle Double -> GridOffset polarOffset r d = GridOffset (r * sin d) (r * cos d) _0 -- | Scale an offset by a scalar. offsetScale :: Dimensionless Double -> GridOffset -> GridOffset offsetScale s off = GridOffset (deltaEast off * s) (deltaNorth off * s) (deltaAltitude off * s) -- | Invert an offset. offsetNegate :: GridOffset -> GridOffset offsetNegate off = GridOffset (negate $ deltaEast off) (negate $ deltaNorth off) (negate $ deltaAltitude off) -- Add an offset on to a point to get another point. applyOffset :: GridOffset -> GridPoint g -> GridPoint g applyOffset off p = GridPoint (eastings p + deltaEast off) (northings p + deltaNorth off) (altitude p + deltaAltitude off) (gridBasis p) -- | The distance represented by an offset. offsetDistance :: GridOffset -> Length Double offsetDistance = sqrt . offsetDistanceSq -- | The square of the distance represented by an offset. offsetDistanceSq :: GridOffset -> Area Double offsetDistanceSq off = deltaEast off ^ pos2 + deltaNorth off ^ pos2 + deltaAltitude off ^ pos2 -- | The direction represented by an offset, as bearing to the right of North. offsetBearing :: GridOffset -> Angle Double offsetBearing off = atan2 (deltaEast off) (deltaNorth off) -- | The offset required to move from p1 to p2. gridOffset :: GridPoint g -> GridPoint g -> GridOffset gridOffset p1 p2 = GridOffset (eastings p2 - eastings p1) (northings p2 - northings p1) (altitude p2 - altitude p1) -- | Coerce a grid point of one type into a grid point of a different type, -- but with the same easting, northing and altitude. This is unsafe because it -- will produce a different position unless the two grids are actually equal. -- -- It should be used only to convert between distinguished grids (e.g. "UkNationalGrid") and -- their equivalent numerical definitions. unsafeGridCoerce :: b -> GridPoint a -> GridPoint b unsafeGridCoerce base p = GridPoint (eastings p) (northings p) (altitude p) base -- | Convert a list of digits to a distance. The first argument is the size of the -- grid square within which these digits specify a position. The first digit is -- in units of one tenth of the grid square, the second one hundredth, and so on. -- The first result is the lower limit of the result, and the second is the size -- of the specified offset. -- So for instance @fromGridDigits (100 *~ kilo meter) "237"@ will return -- -- > Just (23700 meters, 100 meters) -- -- If there are any non-digits in the string then the function returns @Nothing@. fromGridDigits :: Length Double -> String -> Maybe (Length Double, Length Double) fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing where n = length ds d = sum $ zipWith (*) (map ((*~ one) . fromIntegral . digitToInt) ds) (tail $ iterate (/ (10 *~ one)) sq) p = sq / ((10 *~ one) ** (fromIntegral n *~ one)) -- | Convert a distance into a digit string suitable for printing as part -- of a grid reference. The result is the nearest position to the specified -- number of digits, expressed as an integer count of squares and a string of digits. -- If any arguments are invalid then @Nothing@ is returned. toGridDigits :: Length Double -- ^ Size of enclosing grid square. Must be at least 1 km. -> Int -- ^ Number of digits to return. Must be positive. -> Length Double -- ^ Offset to convert into grid. -> Maybe (Integer, String) toGridDigits sq n d = if sq < (1 *~ kilo meter) || n < 0 || d < _0 then Nothing else Just (sqs, pad) where p :: Integer p = 10 P.^ n unit :: Length Double unit = sq / (fromIntegral p *~ one) u = round ((d / unit) /~ one) (sqs, d1) = u `divMod` p s = show d1 pad = if n == 0 then "" else replicate (n P.- length s) '0' ++ s