{-# LANGUAGE 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 Geodetics.Altitude
import Geodetics.Ellipsoids
import Geodetics.Geodetic

-- | 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 {
   forall r. GridPoint r -> Double
eastings, forall r. GridPoint r -> Double
northings, forall r. GridPoint r -> Double
altGP :: Double,
   forall r. GridPoint r -> r
gridBasis :: r
} deriving (GridPoint r -> GridPoint r -> Bool
(GridPoint r -> GridPoint r -> Bool)
-> (GridPoint r -> GridPoint r -> Bool) -> Eq (GridPoint r)
forall r. Eq r => GridPoint r -> GridPoint r -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: forall r. Eq r => GridPoint r -> GridPoint r -> Bool
== :: GridPoint r -> GridPoint r -> Bool
$c/= :: forall r. Eq r => GridPoint r -> GridPoint r -> Bool
/= :: GridPoint r -> GridPoint r -> Bool
Eq, Int -> GridPoint r -> ShowS
[GridPoint r] -> ShowS
GridPoint r -> String
(Int -> GridPoint r -> ShowS)
-> (GridPoint r -> String)
-> ([GridPoint r] -> ShowS)
-> Show (GridPoint r)
forall r. Show r => Int -> GridPoint r -> ShowS
forall r. Show r => [GridPoint r] -> ShowS
forall r. Show r => GridPoint r -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: forall r. Show r => Int -> GridPoint r -> ShowS
showsPrec :: Int -> GridPoint r -> ShowS
$cshow :: forall r. Show r => GridPoint r -> String
show :: GridPoint r -> String
$cshowList :: forall r. Show r => [GridPoint r] -> ShowS
showList :: [GridPoint r] -> ShowS
Show)

instance HasAltitude (GridPoint g) where
   altitude :: GridPoint g -> Double
altitude = GridPoint g -> Double
forall r. GridPoint r -> Double
altGP
   setAltitude :: Double -> GridPoint g -> GridPoint g
setAltitude Double
h GridPoint g
gp = GridPoint g
gp{altGP = h}


-- | A vector relative to a point on a grid. All distances are in meters.
-- 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 {
   GridOffset -> Double
deltaEast, GridOffset -> Double
deltaNorth, GridOffset -> Double
deltaAltitude :: Double
} deriving (GridOffset -> GridOffset -> Bool
(GridOffset -> GridOffset -> Bool)
-> (GridOffset -> GridOffset -> Bool) -> Eq GridOffset
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GridOffset -> GridOffset -> Bool
== :: GridOffset -> GridOffset -> Bool
$c/= :: GridOffset -> GridOffset -> Bool
/= :: GridOffset -> GridOffset -> Bool
Eq, Int -> GridOffset -> ShowS
[GridOffset] -> ShowS
GridOffset -> String
(Int -> GridOffset -> ShowS)
-> (GridOffset -> String)
-> ([GridOffset] -> ShowS)
-> Show GridOffset
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GridOffset -> ShowS
showsPrec :: Int -> GridOffset -> ShowS
$cshow :: GridOffset -> String
show :: GridOffset -> String
$cshowList :: [GridOffset] -> ShowS
showList :: [GridOffset] -> ShowS
Show)

instance Semigroup GridOffset where
  GridOffset
g1 <> :: GridOffset -> GridOffset -> GridOffset
<> GridOffset
g2 = Double -> Double -> Double -> GridOffset
GridOffset (GridOffset -> Double
deltaEast GridOffset
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaEast GridOffset
g2)
                        (GridOffset -> Double
deltaNorth GridOffset
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaNorth GridOffset
g2)
                        (GridOffset -> Double
deltaAltitude GridOffset
g1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaAltitude GridOffset
g2)

instance Monoid GridOffset where
   mempty :: GridOffset
mempty = Double -> Double -> Double -> GridOffset
GridOffset Double
0 Double
0 Double
0
   mappend :: GridOffset -> GridOffset -> GridOffset
mappend = GridOffset -> GridOffset -> GridOffset
forall a. Semigroup a => a -> a -> a
(<>)

-- | An offset defined by a distance (m) and a bearing (radians) 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 :: Double -> Double -> GridOffset
polarOffset :: Double -> Double -> GridOffset
polarOffset Double
r Double
d = Double -> Double -> Double -> GridOffset
GridOffset (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sin Double
d) (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
cos Double
d) Double
0


-- | Scale an offset by a scalar.
offsetScale :: Double -> GridOffset -> GridOffset
offsetScale :: Double -> GridOffset -> GridOffset
offsetScale Double
s GridOffset
off = Double -> Double -> Double -> GridOffset
GridOffset (GridOffset -> Double
deltaEast GridOffset
off Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)
                               (GridOffset -> Double
deltaNorth GridOffset
off Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)
                               (GridOffset -> Double
deltaAltitude GridOffset
off Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s)

-- | Invert an offset.
offsetNegate :: GridOffset -> GridOffset
offsetNegate :: GridOffset -> GridOffset
offsetNegate GridOffset
off = Double -> Double -> Double -> GridOffset
GridOffset (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ GridOffset -> Double
deltaEast GridOffset
off)
                              (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ GridOffset -> Double
deltaNorth GridOffset
off)
                              (Double -> Double
forall a. Num a => a -> a
negate (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ GridOffset -> Double
deltaAltitude GridOffset
off)


-- Add an offset on to a point to get another point.
applyOffset :: GridOffset -> GridPoint g -> GridPoint g
applyOffset :: forall g. GridOffset -> GridPoint g -> GridPoint g
applyOffset GridOffset
off GridPoint g
p = Double -> Double -> Double -> g -> GridPoint g
forall r. Double -> Double -> Double -> r -> GridPoint r
GridPoint (GridPoint g -> Double
forall r. GridPoint r -> Double
eastings GridPoint g
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaEast GridOffset
off)
                           (GridPoint g -> Double
forall r. GridPoint r -> Double
northings GridPoint g
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaNorth GridOffset
off)
                           (GridPoint g -> Double
forall a. HasAltitude a => a -> Double
altitude GridPoint g
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaAltitude GridOffset
off)
                           (GridPoint g -> g
forall r. GridPoint r -> r
gridBasis GridPoint g
p)


-- | The distance represented by an offset.
offsetDistance :: GridOffset -> Double
offsetDistance :: GridOffset -> Double
offsetDistance = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double)
-> (GridOffset -> Double) -> GridOffset -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. GridOffset -> Double
offsetDistanceSq


-- | The square of the distance represented by an offset.
offsetDistanceSq :: GridOffset -> Double
offsetDistanceSq :: GridOffset -> Double
offsetDistanceSq GridOffset
off =
   GridOffset -> Double
deltaEast GridOffset
off Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaNorth GridOffset
off Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ GridOffset -> Double
deltaAltitude GridOffset
off Double -> Int -> Double
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
_2


-- | The direction represented by an offset, as bearing to the right of North.
offsetBearing :: GridOffset -> Double
offsetBearing :: GridOffset -> Double
offsetBearing GridOffset
off = Double -> Double -> Double
forall a. RealFloat a => a -> a -> a
atan2 (GridOffset -> Double
deltaEast GridOffset
off) (GridOffset -> Double
deltaNorth GridOffset
off)


-- | The offset required to move from p1 to p2.
gridOffset :: GridPoint g -> GridPoint g -> GridOffset
gridOffset :: forall g. GridPoint g -> GridPoint g -> GridOffset
gridOffset GridPoint g
p1 GridPoint g
p2 = Double -> Double -> Double -> GridOffset
GridOffset (GridPoint g -> Double
forall r. GridPoint r -> Double
eastings GridPoint g
p2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridPoint g -> Double
forall r. GridPoint r -> Double
eastings GridPoint g
p1)
                              (GridPoint g -> Double
forall r. GridPoint r -> Double
northings GridPoint g
p2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridPoint g -> Double
forall r. GridPoint r -> Double
northings GridPoint g
p1)
                              (GridPoint g -> Double
forall a. HasAltitude a => a -> Double
altitude GridPoint g
p2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- GridPoint g -> Double
forall a. HasAltitude a => a -> Double
altitude GridPoint g
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 :: forall b a. b -> GridPoint a -> GridPoint b
unsafeGridCoerce b
base GridPoint a
p = Double -> Double -> Double -> b -> GridPoint b
forall r. Double -> Double -> Double -> r -> GridPoint r
GridPoint (GridPoint a -> Double
forall r. GridPoint r -> Double
eastings GridPoint a
p) (GridPoint a -> Double
forall r. GridPoint r -> Double
northings GridPoint a
p) (GridPoint a -> Double
forall a. HasAltitude a => a -> Double
altitude GridPoint a
p) b
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 * kilometer) "237"@ will return
--
-- > Just (23700 meters, 100 meters)
--
-- If there are any non-digits in the string then the function returns @Nothing@.
fromGridDigits :: Double -> String -> Maybe (Double, Double)
fromGridDigits :: Double -> String -> Maybe (Double, Double)
fromGridDigits Double
sq String
ds = if (Char -> Bool) -> String -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all Char -> Bool
isDigit String
ds then (Double, Double) -> Maybe (Double, Double)
forall a. a -> Maybe a
Just (Double
d, Double
p) else Maybe (Double, Double)
forall a. Maybe a
Nothing
   where
      n :: Integer
      n :: Integer
n = Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Integer) -> Int -> Integer
forall a b. (a -> b) -> a -> b
$ String -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length String
ds
      d :: Double
d = [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Double] -> Double) -> [Double] -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*)
         ((Char -> Double) -> String -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> (Char -> Int) -> Char -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> Int
digitToInt) String
ds)
         (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop Int
1 ([Double] -> [Double]) -> [Double] -> [Double]
forall a b. (a -> b) -> a -> b
$ (Double -> Double) -> Double -> [Double]
forall a. (a -> a) -> a -> [a]
iterate (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
10) Double
sq)
      p :: Double
p = Double
sq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Integer
10 :: Integer) Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Integer
n)

-- | Convert a distance into a digit string suitable for printing as part
-- of a grid reference. The result is the south or west side of the enclosing grid square,
-- where the size of the square is defined by the number of digits.
-- The result is expressed as an integer count of squares and a string of digits.
-- If any arguments are invalid then @Nothing@ is returned.
toGridDigits ::
   Double           -- ^ Size of enclosing grid square. Must be at least 1000m.
   -> Int           -- ^ Number of digits to return. Must be positive.
   -> Double        -- ^ Offset to convert into grid (m).
   -> Maybe (Integer, String)
toGridDigits :: Double -> Int -> Double -> Maybe (Integer, String)
toGridDigits Double
sq Int
n Double
d =
   if Double
sq Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1000 Bool -> Bool -> Bool
|| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0
   then Maybe (Integer, String)
forall a. Maybe a
Nothing
   else
      (Integer, String) -> Maybe (Integer, String)
forall a. a -> Maybe a
Just (Integer
sqs, String
pad)
   where
      p :: Integer
      p :: Integer
p = Integer
10 Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^ Int
n
      unit :: Double
      unit :: Double
unit = Double
sq Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Integer -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
p
      u :: Integer
u = Double -> Integer
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double
d Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
unit)
      (Integer
sqs, Integer
d1) = Integer
u Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Integer
p
      s :: String
s = Integer -> String
forall a. Show a => a -> String
show Integer
d1
      pad :: String
pad = if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then String
"" else Int -> Char -> String
forall a. Int -> a -> [a]
replicate (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- String -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length String
s) Char
'0' String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
s