diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -11,9 +11,9 @@ between them and the bearing from one to the other. Once you have a geodetic position defined you can project it onto a flat plane, or Grid. -At present Transverse Mercator and Oblique Stereographic grids are provided. More can be -added by defining new instances of the Grid typeclass: see "AddingProjections.txt" for -detais. +At present Transverse Mercator, Polar Stereographic and Oblique Stereographic grids are +provided. More can be added by defining new instances of the Grid typeclass: see +"AddingProjections.txt" for detais. The Paths module defines a path as a parametric function of distance that returns a position and a bearing. Given two paths you can find their intersection using a fast diff --git a/ToDo.txt b/ToDo.txt --- a/ToDo.txt +++ b/ToDo.txt @@ -1,6 +1,6 @@ To Do ===== -* Add Universal Polar grid. +* More national grid systems. -* Add more national grid systems. +* Great-circle paths. diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -39,3 +39,17 @@ `Geodetics.LatLongParser`. * Added Universal Transverse Mercator (UTM) grid, with both UTM and MGRS grid references + +## Version 2.0.0 + +This version has breaking changes relative to version 1.1.0. + +* Added PolarStereographic projection with UPS grid + +* Moved MGRS grid types and functions to its own module and extended it to include UPS gridpoints. + +* Replaced the UtmGridUnit type with a more general GridUnit type in Geodetics.Grid. + +* Provided standalone Parsec parsers for MGRS and UTM grid references. + +* Modified code formatting to have consistently 2 spaces per indent level. diff --git a/geodetics.cabal b/geodetics.cabal --- a/geodetics.cabal +++ b/geodetics.cabal @@ -1,6 +1,6 @@ cabal-version: 3.0 name: geodetics -version: 1.1.0 +version: 2.0.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 @@ -46,7 +46,9 @@ Geodetics.Geodetic, Geodetics.Grid, Geodetics.LatLongParser, + Geodetics.MGRS, Geodetics.Path, + Geodetics.PolarStereographic, Geodetics.Stereographic, Geodetics.TransverseMercator, Geodetics.UK, diff --git a/src/Geodetics/Grid.hs b/src/Geodetics/Grid.hs --- a/src/Geodetics/Grid.hs +++ b/src/Geodetics/Grid.hs @@ -18,7 +18,8 @@ unsafeGridCoerce, -- * Utility functions for grid references fromGridDigits, - toGridDigits + toGridDigits, + GridUnit (..) ) where import Data.Char @@ -67,11 +68,12 @@ mempty = GridOffset 0 0 0 mappend = (<>) + -- | An offset defined by a distance (m) and a bearing (radians) to the right of North. -- -- There is no elevation parameter because we are using a plane to approximate an ellipsoid, -- so elevation would not provide a useful result. If you want to work with elevations --- then "rayPath" will give meaningful results. +-- then 'Geodetics.Path.rayPath' will give meaningful results. polarOffset :: Double -> Double -> GridOffset polarOffset r d = GridOffset (r * sin d) (r * cos d) 0 @@ -124,22 +126,22 @@ -- but with the same easting, northing and altitude. This is unsafe because it -- will produce a different position unless the two grids are actually equal. -- --- It should be used only to convert between distinguished grids (e.g. "UkNationalGrid") and --- their equivalent numerical definitions. +-- This function should be used only to convert between distinguished grids +-- (e.g. "UkNationalGrid") and their equivalent numerical definitions. unsafeGridCoerce :: b -> GridPoint a -> GridPoint b unsafeGridCoerce base p = GridPoint (eastings p) (northings p) (altitude p) base - -- | Convert a list of digits to a distance. The first argument is the size of the -- grid square within which these digits specify a position. The first digit is -- in units of one tenth of the grid square, the second one hundredth, and so on. -- The first result is the lower limit of the result, and the second is the size -- of the specified offset. --- So for instance @fromGridDigits (100 * kilometer) "237"@ will return -- --- > Just (23700 meters, 100 meters) +-- For instance @fromGridDigits (100 * kilometer) "237"@ will return -- +-- > Just (23700, 100) +-- -- If there are any non-digits in the string then the function returns @Nothing@. fromGridDigits :: Double -> String -> Maybe (Double, Double) fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing @@ -151,6 +153,7 @@ (drop 1 $ iterate (/ 10) sq) p = sq / fromIntegral ((10 :: Integer) ^ n) + -- | Convert a distance into a digit string suitable for printing as part -- of a grid reference. The result is the south or west side of the enclosing grid square, -- where the size of the square is defined by the number of digits. @@ -175,3 +178,11 @@ (sqs, d1) = u `divMod` p s = show d1 pad = if n == 0 then "" else replicate (n - length s) '0' ++ s + + +-- | Some grids (notably UTM and UPS) can have optional units in their grid references. +data GridUnit = GridMeters | GridKilometers deriving Eq + +instance Show GridUnit where + show GridMeters = "m" + show GridKilometers = "km" diff --git a/src/Geodetics/MGRS.hs b/src/Geodetics/MGRS.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/MGRS.hs @@ -0,0 +1,409 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE NumericUnderscores #-} + +{- | The Military Grid Reference System (MGRS) + +In MGRS there are two syntaxes for grid references: + + 1. Between 80 South and 84 North a grid reference has a zone number, latitude band letter, a + 2 letter code for the 100km square within the zone, and then northings and eastings within + that square. + + 2. In the polar regions a grid reference has a latitude band letter (A or B for South, Y or Z + for North), a 2 letter code for the 100km square within the polar region, and then northings + and eastings within that square. There is no zone number in the polar regions. +-} +module Geodetics.MGRS ( + -- * MGRS Grid Basis + MgrsGrid (..), + mgrsGrid, + utmToMgrsPoint, + upsToMgrsPoint, + fromMgrsPoint, + -- * Textual Representation + mgrsBandLetterToLatitude, + mgrsLatitudeToBandLetter, + fromMgrsGridReference, + toMgrsGridReference +) where + +import Control.Monad +import Data.Array +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Geodetics.Grid +import Geodetics.PolarStereographic +import Geodetics.UTM +import Text.Parsec +import Text.Parsec.Error +import Text.Printf + + +-- | MGRS grid references can be anywhere on Earth. Hence the position can be either on a UTM or +-- a UPS grid. +data MgrsGrid = + MgrsUtm UtmZone + | MgrsUps UpsGrid + deriving (Eq, Show) + +instance GridClass MgrsGrid WGS84 where + fromGrid p = case gridBasis p of + MgrsUtm zone -> fromGrid $ unsafeGridCoerce zone p + MgrsUps grid -> fromGrid $ unsafeGridCoerce grid p + + toGrid (MgrsUtm zone) = unsafeGridCoerce (MgrsUtm zone) . toGrid zone + toGrid (MgrsUps grid) = unsafeGridCoerce (MgrsUps grid) . toGrid grid + + gridEllipsoid _ = WGS84 + + +-- | Find the most appropriate grid for the given geodetic position +mgrsGrid :: Geodetic WGS84 -> MgrsGrid +mgrsGrid geo = + case utmZone geo of + Just zone -> MgrsUtm zone + Nothing -> MgrsUps $ if latitude geo > 0 then upsNorth else upsSouth + + +-- | Generalise from UTM to MGRS +utmToMgrsPoint :: GridPoint UtmZone -> GridPoint MgrsGrid +utmToMgrsPoint gp = unsafeGridCoerce (MgrsUtm $ gridBasis gp) gp + + +-- | Generalise from UPS to MGRS +upsToMgrsPoint :: GridPoint UpsGrid -> GridPoint MgrsGrid +upsToMgrsPoint gp = unsafeGridCoerce (MgrsUps $ gridBasis gp) gp + + +-- | Convert an MGRS grid point to either UTM or UPS, depending on its basis. +fromMgrsPoint :: GridPoint MgrsGrid -> Either (GridPoint UtmZone) (GridPoint UpsGrid) +fromMgrsPoint gp = case gridBasis gp of + MgrsUtm zone -> Left $ unsafeGridCoerce zone gp + MgrsUps grid -> Right $ unsafeGridCoerce grid gp + + +-- | Convert a list of letters into an inverse lookup function. +letterTable :: [Char] -> Char -> Maybe Int +letterTable cs = \c -> if inRange (bounds arr) c then arr ! c else Nothing + where arr = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip cs [0..]] + -- Second argument in a lambda so that partial application returns a function with a CAF + -- rather than computing the array for every call. + + +-- | The MGRS latitude band code letters, excluding A and B used for Antarctica (south of -80 degrees) +-- and Y and Z used for the Arctic (north of 84 degrees). +mgrsBandLetters :: [Char] +mgrsBandLetters = "CDEFGHJKLMNPQRSTUVWX" + +mgrsBandIx :: Char -> Maybe Int +mgrsBandIx = letterTable mgrsBandLetters + +-- | Polar regions have special bands A, B, Y and Z. +mgrsBandLetterToPole :: Char -> Maybe Pole +mgrsBandLetterToPole 'A' = Just SouthPole +mgrsBandLetterToPole 'B' = Just SouthPole +mgrsBandLetterToPole 'Y' = Just NorthPole +mgrsBandLetterToPole 'Z' = Just NorthPole +mgrsBandLetterToPole _ = Nothing + +polarBandLetters :: [Char] +polarBandLetters = "ABYZ" + +-- | Polar regions use different northings and eastings letters to the rest of the world. +polarEastingLetters :: [Char] +polarEastingLetters = "ABCFGHJKLPQRSTUXYZ" + +polarEastingIx :: Char -> Maybe Int +polarEastingIx = letterTable polarEastingLetters + +polarNorthingLetters :: [Char] +polarNorthingLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ" + +polarNorthingIx :: Char -> Maybe Int +polarNorthingIx = letterTable polarNorthingLetters + + +-- | Convert the three letters of an MGRS grid reference to the south-west corner of a 100km square. +-- +-- If any of the letters are invalid for the positions they are in then an error message is returned. +mgrsLettersToPolar :: + Char -- ^ Latitude band: A, B, Y or Z. + -> Char -- ^ Eastings letter. + -> Char -- ^ Northings letter. + -> Parsec s u (Pole, Double, Double) +mgrsLettersToPolar bandC eastC northC = + case (pole, baseEasting, baseNorthing) of + (Just pole1, Just easting, Just northing) -> pure (pole1, easting, northing + polarNorthingLetterOrigin pole1) + _ -> fail "Invalid polar grid letters" + where + pole = mgrsBandLetterToPole bandC + baseEasting = (bandEasting +) . (100_000 *) . fromIntegral <$> polarEastingIx eastC + baseNorthing = (100_000 *) . fromIntegral <$> polarNorthingIx northC + bandEasting = if bandC `elem` "AY" then 200*kilometer else 2000*kilometer + -- A and Y denote the "eastern" halves of the south and north polar regions respectively. + + +-- | Find the southern boundary of a latitude band letter (excluding poles). +mgrsBandLetterToLatitude :: Char -> Maybe Double +mgrsBandLetterToLatitude band = do + n1 <- mgrsBandIx band + return $ degree * fromIntegral (-80 + n1 * 8) + + +-- | Find the band letter for a latitude, if it is in the range (-80, 84) degrees. +-- (Argument in radians) +mgrsLatitudeToBandLetter :: Double -> Maybe Char +mgrsLatitudeToBandLetter lat = do + guard $ -80 <= dlat && dlat <= 84 + return $ indexMap ! latIdx + where + dlat = lat / degree + ilat :: Int + ilat = floor dlat + latIdx = min 19 $ (ilat + 80) `div` 8 -- Band 19 (X) extends an extra 4 degrees. + indexMap = listArray (0,19) mgrsBandLetters + + + +-- | Letters A-Z except for I and O. +mgrsEastingsLetters :: [Char] +mgrsEastingsLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ" + +mgrsEastingIx :: Char -> Maybe Int +mgrsEastingIx = letterTable mgrsEastingsLetters + +-- | If zone number is in range and the letter is one of the valid Eastings letters for that zone +-- then return the UTM easting in meters. +-- +-- Zone 1 starts with \'A\'. Each zone is 8 characters wide. Hence the letters repeat every 3 zones. +mgrsLetterToEasting :: UtmZoneNumber -> Char -> Maybe Double +mgrsLetterToEasting zn c = do + guard $ 1 <= zn && zn <= 60 + n1 <- mgrsEastingIx c + let n2 = n1 - base + guard $ 0 <= n2 && n2 <= 7 + return $ fromIntegral (n2+1) * 100_000 + where + base = ((zn-1) `mod` 3) * 8 + + +-- | If the zone number is in range and the eastings are between 100,000 and 900,000 then +-- return the Eastings letter. +mgrsEastingToLetter :: UtmZoneNumber -> Double -> Maybe Char +mgrsEastingToLetter zn east = do + guard $ 1 <= zn && zn <= 60 + guard $ 100 * kilometer <= east && east < 900 * kilometer + return $ indexMap ! ix + where + indexMap = listArray (0,23) mgrsEastingsLetters + base = ((zn-1) `mod` 3) * 8 + square = max 0 $ min 7 $ floor $ (east - 100 * kilometer)/(100 * kilometer) -- Clamped in range (0,7). + ix = base + square -- Must be in range (0,23) + + +-- | Letters A-V except for I and O. +mgrsNorthingsLetters :: [Char] +mgrsNorthingsLetters = "ABCDEFGHJKLMNPQRSTUV" + +mgrsNorthingIx :: Char -> Maybe Int +mgrsNorthingIx = letterTable mgrsEastingsLetters + + +-- | MGRS Northings letters have rather complex relationship to the latitude bands. The 20 letters +-- repeat every 2,000km going north and south from the equator, so the latitude band is needed +-- to disambiguate which repetition of the Northings letter is meant. + +-- Unfortunately this repetition of the letters does not neatly coincide with +-- the latitude band boundaries, which are based on degrees of latitude. + +-- The base letter just north of the equator is A in odd-numbered zones and F in even numbered zones. +-- +-- This uses the latitude band to estimate the range of northings that would be valid there, +-- and hence determine which possible grid band is meant by the northings letter. +-- The algorithm used is approximate and deliberately very forgiving: it will accept some grid squares +-- which are north or south of the band given. +mgrsLetterToNorthings :: + UtmZoneNumber + -> Char -- ^ Latitude band letter (@C@ - @X@ excluding @I@ and @O@). + -> Char -- ^ MGRS Northings letter (@A@ - @V@ excluding @I@ and @O@). + -> Maybe Double +mgrsLetterToNorthings zone bandC northingsC = do + guard $ 1 <= zone && zone <= 60 + band <- (/degree) <$> mgrsBandLetterToLatitude bandC + northings0 <- (baseNorthingsOffset +) <$> mgrsNorthingIx northingsC + let bandDist = band * metersPerDegree -- Approx dist from equator to southern edge of band. + bandGridLower = floor $ bandDist / 100_000 - 2 -- Lower limit of band in 100km units + bandGridUpper = ceiling $ if band > 71 * degree -- Upper limit of band in 100km units. + then (bandDist + 12 * metersPerDegree) / 100_000 + 1 -- Band X. + else (bandDist + 8 * metersPerDegree) / 100_000 + 2 -- Other bands. + rep = (bandGridLower - northings0 - 1) `div` 20 -- Lower limit in 2,000,000km units. + grid = (rep+1)*20 + northings0 + guard $ grid >= bandGridLower + guard $ grid <= bandGridUpper + return $ fromIntegral $ grid * 100_000 + where + metersPerDegree = 10_002_000 / 90 -- Equator to north pole. + baseNorthingsOffset :: Int + baseNorthingsOffset = if odd zone then 0 else -5 + + +-- | Find the northings letter of the 100km square containing the given Northings. +-- +-- The input is not range checked. It just assumes that the northings letters repeat forever. +mgrsNorthingToLetter :: UtmZoneNumber -> Double -> Char +mgrsNorthingToLetter zone northings1 = + letters ! ((gridNum + baseNorthingsOffset) `mod` 20) + where + gridNum :: Int + gridNum = floor $ northings1 / (100 * kilometer) + baseNorthingsOffset = if odd zone then 0 else 5 + letters = listArray (0,19) mgrsNorthingsLetters + + +polarEastingsToLetter :: Double -> Char +polarEastingsToLetter eastings1 = + letters ! ((floor (eastings1 / (100 * kilometer)) - 20) `mod` 18 :: Int) + where + letters = listArray (0,17) polarEastingLetters + + +-- | The southern edge of the @A@ northings letter differs between north and south poles. +polarNorthingLetterOrigin :: Pole -> Double +polarNorthingLetterOrigin NorthPole = 1_300 * kilometer +polarNorthingLetterOrigin SouthPole = 800 * kilometer + + +polarNorthingsToLetter :: Pole -> Double -> Maybe Char +polarNorthingsToLetter pole northings1 = do + let i :: Int + i = floor ((northings1 - polarNorthingLetterOrigin pole) / (100 * kilometer)) `mod` 24 + guard $ i >= 0 && i < 24 + pure $ letters ! i + where + letters = listArray (0,23) polarNorthingLetters + +-- | Convert an MGRS grid reference to a UTM @GridPoint@, if the reference is valid. +-- E.g. \"30U XC 99304 10208\" is the grid reference for Nelson's Column in London. +-- +-- If the input contains spaces then these are used to delimit the fields. Multiple +-- spaces are treated as a single space. +-- +-- If the reference is valid this returns the position of the south-west corner of the +-- nominated grid square and an offset to its centre. Altitude is set to zero. +fromMgrsGridReference :: String -> Either [String] (GridPoint MgrsGrid, GridOffset) +fromMgrsGridReference str = case parse parseMgrsGridReference "" str of + Left err -> Left $ filter (not . null) $ lines $ showErrorMessages + "or" "unknown parse error" "expecting" "unexpected" "end of input" + (errorMessages err) + Right r -> Right r + +parseMgrsGridReference :: Parsec String u (GridPoint MgrsGrid, GridOffset) +parseMgrsGridReference = mgrsUtmParser <|> mgrsUpsParser + where + mgrsUtmParser = do + zoneNum <- read <$> many1 digit "UTM zone number" + -- Safe because we can only read digits + spaces + band <- oneOf mgrsBandLetters "latitude band letter" + let (hemi, falseNorthing) = if band >= 'N' + then (UtmNorth, 0) + else (UtmSouth, -10_000_000) + zone <- maybe (fail "Invalid zone") (pure . MgrsUtm) $ mkUtmZone hemi zoneNum + spaces + squareEast <- oneOf mgrsEastingsLetters "eastings letter" + eastingBase <- maybe (fail "Invalid eastings letter") pure $ + mgrsLetterToEasting zoneNum squareEast + spaces + squareNorth <- oneOf mgrsNorthingsLetters "northings letter" + northingBase <- maybe (fail "Invalid northings letter") pure $ + mgrsLetterToNorthings zoneNum band squareNorth + spaces + (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits + when (length eastingChars /= length northingChars) $ + fail "Northings and Eastings must be the same length." + if northingChars == "" && eastingChars == "" + then -- No digits, just return the outer 100km grid square + pure (GridPoint eastingBase (northingBase - falseNorthing) 0 zone, + GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0) + else do + (northing, offset) <- maybe (fail "Invalid northing digits") pure $ + fromGridDigits (100 * kilometer) northingChars + (easting, _) <- maybe (fail "Invalid easting digits") pure $ + fromGridDigits (100 * kilometer) eastingChars + pure (GridPoint (eastingBase + easting) (northingBase + northing - falseNorthing) 0 zone, + GridOffset (offset/2) (offset/2) 0) + + mgrsUpsParser = do + spaces + band <- oneOf polarBandLetters "polar band letter" + spaces + squareEast <- oneOf polarEastingLetters "eastings letter" + spaces + squareNorth <- oneOf polarNorthingLetters "northings letter" + spaces + (pole, eastingBase, northingBase) <- mgrsLettersToPolar band squareEast squareNorth + let grid = case pole of + NorthPole -> MgrsUps upsNorth + SouthPole -> MgrsUps upsSouth + (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits + when (length eastingChars /= length northingChars) $ + fail "Northings and Eastings must be the same length." + if northingChars == "" && eastingChars == "" + then -- No digits, just return the outer 100km grid square + pure (GridPoint eastingBase northingBase 0 grid, + GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0) + else do + (northing, offset) <- maybe (fail "Invalid northing digits") pure $ + fromGridDigits (100 * kilometer) northingChars + (easting, _) <- maybe (fail "Invalid easting digits") pure $ + fromGridDigits (100 * kilometer) eastingChars + pure (GridPoint (eastingBase + easting) (northingBase + northing) 0 grid, + GridOffset (offset/2) (offset/2) 0) + spaced = do + e <- many1 digit + skipMany1 space -- A space is mandatory here. + n <- many1 digit + pure (e,n) + unspaced = do + digits <- many1 digit + let c = length digits + when (odd c) $ fail "Northings and Eastings must be the same length." + pure (splitAt (c `div` 2) digits) + noDigits = do + eof + pure ("", "") + +-- | Convert a UTM or UPS @GridPoint@ to an MGRS grid reference. +toMgrsGridReference :: + Bool -- ^ Include spaces in the output. The standard says no spaces, but they make + -- the output easier to read. + -> Int -- ^ Number of digits of precision in the easting and northing. Must be 0-8. + -> GridPoint MgrsGrid + -> Maybe String +toMgrsGridReference withSpaces precision gp = do + guard $ precision >= 0 && precision <= 8 + case gridBasis gp of + MgrsUtm zone -> do + band <- mgrsLatitudeToBandLetter $ latitude $ fromGrid gp + let + zoneNum = utmZoneNum zone + northLetter = mgrsNorthingToLetter zoneNum $ northings gp + eastLetter <- mgrsEastingToLetter zoneNum $ eastings gp + (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp + (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp + let part1 = printf "%02d" zoneNum <> [band] + part2 = [eastLetter, northLetter] + pure $ if withSpaces + then part1 <> " " <> part2 <> " " <> eastDigits <> " " <> northDigits + else part1 <> part2 <> eastDigits <> northDigits + MgrsUps grid -> do + let zoneLetter = case trueOrigin grid of + NorthPole -> if eastings gp < 2_000 * kilometer then 'Y' else 'Z' + SouthPole -> if eastings gp < 2_000 * kilometer then 'A' else 'B' + eastLetter = polarEastingsToLetter $ eastings gp + northLetter <- polarNorthingsToLetter (trueOrigin grid) $ northings gp + (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp + (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp + pure $ if withSpaces + then zoneLetter : ' ' : eastLetter : northLetter : " " <> eastDigits <> " " <> northDigits + else zoneLetter : eastLetter : northLetter : eastDigits <> northDigits diff --git a/src/Geodetics/Path.hs b/src/Geodetics/Path.hs --- a/src/Geodetics/Path.hs +++ b/src/Geodetics/Path.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE TypeOperators, TypeFamilies, FlexibleContexts #-} +{-# LANGUAGE TypeFamilies, FlexibleContexts #-} -- | The implementation assumes IEEE 754 arithmetic. module Geodetics.Path where @@ -26,7 +26,7 @@ -- ^ Takes a length and returns a position, and direction as heading and elevation angles. pathValidity :: PathValidity } - + -- | Convenience value for paths that are valid for all distances. alwaysValid :: PathValidity alwaysValid = (negate inf, inf) where @@ -46,8 +46,8 @@ -- -- The initial bounds must return one GT or EQ value and one LT or EQ value. If they -- do not then @Nothing@ is returned. -bisect :: - Path e +bisect :: + Path e -> (Geodetic e -> Ordering) -- ^ Evaluation function. -> Double -- ^ Required accuracy in terms of distance along the path. -> Double -> Double -- ^ Initial bounds. @@ -59,7 +59,7 @@ guard $ hasRoot r bisect1 $ sortPair r where - f' d = let (p, _, _) = pathFunc path d in f p + f' d = let (p, _, _) = pathFunc path d in f p pairResults d1 d2 = ((d1, f' d1), (d2, f' d2)) hasRoot (v1, v2) = snd v1 <= EQ && EQ <= snd v2 sortPair (v1, v2) = if snd v1 <= snd v2 then (v1, v2) else (v2, v1) @@ -68,7 +68,7 @@ r3 = f' d3 c1 = ((d1, r1), (d3, r3)) c2 = ((d3, r3), (d2, r2)) - in if abs (d1 - d2) <= t + in if abs (d1 - d2) <= t then return d3 else bisect1 $ if hasRoot c1 then c1 else c2 @@ -96,14 +96,14 @@ | mag < 1e-15 = Nothing | mag3 (nv1 `cross3` nv2) * r <= accuracy = Just (d1, d2) -- Assumes that sin (accuracy/r) == accuracy/r - | otherwise = + | otherwise = if abs d1a + abs d2a < abs d1b + abs d2b then intersect (d1 + d1a) (d2 + d2a) accuracy (pred n) path1 path2 else intersect (d1 + d1b) (d2 + d2b) accuracy (pred n) path1 path2 where (pt1, h1, _) = pathFunc path1 d1 (pt2, h2, _) = pathFunc path2 d2 - vectors :: Double -> Double -> Double + vectors :: Double -> Double -> Double -> (Vec3 Double, Vec3 Double) vectors lat lon b = ( -- Unit vector of normal to surface at (lat,lon) @@ -132,11 +132,11 @@ d1b = gcDist gc1 nv1 nv3b * r d2b = gcDist gc2 nv2 nv3b * r -- Signed angle between v1 and v2, - gcDist norm v1 v2 = - let c = v1 `cross3` v2 - in (if c `dot3` norm < 0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) + gcDist norm v1 v2 = + let c = v1 `cross3` v2 + in (if c `dot3` norm < 0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) r = majorRadius $ ellipsoid pt1 - + {- Note on derivation of "intersect" The algorithm is a variant of the Newton-Raphson method, and shares its advantage @@ -167,7 +167,7 @@ -} -- | A ray from a point heading in a straight line in 3 dimensions. -rayPath :: (Ellipsoid e) => +rayPath :: (Ellipsoid e) => Geodetic e -- ^ Start point. -> Double -- ^ Bearing. -> Double -- ^ Elevation. @@ -181,7 +181,7 @@ (dE,dN,dU) = transform3 (trans3 $ ecefMatrix lat long) delta -- Direction of ray at result point. elevation2 = asin dU bearing2 = if dE == 0 && dN == 0 then bearing else atan2 dE dN -- Allow for vertical elevation. - + ecefMatrix lat long = -- Transform matrix for vectors from (East, North, Up) to (X,Y,Z). ((negate sinLong, negate cosLong*sinLat, cosLong*cosLat), -- East X North X Up X @@ -194,7 +194,7 @@ cosLong = cos long sinLat = sin lat cosLat = cos lat - + direction = (sinB*cosE, cosB*cosE, sinE) -- Direction of ray in ENU delta = transform3 (ecefMatrix (latitude pt1) (longitude pt1)) direction -- Convert to ECEF pt1' = geoToEarth pt1 -- ECEF of origin point. @@ -202,7 +202,7 @@ cosB = cos bearing sinE = sin elevation cosE = cos elevation - + -- | Rhumb line: path following a constant course. Also known as a loxodrome. -- -- The valid range stops a few arc-minutes short of the poles to ensure that the @@ -244,8 +244,8 @@ m0 = meridianRadius (ellipsoid pt) lat0 a = majorRadius $ ellipsoid pt b = minorRadius $ ellipsoid pt - + -- | A path following the line of latitude around the Earth eastwards. -- -- This is equivalent to @rhumbPath pt (pi/2)@ @@ -256,7 +256,7 @@ where line distance = (pt2, pi/2, 0) where - pt2 = Geodetic + pt2 = Geodetic (latitude pt) (longitude pt + distance / r) 0 (ellipsoid pt) r = latitudeRadius (ellipsoid pt) (latitude pt) diff --git a/src/Geodetics/PolarStereographic.hs b/src/Geodetics/PolarStereographic.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/PolarStereographic.hs @@ -0,0 +1,220 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE FlexibleContexts #-} + +module Geodetics.PolarStereographic ( + Pole (..), + PolarStereographic (trueOrigin, falseOrigin, polarEllipsoid, gridScale), + mkGridPolarStereographic, + UpsGrid, + upsNorth, + upsSouth, + fromUpsGridReference, + parseUpsGridReference, + toUpsGridReference +) where + +import Control.Monad +import Data.Char +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Geodetics.Grid +import Text.Parsec +import Text.Parsec.Error +import Text.Read (readMaybe) +import Text.Printf + +-- | Polar stereographic grids are defined for true origins at the north and south poles. +data Pole = NorthPole | SouthPole deriving (Show, Ord, Eq, Enum, Bounded) + + +{- | Polar Stereographic Grids + +Formulae are taken from +/The Universal Grids: Univerersal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS)/ +DMA Technical Manual 8358.2, Defense Mapping Agency, Fairfax, VA. https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf + +When working with polar grids all directions are relative to the grid rather than the actual pole. +So in the Arctic \"North\" on the Universal Polar Stereographic grid means towards the Bering Sea +rather than towards the North Pole. +-} +data PolarStereographic e = PolarStereographic { + trueOrigin :: Pole, + falseOrigin :: GridOffset, + -- ^ The negation of the grid position of the true origin. Used to avoid negative coordinates over the area + -- of interest. The altitude gives a vertical offset from the ellipsoid. + polarEllipsoid :: e, + -- ^ The ellipsoid for the projection. Arguments passed to `toGrid` *must* use this ellipsoid. + -- The type system cannot verify this for `LocalEllipsoid`. + gridScale :: Double, + -- ^ The scaling factor applied at the pole. This balances the distortion between the center + -- and edges of the projection. + + -- Remaining elements are memoised parameters computed from the ellipsoid. + gridA, gridB, gridC, gridD, gridC0 :: !Double +} deriving (Show) + +instance (Eq e) => Eq (PolarStereographic e) where + g1 == g2 = + trueOrigin g1 == trueOrigin g2 && + falseOrigin g1 == falseOrigin g2 && + polarEllipsoid g1 == polarEllipsoid g2 && + gridScale g1 == gridScale g2 + +instance (Ellipsoid e) => GridClass (PolarStereographic e) e where + fromGrid p = Geodetic lat long (altGP p) (polarEllipsoid gb) + where + gridZero = GridPoint 0 0 0 gb + gb = gridBasis p + p' = gridZero `gridOffset` (falseOrigin gb `applyOffset` p) + radius = offsetDistance p' + isoColat = 2 * atan (radius / (gridScale gb * gridC0 gb)) + isoLat = pi/2 - isoColat + lat1 = isoLat + + gridA gb * sin (2*isoLat) + + gridB gb * sin (4*isoLat) + + gridC gb * sin (6*isoLat) + + gridD gb * sin (8*isoLat) + lat = case trueOrigin gb of + NorthPole -> lat1 + SouthPole -> negate lat1 + long = case trueOrigin gb of + NorthPole -> offsetBearing p' { deltaNorth = negate $ deltaNorth p'} + SouthPole -> offsetBearing p' + + toGrid r geo = offsetNegate (falseOrigin r) `applyOffset` GridPoint east north 0 r + where + absLat = abs $ latitude geo + e = sqrt (eccentricity2 $ polarEllipsoid r) + eSinLat = e * sin absLat + tz2 = ((1 + eSinLat)/(1-eSinLat))**(e/2) * tan (pi/4 - absLat / 2) + radius = gridScale r * gridC0 r * tz2 + north = case trueOrigin r of + NorthPole -> negate $ radius * cos (longitude geo) + SouthPole -> radius * cos (longitude geo) + east = radius * sin (longitude geo) + gridEllipsoid = polarEllipsoid + + +mkGridPolarStereographic :: (Ellipsoid e) => + Pole -- ^ True origin at north or south pole. + -> e -- ^ The ellipsoid used for the projection. + -> GridOffset -- ^ Vector from true origin to the false origin. + -> Double -- ^ Scale factor. + -> PolarStereographic e +mkGridPolarStereographic pole ellip offset scale = + PolarStereographic { + trueOrigin = pole, + falseOrigin = offset, + polarEllipsoid = ellip, + gridScale = scale, + gridA = e2/2 + (5/24)*e4 + e6/12 + (13/360)*e8, + gridB = (7/48)*e4 + (29/240)*e6 + (811/11520)*e8, + gridC = (7/120)*e6 + (81/1120)*e8, + gridD = (4279/161280)*e8, + gridC0 = (2 * majorRadius ellip / sqrt (1 - e2)) * ((1-e1)/(1+e1))**(e1/2) + } + where + e1 = sqrt $ eccentricity2 ellip + e2 = eccentricity2 ellip + e4 = e2^_2 + e6 = e2^_3 + e8 = e2^_4 + + +-- | The Universal Polar Stereographic (UPS) grids for north and south poles. +type UpsGrid = PolarStereographic WGS84 + + +-- | UPS grid for the North Pole. +upsNorth :: UpsGrid +upsNorth = mkGridPolarStereographic + NorthPole + WGS84 + (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 }) + 0.994 -- Scale factor + + +-- | UPS grid for the South Pole. +upsSouth :: UpsGrid +upsSouth = mkGridPolarStereographic + SouthPole + WGS84 + (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 }) + 0.994 -- Scale factor + + +-- | Convert a grid reference into a UPS grid location. +-- +-- There doesn't appear to be any conventional representation for polar grid references, +-- so this is an attempt to cover as many bases as possible. It takes an Easting followed +-- by a Northing with spaces in between. Both can have optional units of m or km, +-- and be optionally followed by an \"N\" or \"E\" as appropriate. +-- +-- The choice of pole is provided in an extra argument rather than within the string because humans +-- will normally assume this from the context and so not provide it. +-- +-- If the string cannot be parsed then one or more error messages are returned. +fromUpsGridReference :: Pole -> String -> Either [String] (GridPoint UpsGrid) +fromUpsGridReference pole str = case parse (parseUpsGridReference pole) "" str of + Left err -> Left $ lines $ showErrorMessages + "or" "unknown parse error" "expecting" "unexpected" "end of input" + (errorMessages err) + Right r -> Right r + + +parseUpsGridReference :: Stream s m Char => Pole -> ParsecT s u m (GridPoint UpsGrid) +parseUpsGridReference pole = do + spaces1 + (eastings1, eastUnit) <- readDistance + spaces + optional (oneOf "Ee" "E") + spaces1 + (northings1, northUnit) <- readDistance + unless (eastUnit == northUnit) $ fail "Northings and Eastings units don't match." + spaces1 + optional (oneOf "Nn" "N") + spaces1 + eof + return $ GridPoint eastings1 northings1 0 $ + case pole of + NorthPole -> upsNorth + SouthPole -> upsSouth + where + readDistance = do -- Returns (distance in meters, unit from input) + digits <- many1 (digit <|> char '.' "number") + spaces1 + when (length digits > 20) $ fail "Too many digits." + (multiplier, unit) <- do + unit <- optionMaybe (string1' "m" <|> string1' "km" "units (m or km)") + case unit of + Just "km" -> pure (1000, Just GridKilometers) + Just _ -> pure (1, Just GridMeters) + Nothing -> pure (1, Nothing) + case readMaybe digits of + Just d -> return (d * multiplier, unit) + Nothing -> fail $ "Cannot read number: " <> digits + string1' target = try $ do -- Case-insensitive version of string' + cs <- count (length target) anyToken + if map toLower target == map toLower cs then return cs else unexpected cs + spaces1 = void $ many (char ' ' "space") -- Other white space not permitted. + + +toUpsGridReference :: + Maybe GridUnit -- ^ Include explicit units in the output. @Nothing@ means meters without units. + -> Bool -- ^ Include \"E\" and \"N\" in the output. + -> Int -- ^ Digits of resolution. 0 = 1m resolution, 1 = 10m, 2 = 100m etc. (-2) = 1cm. + -> GridPoint UpsGrid + -> String +toUpsGridReference unit letters res gp = + dist (eastings gp) <> (if letters then "E " else " ") <> + dist (northings gp) <> (if letters then "N" else "") + where + res1 :: Double + res1 = 10 ** fromIntegral res -- Resolution in meters. + floorRes :: Double -> Double + floorRes d = res1 * fromIntegral (floor (d/res1) :: Integer) + dist d = case unit of + Nothing -> printf "%.*f" (-res) $ floorRes d + Just GridMeters -> printf "%.*fm" (-res) $ floorRes d + Just GridKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000 diff --git a/src/Geodetics/TransverseMercator.hs b/src/Geodetics/TransverseMercator.hs --- a/src/Geodetics/TransverseMercator.hs +++ b/src/Geodetics/TransverseMercator.hs @@ -15,6 +15,9 @@ -- a sheet curved around the ellipsoid so that it touches it at one north-south line (hence making it part of -- a slightly elliptical cylinder). -- +-- Arguments passed to `toGrid` *must* use the same ellipsoid as the `trueOrigin`. The type system +-- cannot verify this for 'LocalEllipsoid'. +-- -- The calculations here are based on \"Transverse Mercator Projection: Constants, Formulae and Methods\" -- by the Ordnance Survey, March 1983. -- Retrieved from http://www.threelittlemaids.co.uk/magdec/transverse_mercator_projection.pdf @@ -22,14 +25,14 @@ trueOrigin :: Geodetic e, -- ^ A point on the line where the projection touches the ellipsoid (altitude is ignored). falseOrigin :: GridOffset, - -- ^ The grid position of the true origin. Used to avoid negative coordinates over + -- ^ The negation of the grid position of the true origin. Used to avoid negative coordinates over -- the area of interest. The altitude gives a vertical offset from the ellipsoid. gridScale :: Double, -- ^ A scaling factor that balances the distortion between the east & west edges and the middle -- of the projection. -- Remaining elements are memoised parameters computed from the ellipsoid underlying the true origin. - gridN1, gridN2, gridN3, gridN4 :: Double + gridN1, gridN2, gridN3, gridN4 :: !Double } deriving (Show) diff --git a/src/Geodetics/UTM.hs b/src/Geodetics/UTM.hs --- a/src/Geodetics/UTM.hs +++ b/src/Geodetics/UTM.hs @@ -41,15 +41,11 @@ mkUtmZone, mkUtmZoneUnsafe, fromUtmGridReference, - toUtmGridReference, - mgrsBandLetterToLatitude, - mgrsLatitudeToBandLetter, - fromMgrsGridReference, - toMgrsGridReference + parseUtmGridReference, + toUtmGridReference ) where -import Control.Monad (mplus, guard, void, when, unless) -import Data.Array +import Control.Monad (guard, void, when, unless) import Data.Char import Data.List import Geodetics.Ellipsoids @@ -155,10 +151,6 @@ scale = 0.999_6 --- | Units for UTM grid coordinates. -data UtmGridUnit = UtmMeters | UtmKilometers deriving (Eq, Show) - - -- | Convert a grid reference to a position, if the reference is valid. -- -- The northings and eastings cannot contain more than 20 digits each, @@ -172,13 +164,15 @@ -- -- If the argument cannot be parsed then one or more error messages are returned. fromUtmGridReference :: String -> Either [String] (GridPoint UtmZone) -fromUtmGridReference str = case parse gridP str str of +fromUtmGridReference str = case parse parseUtmGridReference str str of Left err -> Left $ lines $ showErrorMessages "or" "unknown parse error" "expecting" "unexpected" "end of input" (errorMessages err) Right r -> Right r - where - gridP = do + + +parseUtmGridReference :: Stream s m Char => ParsecT s u m (GridPoint UtmZone) +parseUtmGridReference = do spaces1 zone <- readZone "Zone number" hemi <- readHemi "Hemisphere (N or S)" @@ -193,8 +187,9 @@ optional (oneOf "Nn" "N") spaces1 eof - return $ GridPoint eastings1 northings1 0 $ mkUtmZoneUnsafe hemi zone - readZone :: Parsec String () UtmZoneNumber + pure $ GridPoint eastings1 northings1 0 $ mkUtmZoneUnsafe hemi zone + where + readZone :: Stream s m Char => ParsecT s u m UtmZoneNumber readZone = do ds <- many1 digit case readMaybe ds of @@ -202,28 +197,32 @@ Just n -> if n < 1 || n > 60 then fail $ "Zone number " <> show n <> " out of range." - else return n - readHemi :: Parsec String () UtmHemisphere + else pure n + readHemi :: Stream s m Char => ParsecT s u m UtmHemisphere readHemi = do h <- oneOf "NSns" case toUpper h of - 'N' -> return UtmNorth - 'S' -> return UtmSouth + 'N' -> pure UtmNorth + 'S' -> pure UtmSouth _ -> fail $ "Invalid hemisphere: " <> (h : ". Must be N or S.") - readDistance :: Parsec String () (Double, UtmGridUnit) -- (Distance, unit) + readDistance :: Stream s m Char => ParsecT s u m (Double, Maybe GridUnit) -- (Distance, unit) readDistance = do digits <- many1 (digit <|> char '.' "number") spaces1 when (length digits > 20) $ fail "Too many digits." (multiplier, unit) <- do - unit <- option "m" (string1' "m" <|> string1' "km" "units (m or km)") - if unit == "km" then return (1000, UtmKilometers) else return (1, UtmMeters) + unit <- optionMaybe (string1' "m" <|> string1' "km" "units (m or km)") + case unit of + Just "km" -> pure (1000, Just GridKilometers) + Just _ -> pure (1, Just GridMeters) + Nothing -> pure (1, Nothing) case readMaybe digits of - Just d -> return (d * multiplier, unit) + Just d -> pure (d * multiplier, unit) Nothing -> fail $ "Cannot read number: " <> digits string1' target = try $ do -- Case-insensitive version of string' cs <- count (length target) anyToken - if map toLower target == map toLower cs then return cs else unexpected cs + if map toLower target == map toLower cs then pure cs else unexpected cs + spaces1 :: Stream s m Char => ParsecT s u m () spaces1 = void $ many (char ' ' "space") -- Other white space not permitted. @@ -231,7 +230,7 @@ -- The northings and eastings are rounded down to the resolution, so the result is the south-west -- corner of the grid square enclosing the grid point. toUtmGridReference :: - Maybe UtmGridUnit -- ^ Include explicit units in the output. @Nothing@ means meters without units. + Maybe GridUnit -- ^ Include explicit units in the output. @Nothing@ means meters without units. -> Bool -- ^ Include \"E\" and \"N\" in the output. -> Int -- ^ Digits of resolution. 0 = 1m resolution, 1 = 10m, 2 = 100m etc. (-2) = 1cm. -> GridPoint UtmZone @@ -249,216 +248,5 @@ zoneStr = printf "%02d" (utmZoneNum b) <> show (utmHemisphere b) dist d = case unit of Nothing -> printf "%.*f" (-res) $ floorRes d - Just UtmMeters -> printf "%.*fm" (-res) $ floorRes d - Just UtmKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000 - - - --- | The MGRS latitude band code letters, excluding A and B used for Antarctica (south of -80 degrees) --- and Y and Z used for the Arctic (north of 84 degrees). -mgrsBandLetters :: [Char] -mgrsBandLetters = "CDEFGHJKLMNPQRSTUVWX" - - --- | Find the southern boundary of a latitude band letter. -mgrsBandLetterToLatitude :: Char -> Maybe Double -mgrsBandLetterToLatitude band = do - n1 <- ix band - return $ degree * fromIntegral (-80 + n1 * 8) - where - indexMap :: Array Char (Maybe Int) - indexMap = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip mgrsBandLetters [0..]] - ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing - - --- | Find the band letter for a latitude, if it is in the range (-80, 84) degrees. --- (Argument in radians) -mgrsLatitudeToBandLetter :: Double -> Maybe Char -mgrsLatitudeToBandLetter lat = do - guard $ -80 <= dlat && dlat <= 84 - return $ indexMap ! latIdx - where - dlat = lat / degree - ilat :: Int - ilat = floor dlat - latIdx = min 19 $ (ilat + 80) `div` 8 -- Band 19 (X) extends an extra 4 degrees. - indexMap = listArray (0,19) mgrsBandLetters - - - --- | Letters A-Z except for I and O. -mgrsEastingsLetters :: [Char] -mgrsEastingsLetters = "ABCDEFGHJKLMNPQRSTUVWXYZ" - --- | If zone number is in range and the letter is one of the valid Eastings letters for that zone --- then return the UTM easting in meters. --- --- Zone 1 starts with \'A\'. Each zone is 8 characters wide. Hence the letters repeat every 3 zones. -mgrsLetterToEasting :: UtmZoneNumber -> Char -> Maybe Double -mgrsLetterToEasting zn c = do - guard $ 1 <= zn && zn <= 60 - n1 <- ix c - let n2 = n1 - base - guard $ 0 <= n2 && n2 <= 7 - return $ fromIntegral (n2+1) * 100_000 - where - indexMap = accumArray mplus Nothing ('A', 'Z') [(c1, Just n) | (c1, n) <- zip mgrsEastingsLetters [0..]] - base = ((zn-1) `mod` 3) * 8 - ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing - - --- | If the zone number is in range and the eastings are between 100,000 and 900,000 then --- return the Eastings letter. -mgrsEastingToLetter :: UtmZoneNumber -> Double -> Maybe Char -mgrsEastingToLetter zn east = do - guard $ 1 <= zn && zn <= 60 - guard $ 100 * kilometer <= east && east < 900 * kilometer - return $ indexMap ! ix - where - indexMap = listArray (0,23) mgrsEastingsLetters - base = ((zn-1) `mod` 3) * 8 - square = max 0 $ min 7 $ floor $ (east - 100 * kilometer)/(100 * kilometer) -- Clamped in range (0,7). - ix = base + square -- Must be in range (0,23) - - --- | Letters A-V except for I and O. -mgrsNorthingsLetters :: [Char] -mgrsNorthingsLetters = "ABCDEFGHJKLMNPQRSTUV" - - --- | MGRS Northings letters have rather complex relationship to the latitude bands. The 20 letters --- repeat every 2,000km going north and south from the equator, so the latitude band is needed --- to disambiguate which repetition of the Northings letter is meant. - --- Unfortunately this repetition of the letters does not neatly coincide with --- the latitude band boundaries, which are based on degrees of latitude. - --- The base letter just north of the equator is A in odd-numbered zones and F in even numbered zones. --- --- This uses the latitude band to estimate the range of northings that would be valid there, --- and hence determine which possible grid band is meant by the northings letter. --- The algorithm used is approximate and deliberately very forgiving: it will accept some grid squares --- which are north or south of the band given. -mgrsLetterToNorthings :: - UtmZoneNumber - -> Char -- ^ Latitude band letter (@C@ - @X@ excluding @I@ and @O@). - -> Char -- ^ MGRS Northings letter (@A@ - @V@ excluding @I@ and @O@). - -> Maybe Double -mgrsLetterToNorthings zone bandC northingsC = do - guard $ 1 <= zone && zone <= 60 - band <- (/degree) <$> mgrsBandLetterToLatitude bandC - northings0 <- (baseNorthingsOffset +) <$> ix northingsC - let bandDist = band * metersPerDegree -- Approx dist from equator to southern edge of band. - bandGridLower = floor $ bandDist / 100_000 - 2 -- Lower limit of band in 100km units - bandGridUpper = ceiling $ if band > 71 * degree -- Upper limit of band in 100km units. - then (bandDist + 12 * metersPerDegree) / 100_000 + 1 -- Band X. - else (bandDist + 8 * metersPerDegree) / 100_000 + 2 -- Other bands. - rep = (bandGridLower - northings0 - 1) `div` 20 -- Lower limit in 2,000,000km units. - grid = (rep+1)*20 + northings0 - guard $ grid >= bandGridLower - guard $ grid <= bandGridUpper - return $ fromIntegral $ grid * 100_000 - where - metersPerDegree = 10_002_000 / 90 -- Equator to north pole. - baseNorthingsOffset :: Int - baseNorthingsOffset = if odd zone then 0 else -5 - indexMap :: Array Char (Maybe Int) - indexMap = accumArray mplus Nothing ('A', 'Z') [(c, Just n) | (c, n) <- zip mgrsEastingsLetters [0..]] - ix c1 = if inRange (bounds indexMap) c1 then indexMap ! c1 else Nothing - - --- | Find the northings letter of the 100km square containing the given Northings. --- --- The input is not range checked. It just assumes that the northings letters repeat forever. -mgrsNorthingToLetter :: UtmZoneNumber -> Double -> Char -mgrsNorthingToLetter zone northings1 = - letters ! ((gridNum + baseNorthingsOffset) `mod` 20) - where - gridNum :: Int - gridNum = floor $ northings1 / (100 * kilometer) - baseNorthingsOffset = if odd zone then 0 else 5 - letters = listArray (0,19) mgrsNorthingsLetters - - --- | Convert an MGRS grid reference to a UTM @GridPoint@, if the reference is valid. --- E.g. \"30U XC 99304 10208\" is the grid reference for Nelson's Column in London. --- --- If the input contains spaces then these are used to delimit the fields. Any or all spaces --- may be omitted. Multiple spaces are treated as a single space. --- --- If the reference is valid this returns the position of the south-west corner of the --- nominated grid square and an offset to its centre. Altitude is set to zero. -fromMgrsGridReference :: String -> Either [String] (GridPoint UtmZone, GridOffset) -fromMgrsGridReference str = case parse mgrsP str str of - Left err -> Left $ filter (not . null) $ lines $ showErrorMessages - "or" "unknown parse error" "expecting" "unexpected" "end of input" - (errorMessages err) - Right r -> Right r - where - mgrsP = do - zoneNum <- read <$> many1 digit -- Safe because we can only read digits here. - spaces - band <- upper "latitude band letter" - let (hemi, falseNorthing) = if band >= 'N' - then (UtmNorth, 0) - else (UtmSouth, -10_000_000) - zone <- maybe (fail "Invalid zone") return $ mkUtmZone hemi zoneNum - spaces - squareEast <- upper "eastings letter" - eastingBase <- maybe (fail "Invalid eastings letter") return $ - mgrsLetterToEasting zoneNum squareEast - spaces - squareNorth <- upper "northings letter" - northingBase <- maybe (fail "Invalid northings letter") return $ - mgrsLetterToNorthings zoneNum band squareNorth - spaces - (eastingChars, northingChars) <- try spaced <|> try unspaced <|> noDigits - when (length eastingChars /= length northingChars) $ - fail "Northings and Eastings must be the same length." - if northingChars == "" && eastingChars == "" - then -- No digits, just return the outer 100km grid square - return (GridPoint eastingBase (northingBase - falseNorthing) 0 zone, - GridOffset (100 * kilometer / 2) (100 * kilometer / 2) 0) - else do - (northing, offset) <- maybe (fail "Invalid northing digits") return $ - fromGridDigits (100 * kilometer) northingChars - (easting, _) <- maybe (fail "Invalid easting digits") return $ - fromGridDigits (100 * kilometer) eastingChars - return (GridPoint (eastingBase + easting) (northingBase + northing - falseNorthing) 0 zone, - GridOffset (offset/2) (offset/2) 0) - spaced = do - e <- many1 digit - skipMany1 space -- A space is mandatory here. - n <- many1 digit - return (e,n) - unspaced = do - digits <- many1 digit - let c = length digits - when (odd c) $ fail "Northings and Eastings must be the same length." - return (splitAt (c `div` 2) digits) - noDigits = do - eof - return ("", "") - --- | Convert UTM @GridPoint@ to an MGRS grid reference. -toMgrsGridReference :: - Bool -- ^ Include spaces in the output. The standard says no spaces, but they make - -- the output easier to read. - -> Int -- ^ Number of digits of precision in the easting and northing. Must be 0-5. - -> GridPoint UtmZone - -> Maybe String -toMgrsGridReference withSpaces precision gp = do - guard $ precision >= 0 && precision <= 5 - band <- mgrsLatitudeToBandLetter $ latitude $ fromGrid gp - let - zoneNum = utmZoneNum $ gridBasis gp - northLetter = mgrsNorthingToLetter zoneNum $ northings gp - eastLetter <- mgrsEastingToLetter zoneNum $ eastings gp - (_, northDigits) <- toGridDigits (100 * kilometer) precision $ northings gp - (_, eastDigits) <- toGridDigits (100 * kilometer) precision $ eastings gp - let part1 = printf "%02d" zoneNum <> [band] - part2 = [eastLetter, northLetter] - return $ if withSpaces - then part1 <> " " <> part2 <> " " <> eastDigits <> " " <> northDigits - else part1 <> part2 <> eastDigits <> northDigits - + Just GridMeters -> printf "%.*fm" (-res) $ floorRes d + Just GridKilometers -> printf "%.*fkm" (3-res) $ floorRes d / 1000 diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -1,6 +1,7 @@ {-# OPTIONS_GHC -fno-warn-orphans #-} {-# OPTIONS_GHC -Wno-unrecognised-pragmas #-} {-# HLINT ignore "Redundant bracket" #-} +{-# LANGUAGE ViewPatterns #-} module Main where @@ -20,12 +21,14 @@ import Geodetics.Ellipsoids import Geodetics.Geodetic import Geodetics.Grid +import Geodetics.MGRS import Geodetics.Path import Geodetics.Stereographic import Geodetics.TransverseMercator import Geodetics.UK import Geodetics.UTM import LatLongParser (parserTests) +import Geodetics.PolarStereographic as PS main :: IO () main = hspec $ do @@ -58,6 +61,12 @@ prop "MGRS Grid 2" prop_mgrs_gridTest2 describe "MGRS Grid 3" $ mapM_ mgrsGridTest3 utmSampleGrid describe "MGRS Grid 4" $ mapM_ mgrsGridTest4 utmSampleGrid + describe "UPS" $ do + describe "UPS Grid 4" $ mapM_ upsGridTest4 upsSampleGrid + describe "UPS Grid 5" $ mapM_ upsGridTest5 upsSampleGrid + describe "UPS Grid 6" $ mapM_ upsGridTest6 upsSampleGrid + describe "UPS Grid 7" $ mapM_ upsGridTest7 upsSampleGrid + describe "UPS Grid 8" $ mapM_ upsGridTest8 upsSampleGrid describe "Stereographic" $ do it "toGrid north" stereographicToGridN it "fromGrid north" stereographicFromGridN @@ -369,7 +378,7 @@ prop_mgrs_gridTest1 (UtmGridRef str) = case fromUtmGridReference str of Left msg -> assertFailure $ "gridTest1 bogus UTM ref: " <> str <> ". Messages = " <> show msg - Right gp -> + Right (utmToMgrsPoint -> gp) -> fromMgrsGridReference <$> toMgrsGridReference True 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0)) -- | Check that a UTM grid point round-trips to MGRS and back, without spaces. @@ -377,19 +386,73 @@ prop_mgrs_gridTest2 (UtmGridRef str) = case fromUtmGridReference str of Left msg -> assertFailure $ "gridTest2 bogus UTM ref: " <> str <> ". Messages = " <> show msg - Right gp -> + Right (utmToMgrsPoint -> gp) -> fromMgrsGridReference <$> toMgrsGridReference False 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0)) -- | Check that MGRS reference to grid point works for sample points. mgrsGridTest3 :: UtmGridPointTest -mgrsGridTest3 (_, mgrs, gp, _, testName) = +mgrsGridTest3 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) = it testName $ fromMgrsGridReference mgrs `shouldBe` Right (gp, GridOffset 0.5 0.5 0) -- | Check that grid point to MGRS reference works for sample points. mgrsGridTest4 :: UtmGridPointTest -mgrsGridTest4 (_, mgrs, gp, _, testName) = do +mgrsGridTest4 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) = do it (testName <> " with spaces") $ toMgrsGridReference True 5 gp `shouldBe` Just mgrs it (testName <> " without spaces") $ toMgrsGridReference False 5 gp `shouldBe` Just (filter (not . isSpace) mgrs) + + +-- | DMA TM 8358.2 can be found at https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf +-- Oracle for other points is https://geographiclib.sourceforge.io/C++/doc/GeoConvert.1.html +upsSampleGrid :: [(String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String)] +upsSampleGrid = map convert [ + -- MGRS UPS, X, Y, Latitude, Longitude + ("ZAH0000000000", "2000000 2000000", 2000000.00, 2000000.00, 90, 0, "North pole"), + ("BAN0000000000", "2000000 2000000", 2000000.00, 2000000.00, -90, 0, "South pole"), + ("ZAE0000077930", "2000000 1777930", 2000000.00, 1777930.73, 88, 0, "88N 0W"), + ("YXH7793000000", "1777930 2000000", 1777930.73, 2000000.00, 88, -90, "88N, 90W"), + ("BAQ0000022069", "2000000 2222069", 2000000.00, 2222069.27, -88, 0, "88S 0W"), + ("AXN7793000000", "1777930 2000000", 1777930.73, 2000000.00, -88, -90, "88S, 90W"), + ("YTM3012526773", "1530125 2426773", 1530125.78, 2426773.60, 84.28723389, -132.2479892, "DMA TM 8358.2 sample 1"), + ("BCK2297997474", "2222979 1797474", 2222979.47, 1797474.90, -87.28733333, 132.2478619, "DMA TM 8358.2 sample 3") ] + -- DMA TM 8358.2 sample 2 is outside the UPS grid regions and is therefore not a valid test. + where + convert (mgrs, ups, x, y, lat, long, desc) = (mgrs, ups, GridPoint { eastings = x, northings = y, altGP = 0, gridBasis = gb }, geo, desc) + where + geo = Geodetic { latitude = (lat * degree), longitude = (long * degree), geoAlt = 0, ellipsoid = WGS84 } + gb = mkGridPolarStereographic + (if lat < 0 then SouthPole else NorthPole) + WGS84 + (GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 }) + 0.994 + +type UpsGridPointTest = (String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String) -> SpecWith (Arg Expectation) + +-- | Check that the UPS grid point to WGS84 works close enough for sample points. +upsGridTest4 :: UpsGridPointTest +upsGridTest4 (_, _, gp, geo, testName) = + it testName $ geo `closeEnough` fromGrid gp + +-- | Check that WGS84 to UPS grid point works close enough for sample points. +upsGridTest5 :: UpsGridPointTest +upsGridTest5 (_, _, gp, geo, testName) = + it testName $ gp `closeGrid` toGrid (gridBasis gp) geo + +-- | Check that UPS and MGRS parsers agree. +upsGridTest6 :: UpsGridPointTest +upsGridTest6 (mgrs, ups, gp, _, testName) = + it testName $ + fst <$> fromMgrsGridReference mgrs `shouldBe` + upsToMgrsPoint <$> fromUpsGridReference (PS.trueOrigin $ gridBasis gp) ups + +-- | Check MGRS grid generator +upsGridTest7 :: UpsGridPointTest +upsGridTest7 (mgrs, _, gp, _, testName) = + it testName $ toMgrsGridReference False 5 (upsToMgrsPoint gp) `shouldBe` Just mgrs + +-- | Check UPS grid generator +upsGridTest8 :: UpsGridPointTest +upsGridTest8 (_, ups, gp, _, testName) = + it testName $ toUpsGridReference Nothing False 0 gp `shouldBe` ups -- | Standard stereographic grid for point tests in the Northern Hemisphere.