{-# 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