diff --git a/ToDo.txt b/ToDo.txt --- a/ToDo.txt +++ b/ToDo.txt @@ -1,6 +1,6 @@ To Do ===== -* Add UTM and Universal Polar grids. +* Add Universal Polar grid. * Add more national grid systems. diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/geodetics.cabal b/geodetics.cabal --- a/geodetics.cabal +++ b/geodetics.cabal @@ -1,8 +1,9 @@ cabal-version: 3.0 name: geodetics -version: 1.0.0 +version: 1.1.0 build-type: Simple author: Paul Johnson +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 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 diff --git a/src/Geodetics/Ellipsoids.hs b/src/Geodetics/Ellipsoids.hs --- a/src/Geodetics/Ellipsoids.hs +++ b/src/Geodetics/Ellipsoids.hs @@ -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 diff --git a/src/Geodetics/Geodetic.hs b/src/Geodetics/Geodetic.hs --- a/src/Geodetics/Geodetic.hs +++ b/src/Geodetics/Geodetic.hs @@ -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 diff --git a/src/Geodetics/Grid.hs b/src/Geodetics/Grid.hs --- a/src/Geodetics/Grid.hs +++ b/src/Geodetics/Grid.hs @@ -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 diff --git a/src/Geodetics/LatLongParser.hs b/src/Geodetics/LatLongParser.hs --- a/src/Geodetics/LatLongParser.hs +++ b/src/Geodetics/LatLongParser.hs @@ -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 diff --git a/src/Geodetics/Path.hs b/src/Geodetics/Path.hs --- a/src/Geodetics/Path.hs +++ b/src/Geodetics/Path.hs @@ -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) => diff --git a/src/Geodetics/UK.hs b/src/Geodetics/UK.hs --- a/src/Geodetics/UK.hs +++ b/src/Geodetics/UK.hs @@ -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 diff --git a/src/Geodetics/UTM.hs b/src/Geodetics/UTM.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/UTM.hs @@ -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 + diff --git a/test/ArbitraryInstances.hs b/test/ArbitraryInstances.hs --- a/test/ArbitraryInstances.hs +++ b/test/ArbitraryInstances.hs @@ -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 +/- 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 - + diff --git a/test/LatLongParser.hs b/test/LatLongParser.hs new file mode 100644 --- /dev/null +++ b/test/LatLongParser.hs @@ -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)) diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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 -- -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)