geodetics 1.0.0 → 1.1.0

raw patch · 13 files changed

+1175/−288 lines, 13 filesdep +hspecdep +parsecdep −test-frameworkdep −test-framework-hunitdep −test-framework-quickcheck2dep ~HUnitdep ~arrayPVP ok

version bump matches the API change (PVP)

Dependencies added: hspec, parsec

Dependencies removed: test-framework, test-framework-hunit, test-framework-quickcheck2

Dependency ranges changed: HUnit, array

API changes (from Hackage documentation)

Files

• ToDo.txt +1/−1
• changelog.md +18/−2
• geodetics.cabal +24/−12
• src/Geodetics/Ellipsoids.hs +1/−5
• src/Geodetics/Geodetic.hs +113/−31
• src/Geodetics/Grid.hs +6/−12
• src/Geodetics/LatLongParser.hs +26/−32
• src/Geodetics/Path.hs +1/−1
• src/Geodetics/UK.hs +6/−6
• src/Geodetics/UTM.hs +464/−0
• test/ArbitraryInstances.hs +71/−52
• test/LatLongParser.hs +209/−0
• test/Main.hs +235/−134

@@ -1,6 +1,6 @@ To Do ===== -* Add UTM and Universal Polar grids.+* Add Universal Polar grid.  * Add more national grid systems.

@@ -21,5 +21,21 @@ Version 0.1.2: Fixed bugs #16 and #17: Unicode PRIME and DOUBLE PRIME now allowed in    position strings, and the degree symbol is allowed for decimal degrees. -Version 1.0.0: Removed dependency on Dimensional library. This is a breaking change:-   hence the major version bump. Also fixed bug #18 (and #19).+## Version 1.0.0++* Removed dependency on Dimensional library. This is a breaking change:+   hence the major version bump.++* Fixed bug #18 (and #19).++## Version 1.1.0++* Dropped a redundant `Ellipsoid` constraint on `antipode`.++* Added functions `showGeodeticLatLong`, `showGeodeticSignedDecimal`,+  `showGeodeticNSEWDecimal`, and `showGeodeticDDDMMSS`.++* Fixed a few minor bugs in the `readGroundPosition` and the+  `Geodetics.LatLongParser`.++* Added Universal Transverse Mercator (UTM) grid, with both UTM and MGRS grid references

@@ -1,8 +1,9 @@ cabal-version:  3.0 name:           geodetics-version:        1.0.0+version:        1.1.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 extra-doc-files:                 AddingProjections.txt,                 LICENSE,@@ -10,31 +11,35 @@                 changelog.md,                 ToDo.txt license:        BSD-3-Clause-copyright:      Paul Johnson 2018,2024+copyright:      Paul Johnson 2018, 2025 synopsis:       Terrestrial coordinate systems and geodetic calculations. description:    Precise geographical coordinates (latitude & longitude), with conversion between                 different reference frames and projections.                                  Certain distinguished reference frames and grids are given distinct                 types so that coordinates expressed within them cannot be confused with-                from coordinates in other frames.+                coordinates in other frames. license-file:   LICENSE maintainer:     Paul Johnson <paul@cogito.org.uk> homepage:       https://github.com/PaulJohnson/geodetics category:       Geography-tested-with:    GHC==9.10.1  source-repository head   type:     git   location: https://github.com/PaulJohnson/geodetics +flag dev+    description: Turn on development settings.+    default:     False+    manual:      True+ library   hs-source-dirs:  src   build-depends:                    base >= 4.17 && < 5,                    array >= 0.1 && < 0.6,+                   parsec >= 3.1 && < 3.2,                    Stream >= 0.4.6 && < 0.5-  ghc-options:     -Wall   exposed-modules:                    Geodetics.Altitude,                    Geodetics.Ellipsoids,@@ -44,8 +49,13 @@                    Geodetics.Path,                    Geodetics.Stereographic,                    Geodetics.TransverseMercator,-                   Geodetics.UK+                   Geodetics.UK,+                   Geodetics.UTM   Default-Language: Haskell2010+  if flag(dev)+      ghc-options:+          -Wall -Werror -Wredundant-constraints -Wpartial-fields+          -Wunused-packages  test-suite GeodeticTest   type:            exitcode-stdio-1.0@@ -53,16 +63,18 @@   x-uses-tf:       true   build-depends:   geodetics,                    base >= 4.6 && < 5,-                   HUnit >= 1.2,+                   hspec >= 2.11,+                   HUnit >= 1.6,                    QuickCheck >= 2.4,-                   test-framework >= 0.4.1,-                   test-framework-quickcheck2,-                   test-framework-hunit,-                   array,                    checkers   hs-source-dirs:                    test-  ghc-options:     -Wall -rtsopts+  ghc-options:     -rtsopts   other-modules:                    ArbitraryInstances+                   LatLongParser   Default-Language: Haskell2010+  if flag(dev)+      ghc-options:+          -Wall -Werror -Wredundant-constraints -Wpartial-fields+          -Wunused-packages

@@ -207,11 +207,7 @@ -- The WGS84 has a special place in this library as the standard Ellipsoid against -- which all others are defined. data WGS84 = WGS84--instance Eq WGS84 where _ == _ = True--instance Show WGS84 where-   show _ = "WGS84"+  deriving (Eq, Show)  instance Ellipsoid WGS84 where    majorRadius _ = 6378137.0

@@ -10,6 +10,10 @@    groundDistance,    properAngle,    showAngle,+   showGeodeticLatLong,+   showGeodeticSignedDecimal,+   showGeodeticNSEWDecimal,+   showGeodeticDDDMMSS,    -- * Earth Centred Earth Fixed Coordinates    ECEF,    geoToEarth,@@ -18,13 +22,13 @@    WGS84 (..) ) where --import Data.Char (chr)+import Data.Bool (bool) import Data.Maybe import Geodetics.Altitude import Geodetics.Ellipsoids import Geodetics.LatLongParser import Text.ParserCombinators.ReadP+import Text.Printf  -- | Defines a three-D position on or around the Earth using latitude, -- longitude and altitude with respect to a specified ellipsoid, with@@ -33,7 +37,7 @@ -- places, which is a -- resolution of about 1m on the Earth's surface. Internally latitude -- and longitude are stored as double precision radians. Convert to--- degrees using e.g.  @latitude g /~ degree@.+-- degrees using e.g.  @latitude g / degree@. -- -- The functions here deal with altitude by assuming that the local -- height datum is always co-incident with the ellipsoid in use,@@ -62,12 +66,8 @@ }  instance (Ellipsoid e) => Show (Geodetic e) where-   show g = concat [-      showAngle (abs $ latitude g),  " ", letter "SN" (latitude g),  ", ",-      showAngle (abs $ longitude g), " ", letter "WE" (longitude g), ", ",-      show (altitude g), " ", show (ellipsoid g)]-      where letter s n = [s !! (if n < 0 then 0 else 1)]-+   show g =+     showGeodeticLatLong g ++ ", " ++ show (altitude g) ++ " " ++ show (ellipsoid g)  -- | Read the latitude and longitude of a ground position and -- return a Geodetic position on the specified ellipsoid.@@ -80,7 +80,7 @@ -- -- * Decimal degrees NSEW: 34.52327N, 46.23234W ----- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92'+-- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92' W -- -- * Degrees, minutes and seconds (units optional): 34° 31' 23.52\" N, 46° 13' 56.43\" W --@@ -88,28 +88,110 @@ readGroundPosition :: (Ellipsoid e) => e -> String -> Maybe (Geodetic e) readGroundPosition e str =    case map fst $ filter (null . snd) $ readP_to_S latLong str of-      [] -> Nothing-      (lat,long) : _ -> Just $ groundPosition $ Geodetic (lat * degree) (long * degree) undefined e-+      [(lat,long)] -> Just $ groundPosition $ Geodetic+        { latitude = lat * degree,+          longitude = long * degree,+          geoAlt = 0.0,+          ellipsoid = e+        }+      -- It appears incorrect to accept more than 1 interpretation since+      -- that would mean a pathological ambiguity in the parser. As stated,+      -- the accepted grammars do not overlap. Furthermore, should there+      -- exist several different ways to parse the input, it appears+      -- dangerous to arbitrarily choose one of them while disregarding the+      -- others, hence that case should be treated as a failure.+      _ -> Nothing  -- | Show an angle as degrees, minutes and seconds to two decimal places. showAngle :: Double -> String showAngle a-   | isNaN a        = "NaN"  -- Not a Nangle-   | isInfinite a   = sgn ++ "Infinity"-   | otherwise      = concat [sgn, show d, [chr 0xB0, ' '],-                              show m, "\8242 ",-                              show s, ".", dstr, "\8243" ]+   | isNaN a        = "NaN"  -- Not an angle+   | isInfinite a   = sign ++ "Infinity"+   | otherwise      = printf "%s%d° %d′ %.2f″" sign  d m s    where-      sgn = if a < 0 then "-" else ""-      centisecs :: Integer-      centisecs = abs $ round $ (a / (arcsecond / 100))-      (d, m1) = centisecs `divMod` 360000-      (m, s1) = m1 `divMod` 6000   -- hundredths of arcsec per arcmin-      (s, ds) = s1 `divMod` 100-      dstr = reverse $ take 2 $ reverse (show ds) ++ "00" -- Decimal fraction with zero padding.+     sign = if isPositive then "" else "-"+     (d, m, s, isPositive) = radianToDegrees a +-- | Show 'Geodetic' as a pair of degrees, minutes, and seconds. This is+-- similar to the 'Show' instance for 'Geodetic', except it does not include+-- the altitude and the ellipsoid.+--+-- @since 1.1.0+showGeodeticLatLong :: Geodetic e -> String+showGeodeticLatLong x =+  printf+    "%d° %d′ %.2f″ %c, %d° %d′ %.2f″ %c"+    latD latM latS (bool 'S' 'N' isNorth)+    longD longM longS (bool 'W' 'E' isEast)+  where+    (latD, latM, latS, isNorth) = radianToDegrees (latitude x)+    (longD, longM, longS, isEast) = radianToDegrees (longitude x) +-- | Show 'Geodetic' as a pair of signed decimal degrees (5 decimal places+-- of precision), e.g. 34.52327, -46.23234.+--+-- @since 1.1.0+showGeodeticSignedDecimal :: Geodetic e -> String+showGeodeticSignedDecimal x =+  printf+    "%.5f, %.5f"+    (toDegrees (latitude x))+    (toDegrees (longitude x))+  where+    toDegrees r = r * 180 / pi++-- | Show 'Geodetic' as a pair of decimal degrees NSEW, e.g. 34.52327N,+-- 46.23234W.+--+-- @since 1.1.0+showGeodeticNSEWDecimal :: Geodetic e -> String+showGeodeticNSEWDecimal x =+  printf+    "%.5f%c, %.5f%c"+    (toAbsDegrees (latitude x))+    (c 'N' 'S' (latitude x))+    (toAbsDegrees (longitude x))+    (c 'E' 'W' (longitude x))+  where+    toAbsDegrees r = abs (r * 180 / pi)+    c forPositive forNegative r =+      if r < 0+        then forNegative+        else forPositive++-- | Show 'Geodetic' as a pair of angles in the DDDMMSS format, e.g.+-- 343123.52N, 461356.43W.+--+-- @since 1.1.0+showGeodeticDDDMMSS ::+  -- | Use leading zeros+  Bool ->+  Geodetic e ->+  String+showGeodeticDDDMMSS useLeadingZeros x =+  printf+    fmt+    latD latM latS (bool 'S' 'N' isNorth)+    longD longM longS (bool 'W' 'E' isEast)+  where+    fmt =+      if useLeadingZeros+        then "%03d%02d%05.2f%c, %03d%02d%05.2f%c"+        else "%d%02d%05.2f%c, %d%02d%05.2f%c"+    (latD, latM, latS, isNorth) = radianToDegrees (latitude x)+    (longD, longM, longS, isEast) = radianToDegrees (longitude x)++-- | An internal helper function for 'showAngle' and 'showGeodeticDDDMMSS'.+radianToDegrees :: Double -> (Integer, Integer, Double, Bool)+radianToDegrees a = (d, m, s, a >= 0)+  where+    centisecs :: Integer+    centisecs = abs $ round (a / (arcsecond / 100))+    (d, m1) = centisecs `divMod` 360000+    (m, _) = m1 `divMod` 6000   -- hundredths of arcsec per arcmin+    s = abs $+      (abs a - fromIntegral d * degree - fromIntegral m * arcminute) / arcsecond+ instance (Ellipsoid e) => HasAltitude (Geodetic e) where    altitude = geoAlt    setAltitude h g = g{geoAlt = h}@@ -118,7 +200,7 @@  -- | The point on the Earth diametrically opposite the argument, with -- the same altitude.-antipode :: (Ellipsoid e) => Geodetic e -> Geodetic e+antipode :: Geodetic e -> Geodetic e antipode g = Geodetic lat long (geoAlt g) (ellipsoid g)    where       lat = negate $ latitude g@@ -230,16 +312,16 @@        listToMaybe $ dropWhile converging $ take 100 $ zip lambdas $ drop 1 lambdas      let        uSq = cos2Alpha * (a^ _2 - b^ _2) / b^ _2-       bigA = 1 + uSq/16384 * (4096 + uSq * ((-768) + uSq * ((320 - 175*uSq))))-       bigB =     uSq/1024  * (256  + uSq * ((-128) + uSq * ((74 -  47* uSq))))+       bigA = 1 + uSq/16384 * (4096 + uSq * ((-768) + uSq * (320 - 175*uSq)))+       bigB =     uSq/1024  * (256  + uSq * ((-128) + uSq * (74 -  47* uSq)))        deltaDelta =          bigB * sinDelta * (cos2DeltaM +                              bigB/4 * (cosDelta * (2 * cos2DeltaM^ _2 - 1)                                        - bigB/6 * cos2DeltaM * (4 * sinDelta^ _2 - 3)                                           * (4 * cos2DeltaM - 3)))        s = b * bigA * (delta - deltaDelta)-       alpha1 = atan2(cosU2 * sin lambda) (cosU1 * sinU2 - sinU1 * cosU2 * cos lambda)-       alpha2 = atan2(cosU1 * sin lambda) (cosU1 * sinU2 * cos lambda - sinU1 * cosU2)+       alpha1 = atan2 (cosU2 * sin lambda) (cosU1 * sinU2 - sinU1 * cosU2 * cos lambda)+       alpha2 = atan2 (cosU1 * sin lambda) (cosU1 * sinU2 * cos lambda - sinU1 * cosU2)      return (s, alpha1, alpha2)   where     f = flattening $ ellipsoid p1@@ -257,7 +339,7 @@       where         sinLambda = sin lambda         cosLambda = cos lambda-        sinDelta = sqrt((cosU2 * sinLambda) ^ _2 ++        sinDelta = sqrt ((cosU2 * sinLambda) ^ _2 +                         (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ _2)         cosDelta = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda         delta = atan2 sinDelta cosDelta

@@ -42,14 +42,7 @@ data GridPoint r = GridPoint {    eastings, northings, altGP :: Double,    gridBasis :: r-} deriving (Show)---instance Eq (GridPoint r) where-   p1 == p2  =-      eastings p1 == eastings p2 &&-      northings p1 == northings p2 &&-      altGP p1 == altGP p2+} deriving (Eq, Show)  instance HasAltitude (GridPoint g) where    altitude = altGP@@ -156,11 +149,12 @@       d = sum $ zipWith (*)          (map (fromIntegral . digitToInt) ds)          (drop 1 $ iterate (/ 10) sq)-      p = sq / fromIntegral ((10 :: Int) ^ n)+      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 nearest position to the specified--- number of digits, expressed as an integer count of squares and a string of digits.+-- 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.@@ -177,7 +171,7 @@       p = 10 ^ n       unit :: Double       unit = sq / fromIntegral p-      u = round (d / unit)+      u = floor (d / unit)       (sqs, d1) = u `divMod` p       s = show d1       pad = if n == 0 then "" else replicate (n - length s) '0' ++ s

@@ -43,6 +43,10 @@ decimal :: ReadP Double decimal = do    str1 <- munch1 isDigit+   -- In order to avoid ambiguity where 'decimal' and 'dms7' both match+   -- (which can happen in cases with several leading zeros for very small+   -- angles), stipulate that str1 is sufficiently short.+   guard (length str1 < 5)    option (read str1) $ do       str2 <- char '.' *> munch1 isDigit       return $ read $ str1 ++ '.' : str2@@ -69,33 +73,26 @@    d <- fromIntegral <$> natural    guard $ d <= 360    skipSpaces-   ms <- option 0 $ do-      m <- fromIntegral <$> natural-      guard $ m < 60-      skipSpaces-      s <- option 0 decimal-      guard $ s < 60-      return $ m / 60 + s / 3600-   return $ d + ms-+   m <- fromIntegral <$> natural+   guard $ m < 60+   skipSpaces+   s <- decimal+   guard $ s < 60+   return $ d + m / 60 + s / 3600  -- | Parse an unsigned angle written using degrees, minutes and seconds with units (° ' \"). -- At least one component must be specified. degreesMinutesSecondsUnits :: ReadP Double degreesMinutesSecondsUnits = do-   (s, a) <- gather $ do-      d <- fromIntegral <$> option 0 (natural <* char '°')-      guard $ d <= 360-      skipSpaces-      m <- fromIntegral <$> option 0 (natural <* minuteTick)-      guard $ m < 60-      skipSpaces-      s <- option 0 (decimal <* secondTick)-      guard $ s < 60-      return $ d + m / 60 + s / 3600-   guard $ not $ null s  -- Must specify at least one component.-   return a-+   d <- fromIntegral <$> natural <* char '°'+   guard $ d <= 360+   skipSpaces+   m <- fromIntegral <$> natural <* minuteTick+   guard $ m < 60+   skipSpaces+   s <- decimal <* secondTick+   guard $ s < 60+   return $ d + m / 60 + s / 3600  -- | Parse an unsigned angle written using degrees and decimal minutes. degreesDecimalMinutes :: ReadP Double@@ -103,7 +100,7 @@    d <- fromIntegral <$> natural    skipSpaces    guard $ d <= 360   -- Difference from degreesMinutesSeconds just to shut style checker up.-   m <- option 0 decimal+   m <- decimal    guard $ m < 60    return $ d + m/60 @@ -111,15 +108,12 @@ -- | Parse an unsigned angle written using degrees and decimal minutes with units (° ') degreesDecimalMinutesUnits :: ReadP Double degreesDecimalMinutesUnits = do-   (s, a) <- gather $ do-      d <- fromIntegral <$> option 0 (natural <* char '°')-      guard $ d <= 360-      m <- option 0 (decimal <* minuteTick)-      guard $ m < 60-      return $ d + m / 60-   guard $ not $ null s  -- Must specify at least one component.-   return a-+   d <- fromIntegral <$> natural <* char '°'+   guard $ d <= 360+   skipSpaces+   m <- decimal <* minuteTick+   guard $ m < 60+   return $ d + m / 60  -- | Parse an unsigned angle written in DDDMMSS.ss format. -- Leading zeros on the degrees and decimal places on the seconds are optional

@@ -209,7 +209,7 @@ -- polar singularities are not included. Anyone using a rhumb line that close to a pole -- must be going round the twist anyway. ----- Based on *Practical Sailing Formulas for Rhumb-Line Tracks on an Oblate Earth* +-- Based on /Practical Sailing Formulas for Rhumb-Line Tracks on an Oblate Earth/ -- by G.H. Kaplan, U.S. Naval Observatory. Except for points close to the poles  -- the approximation is accurate to within a few meters over 1000km. rhumbPath :: (Ellipsoid e) =>

@@ -33,8 +33,8 @@    majorRadius _ = 6377563.396    flatR _ = 299.3249646    helmert _ = Helmert {-      cX = 446.448, cY = (-125.157), cZ = 542.06,-      helmertScale = (-20.4894),+      cX = 446.448, cY = -125.157, cZ = 542.06,+      helmertScale = -20.4894,       rX = 0.1502 * arcsecond, rY = 0.247 * arcsecond, rZ = 0.8421 * arcsecond }  -- | The UK National Grid is a Transverse Mercator projection with a true origin at@@ -100,10 +100,10 @@ fromUkGridLetters c1 c2 = applyOffset <$> (mappend <$> g1 <*> g2) <*> letterOrigin    where       letterOrigin = Just $ GridPoint ((-1000) * kilometer) ((-500) * kilometer) m0 UkNationalGrid-      gridIndex c =-         if inRange ('A', 'H') c then Just $ ord c - ord 'A'  -- 'I' is not used.-         else if inRange ('J', 'Z') c then Just $ ord c - ord 'B'-         else Nothing+      gridIndex c+         | inRange ('A', 'H') c  = Just $ ord c - ord 'A'  -- 'I' is not used.+         | inRange ('J', 'Z') c  = Just $ ord c - ord 'B'+         | otherwise             = Nothing       gridSquare c = do -- Maybe monad          g <- gridIndex c          let (y,x) = g `divMod` 5

@@ -0,0 +1,464 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE NumericUnderscores    #-}+{-# LANGUAGE FlexibleContexts #-}++{- | Universal Transverse Mercator (UTM)++The UTM grid system covers the whole world between 84°N and 80°S. It divides the world into +grid zones of 6° longitude by 8° latitude. Each zone has a 2 digit number for longitude and+a letter for latitude. This regular system has two exceptions:++* North of Norway the zones 32X, 34X and 36X are not used, with 31X, 33X, 35X and 37X being+  wider instead.++* Zone 32V is widened to cover the south-western end of Norway.++There are two notations for writing UTM grid positions:++* The UTM standard: Zone number, N or S for hemisphere, and then northings and eastings+  relative to the equator.++* The Military Grid Reference System (MGRS): Zone number, latitude band letter, a+  2 letter code for the 100km square within the zone, and then northings and eastings within+  that square.++In this library each UTM longitude zone has two grids, one for the northern hemisphere and+one for the south.++For more details see++* https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system.++* THE UNIVERSAL GRIDS: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS).+  DMA Technical Manual. AD-A226497. https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf+-}+module Geodetics.UTM (+  UtmHemisphere (..),+  UtmZoneNumber,+  utmZoneNumber,+  UtmZone (utmHemisphere, utmZoneNum, utmProjection),+  utmZone,+  mkUtmZone,+  mkUtmZoneUnsafe,+  fromUtmGridReference,+  toUtmGridReference,+  mgrsBandLetterToLatitude,+  mgrsLatitudeToBandLetter,+  fromMgrsGridReference,+  toMgrsGridReference+) where++import Control.Monad (mplus, guard, void, when, unless)+import Data.Array+import Data.Char+import Data.List+import Geodetics.Ellipsoids+import Geodetics.Geodetic+import Geodetics.Grid+import Geodetics.TransverseMercator+import Text.Parsec+import Text.Parsec.Error+import Text.Printf+import Text.Read++++-- | In UTM the northern and southern hemispheres have different false origins.+data UtmHemisphere = UtmNorth | UtmSouth deriving Eq++instance Show UtmHemisphere where+  show UtmNorth = "N"+  show UtmSouth = "S"+++-- | A UTM Zone number. Must be between 1 and 60.+type UtmZoneNumber = Int+++-- | A UTM Zone, representing a band of typically 6 degrees of latitude between the equator and one of+-- the poles. The projection *must* match the hemisphere and zone.+data UtmZone = UtmZone {+  utmHemisphere :: UtmHemisphere,+  utmZoneNum :: UtmZoneNumber,+  utmProjection :: GridTM WGS84+} deriving (Show)++instance Eq UtmZone where+  z1 == z2  = utmHemisphere z1 == utmHemisphere z2 && utmZoneNum z1 == utmZoneNum z2++instance GridClass UtmZone WGS84 where+  fromGrid p = fromGrid $ unsafeGridCoerce (utmProjection $ gridBasis p) p+  toGrid grid = unsafeGridCoerce grid . toGrid (utmProjection grid)+  gridEllipsoid _ = WGS84+++-- Internal data type representing a "rectangle" of latitude/longitude with an exceptional zone number.+data UtmException = UtmE {+  uteSW :: (Int, Int),  -- South west corner in integer degrees (lat, long), inclusive.+  uteNE :: (Int, Int),  -- North east corner in integer degrees (lat, long), exclusive.+  uteActual :: UtmZoneNumber+} deriving Show+++-- | Determine if the integer latitude and longitude are within the exception area.+inException :: Int -> Int -> UtmException -> Bool+inException lat long e =+    inR lat  (fst $ uteSW e) (fst $ uteNE e) &&+    inR long (snd $ uteSW e) (snd $ uteNE e)+  where+    inR v v1 v2 = v1 <= v && v < v2+++-- The UTM zone that encloses a given geodetic position. For most of the world this is based on+-- @longitude/6@, but there are exceptions around Norway and Svalbard.+utmZoneNumber :: Geodetic a -> Maybe UtmZoneNumber+utmZoneNumber geo = do+    guard $ lat1 >= (-80) && lat1 < 84+    return $ maybe zone1 uteActual exception+  where+    lat1 = floor $ latitude geo / degree+    long1 = floor $ longitude geo / degree+    zone1 = (long1 `div` 6 + 30) `mod` 60 + 1+    exception = find (inException lat1 long1)+      [+        UtmE (56,03) (64,12) 32,  -- Southwestern end of Norway around Bergen.+        UtmE (72,00) (84,09) 31,+        UtmE (72,09) (84,21) 33,  -- Svalbard.+        UtmE (72,21) (84,33) 35,+        UtmE (72,33) (84,42) 37+      ]+++-- | The UTM Zone for the given location, if it exists.+utmZone :: Geodetic a -> Maybe UtmZone+utmZone geo = do+  let hemi = if latitude geo >= 0 then UtmNorth else UtmSouth+  zn <- utmZoneNumber geo+  mkUtmZone hemi zn+++-- | Construct a UTM Zone value. Returns @Nothing@ if the zone number is out of range.+mkUtmZone :: UtmHemisphere -> UtmZoneNumber -> Maybe UtmZone+mkUtmZone h n = do+    guard $ n >= 1 && n <= 60+    return $ mkUtmZoneUnsafe h n+++-- | Construct a UTM Zone value without checking whether the zone number is valid.+mkUtmZoneUnsafe :: UtmHemisphere -> UtmZoneNumber -> UtmZone+mkUtmZoneUnsafe h n = UtmZone h n $ mkGridTM trueO falseO scale+  where+    trueO = Geodetic 0 (degree * fromIntegral (n * 6 - 183)) 0 WGS84+    falseO = case h of+      UtmNorth -> GridOffset (-500_000) 0 0+      UtmSouth -> GridOffset (-500_000) (-10_000_000) 0+    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,+-- including an optional decimal point. Negative values are not permitted.+--+-- Northings and eastings can each be followed by an optional unit. The unit+-- must be either \"m\" or \"km\". The units for both+-- must be the same because otherwise its probably an error. The default is meters.+--+-- Northings may be followed by an \"N\" and Eastings may be followed by an \"E\".+--+-- 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+    Left err -> Left $ lines $ showErrorMessages+      "or" "unknown parse error" "expecting" "unexpected" "end of input"+      (errorMessages err)+    Right r -> Right r+  where+    gridP = do+      spaces1+      zone <- readZone <?> "Zone number"+      hemi <- readHemi <?> "Hemisphere (N or S)"+      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 $ mkUtmZoneUnsafe hemi zone+    readZone :: Parsec String () UtmZoneNumber+    readZone = do+      ds <- many1 digit+      case readMaybe ds of+        Nothing -> fail "Zone number not found."+        Just n ->+          if n < 1 || n > 60+            then fail $ "Zone number " <> show n <> " out of range."+            else return n+    readHemi :: Parsec String () UtmHemisphere+    readHemi = do+      h <- oneOf "NSns"+      case toUpper h of+        'N' -> return UtmNorth+        'S' -> return UtmSouth+        _ -> fail $ "Invalid hemisphere: " <> (h : ". Must be N or S.")+    readDistance :: Parsec String () (Double, UtmGridUnit)  -- (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)+      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.+++-- | Convert a grid point to a UTM grid reference.+-- 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.+  -> Bool -- ^ Include \"E\" and \"N\" in the output.+  -> Int  -- ^ Digits of resolution. 0 = 1m resolution, 1 = 10m, 2 = 100m etc. (-2) = 1cm.+  -> GridPoint UtmZone+  -> String+toUtmGridReference unit letters res gp =+    zoneStr <> " " <>+    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)+    b = gridBasis gp+    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+

@@ -1,5 +1,9 @@-{-# LANGUAGE FlexibleInstances, RankNTypes, KindSignatures, DataKinds #-}+{-# LANGUAGE FlexibleInstances, RankNTypes, KindSignatures, DataKinds, NumericUnderscores #-} {-# OPTIONS_GHC -fno-warn-orphans #-}+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}+{-# HLINT ignore "Redundant <$>" #-}+{-# HLINT ignore "Use underscore" #-}+{-# HLINT ignore "Redundant bracket" #-}  -- | Orphan "Arbitrary" and related instances for testing purposes.  @@ -14,6 +18,7 @@ import Geodetics.Stereographic as SG import Geodetics.TransverseMercator as TM import Test.QuickCheck+import Text.Printf   -- | Shrink an angle so that shrunk values are round numbers of degrees.@@ -31,40 +36,40 @@    show (Bearing b) = "Bearing " ++ showAngle b  instance Arbitrary Bearing where-   arbitrary = Bearing <$> (* degree) <$> choose (-180,180)+   arbitrary = Bearing . (* degree) <$> choose (-180,180)    shrink (Bearing b) = Bearing <$> shrinkAngle b-   -   ++ newtype Azimuth = Azimuth Double  instance Show Azimuth where    show (Azimuth a) = "Azimuth " ++ showAngle a-   + instance Arbitrary Azimuth where-   arbitrary = Azimuth <$> (* degree) <$> choose (0,90)+   arbitrary = Azimuth . (* degree) <$> choose (0,90)    shrink (Azimuth a) = Azimuth <$> shrinkAngle a-   -   ++ -- | Wrapper for arbitrary distances up to 10,000 km newtype Distance = Distance Double deriving (Show)  instance Arbitrary Distance where-   arbitrary = Distance <$> (* kilometer) <$> choose (0,10000)+   arbitrary = Distance . (* kilometer) <$> choose (0,10_000)    shrink (Distance d) = Distance <$> shrinkDistance d-    + -- | Wrapper for arbitrary distances up to 1,000 km newtype Distance2 = Distance2 Double deriving (Show)  instance Arbitrary Distance2 where-   arbitrary = Distance2 <$> (* kilometer) <$> choose (0,1000)+   arbitrary = Distance2 . (* kilometer) <$> choose (0,1_000)    shrink (Distance2 d) = Distance2 <$> shrinkDistance d  -- | Wrapper for arbitrary altitudes up to 10 km newtype Altitude = Altitude Double deriving (Show)  instance Arbitrary Altitude where-   arbitrary = Altitude <$> (* kilometer) <$> choose (0,10)+   arbitrary = Altitude . (* kilometer) <$> choose (0,10)    shrink (Altitude h) = Altitude <$> shrinkDistance h  @@ -76,20 +81,34 @@    shrink (Scalar s) = Scalar <$> shrink s  --- | Wrapper for arbitrary grid references.-newtype GridRef = GridRef String deriving Show+-- | Wrapper for arbitrary UK grid references.+newtype UkGridRef = UkGridRef String deriving Show -instance Arbitrary GridRef where+instance Arbitrary UkGridRef where    arbitrary = do       n <- choose (0,4)       c1 <- elements "HJNOST" -- General vicinity of UK       c2 <- elements $ ['A'..'H'] ++ ['J'..'Z']       dx <- vectorOf n $ choose ('0','9')       dy <- vectorOf n $ choose ('0','9')-      return $ GridRef $ c1 : c2 : (dx ++ dy)+      return $ UkGridRef $ c1 : c2 : (dx ++ dy)    shrink = shrinkNothing  +-- | Wrapper for arbitrary UTM grid references.+newtype UtmGridRef = UtmGridRef String deriving Show++instance Arbitrary UtmGridRef where+   arbitrary = do+      zone <- choose (1,60 :: Int)+      hemi <- elements "NS"+      e <- choose (100_000, 899_999 :: Integer)+      n <- case hemi of+         'N' -> choose (0, 9_328_195 :: Integer) -- Equator to 84 degrees north.+         _   -> choose (1_118_248, 10_000_000)  -- 80 degrees south to Equator.+      return $ UtmGridRef $ printf "%02d" zone <> [hemi] <> " " <> show e <> "E " <> show n <> "N"++ -- | Generate in range +/- <arg> m. genOffset :: Double -> Gen Double genOffset d = choose (-d, d)@@ -106,8 +125,8 @@  genSeconds :: Gen Double genSeconds = (* arcsecond) <$> choose (-10,10)-     + -- | Shrinking with the original value preserved. Used for shrinking records.  See  -- http://stackoverflow.com/questions/14006005/idiomatic-way-to-shrink-a-record-in-quickcheck for details. shrink' :: (Arbitrary a) => a -> [a]@@ -118,70 +137,70 @@   instance Arbitrary Helmert where-   arbitrary = -      Helmert <$> genOffset 300 <*> genOffset 300 <*> genOffset 300 <*> +   arbitrary =+      Helmert <$> genOffset 300 <*> genOffset 300 <*> genOffset 300 <*>          (choose (-5,10)) <*> genSeconds <*> genSeconds <*> genSeconds-   shrink h = +   shrink h =       drop 1 $ Helmert <$> shrink' (cX h) <*> shrink' (cY h) <*> shrink' (cZ h) <*>          shrink' (helmertScale h) <*>-         shrink' (rX h) <*> shrink' (rY h) <*> shrink' (rZ h)      +         shrink' (rX h) <*> shrink' (rY h) <*> shrink' (rZ h)   instance Arbitrary WGS84 where    arbitrary = return WGS84    shrink = shrinkNothing-    + instance Arbitrary LocalEllipsoid where    arbitrary =-      LocalEllipsoid <$> (("Local_" ++) <$> replicateM 3 (choose ('A','Z'))) <*>  -- name+      (LocalEllipsoid . ("Local_" ++) <$> replicateM 3 (choose ('A','Z'))) <*>  -- name          (choose (6378100, 6378400)) <*>                    -- majorRadius          (choose (297,300)) <*>                             -- flatR          arbitrary                                          -- helmert    shrink e = drop 1 $ LocalEllipsoid (nameLocal e) (majorRadius e) (flatR e) <$> shrink' (helmert e) -        + instance (Ellipsoid e, Arbitrary e) => Arbitrary (Geodetic e) where-   arbitrary = +   arbitrary =       Geodetic <$> genLatitude <*> genLongitude <*> genOffset 1 <*> arbitrary-   shrink g = -      drop 1 $ Geodetic <$> shrinkAngle' (latitude g) <*> shrinkAngle' (longitude g) <*> +   shrink g =+      drop 1 $ Geodetic <$> shrinkAngle' (latitude g) <*> shrinkAngle' (longitude g) <*>          shrink' (altitude g) <*> shrink' (ellipsoid g)  instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridTM e)) where    arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary-   shrink p = drop 1 $ GridPoint <$> -      shrink' (eastings p) <*> -      shrink' (northings p) <*> -      shrink' (altitude p) <*> +   shrink p = drop 1 $ GridPoint <$>+      shrink' (eastings p) <*>+      shrink' (northings p) <*>+      shrink' (altitude p) <*>       shrink' (gridBasis p)   instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridStereo e)) where    arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary-   shrink p = drop 1 $ GridPoint <$> -      shrink' (eastings p) <*> -      shrink' (northings p) <*> -      shrink' (altitude p) <*> +   shrink p = drop 1 $ GridPoint <$>+      shrink' (eastings p) <*>+      shrink' (northings p) <*>+      shrink' (altitude p) <*>       shrink' (gridBasis p)   instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridTM e) where    arbitrary = mkGridTM <$> arbitrary <*> arbitrary <*> choose (0.95,1.0)    shrink tm = drop 1 $ mkGridTM <$> shrink' (trueOrigin tm) <*> shrink' (falseOrigin tm) <*> [TM.gridScale tm]-   -   ++ instance Arbitrary GridOffset where    arbitrary = GridOffset <$> genOffset 100000 <*> genOffset 100000 <*> genAlt-   shrink d = drop 1 $ GridOffset <$> +   shrink d = drop 1 $ GridOffset <$>       shrink' (deltaEast d) <*> shrink' (deltaNorth d) <*> shrink' (deltaAltitude d)   instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridStereo e) where    arbitrary = mkGridStereo <$> arbitrary <*> arbitrary <*> choose (0.95,1.0)    shrink sg = drop 1 $ mkGridStereo <$> shrink' (gridTangent sg) <*> shrink' (gridOrigin sg) <*> [SG.gridScale sg]-    + -- | Wrapper for arbitrary rays, along with creation parameters for printing and shrinking. data Ray e = Ray (Geodetic e) Double Double @@ -202,8 +221,8 @@       b' <- shrinkAngle' b       e' <- shrinkAngle' e       return $ Ray p0' b' e'-      -     ++ -- | Two rhumb paths starting not more than 1000 km apart. data RhumbPaths2 = RP2 {       rp2Point0 :: Geodetic WGS84,@@ -212,24 +231,24 @@       rp2Bearing1 :: Bearing,       rp2Bearing2 :: Bearing    }-   + instance Show RhumbPaths2 where    show rp2 = show (pt1, Bearing b1) ++ show (pt2, Bearing b2)-      where +      where          (p1, p2) = mk2RhumbPaths rp2          (pt1, b1, _) = pathFunc p1 0          (pt2, b2, _) = pathFunc p2 0-          + instance Arbitrary RhumbPaths2 where-   arbitrary = RP2 +   arbitrary = RP2       <$> arbitrary `suchThat` ((< (70 * degree)) . abs . latitude)       <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary-   shrink rp = -      drop 1 $ RP2 <$> -         shrink' (rp2Point0 rp) <*> -         shrink' (rp2Bearing0 rp) <*> -         shrink' (rp2Distance rp) <*> -         shrink' (rp2Bearing1 rp) <*> +   shrink rp =+      drop 1 $ RP2 <$>+         shrink' (rp2Point0 rp) <*>+         shrink' (rp2Bearing0 rp) <*>+         shrink' (rp2Distance rp) <*>+         shrink' (rp2Bearing1 rp) <*>          shrink' (rp2Bearing2 rp)  mk2RhumbPaths :: RhumbPaths2 -> (Path WGS84, Path WGS84)@@ -240,4 +259,4 @@       path1 = rhumbPath pt0 b1       (pt2, _, _) = pathFunc path0 d       path2 = rhumbPath pt2 b2- +

@@ -0,0 +1,209 @@+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE ScopedTypeVariables #-}++module LatLongParser (parserTests) where++import ArbitraryInstances ()+import Geodetics.Geodetic+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.QuickCheck hiding ((===))++-- | The 'Coordinates' that certain test cases must return.+expected0, expected1, expected2, expected3, expected4, expected5 :: Coordinates+expected0 = Coordinates 0.6025447300549828 (-0.8069065539014753)+expected1 = Coordinates 0.6025545620764358 (-0.8069006197820184)+expected2 = Coordinates 0.6024294801467094 (-0.806633002630046)+expected3 = Coordinates 0.6025435083245063 (-0.8069065829902962)+expected4 = Coordinates 0.6025409872933646 (-0.8069044982914674)+expected5 = Coordinates 0.3139836658436815 (-5.259743626357356e-4)++parserTests :: SpecWith ()+parserTests =+  describe "LatLongParser" $ do+    testTree+      "Signed decimal degrees"+      ("34.52327", "-46.23234")+      (commaOrNot (produces expected0))+    testTree+      "Decimal degrees NSEW"+      ("34.52327N", "46.23234W")+      (abOrBa (commaOrNot (produces expected0)))+    testTree+      "Degrees and decimal minutes, with decimal point"+      ("34° 31.43' N", "46° 13.92' W")+      (allVariations expected1)+    testTree+      "Degrees and decimal minutes, without decimal point"+      ("34° 31' N", "46° 13' W")+      (allVariations expected2)+    testTree+      "Degrees, minutes and seconds, with decimal point"+      ("34° 31' 23.52\" N", "46° 13' 56.43\" W")+      (allVariations expected3)+    testTree+      "Degrees, minutes and seconds, without decimal point"+      ("34° 31' 23\" N", "46° 13' 56\" W")+      (allVariations expected4)+    testTree+      "DDDMMSS format, no leading zeros, with decimal point"+      ("343123.52N", "461356.43W")+      (abOrBa (commaOrNot (produces expected3)))+    testTree+      "DDDMMSS format, with leading zeros, with decimal point"+      ("343123.52N", "0461356.43W")+      (abOrBa (commaOrNot (produces expected3)))+    testTree+      "DDDMMSS format, no decimal point"+      ("343123N", "0461356W")+      (abOrBa (commaOrNot (produces expected4)))+    testTree+      "Very small angles in the DDDMMSS format are unambiguous"+      ("175923.78N", "00148.49W")+      (abOrBa (commaOrNot (produces expected5)))+    prop "degrees, minutes, seconds roundtrip" dmsRoundtrip+    prop "signed decimal roundtrip" signedDecimalRoundtrip+    prop "decimal degrees NSEW roundtrip" decimalDegreesNSEWRoundtrip+    prop "DDDMMSS roundtrip" dddmmssRoundtrip+++----------------------------------------------------------------------------+-- The Coordinates type++-- | A pair of coordinates: latitude and longitude expressed in radians.+data Coordinates = Coordinates Double Double+  deriving (Show, Eq)++instance Arbitrary Coordinates where+  arbitrary = do+    x :: Geodetic WGS84 <- arbitrary+    return (Coordinates (latitude x) (longitude x))++-- | Convert 'Coordinates' to a 'Geodetic'.+coordinatesToGeodetic :: Coordinates -> Geodetic WGS84+coordinatesToGeodetic (Coordinates lat long) =+  Geodetic+    { latitude = lat,+      longitude = long,+      geoAlt = 0.0,+      ellipsoid = WGS84+    }++-- | Convert a 'Geodetic' to 'Coordinates'.+geodeticToCoordinates :: Geodetic WGS84 -> Coordinates+geodeticToCoordinates x = Coordinates (latitude x) (longitude x)++----------------------------------------------------------------------------+-- The machinery for unit tests++-- | A tree-like structure of test cases.+testTree ::+  -- | Name of the test tree+  String ->+  -- | The input+  input ->+  -- | The test generator+  (input -> SpecWith ()) ->+  -- | The resulting test+  SpecWith ()+testTree testTreeName input generator =+  describe testTreeName (generator input)+++-- | Stipulate that parsing the input produces the expected result.+produces ::+  -- | The value that should be produced by the parser+  Coordinates ->+  -- | The input+  String ->+  -- | The resulting collection of tests+  SpecWith (Arg Expectation)+produces expected input =+  it caseName $ Just expected `shouldBe` mcoordinates+  where+    caseName = "Input: \"" ++ input ++ "\" -> " ++ show expected+    mcoordinates = geodeticToCoordinates <$> readGroundPosition WGS84 input+++-- | The two 'String's will be examined in various orders.+abOrBa :: ((String, String) -> SpecWith (Arg Expectation)) -> (String, String) -> SpecWith (Arg Expectation)+abOrBa f (a, b) = do+  describe "straight order" (f (a, b))+  describe "reverse order" (f (b, a))+++-- | Two two 'String's will be concatenated with a space between them and+-- with a comma and a space between them.+commaOrNot :: (String -> SpecWith (Arg Expectation)) -> (String, String) -> SpecWith (Arg Expectation)+commaOrNot f (a, b) = do+  describe "with comma" (f (a ++ ", " ++ b))+  describe "without comma" (f (a ++ " " ++ b))+++-- | Test with and without units both in normal and Unicode version.+variousUnits :: (String -> SpecWith ()) -> String -> SpecWith ()+variousUnits f input = do+  describe "without Unicode" (f input)+  describe "with Unicode" (f unicode)+  describe "without units" (f noUnits)+  where+    unicode = g <$> input+    g '\'' = '′'+    g '\"' = '″'+    g x = x+    noUnits = filter (not . isUnit) input+    isUnit = \case+      '°' -> True+      '\'' -> True+      '\"' -> True+      '′' -> True+      '″' -> True+      _ -> False+++-- | Tests with all possible variations.+allVariations :: Coordinates -> (String, String) -> SpecWith (Arg Expectation)+allVariations expected =+  abOrBa (commaOrNot (variousUnits (produces expected)))+++----------------------------------------------------------------------------+-- The machinery for roundtrip property tests++-- | A lax version of 'Test.QuickCheck.===' that is forgiving with respect+-- to precision lost due to printing.+infix 4 ===++(===) :: Coordinates -> Coordinates -> Property+Coordinates lata longa === Coordinates latb longb =+  consider "latitude" lata latb+    .&&. consider "longitude" longa longb+  where+    consider str x y =+      counterexample+        (str ++ " " ++ show x ++ " vs " ++ show y)+        (abs (x - y) < epsilon)+    epsilon = 0.0000001++roundtrip :: Coordinates -> String -> Property+roundtrip x rendered =+  counterexample rendered $+    case geodeticToCoordinates <$> readGroundPosition WGS84 rendered of+      Nothing -> property False+      Just y -> y === x++dmsRoundtrip :: Coordinates -> Property+dmsRoundtrip x =+  roundtrip x (showGeodeticLatLong (coordinatesToGeodetic x))++signedDecimalRoundtrip :: Coordinates -> Property+signedDecimalRoundtrip x =+  roundtrip x (showGeodeticSignedDecimal (coordinatesToGeodetic x))++decimalDegreesNSEWRoundtrip :: Coordinates -> Property+decimalDegreesNSEWRoundtrip x =+  roundtrip x (showGeodeticNSEWDecimal (coordinatesToGeodetic x))++dddmmssRoundtrip :: Bool -> Coordinates -> Property+dddmmssRoundtrip useLeadingZeros x =+  roundtrip x (showGeodeticDDDMMSS useLeadingZeros (coordinatesToGeodetic x))

@@ -1,14 +1,16 @@ {-# OPTIONS_GHC -fno-warn-orphans #-}+{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}+{-# HLINT ignore "Redundant bracket" #-}  module Main where +import Control.Monad+import Data.Char+import Data.Either import Data.Maybe-import Test.Framework (Test, defaultMainWithOpts, testGroup)-import Test.Framework.Options (TestOptions, TestOptions'(..))-import Test.Framework.Runners.Options (RunnerOptions, RunnerOptions'(..))-import Test.Framework.Providers.HUnit-import Test.Framework.Providers.QuickCheck2 (testProperty)-import qualified Test.HUnit as HU+import Test.Hspec+import Test.Hspec.QuickCheck+import Test.HUnit (assertFailure) import Test.QuickCheck import Test.QuickCheck.Checkers (EqProp, eq, (=-=), unbatch) import Test.QuickCheck.Classes (monoid)@@ -22,21 +24,54 @@ import Geodetics.Stereographic import Geodetics.TransverseMercator import Geodetics.UK-+import Geodetics.UTM+import LatLongParser (parserTests)  main :: IO ()-main = do-   let empty_test_opts = mempty :: TestOptions-   let my_test_opts = empty_test_opts {-     topt_maximum_generated_tests = Just 1000-   }--   let empty_runner_opts = mempty :: RunnerOptions-   let my_runner_opts = empty_runner_opts {-     ropt_test_options = Just my_test_opts-   }--   defaultMainWithOpts tests my_runner_opts+main = hspec $ do+   describe "Geodetic" $ do+      prop "WGS84 and back" prop_WGS84_and_back+      prop "Zero ground distance" prop_zero_ground+      describe "UK Points" $ mapM_ pointTest ukPoints+      describe "World line" $ mapM_ worldLineTests worldLines+   parserTests+   describe "Grid" $ do+      prop "Grid Offset 1" prop_offset1+      prop "Grid Offset 2" prop_offset2+      prop "Grid Offset 3" prop_offset3+      prop "Grid 1" prop_grid1+   describe "TransverseMercator" $ prop "fromGrid . toGrid == id" prop_tmGridInverse+   describe "UK" $ do+      prop "UK Grid 1" prop_ukGrid1+      describe "UK Grid 2" $ mapM_ ukGridTest2 ukSampleGrid+      describe "UK Grid 3" $ mapM_ ukGridTest3 ukSampleGrid+      describe "UK Grid 4" $ mapM_ ukGridTest4 ukSampleGrid+      describe "UK Grid 5" $ mapM_ ukGridTest5 ukSampleGrid+   describe "UTM" $ do+      prop "UTM Grid 1" prop_utmGridTest1+      describe "UTM Grid 2" $ mapM_ utmGridTest2 utmSampleGrid+      describe "UTM Grid 3" $ mapM_ utmGridTest3 utmSampleGrid+      describe "UTM Grid 4" $ mapM_ utmGridTest4 utmSampleGrid+      describe "UTM Grid 5" $ mapM_ utmGridTest5 utmSampleGrid+   describe "MGRS" $ do+      prop "MGRS Grid 1" prop_mgrs_gridTest1+      prop "MGRS Grid 2" prop_mgrs_gridTest2+      describe "MGRS Grid 3" $ mapM_ mgrsGridTest3 utmSampleGrid+      describe "MGRS Grid 4" $ mapM_ mgrsGridTest4 utmSampleGrid+   describe "Stereographic" $ do+      it "toGrid north" stereographicToGridN+      it "fromGrid north" stereographicFromGridN+      it "toGrid south" stereographicToGridS+      it "fromGrid south" stereographicFromGridS+      prop "Stereographic round trip" prop_stereographic+   describe "Paths" $ do+      prop "Ray Path 1" prop_rayPath1+      prop "Ray Continuity" prop_rayContinuity+      prop "Ray Bisection" prop_rayBisect+      prop "Rhumb Continuity" prop_rhumbContinuity+      prop "Rhumb Intersection" prop_rhumbIntersect+   describe "GridOffset monoid" $ mapM_ (uncurry prop) $ unbatch $ monoid (mempty :: GridOffset)+   describe "Helmert monoid" $ mapM_ (uncurry prop) $ unbatch $ monoid (mempty :: Helmert)  instance EqProp GridOffset where   (GridOffset a b c) =-= (GridOffset a' b' c') =@@ -50,86 +85,65 @@                    rX' ≈- rX'', rY' ≈- rY'', rZ' ≈- rZ'']      where x ≈ y = abs (x - y) < 0.00001-          x ≈- y = abs (x - y) < (1 / (5 * 2) ^ _5)--tests :: [Test]-tests = [-   testGroup "Geodetic" [-      testProperty "WGS84 and back" prop_WGS84_and_back,-      testProperty "Zero ground distance" prop_zero_ground,-      testGroup "UK Points" $ map pointTest ukPoints],-      testGroup "World lines" $ map worldLineTests worldLines,-   testGroup "Grid" [-      testProperty "Grid Offset 1" prop_offset1,-      testProperty "Grid Offset 2" prop_offset2,-      testProperty "Grid Offset 3" prop_offset3,-      testProperty "Grid 1" prop_grid1 ],-   testGroup "TransverseMercator" [-      testCase "fromGrid . toGrid == id" $ HU.assertBool "" prop_tmGridInverse-      ],-   testGroup "UK" [-      testProperty "UK Grid 1" prop_ukGrid1,-      testGroup "UK Grid 2" $ map ukGridTest2 ukSampleGrid,-      testGroup "UK Grid 3" $ map ukGridTest3 ukSampleGrid,-      testGroup "UK Grid 4" $ map ukGridTest4 ukSampleGrid,-      testGroup "UK Grid 5" $ map ukGridTest5 ukSampleGrid-      ],-   testGroup "Stereographic" [-      testCase "toGrid north" $ HU.assertBool "" stereographicToGridN,-      testCase "fromGrid north" $ HU.assertBool "" stereographicFromGridN,-      testCase "toGrid south" $ HU.assertBool "" stereographicToGridS,-      testCase "fromGrid south" $ HU.assertBool "" stereographicFromGridS,-      testProperty "Stereographic round trip" prop_stereographic-      ],-   testGroup "Paths" [-      testProperty "Ray Path 1" prop_rayPath1,-      testProperty "Ray Continuity" prop_rayContinuity,-      testProperty "Ray Bisection" prop_rayBisect,-      testProperty "Rhumb Continuity" prop_rhumbContinuity,-      testProperty "Rhumb Intersection" prop_rhumbIntersect-      ],-   testGroup "GridOffset" $ map (uncurry testProperty) $ unbatch $ monoid (mempty :: GridOffset),-   testGroup "Helmert" $ map (uncurry testProperty) $ unbatch $ monoid (mempty :: Helmert)-   ]+          x ≈- y = abs (x - y) < 1 / (5 * 2) ^ _5   -- | The positions are within 30 cm.-samePlace :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool-samePlace p1 p2 = geometricalDistance p1 p2 < 0.3+samePlace :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Expectation+samePlace p1 p2 = expectTrue msg $ geometricalDistance p1 p2 < 0.3+   where+      msg = "location " <> show p2 <> " is > 30cm from expected " <> show p1 +samePlace' :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool+samePlace' p1 p2 = geometricalDistance p1 p2 < 0.3  -- | The positions are within 10 m.-closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool-closeEnough p1 p2 = geometricalDistance p1 p2 < 10+closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Expectation+closeEnough p1 p2 = expectTrue msg $ geometricalDistance p1 p2 < 10+   where+      msg = "location " <> show p2 <> " is > 10m from expected " <> show p1 +closeEnough' :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool+closeEnough' p1 p2 = geometricalDistance p1 p2 < 10  -- | The angles are within 0.01 arcsec-sameAngle :: Double -> Double -> Bool-sameAngle v1 v2 = abs (properAngle (v1 - v2)) < 0.01 * arcsecond+sameAngle :: Double -> Double -> Expectation+sameAngle v1 v2 = expectTrue msg $ abs (properAngle (v1 - v2)) < 0.01 * arcsecond+   where+      msg = "expected angle " <> show (v1 / degree) <> ", got " <> show (v2 / degree) --- | The grid positions are within 1mm-sameGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool-sameGrid p1 p2 = check eastings && check northings && check altitude-   where check f = f p1 - f p2 < 1e-3+sameAngle' :: Double -> Double -> Bool+sameAngle' v1 v2 = abs (properAngle (v1 - v2)) < 0.01 * arcsecond +-- | The grid positions are within 1mm+sameGrid :: (Show r) => GridPoint r -> GridPoint r -> Expectation+sameGrid p1 p2 = expectTrue msg $ check eastings && check northings && check altitude+   where+      msg = "expected " <> show p1 <> ", got " <> show p2+      check f = f p1 - f p2 < 1e-3  -- | Grid offsets are within 1mm.-sameOffset :: GridOffset -> GridOffset -> Bool-sameOffset go1 go2 = check deltaNorth && check deltaEast && check deltaAltitude-   where check f = f go1 - f go2 < 1e-3+sameOffset :: GridOffset -> GridOffset -> Expectation+sameOffset go1 go2 = expectTrue msg $ check deltaNorth && check deltaEast && check deltaAltitude+   where+      msg = "expected " <> show go1 <> ", got " <> show go2+      check f = f go1 - f go2 < 1e-3   -- | The grid X and Y are both within 1 meter-closeGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool-closeGrid p1 p2 = check eastings && check northings && check altitude-   where check f = f p1 - f p2 < 1+closeGrid :: (Show r) => GridPoint r -> GridPoint r -> Expectation+closeGrid p1 p2 = expectTrue msg $ check eastings && check northings && check altitude+   where+      msg = "expected " <> show p1 <> ", got " <> show p2+      check f = f p1 - f p2 < 1 + -- | Degrees, minutes and seconds into radians. dms :: Int -> Int -> Double -> Double dms d m s = fromIntegral d * degree + fromIntegral m * arcminute + s * arcsecond  -- | Round-trip from local to WGS84 and back is identity (approximately)-prop_WGS84_and_back :: Geodetic LocalEllipsoid -> Bool+prop_WGS84_and_back :: Geodetic LocalEllipsoid -> Expectation prop_WGS84_and_back p = samePlace p $ toLocal (ellipsoid p) $ toWGS84 p  @@ -154,8 +168,8 @@       10001965.729, 0 * degree, 180 * degree)]  -worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Double, Double, Double) -> Test-worldLineTests (str, g1, g2, d, a, b) = testCase str $ HU.assertBool "" $ ok $ groundDistance g1 g2+worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Double, Double, Double) -> SpecWith (Arg Expectation)+worldLineTests (str, g1, g2, d, a, b) = it str $ ok $ groundDistance g1 g2    where       ok Nothing = False       ok (Just (d1, a1, b1)) =@@ -182,37 +196,37 @@   -- Convert a named point into a test-pointTest :: (Ellipsoid e2) => (String, Geodetic WGS84, Geodetic e2) -> Test-pointTest (testName, wgs84, local) =  testCase testName $ HU.assertBool "" $ samePlace wgs84 (toWGS84 local)+pointTest :: (Ellipsoid e2) => (String, Geodetic WGS84, Geodetic e2) -> SpecWith (Arg Expectation)+pointTest (testName, wgs84, local) =  it testName $ wgs84 `samePlace` toWGS84 local   -- The negation of the sum of a list of offsets is equal to the sum of the negated items.-prop_offset1 :: [GridOffset] -> Bool+prop_offset1 :: [GridOffset] -> Expectation prop_offset1 offsets = sameOffset (offsetNegate $ mconcat offsets) (mconcat $ map offsetNegate offsets)  -- A polar offset multiplied by a scalar is equal to an offset in the same direction with the length multiplied.-prop_offset2 :: Distance -> Bearing -> Scalar -> Bool+prop_offset2 :: Distance -> Bearing -> Scalar -> Expectation prop_offset2 (Distance d) (Bearing h) (Scalar s) = sameOffset go1 go2    where       go1 = offsetScale s $ polarOffset d h       go2 = polarOffset (d * s) h  -- | A polar offset has the offset distance and bearing of its arguments.-prop_offset3 :: GridOffset -> Bool+prop_offset3 :: GridOffset -> Expectation prop_offset3 delta = sameOffset delta0                                 (polarOffset (offsetDistance delta0) (offsetBearing delta))    where delta0 = delta {deltaAltitude = 0}  -- | Given a grid point and an offset, applying the offset to the point gives a new point which -- is offset from the first point by the argument offset.-prop_grid1 :: GridPoint (GridTM LocalEllipsoid) -> GridOffset -> Bool+prop_grid1 :: GridPoint (GridTM LocalEllipsoid) -> GridOffset -> Expectation prop_grid1 p d = sameOffset d $ p `gridOffset` applyOffset d p  -- | Check that using toGrid/fromGrid for TransverseMercator projection are inverses -- | for negative latitudes near the coordinates 0,0-prop_tmGridInverse :: Bool-prop_tmGridInverse = -   let origin = Geodetic +prop_tmGridInverse :: Expectation+prop_tmGridInverse =+   let origin = Geodetic          { latitude = 0 * degree          , longitude = 0 * degree          , geoAlt = 0@@ -223,12 +237,12 @@        tp1 = toGrid g testPoint        tp2 = fromGrid tp1    in tp2 `closeEnough` testPoint-   + -- | Converting a UK grid reference to a GridPoint and back is a null operation.-prop_ukGrid1 :: GridRef -> Bool-prop_ukGrid1 (GridRef str) =-   str ==-   (fromJust $ toUkGridReference ((length str - 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str)+prop_ukGrid1 :: UkGridRef -> Expectation+prop_ukGrid1 (UkGridRef str) =+   str `shouldBe`+   fromJust (toUkGridReference ((length str - 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str)  -- | UK Grid Reference points. The oracle for these points was the -- UK Grid Reference Finder (gridreferencefinder.com), retrieved on 26 Jan 2013.@@ -249,35 +263,35 @@    -- Caister and Framingham are taken from Ordnance Survey worked examples.    where       convert (grid, x, y, lat, long, desc) =-         (grid, GridPoint (x) (y) (0) UkNationalGrid,-          Geodetic (lat * degree) (long * degree) (0) WGS84, desc)+         (grid, GridPoint x y 0 UkNationalGrid,+          Geodetic (lat * degree) (long * degree) 0 WGS84, desc) -type GridPointTest = (String, GridPoint UkNationalGrid, Geodetic WGS84, String) -> Test+type UkGridPointTest = (String, GridPoint UkNationalGrid, Geodetic WGS84, String) -> SpecWith (Arg Expectation)  -- | Check that grid reference to grid point works for sample points.-ukGridTest2 :: GridPointTest-ukGridTest2 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool ""-   $ (fst $ fromJust $ fromUkGridReference gridRef) == gp+ukGridTest2 :: UkGridPointTest+ukGridTest2 (gridRef, gp, _, testName) =+   it testName $ fst (fromJust $ fromUkGridReference gridRef) `shouldBe` gp  -- | Check that grid point to grid reference works for sample points.-ukGridTest3 :: GridPointTest-ukGridTest3 (gridRef, gp, _, testName) = testCase testName $ HU.assertBool ""-   $ toUkGridReference 5 gp == Just gridRef+ukGridTest3 :: UkGridPointTest+ukGridTest3 (gridRef, gp, _, testName) =+   it testName $ toUkGridReference 5 gp `shouldBe` Just gridRef  -- | Check that grid point to WGS84 works close enough for sample points.-ukGridTest4 :: GridPointTest-ukGridTest4 (_, gp, geo, testName) = testCase testName $ HU.assertBool ""-   $ closeEnough geo $ toWGS84 $ fromGrid gp+ukGridTest4 :: UkGridPointTest+ukGridTest4 (_, gp, geo, testName) =+   it testName $ geo `closeEnough` toWGS84 (fromGrid gp)  -- | Check that WGS84 to grid point works close enough for sample points.-ukGridTest5 :: GridPointTest-ukGridTest5 (_, gp, geo, testName) = testCase testName $ HU.assertBool ""-   $ offsetDistance (gridOffset gp $ toGrid UkNationalGrid $ toLocal OSGB36 geo) < 1+ukGridTest5 :: UkGridPointTest+ukGridTest5 (_, gp, geo, testName) =+   it testName $ gp `closeGrid` toGrid UkNationalGrid (toLocal OSGB36 geo)   -- | Worked example for UK Geodetic to GridPoint, taken from "A Guide to Coordinate Systems in Great Britain" [1] ukTest :: Geodetic OSGB36-ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) (0) OSGB36+ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) 0 OSGB36  {-    v = 6.3885023333E+06@@ -296,56 +310,138 @@ -}  +-- | Check that a UTM grid reference round-trips to a gridpoint and back.+prop_utmGridTest1 :: UtmGridRef -> Expectation+prop_utmGridTest1 (UtmGridRef str) =+   str `shouldBe`+   toUtmGridReference Nothing True 0 (fromRight (error str) $ fromUtmGridReference str)++-- | Sample points for UTM, in both UTM and MGRS formats. The oracle for these points was the Montana+-- State University converter. http://rcn.montana.edu/resources/Converter.aspx.+-- Retrieved on 1st Feb 2025.+utmSampleGrid :: [(String, String, GridPoint UtmZone, Geodetic WGS84, String)]+utmSampleGrid = map convert [+ -- UTM Reference,        MGRS Reference,       X,      Y,        Latitude,    Longitude,   Description  +   ("30N 699304 5710208", "30U XC 99304 10208", 699304, 5710208,   51.5078064,  -0.1279388, "Nelson's Column"),+   ("35S 379101 8017747", "35K LA 79101 17747", 379101, 8017747,  -17.9249501,  25.8585071, "Victoria Falls"),+   ("27N 454366 7111715", "27W VM 54366 11715", 454366, 7111715,   64.1289075, -21.9373288, "Reykjavik Airport"),+   ("32N 297697 6700532", "32V KN 97697 00532", 297697, 6700532,   60.3904298,   5.3284215, "Bergen"),+   ("23S 683473 7460697", "23K PQ 83473 60697", 683473, 7460697,  -22.9518122, -43.2105383, "Christ the Redeemer"),+   ("21S 439699 4272868", "21F VC 39699 72868", 439699, 4272868,  -51.6919073, -57.8724006, "Port Stanley"),+   ("59S 461474 4919822", "59G MK 61474 19822", 461474, 4919822,  -45.8740880, 170.5035807, "Dunedin"),+   ("31N 166022 0"      , "31N AA 66022 00000", 166022,       0,    0.0      ,   0.0      , "The Origin"),+   ("31S 166022 9999999", "31M AV 66022 99999", 166022, 9999999,  -0.00000903,   0.0      , "1 meter south")]+   where+      convert (utm, mgrs, x, y, lat, long, desc) =+         (utm, mgrs, GridPoint x y 0 zone, geo, desc)+         where+            geo = Geodetic (lat * degree) (long * degree) 0 WGS84+            znum = fromJust $ utmZoneNumber geo+            hemi = if lat < 0 then UtmSouth else UtmNorth+            zone = fromJust $ mkUtmZone hemi znum++type UtmGridPointTest = (String, String, GridPoint UtmZone, Geodetic WGS84, String) -> SpecWith (Arg Expectation)+++-- | Check that UTM reference to grid point works for sample points.+utmGridTest2 :: UtmGridPointTest+utmGridTest2 (gridRef, _, gp, _, testName) =+   it testName $ fromUtmGridReference gridRef `shouldBe` Right gp++-- | Check that Grid point to UTM reference works for sample points.+utmGridTest3 :: UtmGridPointTest+utmGridTest3 (gridRef, _, gp, _, testName) =+   it testName $ toUtmGridReference Nothing False 0 gp `shouldBe` gridRef++-- | Check that grid point to WGS84 works close enough for sample points.+utmGridTest4 :: UtmGridPointTest+utmGridTest4 (_, _, gp, geo, testName) =+   it testName $ geo `closeEnough` fromGrid gp++-- | Check that WGS84 to grid point works close enough for sample points.+utmGridTest5 :: UtmGridPointTest+utmGridTest5 (_, _, gp, geo, testName) =+   it testName $ gp `closeGrid` toGrid (fromJust $ utmZone geo) geo+++-- | Check that a UTM grid point round-trips to MGRS and back, with spaces.+prop_mgrs_gridTest1 :: UtmGridRef -> Expectation+prop_mgrs_gridTest1 (UtmGridRef str) =+   case fromUtmGridReference str of+      Left msg -> assertFailure $ "gridTest1 bogus UTM ref: " <> str <> ". Messages = " <> show msg+      Right 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.+prop_mgrs_gridTest2 :: UtmGridRef -> Expectation+prop_mgrs_gridTest2 (UtmGridRef str) =+   case fromUtmGridReference str of+      Left msg -> assertFailure $ "gridTest2 bogus UTM ref: " <> str <> ". Messages = " <> show msg+      Right 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) =+   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+   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)++ -- | Standard stereographic grid for point tests in the Northern Hemisphere. stereoGridN :: GridStereo LocalEllipsoid-stereoGridN = mkGridStereo tangent origin (0.9999079)+stereoGridN = mkGridStereo tangent origin 0.9999079    where-      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155) (299.15281) mempty-      tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) (0) ellipse-      origin = GridOffset (155000) (463000) (0)+      ellipse = LocalEllipsoid "Bessel 1841" 6377397.155 299.15281 mempty+      tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) 0 ellipse+      origin = GridOffset 155000 463000 0   -- | Standard steregraphic grid for point tests in the Southern Hemisphere. -- -- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated. stereoGridS :: GridStereo LocalEllipsoid-stereoGridS = mkGridStereo tangent origin (0.9999079)+stereoGridS = mkGridStereo tangent origin 0.9999079    where-      ellipse = LocalEllipsoid "Bessel 1841" (6377397.155) (299.15281) mempty-      tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) (0) ellipse-      origin = GridOffset ((-155000)) (463000) (0)+      ellipse = LocalEllipsoid "Bessel 1841" 6377397.155 299.15281 mempty+      tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) 0 ellipse+      origin = GridOffset (-155000) 463000 0   -- | Data for the stereographic tests taken from -- <http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf>-stereographicToGridN :: Bool+stereographicToGridN :: Expectation stereographicToGridN = sameGrid g1 g1'    where-      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridN-      g1 = GridPoint (196105.283) (557057.739) (0) stereoGridN+      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridN+      g1 = GridPoint 196105.283 557057.739 0 stereoGridN       g1' = toGrid stereoGridN p1 -stereographicFromGridN :: Bool+stereographicFromGridN :: Expectation stereographicFromGridN = samePlace p1 p1'    where-      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridN-      g1 = GridPoint (196105.283) (557057.739) (0) stereoGridN+      p1 = Geodetic (dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridN+      g1 = GridPoint 196105.283 557057.739 0 stereoGridN       p1' = fromGrid g1  -stereographicToGridS :: Bool+stereographicToGridS :: Expectation stereographicToGridS = sameGrid g1 g1'    where-      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridS-      g1 = GridPoint ((-196105.283)) (557057.739) (0) stereoGridS+      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridS+      g1 = GridPoint (-196105.283) 557057.739 0 stereoGridS       g1' = toGrid stereoGridS p1  -stereographicFromGridS :: Bool+stereographicFromGridS :: Expectation stereographicFromGridS = samePlace p1 p1'    where-      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0) $ gridEllipsoid stereoGridS-      g1 = GridPoint ((-196105.283)) (557057.739) (0) stereoGridS+      p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridS+      g1 = GridPoint (-196105.283) 557057.739 0 stereoGridS       p1' = fromGrid g1  @@ -362,7 +458,7 @@ -- | A ray at distance zero returns its original arguments. prop_rayPath1 :: Ray WGS84 -> Bool prop_rayPath1 r@(Ray pt b e) =-      samePlace pt pt1 && sameAngle b b1 && sameAngle e e1+      samePlace' pt pt1 && sameAngle' b b1 && sameAngle' e e1    where (pt1,b1,e1) = pathFunc (getRay r) 0  @@ -378,7 +474,7 @@ prop_pathContinuity pf pt0 (Bearing b0) (Azimuth a0) (Distance d1) (Distance d2) =    counterexample (show ((pt2, Bearing b2, Azimuth a2), (pt3, Bearing b3, Azimuth a3))) $       pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>-      closeEnough pt2 pt3 && sameAngle b2 b3 && sameAngle a2 a3+      closeEnough' pt2 pt3 && sameAngle' b2 b3 && sameAngle' a2 a3    where       path0 = pf pt0 b0 a0       (pt1, b1, a1) = pathFunc path0 d1@@ -393,7 +489,7 @@ prop_pathContinuity1 pf pt0 (Bearing b0) (Distance2 d1) (Distance2 d2) =    counterexample (show ((pt2, Bearing b2), (pt3, Bearing b3))) $       pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>-      closeEnough pt2 pt3 && sameAngle b2 b3+      closeEnough' pt2 pt3 && sameAngle' b2 b3    where       path0 = pf pt0 b0       (pt1, b1, _) = pathFunc path0 d1@@ -411,7 +507,7 @@ -- This is a test of bisection rather than rays. prop_rayBisect :: Ray WGS84 -> Altitude -> Bool prop_rayBisect r (Altitude height) =-   case bisect ray0 f (1e-2) (0) (1000 * kilometer) of+   case bisect ray0 f 1e-2 0 (1000 * kilometer) of       Nothing -> False       Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1e-2    where@@ -427,7 +523,7 @@ -- | Two rhumb paths intersect at the same place. prop_rhumbIntersect :: RhumbPaths2 -> Property prop_rhumbIntersect rp =-   case intersect 0 0 (0.1) 100 path1 path2 of+   case intersect 0 0 0.1 100 path1 path2 of       Just (d1, d2) ->          let (pt1, _, _) = pathFunc path1 d1              (pt2, _, _) = pathFunc path2 d2@@ -435,3 +531,8 @@       Nothing -> label "No intersection" True    where       (path1, path2) = mk2RhumbPaths rp+++-- | Copied from `Test.Hspec.Expectations` source. It ought to be exported from there.+expectTrue :: HasCallStack => String -> Bool -> Expectation+expectTrue msg b = unless b (expectationFailure msg)