geodetics 1.1.0 → 2.0.0

raw patch · 11 files changed

+786/−276 lines, 11 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• README.md +3/−3
• ToDo.txt +2/−2
• changelog.md +14/−0
• geodetics.cabal +3/−1
• src/Geodetics/Grid.hs +18/−7
• src/Geodetics/MGRS.hs +409/−0
• src/Geodetics/Path.hs +18/−18
• src/Geodetics/PolarStereographic.hs +220/−0
• src/Geodetics/TransverseMercator.hs +5/−2
• src/Geodetics/UTM.hs +27/−239
• test/Main.hs +67/−4

@@ -11,9 +11,9 @@ between them and the bearing from one to the other.  Once you have a geodetic position defined you can project it onto a flat plane, or Grid.-At present Transverse Mercator and Oblique Stereographic grids are provided. More can be-added by defining new instances of the Grid typeclass: see "AddingProjections.txt" for-detais.+At present Transverse Mercator, Polar Stereographic and Oblique Stereographic grids are+provided. More can be added by defining new instances of the Grid typeclass: see+"AddingProjections.txt" for detais.  The Paths module defines a path as a parametric function of distance that returns a position and a bearing. Given two paths you can find their intersection using a fast

@@ -1,6 +1,6 @@ To Do ===== -* Add Universal Polar grid.+* More national grid systems. -* Add more national grid systems.+* Great-circle paths.

@@ -39,3 +39,17 @@   `Geodetics.LatLongParser`.  * Added Universal Transverse Mercator (UTM) grid, with both UTM and MGRS grid references++## Version 2.0.0++This version has breaking changes relative to version 1.1.0.++* Added PolarStereographic projection with UPS grid++* Moved MGRS grid types and functions to its own module and extended it to include UPS gridpoints.++* Replaced the UtmGridUnit type with a more general GridUnit type in Geodetics.Grid.++* Provided standalone Parsec parsers for MGRS and UTM grid references.++* Modified code formatting to have consistently 2 spaces per indent level.

@@ -1,6 +1,6 @@ cabal-version:  3.0 name:           geodetics-version:        1.1.0+version:        2.0.0 build-type:     Simple author:         Paul Johnson <paul@cogito.org.uk> tested-with:    ghc ==9.6.6 ghc ==9.8.4 ghc ==9.10.1 ghc ==9.12.1@@ -46,7 +46,9 @@                    Geodetics.Geodetic,                    Geodetics.Grid,                    Geodetics.LatLongParser,+                   Geodetics.MGRS,                    Geodetics.Path,+                   Geodetics.PolarStereographic,                    Geodetics.Stereographic,                    Geodetics.TransverseMercator,                    Geodetics.UK,

@@ -18,7 +18,8 @@    unsafeGridCoerce,    -- * Utility functions for grid references    fromGridDigits,-   toGridDigits+   toGridDigits,+   GridUnit (..) ) where  import Data.Char@@ -67,11 +68,12 @@    mempty = GridOffset 0 0 0    mappend = (<>) + -- | 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.+-- then 'Geodetics.Path.rayPath' will give meaningful results. polarOffset :: Double -> Double -> GridOffset polarOffset r d = GridOffset (r * sin d) (r * cos d) 0 @@ -124,22 +126,22 @@ -- 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.+-- This function 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 * kilometer) "237"@ will return ----- > Just (23700 meters, 100 meters)+-- For instance @fromGridDigits (100 * kilometer) "237"@ will return --+-- > Just (23700, 100)+-- -- If there are any non-digits in the string then the function returns @Nothing@. fromGridDigits :: Double -> String -> Maybe (Double, Double) fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing@@ -151,6 +153,7 @@          (drop 1 $ iterate (/ 10) sq)       p = sq / fromIntegral ((10 :: 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.@@ -175,3 +178,11 @@       (sqs, d1) = u `divMod` p       s = show d1       pad = if n == 0 then "" else replicate (n - length s) '0' ++ s+++-- | Some grids (notably UTM and UPS) can have optional units in their grid references.+data GridUnit = GridMeters | GridKilometers deriving Eq++instance Show GridUnit where+   show GridMeters = "m"+   show GridKilometers = "km"

@@ -0,0 +1,409 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NumericUnderscores #-}++{- | The Military Grid Reference System (MGRS)++In MGRS there are two syntaxes for grid references:++  1. Between 80 South and 84 North a grid reference has a zone number, latitude band letter, a+    2 letter code for the 100km square within the zone, and then northings and eastings within+    that square.++  2. In the polar regions a grid reference has a latitude band letter (A or B for South, Y or Z +    for North), a 2 letter code for the 100km square within the polar region, and then northings+    and eastings within that square. There is no zone number in the polar regions.+-}+module Geodetics.MGRS (+  -- * MGRS Grid Basis+  MgrsGrid (..),+  mgrsGrid,+  utmToMgrsPoint,+  upsToMgrsPoint,+  fromMgrsPoint,+  -- * Textual Representation+  mgrsBandLetterToLatitude,+  mgrsLatitudeToBandLetter,+  fromMgrsGridReference,+  toMgrsGridReference+) where++import Control.Monad+import Data.Array+import Geodetics.Ellipsoids+import Geodetics.Geodetic+import Geodetics.Grid+import Geodetics.PolarStereographic+import Geodetics.UTM+import Text.Parsec+import Text.Parsec.Error+import Text.Printf+++-- | MGRS grid references can be anywhere on Earth. Hence the position can be either on a UTM or+-- a UPS grid.+data MgrsGrid =+  MgrsUtm UtmZone+  | MgrsUps UpsGrid+  deriving (Eq, Show)++instance GridClass MgrsGrid WGS84 where+  fromGrid p = case gridBasis p of+    MgrsUtm zone -> fromGrid $ unsafeGridCoerce zone p+    MgrsUps grid -> fromGrid $ unsafeGridCoerce grid p++  toGrid (MgrsUtm zone) = unsafeGridCoerce (MgrsUtm zone) . toGrid zone+  toGrid (MgrsUps grid) = unsafeGridCoerce (MgrsUps grid) . toGrid grid++  gridEllipsoid _ = WGS84+++-- | Find the most appropriate grid for the given geodetic position+mgrsGrid :: Geodetic WGS84 -> MgrsGrid+mgrsGrid geo =+  case utmZone geo of+    Just zone -> MgrsUtm zone+    Nothing -> MgrsUps $ if latitude geo > 0 then upsNorth else upsSouth+++-- | Generalise from UTM to MGRS+utmToMgrsPoint :: GridPoint UtmZone -> GridPoint MgrsGrid+utmToMgrsPoint gp = unsafeGridCoerce (MgrsUtm $ gridBasis gp) gp+++-- | Generalise from UPS to MGRS+upsToMgrsPoint :: GridPoint UpsGrid -> GridPoint MgrsGrid+upsToMgrsPoint gp = unsafeGridCoerce (MgrsUps $ gridBasis gp) gp+++-- | Convert an MGRS grid point to either UTM or UPS, depending on its basis.+fromMgrsPoint :: GridPoint MgrsGrid -> Either (GridPoint UtmZone) (GridPoint UpsGrid)+fromMgrsPoint gp = case gridBasis gp of+  MgrsUtm zone -> Left $ unsafeGridCoerce zone gp+  MgrsUps grid -> Right $ unsafeGridCoerce grid gp+++-- | Convert a list of letters into an inverse lookup function.+letterTable :: [Char] -> Char -> Maybe Int+letterTable cs = \c -> if inRange (bounds arr) c then arr ! c else Nothing+  where arr = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip cs [0..]]+  -- Second argument in a lambda so that partial application returns a function with a CAF+  -- rather than computing the array for every call.+++-- | The MGRS latitude band code letters, excluding A and B used for Antarctica (south of -80 degrees)+-- and Y and Z used for the Arctic (north of 84 degrees).+mgrsBandLetters :: [Char]+mgrsBandLetters = "CDEFGHJKLMNPQRSTUVWX"++mgrsBandIx :: Char -> Maybe Int+mgrsBandIx = letterTable mgrsBandLetters++-- | Polar regions have special bands A, B, Y and Z.+mgrsBandLetterToPole :: Char -> Maybe Pole+mgrsBandLetterToPole 'A' = Just SouthPole+mgrsBandLetterToPole 'B' = Just SouthPole+mgrsBandLetterToPole 'Y' = Just NorthPole+mgrsBandLetterToPole 'Z' = Just NorthPole+mgrsBandLetterToPole _ = Nothing++polarBandLetters :: [Char]+polarBandLetters = "ABYZ"++-- | Polar regions use different northings and eastings letters to the rest of the world.+polarEastingLetters :: [Char]+polarEastingLetters = "ABCFGHJKLPQRSTUXYZ"++polarEastingIx :: Char -> Maybe Int+polarEastingIx = letterTable polarEastingLetters++polarNorthingLetters :: [Char]+polarNorthingLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ"++polarNorthingIx :: Char -> Maybe Int+polarNorthingIx = letterTable polarNorthingLetters+++-- | Convert the three letters of an MGRS grid reference to the south-west corner of a 100km square.+--+-- If any of the letters are invalid for the positions they are in then an error message is returned.+mgrsLettersToPolar ::+  Char   -- ^ Latitude band: A, B, Y or Z.+  -> Char  -- ^ Eastings letter.+  -> Char  -- ^ Northings letter.+  -> Parsec s u (Pole, Double, Double)+mgrsLettersToPolar bandC eastC northC =+  case (pole, baseEasting, baseNorthing) of+    (Just pole1, Just easting, Just northing) -> pure (pole1, easting, northing + polarNorthingLetterOrigin pole1)+    _ -> fail "Invalid polar grid letters"+  where+    pole = mgrsBandLetterToPole bandC+    baseEasting = (bandEasting +) . (100_000 *) . fromIntegral <$> polarEastingIx eastC+    baseNorthing = (100_000 *) . fromIntegral <$> polarNorthingIx northC+    bandEasting = if bandC `elem` "AY" then 200*kilometer else 2000*kilometer+      -- A and Y denote the "eastern" halves of the south and north polar regions respectively.+++-- | Find the southern boundary of a latitude band letter (excluding poles).+mgrsBandLetterToLatitude :: Char -> Maybe Double+mgrsBandLetterToLatitude band = do+    n1 <- mgrsBandIx band+    return $ degree * fromIntegral (-80 + n1 * 8)+++-- | Find the band letter for a latitude, if it is in the range (-80, 84) degrees.+-- (Argument in radians)+mgrsLatitudeToBandLetter :: Double -> Maybe Char+mgrsLatitudeToBandLetter lat = do+    guard $ -80 <= dlat && dlat <= 84+    return $ indexMap ! latIdx+  where+    dlat = lat / degree+    ilat :: Int+    ilat = floor dlat+    latIdx = min 19 $ (ilat + 80) `div` 8  -- Band 19 (X) extends an extra 4 degrees.+    indexMap = listArray (0,19) mgrsBandLetters++++-- | Letters A-Z except for I and O.+mgrsEastingsLetters :: [Char]+mgrsEastingsLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ"++mgrsEastingIx :: Char -> Maybe Int+mgrsEastingIx = letterTable mgrsEastingsLetters++-- | If zone number is in range and the letter is one of the valid Eastings letters for that zone+-- then return the UTM easting in meters.+--+-- Zone 1 starts with \'A\'. Each zone is 8 characters wide. Hence the letters repeat every 3 zones.+mgrsLetterToEasting :: UtmZoneNumber -> Char -> Maybe Double+mgrsLetterToEasting zn c = do+    guard $ 1 <= zn && zn <= 60+    n1 <- mgrsEastingIx c+    let n2 = n1 - base+    guard $ 0 <= n2 && n2 <= 7+    return $ fromIntegral (n2+1) * 100_000+  where+    base = ((zn-1) `mod` 3) * 8+++-- | If the zone number is in range and the eastings are between 100,000 and 900,000 then+-- return the Eastings letter.+mgrsEastingToLetter :: UtmZoneNumber -> Double -> Maybe Char+mgrsEastingToLetter zn east = do+    guard $ 1 <= zn && zn <= 60+    guard $ 100 * kilometer <= east && east < 900 * kilometer+    return $ indexMap ! ix+  where+    indexMap = listArray (0,23) mgrsEastingsLetters+    base = ((zn-1) `mod` 3) * 8+    square = max 0 $ min 7 $ floor $ (east - 100 * kilometer)/(100 * kilometer)  -- Clamped in range (0,7).+    ix = base + square  -- Must be in range (0,23)+++-- | Letters A-V except for I and O.+mgrsNorthingsLetters :: [Char]+mgrsNorthingsLetters = "ABCDEFGHJKLMNPQRSTUV"++mgrsNorthingIx :: Char -> Maybe Int+mgrsNorthingIx = letterTable mgrsEastingsLetters+++-- | MGRS Northings letters have rather complex relationship to the latitude bands. The 20 letters+-- repeat every 2,000km going north and south from the equator, so the latitude band is needed+-- to disambiguate which repetition of the Northings letter is meant.++-- Unfortunately this repetition of the letters does not neatly coincide with+-- the latitude band boundaries, which are based on degrees of latitude.++-- The base letter just north of the equator is A in odd-numbered zones and F in even numbered zones.+--+-- This uses the latitude band to estimate the range of northings that would be valid there,+-- and hence determine which possible grid band is meant by the northings letter.+-- The algorithm used is approximate and deliberately very forgiving: it will accept some grid squares+-- which are north or south of the band given.+mgrsLetterToNorthings ::+  UtmZoneNumber+  -> Char  -- ^ Latitude band letter (@C@ - @X@ excluding @I@ and @O@).+  -> Char  -- ^ MGRS Northings letter (@A@ - @V@ excluding @I@ and @O@).+  -> Maybe Double+mgrsLetterToNorthings zone bandC northingsC = do+  guard $ 1 <= zone && zone <= 60+  band <- (/degree) <$> mgrsBandLetterToLatitude bandC+  northings0 <- (baseNorthingsOffset +) <$> mgrsNorthingIx northingsC+  let bandDist = band * metersPerDegree -- Approx dist from equator to southern edge of band.+      bandGridLower = floor $ bandDist / 100_000 - 2  -- Lower limit of band in 100km units+      bandGridUpper = ceiling $ if band > 71 * degree  -- Upper limit of band in 100km units.+        then (bandDist + 12 * metersPerDegree) / 100_000 + 1  -- Band X.+        else (bandDist + 8 * metersPerDegree) / 100_000 + 2  -- Other bands.+      rep = (bandGridLower - northings0 - 1) `div` 20  -- Lower limit in 2,000,000km units.+      grid = (rep+1)*20 + northings0+  guard $ grid >= bandGridLower+  guard $ grid <= bandGridUpper+  return $ fromIntegral $ grid * 100_000+  where+    metersPerDegree = 10_002_000 / 90  -- Equator to north pole.+    baseNorthingsOffset :: Int+    baseNorthingsOffset = if odd zone then 0 else -5+++-- | Find the northings letter of the 100km square containing the given Northings.+--+-- The input is not range checked. It just assumes that the northings letters repeat forever.+mgrsNorthingToLetter :: UtmZoneNumber -> Double -> Char+mgrsNorthingToLetter zone northings1 =+  letters ! ((gridNum + baseNorthingsOffset) `mod` 20)+  where+    gridNum :: Int+    gridNum = floor $ northings1 / (100 * kilometer)+    baseNorthingsOffset = if odd zone then 0 else 5+    letters = listArray (0,19) mgrsNorthingsLetters+++polarEastingsToLetter :: Double -> Char+polarEastingsToLetter eastings1 =+  letters ! ((floor (eastings1 / (100 * kilometer)) - 20) `mod` 18 :: Int)+  where+    letters = listArray (0,17) polarEastingLetters+++-- | The southern edge of the @A@ northings letter differs between north and south poles.+polarNorthingLetterOrigin :: Pole -> Double+polarNorthingLetterOrigin NorthPole = 1_300 * kilometer+polarNorthingLetterOrigin SouthPole =   800 * kilometer+++polarNorthingsToLetter :: Pole -> Double -> Maybe Char+polarNorthingsToLetter pole northings1 = do+  let i :: Int+      i = floor ((northings1 - polarNorthingLetterOrigin pole) / (100 * kilometer)) `mod` 24+  guard $ i >= 0 && i < 24+  pure $ letters ! i+  where+    letters = listArray (0,23) polarNorthingLetters++-- | Convert an MGRS grid reference to a UTM @GridPoint@, if the reference is valid.+-- E.g. \"30U XC 99304 10208\" is the grid reference for Nelson's Column in London.+-- +-- If the input contains spaces then these are used to delimit the fields. Multiple+-- spaces are treated as a single space.+--+-- If the reference is valid this returns the position of the south-west corner of the+-- nominated grid square and an offset to its centre. Altitude is set to zero.+fromMgrsGridReference :: String -> Either [String] (GridPoint MgrsGrid, GridOffset)+fromMgrsGridReference str = case parse parseMgrsGridReference "" str of+    Left err -> Left $ filter (not . null) $ lines $ showErrorMessages+      "or" "unknown parse error" "expecting" "unexpected" "end of input"+      (errorMessages err)+    Right r -> Right r++parseMgrsGridReference :: Parsec String u (GridPoint MgrsGrid, GridOffset)+parseMgrsGridReference = mgrsUtmParser <|> mgrsUpsParser+  where+    mgrsUtmParser = do+      zoneNum <- read <$> many1 digit <?> "UTM zone number"+          -- Safe because we can only read digits+      spaces+      band <- oneOf mgrsBandLetters <?> "latitude band letter"+      let (hemi, falseNorthing) = if band >= 'N'+            then (UtmNorth, 0)+            else (UtmSouth, -10_000_000)+      zone <- maybe (fail "Invalid zone") (pure . MgrsUtm) $ mkUtmZone hemi zoneNum+      spaces+      squareEast <- oneOf mgrsEastingsLetters <?> "eastings letter"+      eastingBase <- maybe (fail "Invalid eastings letter") pure $+        mgrsLetterToEasting zoneNum squareEast+      spaces+      squareNorth <- oneOf mgrsNorthingsLetters <?> "northings letter"+      northingBase <- maybe (fail "Invalid northings letter") pure $+        mgrsLetterToNorthings zoneNum band squareNorth+      spaces+      (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits+      when (length eastingChars /= length northingChars) $+        fail "Northings and Eastings must be the same length."+      if northingChars == "" && eastingChars == ""+        then -- No digits, just return the outer 100km grid square+          pure (GridPoint eastingBase (northingBase - falseNorthing) 0 zone,+                  GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0)+        else do+          (northing, offset) <- maybe (fail "Invalid northing digits") pure $+            fromGridDigits (100 * kilometer) northingChars+          (easting, _) <- maybe (fail "Invalid easting digits") pure $+            fromGridDigits (100 * kilometer) eastingChars+          pure (GridPoint (eastingBase + easting) (northingBase + northing - falseNorthing) 0 zone,+                  GridOffset (offset/2) (offset/2) 0)++    mgrsUpsParser = do+      spaces+      band <- oneOf polarBandLetters <?> "polar band letter"+      spaces+      squareEast <- oneOf polarEastingLetters <?> "eastings letter"+      spaces+      squareNorth <- oneOf polarNorthingLetters <?> "northings letter"+      spaces+      (pole, eastingBase, northingBase) <- mgrsLettersToPolar band squareEast squareNorth+      let grid = case pole of+            NorthPole -> MgrsUps upsNorth+            SouthPole -> MgrsUps upsSouth+      (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits+      when (length eastingChars /= length northingChars) $+        fail "Northings and Eastings must be the same length."+      if northingChars == "" && eastingChars == ""+        then -- No digits, just return the outer 100km grid square+          pure (GridPoint eastingBase northingBase 0 grid,+                  GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0)+        else do+          (northing, offset) <- maybe (fail "Invalid northing digits") pure $+            fromGridDigits (100 * kilometer) northingChars+          (easting, _) <- maybe (fail "Invalid easting digits") pure $+            fromGridDigits (100 * kilometer) eastingChars+          pure (GridPoint (eastingBase + easting) (northingBase + northing) 0 grid,+                  GridOffset (offset/2) (offset/2) 0)+    spaced = do+      e <- many1 digit+      skipMany1 space  -- A space is mandatory here.+      n <- many1 digit+      pure (e,n)+    unspaced = do+      digits <- many1 digit+      let c = length digits+      when (odd c) $ fail "Northings and Eastings must be the same length."+      pure (splitAt (c `div` 2) digits)+    noDigits = do+      eof+      pure ("", "")++-- | Convert a UTM or UPS @GridPoint@ to an MGRS grid reference.+toMgrsGridReference ::+  Bool  -- ^ Include spaces in the output. The standard says no spaces, but they make +        -- the output easier to read.+  -> Int  -- ^ Number of digits of precision in the easting and northing. Must be 0-8.+  -> GridPoint MgrsGrid+  -> Maybe String+toMgrsGridReference withSpaces precision gp = do+  guard $ precision >= 0 && precision <= 8+  case gridBasis gp of+    MgrsUtm zone -> do+      band <- mgrsLatitudeToBandLetter $ latitude $ fromGrid gp+      let+        zoneNum = utmZoneNum zone+        northLetter = mgrsNorthingToLetter zoneNum $ northings gp+      eastLetter <- mgrsEastingToLetter zoneNum $ eastings gp+      (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp+      (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp+      let part1 = printf "%02d" zoneNum <> [band]+          part2 = [eastLetter, northLetter]+      pure $ if withSpaces+        then part1 <> " " <> part2 <> " " <> eastDigits <> " " <> northDigits+        else part1 <> part2 <> eastDigits <> northDigits+    MgrsUps grid -> do+      let zoneLetter = case trueOrigin grid of+            NorthPole -> if eastings gp < 2_000 * kilometer then 'Y' else 'Z'+            SouthPole -> if eastings gp < 2_000 * kilometer then 'A' else 'B'+          eastLetter = polarEastingsToLetter $ eastings gp+      northLetter <- polarNorthingsToLetter (trueOrigin grid) $ northings gp+      (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp+      (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp+      pure $ if withSpaces+        then zoneLetter : ' ' : eastLetter : northLetter : " " <> eastDigits <> " " <> northDigits+        else zoneLetter : eastLetter : northLetter : eastDigits <> northDigits

@@ -1,4 +1,4 @@-{-# LANGUAGE TypeOperators, TypeFamilies, FlexibleContexts #-}+{-# LANGUAGE TypeFamilies, FlexibleContexts #-} -- | The implementation assumes IEEE 754 arithmetic.  module Geodetics.Path where@@ -26,7 +26,7 @@          -- ^ Takes a length and returns a position, and direction as heading and elevation angles.       pathValidity :: PathValidity    }-   + -- | Convenience value for paths that are valid for all distances. alwaysValid :: PathValidity alwaysValid = (negate inf, inf) where@@ -46,8 +46,8 @@ --  -- The initial bounds must return one GT or EQ value and one LT or EQ value. If they -- do not then @Nothing@ is returned.-bisect :: -   Path e +bisect ::+   Path e    -> (Geodetic e -> Ordering)        -- ^ Evaluation function.    -> Double                          -- ^ Required accuracy in terms of distance along the path.    -> Double -> Double                -- ^ Initial bounds.@@ -59,7 +59,7 @@       guard $ hasRoot r       bisect1 $ sortPair r    where-      f' d = let (p, _, _) = pathFunc path d in f p +      f' d = let (p, _, _) = pathFunc path d in f p       pairResults d1 d2 = ((d1, f' d1), (d2, f' d2))       hasRoot (v1, v2) = snd v1 <= EQ && EQ <= snd v2       sortPair (v1, v2) = if snd v1 <= snd v2 then (v1, v2) else (v2, v1)@@ -68,7 +68,7 @@              r3 = f' d3              c1 = ((d1, r1), (d3, r3))              c2 = ((d3, r3), (d2, r2))-         in if abs (d1 - d2) <= t +         in if abs (d1 - d2) <= t             then return d3             else bisect1 $ if hasRoot c1 then c1 else c2 @@ -96,14 +96,14 @@    | mag < 1e-15                    = Nothing    | mag3 (nv1 `cross3` nv2) * r <= accuracy = Just (d1, d2)        -- Assumes that sin (accuracy/r) == accuracy/r-   | otherwise = +   | otherwise =       if abs d1a + abs d2a < abs d1b + abs d2b          then intersect (d1 + d1a) (d2 + d2a) accuracy (pred n) path1 path2          else intersect (d1 + d1b) (d2 + d2b) accuracy (pred n) path1 path2    where       (pt1, h1, _) = pathFunc path1 d1       (pt2, h2, _) = pathFunc path2 d2-      vectors :: Double -> Double -> Double +      vectors :: Double -> Double -> Double                  -> (Vec3 Double, Vec3 Double)       vectors lat lon b = (           -- Unit vector of normal to surface at (lat,lon)@@ -132,11 +132,11 @@       d1b = gcDist gc1 nv1 nv3b * r       d2b = gcDist gc2 nv2 nv3b * r       -- Signed angle between v1 and v2, -      gcDist norm v1 v2 = -         let c = v1 `cross3` v2 -         in (if c `dot3` norm < 0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) +      gcDist norm v1 v2 =+         let c = v1 `cross3` v2+         in (if c `dot3` norm < 0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2)       r = majorRadius $ ellipsoid pt1-          + {- Note on derivation of "intersect"  The algorithm is a variant of the Newton-Raphson method, and shares its advantage@@ -167,7 +167,7 @@ -}  -- | A ray from a point heading in a straight line in 3 dimensions. -rayPath :: (Ellipsoid e) => +rayPath :: (Ellipsoid e) =>    Geodetic e          -- ^ Start point.    -> Double           -- ^ Bearing.    -> Double           -- ^ Elevation.@@ -181,7 +181,7 @@             (dE,dN,dU) = transform3 (trans3 $ ecefMatrix lat long) delta  -- Direction of ray at result point.             elevation2 = asin dU             bearing2 = if dE == 0 && dN == 0 then bearing else atan2 dE dN  -- Allow for vertical elevation.-            +       ecefMatrix lat long =   -- Transform matrix for vectors from (East, North, Up) to (X,Y,Z).          ((negate sinLong, negate cosLong*sinLat, cosLong*cosLat),               --    East X      North X               Up X@@ -194,7 +194,7 @@             cosLong = cos long             sinLat = sin lat             cosLat = cos lat-      +       direction = (sinB*cosE, cosB*cosE, sinE)  -- Direction of ray in ENU       delta = transform3 (ecefMatrix (latitude pt1) (longitude pt1)) direction  -- Convert to ECEF       pt1' = geoToEarth pt1    -- ECEF of origin point.@@ -202,7 +202,7 @@       cosB = cos bearing       sinE = sin elevation       cosE = cos elevation-      + -- | Rhumb line: path following a constant course. Also known as a loxodrome. -- -- The valid range stops a few arc-minutes short of the poles to ensure that the @@ -244,8 +244,8 @@       m0 = meridianRadius (ellipsoid pt) lat0       a = majorRadius $ ellipsoid pt       b = minorRadius $ ellipsoid pt-    + -- | A path following the line of latitude around the Earth eastwards. -- -- This is equivalent to @rhumbPath pt (pi/2)@@@ -256,7 +256,7 @@    where       line distance = (pt2, pi/2, 0)          where-            pt2 = Geodetic +            pt2 = Geodetic                (latitude pt) (longitude pt + distance / r)                0 (ellipsoid pt)       r = latitudeRadius (ellipsoid pt) (latitude pt)

@@ -0,0 +1,220 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FlexibleContexts #-}++module Geodetics.PolarStereographic (+  Pole (..),+  PolarStereographic (trueOrigin, falseOrigin, polarEllipsoid, gridScale),+  mkGridPolarStereographic,+  UpsGrid,+  upsNorth,+  upsSouth,+  fromUpsGridReference,+  parseUpsGridReference,+  toUpsGridReference+) where++import Control.Monad+import Data.Char+import Geodetics.Ellipsoids+import Geodetics.Geodetic+import Geodetics.Grid+import Text.Parsec+import Text.Parsec.Error+import Text.Read (readMaybe)+import Text.Printf++-- | Polar stereographic grids are defined for true origins at the north and south poles.+data Pole = NorthPole | SouthPole deriving (Show, Ord, Eq, Enum, Bounded)+++{- | Polar Stereographic Grids++Formulae are taken from+/The Universal Grids: Univerersal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)/+DMA Technical Manual 8358.2, Defense Mapping Agency, Fairfax, VA. https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf++When working with polar grids all directions are relative to the grid rather than the actual pole.+So in the Arctic \"North\" on the Universal Polar Stereographic grid means towards the Bering Sea+rather than towards the North Pole.+-}+data PolarStereographic e = PolarStereographic {+  trueOrigin :: Pole,+  falseOrigin :: GridOffset,+    -- ^ The negation of 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.+  polarEllipsoid :: e,+    -- ^ The ellipsoid for the projection. Arguments passed to `toGrid` *must* use this ellipsoid.+    -- The type system cannot verify this for `LocalEllipsoid`.+  gridScale :: Double,+    -- ^ The scaling factor applied at the pole. This balances the distortion between the center+    -- and edges of the projection.++  -- Remaining elements are memoised parameters computed from the ellipsoid.+  gridA, gridB, gridC, gridD, gridC0 :: !Double+} deriving (Show)++instance (Eq e) => Eq (PolarStereographic e) where+  g1 == g2 =+    trueOrigin g1 == trueOrigin g2 &&+    falseOrigin g1 == falseOrigin g2 &&+    polarEllipsoid g1 == polarEllipsoid g2 &&+    gridScale g1 == gridScale g2++instance (Ellipsoid e) => GridClass (PolarStereographic e) e where+  fromGrid p = Geodetic lat long (altGP p) (polarEllipsoid gb)+    where+      gridZero = GridPoint 0 0 0 gb+      gb = gridBasis p+      p' = gridZero `gridOffset` (falseOrigin gb `applyOffset` p)+      radius = offsetDistance p'+      isoColat = 2 * atan (radius / (gridScale gb * gridC0 gb))+      isoLat = pi/2 - isoColat+      lat1 = isoLat ++        gridA gb * sin (2*isoLat) ++        gridB gb * sin (4*isoLat) ++        gridC gb * sin (6*isoLat) ++        gridD gb * sin (8*isoLat)+      lat = case trueOrigin gb of+        NorthPole -> lat1+        SouthPole -> negate lat1+      long = case trueOrigin gb of+        NorthPole -> offsetBearing p' { deltaNorth = negate $ deltaNorth p'}+        SouthPole -> offsetBearing p'++  toGrid r geo = offsetNegate (falseOrigin r) `applyOffset` GridPoint east north 0 r+    where+      absLat = abs $ latitude geo+      e = sqrt (eccentricity2 $ polarEllipsoid r)+      eSinLat = e * sin absLat+      tz2 = ((1 + eSinLat)/(1-eSinLat))**(e/2) * tan (pi/4 - absLat / 2)+      radius = gridScale r * gridC0 r * tz2+      north = case trueOrigin r of+        NorthPole -> negate $ radius * cos (longitude geo)+        SouthPole -> radius * cos (longitude geo)+      east = radius * sin (longitude geo)+  gridEllipsoid = polarEllipsoid+++mkGridPolarStereographic :: (Ellipsoid e) =>+  Pole          -- ^ True origin at north or south pole.+  -> e          -- ^ The ellipsoid used for the projection.+  -> GridOffset -- ^ Vector from true origin to the false origin.+  -> Double     -- ^ Scale factor.+  -> PolarStereographic e+mkGridPolarStereographic pole ellip offset scale =+  PolarStereographic {+    trueOrigin = pole,+    falseOrigin = offset,+    polarEllipsoid = ellip,+    gridScale = scale,+    gridA = e2/2 + (5/24)*e4 + e6/12 + (13/360)*e8,+    gridB = (7/48)*e4 + (29/240)*e6 + (811/11520)*e8,+    gridC = (7/120)*e6 + (81/1120)*e8,+    gridD = (4279/161280)*e8,+    gridC0 = (2 * majorRadius ellip / sqrt (1 - e2)) * ((1-e1)/(1+e1))**(e1/2)+  }+  where+    e1 = sqrt $ eccentricity2 ellip+    e2 = eccentricity2 ellip+    e4 = e2^_2+    e6 = e2^_3+    e8 = e2^_4+++-- | The Universal Polar Stereographic (UPS) grids for north and south poles.+type UpsGrid = PolarStereographic WGS84+++-- | UPS grid for the North Pole.+upsNorth :: UpsGrid+upsNorth = mkGridPolarStereographic+  NorthPole+  WGS84+  (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 })+  0.994  -- Scale factor+++-- | UPS grid for the South Pole.+upsSouth :: UpsGrid+upsSouth = mkGridPolarStereographic+  SouthPole+  WGS84+  (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 })+  0.994   -- Scale factor+++-- | Convert a grid reference into a UPS grid location.+--+-- There doesn't appear to be any conventional representation for polar grid references,+-- so this is an attempt to cover as many bases as possible. It takes an Easting followed+-- by a Northing with spaces in between. Both can have optional units of m or km,+-- and be optionally followed by an \"N\" or \"E\" as appropriate.+--+-- The choice of pole is provided in an extra argument rather than within the string because humans+-- will normally assume this from the context and so not provide it.+--+-- If the string cannot be parsed then one or more error messages are returned.+fromUpsGridReference :: Pole -> String -> Either [String] (GridPoint UpsGrid)+fromUpsGridReference pole str = case parse (parseUpsGridReference pole) "" str of+    Left err -> Left $ lines $ showErrorMessages+      "or" "unknown parse error" "expecting" "unexpected" "end of input"+      (errorMessages err)+    Right r -> Right r+++parseUpsGridReference :: Stream s m Char => Pole -> ParsecT s u m (GridPoint UpsGrid)+parseUpsGridReference pole = do+      spaces1+      (eastings1, eastUnit) <- readDistance+      spaces+      optional (oneOf "Ee" <?> "E")+      spaces1+      (northings1, northUnit) <- readDistance+      unless (eastUnit == northUnit) $ fail "Northings and Eastings units don't match."+      spaces1+      optional (oneOf "Nn" <?> "N")+      spaces1+      eof+      return $ GridPoint eastings1 northings1 0 $+        case pole of+          NorthPole -> upsNorth+          SouthPole -> upsSouth+  where+    readDistance = do  -- Returns (distance in meters, unit from input)+      digits <- many1 (digit <|> char '.' <?> "number")+      spaces1+      when (length digits > 20) $ fail "Too many digits."+      (multiplier, unit) <- do+        unit <- optionMaybe (string1' "m" <|> string1' "km" <?> "units (m or km)")+        case unit of+          Just "km" -> pure (1000, Just GridKilometers)+          Just _    -> pure (1, Just GridMeters)+          Nothing   -> pure (1, Nothing)+      case readMaybe digits of+        Just d -> return (d * multiplier, unit)+        Nothing -> fail $ "Cannot read number: " <> digits+    string1' target = try $ do  -- Case-insensitive version of string'+      cs <- count (length target) anyToken+      if map toLower target == map toLower cs then return cs else unexpected cs+    spaces1 = void $ many (char ' ' <?> "space")  -- Other white space not permitted.+++toUpsGridReference :: +  Maybe GridUnit  -- ^ Include explicit units in the output. @Nothing@ means meters without units.+  -> Bool  -- ^ Include \"E\" and \"N\" in the output.+  -> Int  -- ^ Digits of resolution. 0 = 1m resolution, 1 = 10m, 2 = 100m etc. (-2) = 1cm.+  -> GridPoint UpsGrid+  -> String+toUpsGridReference unit letters res gp =+  dist (eastings gp)  <> (if letters then "E " else " ") <>+    dist (northings gp) <> (if letters then "N" else "")+  where+    res1 :: Double+    res1 = 10 ** fromIntegral res   -- Resolution in meters.+    floorRes :: Double -> Double+    floorRes d = res1 * fromIntegral (floor (d/res1) :: Integer)+    dist d = case unit of+      Nothing            -> printf "%.*f" (-res) $ floorRes d+      Just GridMeters     -> printf "%.*fm" (-res) $ floorRes d+      Just GridKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000

@@ -15,6 +15,9 @@ -- 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). --+-- Arguments passed to `toGrid` *must* use the same ellipsoid as the `trueOrigin`. The type system+-- cannot verify this for 'LocalEllipsoid'.+-- -- The calculations here are based on \"Transverse Mercator Projection: Constants, Formulae and Methods\" -- by the Ordnance Survey, March 1983. -- Retrieved from http://www.threelittlemaids.co.uk/magdec/transverse_mercator_projection.pdf@@ -22,14 +25,14 @@    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 negation of 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 :: 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 :: Double+   gridN1, gridN2, gridN3, gridN4 :: !Double } deriving (Show)  

@@ -41,15 +41,11 @@   mkUtmZone,   mkUtmZoneUnsafe,   fromUtmGridReference,-  toUtmGridReference,-  mgrsBandLetterToLatitude,-  mgrsLatitudeToBandLetter,-  fromMgrsGridReference,-  toMgrsGridReference+  parseUtmGridReference,+  toUtmGridReference ) where -import Control.Monad (mplus, guard, void, when, unless)-import Data.Array+import Control.Monad (guard, void, when, unless) import Data.Char import Data.List import Geodetics.Ellipsoids@@ -155,10 +151,6 @@     scale = 0.999_6  --- | Units for UTM grid coordinates.-data UtmGridUnit = UtmMeters | UtmKilometers deriving (Eq, Show)-- -- | Convert a grid reference to a position, if the reference is valid. -- -- The northings and eastings cannot contain more than 20 digits each,@@ -172,13 +164,15 @@ -- -- If the argument cannot be parsed then one or more error messages are returned. fromUtmGridReference :: String -> Either [String] (GridPoint UtmZone)-fromUtmGridReference str = case parse gridP str str of+fromUtmGridReference str = case parse parseUtmGridReference str str of     Left err -> Left $ lines $ showErrorMessages       "or" "unknown parse error" "expecting" "unexpected" "end of input"       (errorMessages err)     Right r -> Right r-  where-    gridP = do+++parseUtmGridReference :: Stream s m Char => ParsecT s u m (GridPoint UtmZone)+parseUtmGridReference = do       spaces1       zone <- readZone <?> "Zone number"       hemi <- readHemi <?> "Hemisphere (N or S)"@@ -193,8 +187,9 @@       optional (oneOf "Nn" <?> "N")       spaces1       eof-      return $ GridPoint eastings1 northings1 0 $ mkUtmZoneUnsafe hemi zone-    readZone :: Parsec String () UtmZoneNumber+      pure $ GridPoint eastings1 northings1 0 $ mkUtmZoneUnsafe hemi zone+  where+    readZone :: Stream s m Char => ParsecT s u m UtmZoneNumber     readZone = do       ds <- many1 digit       case readMaybe ds of@@ -202,28 +197,32 @@         Just n ->           if n < 1 || n > 60             then fail $ "Zone number " <> show n <> " out of range."-            else return n-    readHemi :: Parsec String () UtmHemisphere+            else pure n+    readHemi :: Stream s m Char => ParsecT s u m UtmHemisphere     readHemi = do       h <- oneOf "NSns"       case toUpper h of-        'N' -> return UtmNorth-        'S' -> return UtmSouth+        'N' -> pure UtmNorth+        'S' -> pure UtmSouth         _ -> fail $ "Invalid hemisphere: " <> (h : ". Must be N or S.")-    readDistance :: Parsec String () (Double, UtmGridUnit)  -- (Distance, unit)+    readDistance :: Stream s m Char => ParsecT s u m (Double, Maybe GridUnit)  -- (Distance, unit)     readDistance = do       digits <- many1 (digit <|> char '.' <?> "number")       spaces1       when (length digits > 20) $ fail "Too many digits."       (multiplier, unit) <- do-        unit <- option "m" (string1' "m" <|> string1' "km" <?> "units (m or km)")-        if unit == "km" then return (1000, UtmKilometers) else return (1, UtmMeters)+        unit <- optionMaybe (string1' "m" <|> string1' "km" <?> "units (m or km)")+        case unit of+          Just "km" -> pure (1000, Just GridKilometers)+          Just _    -> pure (1, Just GridMeters)+          Nothing   -> pure (1, Nothing)       case readMaybe digits of-        Just d -> return (d * multiplier, unit)+        Just d -> pure (d * multiplier, unit)         Nothing -> fail $ "Cannot read number: " <> digits     string1' target = try $ do  -- Case-insensitive version of string'       cs <- count (length target) anyToken-      if map toLower target == map toLower cs then return cs else unexpected cs+      if map toLower target == map toLower cs then pure cs else unexpected cs+    spaces1 :: Stream s m Char => ParsecT s u m ()     spaces1 = void $ many (char ' ' <?> "space")  -- Other white space not permitted.  @@ -231,7 +230,7 @@ -- The northings and eastings are rounded down to the resolution, so the result is the south-west -- corner of the grid square enclosing the grid point. toUtmGridReference ::-  Maybe UtmGridUnit  -- ^ Include explicit units in the output. @Nothing@ means meters without units.+  Maybe GridUnit  -- ^ Include explicit units in the output. @Nothing@ means meters without units.   -> Bool -- ^ Include \"E\" and \"N\" in the output.   -> Int  -- ^ Digits of resolution. 0 = 1m resolution, 1 = 10m, 2 = 100m etc. (-2) = 1cm.   -> GridPoint UtmZone@@ -249,216 +248,5 @@     zoneStr = printf "%02d" (utmZoneNum b) <> show (utmHemisphere b)     dist d = case unit of       Nothing            -> printf "%.*f" (-res) $ floorRes d-      Just UtmMeters     -> printf "%.*fm" (-res) $ floorRes d-      Just UtmKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000------ | The MGRS latitude band code letters, excluding A and B used for Antarctica (south of -80 degrees)--- and Y and Z used for the Arctic (north of 84 degrees).-mgrsBandLetters :: [Char]-mgrsBandLetters = "CDEFGHJKLMNPQRSTUVWX"----- | Find the southern boundary of a latitude band letter.-mgrsBandLetterToLatitude :: Char -> Maybe Double-mgrsBandLetterToLatitude band = do-    n1 <- ix band-    return $ degree * fromIntegral (-80 + n1 * 8)-  where-    indexMap :: Array Char (Maybe Int)-    indexMap = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip mgrsBandLetters [0..]]-    ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing----- | Find the band letter for a latitude, if it is in the range (-80, 84) degrees.--- (Argument in radians)-mgrsLatitudeToBandLetter :: Double -> Maybe Char-mgrsLatitudeToBandLetter lat = do-    guard $ -80 <= dlat && dlat <= 84-    return $ indexMap ! latIdx-  where-    dlat = lat / degree-    ilat :: Int-    ilat = floor dlat-    latIdx = min 19 $ (ilat + 80) `div` 8  -- Band 19 (X) extends an extra 4 degrees.-    indexMap = listArray (0,19) mgrsBandLetters------ | Letters A-Z except for I and O.-mgrsEastingsLetters :: [Char]-mgrsEastingsLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ"---- | If zone number is in range and the letter is one of the valid Eastings letters for that zone--- then return the UTM easting in meters.------ Zone 1 starts with \'A\'. Each zone is 8 characters wide. Hence the letters repeat every 3 zones.-mgrsLetterToEasting :: UtmZoneNumber -> Char -> Maybe Double-mgrsLetterToEasting zn c = do-    guard $ 1 <= zn && zn <= 60-    n1 <- ix c-    let n2 = n1 - base-    guard $ 0 <= n2 && n2 <= 7-    return $ fromIntegral (n2+1) * 100_000-  where-    indexMap = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip mgrsEastingsLetters [0..]]-    base = ((zn-1) `mod` 3) * 8-    ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing----- | If the zone number is in range and the eastings are between 100,000 and 900,000 then--- return the Eastings letter.-mgrsEastingToLetter :: UtmZoneNumber -> Double -> Maybe Char-mgrsEastingToLetter zn east = do-    guard $ 1 <= zn && zn <= 60-    guard $ 100 * kilometer <= east && east < 900 * kilometer-    return $ indexMap ! ix-  where-    indexMap = listArray (0,23) mgrsEastingsLetters-    base = ((zn-1) `mod` 3) * 8-    square = max 0 $ min 7 $ floor $ (east - 100 * kilometer)/(100 * kilometer)  -- Clamped in range (0,7).-    ix = base + square  -- Must be in range (0,23)----- | Letters A-V except for I and O.-mgrsNorthingsLetters :: [Char]-mgrsNorthingsLetters = "ABCDEFGHJKLMNPQRSTUV"----- | MGRS Northings letters have rather complex relationship to the latitude bands. The 20 letters--- repeat every 2,000km going north and south from the equator, so the latitude band is needed--- to disambiguate which repetition of the Northings letter is meant.---- Unfortunately this repetition of the letters does not neatly coincide with--- the latitude band boundaries, which are based on degrees of latitude.---- The base letter just north of the equator is A in odd-numbered zones and F in even numbered zones.------ This uses the latitude band to estimate the range of northings that would be valid there,--- and hence determine which possible grid band is meant by the northings letter.--- The algorithm used is approximate and deliberately very forgiving: it will accept some grid squares--- which are north or south of the band given.-mgrsLetterToNorthings ::-  UtmZoneNumber-  -> Char  -- ^ Latitude band letter (@C@ - @X@ excluding @I@ and @O@).-  -> Char  -- ^ MGRS Northings letter (@A@ - @V@ excluding @I@ and @O@).-  -> Maybe Double-mgrsLetterToNorthings zone bandC northingsC = do-  guard $ 1 <= zone && zone <= 60-  band <- (/degree) <$> mgrsBandLetterToLatitude bandC-  northings0 <- (baseNorthingsOffset +) <$> ix northingsC-  let bandDist = band * metersPerDegree -- Approx dist from equator to southern edge of band.-      bandGridLower = floor $ bandDist / 100_000 - 2  -- Lower limit of band in 100km units-      bandGridUpper = ceiling $ if band > 71 * degree  -- Upper limit of band in 100km units.-        then (bandDist + 12 * metersPerDegree) / 100_000 + 1  -- Band X.-        else (bandDist + 8 * metersPerDegree) / 100_000 + 2  -- Other bands.-      rep = (bandGridLower - northings0 - 1) `div` 20  -- Lower limit in 2,000,000km units.-      grid = (rep+1)*20 + northings0-  guard $ grid >= bandGridLower-  guard $ grid <= bandGridUpper-  return $ fromIntegral $ grid * 100_000-  where-    metersPerDegree = 10_002_000 / 90  -- Equator to north pole.-    baseNorthingsOffset :: Int-    baseNorthingsOffset = if odd zone then 0 else -5-    indexMap :: Array Char (Maybe Int)-    indexMap = accumArray mplus Nothing ('A', 'Z') [(c, Just n) | (c, n) <- zip mgrsEastingsLetters [0..]]-    ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing----- | Find the northings letter of the 100km square containing the given Northings.------ The input is not range checked. It just assumes that the northings letters repeat forever.-mgrsNorthingToLetter :: UtmZoneNumber -> Double -> Char-mgrsNorthingToLetter zone northings1 =-  letters ! ((gridNum + baseNorthingsOffset) `mod` 20)-  where-    gridNum :: Int-    gridNum = floor $ northings1 / (100 * kilometer)-    baseNorthingsOffset = if odd zone then 0 else 5-    letters = listArray (0,19) mgrsNorthingsLetters----- | Convert an MGRS grid reference to a UTM @GridPoint@, if the reference is valid.--- E.g. \"30U XC 99304 10208\" is the grid reference for Nelson's Column in London.--- --- If the input contains spaces then these are used to delimit the fields. Any or all spaces--- may be omitted. Multiple spaces are treated as a single space.------ If the reference is valid this returns the position of the south-west corner of the--- nominated grid square and an offset to its centre. Altitude is set to zero.-fromMgrsGridReference :: String -> Either [String] (GridPoint UtmZone, GridOffset)-fromMgrsGridReference str = case parse mgrsP str str of-    Left err -> Left $ filter (not . null) $ lines $ showErrorMessages-      "or" "unknown parse error" "expecting" "unexpected" "end of input"-      (errorMessages err)-    Right r -> Right r-  where-    mgrsP = do-      zoneNum <- read <$> many1 digit  -- Safe because we can only read digits here.-      spaces-      band <- upper <?> "latitude band letter"-      let (hemi, falseNorthing) = if band >= 'N'-            then (UtmNorth, 0)-            else (UtmSouth, -10_000_000)-      zone <- maybe (fail "Invalid zone") return $ mkUtmZone hemi zoneNum-      spaces-      squareEast <- upper <?> "eastings letter"-      eastingBase <- maybe (fail "Invalid eastings letter") return $-        mgrsLetterToEasting zoneNum squareEast-      spaces-      squareNorth <- upper <?> "northings letter"-      northingBase <- maybe (fail "Invalid northings letter") return $-        mgrsLetterToNorthings zoneNum band squareNorth-      spaces-      (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits-      when (length eastingChars /= length northingChars) $-        fail "Northings and Eastings must be the same length."-      if northingChars == "" && eastingChars == ""-        then -- No digits, just return the outer 100km grid square-          return (GridPoint eastingBase (northingBase - falseNorthing) 0 zone,-                  GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0)-        else do-          (northing, offset) <- maybe (fail "Invalid northing digits") return $-            fromGridDigits (100 * kilometer) northingChars-          (easting, _) <- maybe (fail "Invalid easting digits") return $-            fromGridDigits (100 * kilometer) eastingChars-          return (GridPoint (eastingBase + easting) (northingBase + northing - falseNorthing) 0 zone,-                  GridOffset (offset/2) (offset/2) 0)-    spaced = do-      e <- many1 digit-      skipMany1 space  -- A space is mandatory here.-      n <- many1 digit-      return (e,n)-    unspaced = do-      digits <- many1 digit-      let c = length digits-      when (odd c) $ fail "Northings and Eastings must be the same length."-      return (splitAt (c `div` 2) digits)-    noDigits = do-      eof-      return ("", "")---- | Convert UTM @GridPoint@ to an MGRS grid reference.-toMgrsGridReference ::-  Bool  -- ^ Include spaces in the output. The standard says no spaces, but they make -        -- the output easier to read.-  -> Int  -- ^ Number of digits of precision in the easting and northing. Must be 0-5.-  -> GridPoint UtmZone-  -> Maybe String-toMgrsGridReference withSpaces precision gp = do-  guard $ precision >= 0 && precision <= 5-  band <- mgrsLatitudeToBandLetter $ latitude $ fromGrid gp-  let-    zoneNum = utmZoneNum $ gridBasis gp-    northLetter = mgrsNorthingToLetter zoneNum $ northings gp-  eastLetter <- mgrsEastingToLetter zoneNum $ eastings gp-  (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp-  (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp-  let part1 = printf "%02d" zoneNum <> [band]-      part2 = [eastLetter, northLetter]-  return $ if withSpaces-    then part1 <> " " <> part2 <> " " <> eastDigits <> " " <> northDigits-    else part1 <> part2 <> eastDigits <> northDigits-+      Just GridMeters     -> printf "%.*fm" (-res) $ floorRes d+      Just GridKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000

@@ -1,6 +1,7 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} {-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} {-# HLINT ignore "Redundant bracket" #-}+{-# LANGUAGE ViewPatterns #-}  module Main where @@ -20,12 +21,14 @@ import Geodetics.Ellipsoids import Geodetics.Geodetic import Geodetics.Grid+import Geodetics.MGRS import Geodetics.Path import Geodetics.Stereographic import Geodetics.TransverseMercator import Geodetics.UK import Geodetics.UTM import LatLongParser (parserTests)+import Geodetics.PolarStereographic as PS  main :: IO () main = hspec $ do@@ -58,6 +61,12 @@       prop "MGRS Grid 2" prop_mgrs_gridTest2       describe "MGRS Grid 3" $ mapM_ mgrsGridTest3 utmSampleGrid       describe "MGRS Grid 4" $ mapM_ mgrsGridTest4 utmSampleGrid+   describe "UPS" $ do+      describe "UPS Grid 4" $ mapM_ upsGridTest4 upsSampleGrid+      describe "UPS Grid 5" $ mapM_ upsGridTest5 upsSampleGrid+      describe "UPS Grid 6" $ mapM_ upsGridTest6 upsSampleGrid+      describe "UPS Grid 7" $ mapM_ upsGridTest7 upsSampleGrid+      describe "UPS Grid 8" $ mapM_ upsGridTest8 upsSampleGrid    describe "Stereographic" $ do       it "toGrid north" stereographicToGridN       it "fromGrid north" stereographicFromGridN@@ -369,7 +378,7 @@ prop_mgrs_gridTest1 (UtmGridRef str) =    case fromUtmGridReference str of       Left msg -> assertFailure $ "gridTest1 bogus UTM ref: " <> str <> ". Messages = " <> show msg-      Right gp ->+      Right (utmToMgrsPoint -> gp) ->          fromMgrsGridReference <$> toMgrsGridReference True 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0))  -- | Check that a UTM grid point round-trips to MGRS and back, without spaces.@@ -377,19 +386,73 @@ prop_mgrs_gridTest2 (UtmGridRef str) =    case fromUtmGridReference str of       Left msg -> assertFailure $ "gridTest2 bogus UTM ref: " <> str <> ". Messages = " <> show msg-      Right gp ->+      Right (utmToMgrsPoint -> gp) ->          fromMgrsGridReference <$> toMgrsGridReference False 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0))  -- | Check that MGRS reference to grid point works for sample points. mgrsGridTest3 :: UtmGridPointTest-mgrsGridTest3 (_, mgrs, gp, _, testName) =+mgrsGridTest3 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) =    it testName $ fromMgrsGridReference mgrs `shouldBe` Right (gp, GridOffset 0.5 0.5 0)  -- | Check that grid point to MGRS reference works for sample points. mgrsGridTest4 :: UtmGridPointTest-mgrsGridTest4 (_, mgrs, gp, _, testName) = do+mgrsGridTest4 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) = do    it (testName <> " with spaces") $ toMgrsGridReference True 5 gp `shouldBe` Just mgrs    it (testName <> " without spaces") $ toMgrsGridReference False 5 gp `shouldBe` Just (filter (not . isSpace) mgrs)+++-- | DMA TM 8358.2 can be found at https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf+-- Oracle for other points is https://geographiclib.sourceforge.io/C++/doc/GeoConvert.1.html+upsSampleGrid :: [(String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String)]+upsSampleGrid = map convert [+    -- MGRS             UPS,               X,          Y,           Latitude,    Longitude+      ("ZAH0000000000", "2000000 2000000", 2000000.00, 2000000.00,  90,             0,         "North pole"),+      ("BAN0000000000", "2000000 2000000", 2000000.00, 2000000.00, -90,             0,         "South pole"),+      ("ZAE0000077930", "2000000 1777930", 2000000.00, 1777930.73,  88,             0,         "88N 0W"),+      ("YXH7793000000", "1777930 2000000", 1777930.73, 2000000.00,  88,           -90,         "88N, 90W"),+      ("BAQ0000022069", "2000000 2222069", 2000000.00, 2222069.27, -88,             0,         "88S 0W"),+      ("AXN7793000000", "1777930 2000000", 1777930.73, 2000000.00, -88,           -90,         "88S, 90W"),+      ("YTM3012526773", "1530125 2426773", 1530125.78, 2426773.60,  84.28723389, -132.2479892, "DMA TM 8358.2 sample 1"),+      ("BCK2297997474", "2222979 1797474", 2222979.47, 1797474.90, -87.28733333,  132.2478619, "DMA TM 8358.2 sample 3") ]+      -- DMA TM 8358.2 sample 2 is outside the UPS grid regions and is therefore not a valid test.+   where+      convert (mgrs, ups, x, y, lat, long, desc) = (mgrs, ups, GridPoint { eastings = x, northings = y, altGP = 0, gridBasis = gb }, geo, desc)+         where+            geo = Geodetic { latitude = (lat * degree), longitude = (long * degree), geoAlt = 0, ellipsoid = WGS84 }+            gb = mkGridPolarStereographic+               (if lat < 0 then SouthPole else NorthPole)+               WGS84+               (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 })+               0.994++type UpsGridPointTest = (String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String) -> SpecWith (Arg Expectation)++-- | Check that the UPS grid point to WGS84 works close enough for sample points.+upsGridTest4 :: UpsGridPointTest+upsGridTest4 (_, _, gp, geo, testName) =+   it testName $ geo `closeEnough` fromGrid gp++-- | Check that WGS84 to UPS grid point works close enough for sample points.+upsGridTest5 :: UpsGridPointTest+upsGridTest5 (_, _, gp, geo, testName) =+   it testName $ gp `closeGrid` toGrid (gridBasis gp) geo++-- | Check that UPS and MGRS parsers agree.+upsGridTest6 :: UpsGridPointTest+upsGridTest6 (mgrs, ups, gp, _, testName) =+   it testName $+         fst <$> fromMgrsGridReference mgrs `shouldBe`+         upsToMgrsPoint <$> fromUpsGridReference (PS.trueOrigin $ gridBasis gp) ups++-- | Check MGRS grid generator+upsGridTest7 :: UpsGridPointTest+upsGridTest7 (mgrs, _, gp, _, testName) =+   it testName $ toMgrsGridReference False 5 (upsToMgrsPoint gp) `shouldBe` Just mgrs++-- | Check UPS grid generator+upsGridTest8 :: UpsGridPointTest+upsGridTest8 (_, ups, gp, _, testName) =+   it testName $ toUpsGridReference Nothing False 0 gp `shouldBe` ups   -- | Standard stereographic grid for point tests in the Northern Hemisphere.