diff --git a/ChangeLog.md b/ChangeLog.md --- a/ChangeLog.md +++ b/ChangeLog.md @@ -1,3 +1,9 @@ -### 0.1.0.0 - -- Initial version +### 0.2.0.0 + +- GeoPos -> LatLong +- Split Position from GreatCircle +- require base >= 4.9 + +### 0.1.0.0 + +- Initial version diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,15 +1,14 @@ -# Jord +# Jord - Geographical Position Calculations [![travis build status](https://img.shields.io/travis/ofmooseandmen/jord/master.svg?label=travis+build)](https://travis-ci.org/ofmooseandmen/jord) +[![Hackage](https://img.shields.io/hackage/v/jord.svg)](http://hackage.haskell.org/package/jord) [![license](https://img.shields.io/badge/license-BSD3-lightgray.svg)](https://opensource.org/licenses/BSD-3-Clause) > __Jord__ [_Swedish_] is __Earth__ [_English_] -Geographic position calculations on great circles. - ## What is this? -Jord is a [Haskell](https://www.haskell.org) library that implements various geographical position calculations on great circles using the algorithms described in [Gade, K. (2010). A Non-singular Horizontal Position Representation](http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf). +Jord is a [Haskell](https://www.haskell.org) library that implements various geographical position calculations using the algorithms described in [Gade, K. (2010). A Non-singular Horizontal Position Representation](http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf). ## How do I build it? @@ -20,6 +19,8 @@ ``` ## How do I use it? + +[See documentation on Hackage](http://hackage.haskell.org/package/jord/docs/Data-Geo-Jord.html) ```haskell import Data.Geo.Jord diff --git a/app/Main.hs b/app/Main.hs --- a/app/Main.hs +++ b/app/Main.hs @@ -8,162 +8,162 @@ -- -- REPL around "Jord". -- -module Main where - -import Data.Geo.Jord -import Data.List ((\\), dropWhileEnd, intercalate, isPrefixOf) -import Prelude hiding (lookup) -import System.Console.Haskeline - -search :: String -> [Completion] -search s = map simpleCompletion $ filterFunc s - -filterFunc :: String -> [String] -filterFunc s = map (\f -> pref ++ f) filtered - where +module Main where + +import Data.Geo.Jord +import Data.List ((\\), dropWhileEnd, intercalate, isPrefixOf) +import Prelude hiding (lookup) +import System.Console.Haskeline + +search :: String -> [Completion] +search s = map simpleCompletion $ filterFunc s + +filterFunc :: String -> [String] +filterFunc s = map (\f -> pref ++ f) filtered + where pref = dropWhileEnd (/= '(') s -- everything before the last '(' inclusive func = (\\) s pref -- everything after the last '(' - filtered = filter (\f -> func `isPrefixOf` f) functions - -mySettings :: Settings IO -mySettings = - Settings - { complete = completeWord Nothing " \t" $ return . search - , historyFile = Nothing - , autoAddHistory = True - } - -main :: IO () -main = do - putStrLn - ("jord interpreter, version " ++ - jordVersion ++ ": https://github.com/ofmooseandmen/jord :? for help") - runInputT mySettings $ withInterrupt $ loop emptyVault - where - loop state = do - input <- handleInterrupt (return (Just "")) $ getInputLine "jord> " - case input of - Nothing -> return () - Just ":quit" -> return () - Just ":q" -> return () - Just i -> do - let (result, newState) = evalS i state - printS result - loop newState - -printS :: Either String String -> InputT IO () -printS (Left err) = outputStrLn ("jord> " ++ err) -printS (Right "") = return () -printS (Right r) = outputStrLn ("jord> " ++ r) - -evalS :: String -> Vault -> (Either String String, Vault) -evalS s vault - | null s = (Right "", vault) - | head s == ':' = evalC w vault - | (v:"=":e) <- w = - let r = eval (unwords e) vault - vault' = save r v vault - in (showR r, vault') - | otherwise = (showR (eval s vault), vault) - where - w = words s - -evalC :: [String] -> Vault -> (Either String String, Vault) -evalC [":show", v] vault = (evalShow v vault, vault) -evalC [":delete", v] vault = evalDel (Just v) vault -evalC [":clear"] vault = evalDel Nothing vault -evalC [":help"] vault = (Right help, vault) -evalC [":?"] vault = (Right help, vault) -evalC c vault = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), vault) - -evalShow :: String -> Vault -> Either String String -evalShow n vault = maybe (Left ("Unbound variable: " ++ n)) (Right . showVar n) (lookup n vault) - -evalDel :: Maybe String -> Vault -> (Either String String, Vault) -evalDel (Just n) vault = (Right ("deleted var: " ++ n), delete n vault) -evalDel Nothing _ = (Right "deleted all variable ", emptyVault) - -help :: String -help = - "\njord interpreter, version " ++ - jordVersion ++ - "\n" ++ - "\n Commands available from the prompt:\n\n" ++ - " :help, :? display this list of commands\n" ++ - " :quit, :q quit jord\n" ++ - " :show {var} shows {var}\n" ++ - " :delete {var} deletes {var}\n" ++ - " :clear deletes all variable(s)\n" ++ - "\n Jord expressions:\n\n" ++ - " (f x y) where f is one of function described below and x and y\n" ++ - " are either parameters in one of the format described below or\n" ++ - " a call to another function\n" ++ - "\n" ++ - " (finalBearing (destination (antipode 54°N,154°E) 54° 1000m) (readGeoPos 54°N,154°E))\n" ++ - "\n" ++ - " Top level () can be ommitted: antipode 54N028E\n" ++ - "\n Position calculations:\n\n" ++ - " antipode pos antipodal point of pos\n" ++ - " crossTrackDistance pos gc signed distance from pos to great circle gc\n" ++ - " distance pos1 pos2 surface distance between pos1 and pos2\n" ++ - " destination pos len ang destination position from pos having travelled len\n" ++ - " on initial bearing ang\n" ++ - " finalBearing pos1 pos2 initial bearing from pos1 to pos2\n" ++ - " initialBearing pos1 pos2 bearing arriving at pos2 from pos1\n" ++ - " interpolate pos1 pos2 (0..1) position at fraction between pos1 and pos2\n" ++ - " intersections gc1 gc2 intersections between great circles gc1 and gc2\n" ++ - " exactly 0 or 2 intersections\n" ++ - " isInside pos [pos] is p inside polygon?\n" ++ - " mean [pos] geographical mean position of [pos]\n" ++ - "\n Constructors and conversions:\n\n" ++ - " decimalDegrees double angle from decimal degrees\n" ++ - " greatCircle pos1 pos2 great circle passing by pos1 and pos2\n" ++ - " greatCircle pos ang great circle passing by pos and heading on bearing ang\n" ++ - " latLong ang ang geographic position from latitude & longitude\n" ++ - " latLongDecimal double double geographic position from latitude & longitude (DD)\n" ++ - " readGeoPos string geographic position from string\n" ++ - " toDecimalDegrees pos latitude and longitude of pos in decimal degrees\n" ++ - " toDecimalDegrees ang decimal degrees of ang\n" ++ - " toKilometres len length to kilometres\n" ++ - " toMetres len length to metres\n" ++ - " toNauticalMiles len length to nautical miles\n" ++ - " toNVector pos n-vector corresponding to pos\n" ++ - "\n Supported Position formats:\n\n" ++ - " DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - 553621N0130209E\n" ++ - " d°m's\"[N|S],d°m's\"[E|W] - 55°36'21\"N,13°2'9\"E\n" ++ - " ^ zeroes can be ommitted and separtors can also be d, m, s\n" ++ - " decimal°[N|S],decimal°[E|W] - 51.885°N,13,1°E\n" ++ - "\n Supported Angle formats:\n\n" ++ - " d°m's - 55°36'21.154\n" ++ - " decimal° - 51.885°\n" ++ - "\n Supported Length formats: {l}m, {l}km, {l}Nm\n\n" ++ - "\n Every evaluated result can be saved by prefixing the expression with \"{var} = \"\n" ++ - " Saved results can subsequently be used when calling a function\n" ++ - " jord> a = antipode 54N028E\n" ++ " jord> antipode a\n" - -save :: Result -> String -> Vault -> Vault -save (Right v) k vault = insert k v vault -save _ _ vault = vault - -showR :: Result -> Either String String -showR (Left err) = Left err -showR (Right v) = Right (showV v) - -showV :: Value -> String -showV (Ang a) = "angle: " ++ show a -showV (AngDec a) = "angle (dd): " ++ show a -showV (Bool b) = show b -showV (Double d) = show d -showV (Len l) = "length: " ++ show l -showV (Geo g) = "geographic position: " ++ show g -showV (Geos gs) = "geographic position: " ++ intercalate "; " (map show gs) -showV (GeoDec ll) = "latitude, longitude (dd): " ++ show (fst ll) ++ ", " ++ show (snd ll) -showV (GeosDec lls) = - "latitudes, longitudes (dd): " ++ - intercalate "; " (map (\ll -> show (fst ll) ++ ", " ++ show (snd ll)) lls) -showV (Vec v) = "n-vector: " ++ show v -showV (Vecs vs) = "n-vectors: " ++ intercalate "; " (map show vs) -showV (Gc gc) = "great circle: " ++ show gc - -showVar :: String -> Value -> String -showVar n v = n ++ "=" ++ showV v + filtered = filter (\f -> func `isPrefixOf` f) functions + +mySettings :: Settings IO +mySettings = + Settings + { complete = completeWord Nothing " \t" $ return . search + , historyFile = Nothing + , autoAddHistory = True + } + +main :: IO () +main = do + putStrLn + ("jord interpreter, version " ++ + jordVersion ++ ": https://github.com/ofmooseandmen/jord :? for help") + runInputT mySettings $ withInterrupt $ loop emptyVault + where + loop state = do + input <- handleInterrupt (return (Just "")) $ getInputLine "jord> " + case input of + Nothing -> return () + Just ":quit" -> return () + Just ":q" -> return () + Just i -> do + let (result, newState) = evalS i state + printS result + loop newState + +printS :: Either String String -> InputT IO () +printS (Left err) = outputStrLn ("jord> " ++ err) +printS (Right "") = return () +printS (Right r) = outputStrLn ("jord> " ++ r) + +evalS :: String -> Vault -> (Either String String, Vault) +evalS s vault + | null s = (Right "", vault) + | head s == ':' = evalC w vault + | (v:"=":e) <- w = + let r = eval (unwords e) vault + vault' = save r v vault + in (showR r, vault') + | otherwise = (showR (eval s vault), vault) + where + w = words s + +evalC :: [String] -> Vault -> (Either String String, Vault) +evalC [":show", v] vault = (evalShow v vault, vault) +evalC [":delete", v] vault = evalDel (Just v) vault +evalC [":clear"] vault = evalDel Nothing vault +evalC [":help"] vault = (Right help, vault) +evalC [":?"] vault = (Right help, vault) +evalC c vault = (Left ("Unsupported command " ++ unwords c ++ "; :? for help"), vault) + +evalShow :: String -> Vault -> Either String String +evalShow n vault = maybe (Left ("Unbound variable: " ++ n)) (Right . showVar n) (lookup n vault) + +evalDel :: Maybe String -> Vault -> (Either String String, Vault) +evalDel (Just n) vault = (Right ("deleted var: " ++ n), delete n vault) +evalDel Nothing _ = (Right "deleted all variable ", emptyVault) + +help :: String +help = + "\njord interpreter, version " ++ + jordVersion ++ + "\n" ++ + "\n Commands available from the prompt:\n\n" ++ + " :help, :? display this list of commands\n" ++ + " :quit, :q quit jord\n" ++ + " :show {var} shows {var}\n" ++ + " :delete {var} deletes {var}\n" ++ + " :clear deletes all variable(s)\n" ++ + "\n Jord expressions:\n\n" ++ + " (f x y) where f is one of function described below and x and y\n" ++ + " are either parameters in one of the format described below or\n" ++ + " a call to another function\n" ++ + "\n" ++ + " (finalBearing (destination (antipode 54°N,154°E) 54° 1000m) (readGeoPos 54°N,154°E))\n" ++ + "\n" ++ + " Top level () can be ommitted: antipode 54N028E\n" ++ + "\n Position calculations:\n\n" ++ + " antipode pos antipodal point of pos\n" ++ + " crossTrackDistance pos gc signed distance from pos to great circle gc\n" ++ + " distance pos1 pos2 surface distance between pos1 and pos2\n" ++ + " destination pos len ang destination position from pos having travelled len\n" ++ + " on initial bearing ang\n" ++ + " finalBearing pos1 pos2 initial bearing from pos1 to pos2\n" ++ + " initialBearing pos1 pos2 bearing arriving at pos2 from pos1\n" ++ + " interpolate pos1 pos2 (0..1) position at fraction between pos1 and pos2\n" ++ + " intersections gc1 gc2 intersections between great circles gc1 and gc2\n" ++ + " exactly 0 or 2 intersections\n" ++ + " isInside pos [pos] is p inside polygon?\n" ++ + " mean [pos] geographical mean position of [pos]\n" ++ + "\n Constructors and conversions:\n\n" ++ + " decimalDegrees double angle from decimal degrees\n" ++ + " greatCircle pos1 pos2 great circle passing by pos1 and pos2\n" ++ + " greatCircle pos ang great circle passing by pos and heading on bearing ang\n" ++ + " latLong ang ang geographic position from latitude & longitude\n" ++ + " latLongDecimal double double geographic position from latitude & longitude (DD)\n" ++ + " readLatLong string geographic position from string\n" ++ + " toDecimalDegrees pos latitude and longitude of pos in decimal degrees\n" ++ + " toDecimalDegrees ang decimal degrees of ang\n" ++ + " toKilometres len length to kilometres\n" ++ + " toMetres len length to metres\n" ++ + " toNauticalMiles len length to nautical miles\n" ++ + " toNVector pos n-vector corresponding to pos\n" ++ + "\n Supported Lat/Long formats:\n\n" ++ + " DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - 553621N0130209E\n" ++ + " d°m's\"[N|S],d°m's\"[E|W] - 55°36'21\"N,13°2'9\"E\n" ++ + " ^ zeroes can be ommitted and separtors can also be d, m, s\n" ++ + " decimal°[N|S],decimal°[E|W] - 51.885°N,13,1°E\n" ++ + "\n Supported Angle formats:\n\n" ++ + " d°m's - 55°36'21.154\n" ++ + " decimal° - 51.885°\n" ++ + "\n Supported Length formats: {l}m, {l}km, {l}Nm\n\n" ++ + "\n Every evaluated result can be saved by prefixing the expression with \"{var} = \"\n" ++ + " Saved results can subsequently be used when calling a function\n" ++ + " jord> a = antipode 54N028E\n" ++ " jord> antipode a\n" + +save :: Result -> String -> Vault -> Vault +save (Right v) k vault = insert k v vault +save _ _ vault = vault + +showR :: Result -> Either String String +showR (Left err) = Left err +showR (Right v) = Right (showV v) + +showV :: Value -> String +showV (Ang a) = "angle: " ++ show a +showV (AngDec a) = "angle (dd): " ++ show a +showV (Bool b) = show b +showV (Double d) = show d +showV (Len l) = "length: " ++ show l +showV (Ll g) = "geographic position: " ++ show g +showV (Lls gs) = "geographic position: " ++ intercalate "; " (map show gs) +showV (LlDec ll) = "latitude, longitude (dd): " ++ show (fst ll) ++ ", " ++ show (snd ll) +showV (LlsDec lls) = + "latitudes, longitudes (dd): " ++ + intercalate "; " (map (\ll -> show (fst ll) ++ ", " ++ show (snd ll)) lls) +showV (Vec v) = "n-vector: " ++ show v +showV (Vecs vs) = "n-vectors: " ++ intercalate "; " (map show vs) +showV (Gc gc) = "great circle: " ++ show gc + +showVar :: String -> Value -> String +showVar n v = n ++ "=" ++ showV v diff --git a/jord.cabal b/jord.cabal --- a/jord.cabal +++ b/jord.cabal @@ -2,11 +2,11 @@ -- -- see: https://github.com/sol/hpack -- --- hash: c02feb8438ae6672aeef69101f3c30532d6cf7132fece36f1fb9b67b207c9332 +-- hash: c1197ea6dc351c977d808ecb2b0eef21c2993f1a1c14e40471ff7a0b18fcd501 name: jord -version: 0.1.0.0 -synopsis: Geographic position calculations on Great Circles +version: 0.2.0.0 +synopsis: Geographical Position Calculations description: Please see the README on GitHub at category: Geography stability: experimental @@ -33,8 +33,9 @@ Data.Geo.Jord.Angle Data.Geo.Jord.Eval Data.Geo.Jord.Length - Data.Geo.Jord.GeoPos + Data.Geo.Jord.LatLong Data.Geo.Jord.GreatCircle + Data.Geo.Jord.Position Data.Geo.Jord.Quantity Data.Geo.Jord.NVector other-modules: @@ -43,7 +44,7 @@ src ghc-options: -Wall build-depends: - base >=4.7 && <5 + base >=4.9 && <5 default-language: Haskell2010 executable jord-exe @@ -54,7 +55,7 @@ app ghc-options: -Wall build-depends: - base >=4.7 && <5 + base >=4.9 && <5 , haskeline >=0.7 && <0.8 , jord default-language: Haskell2010 @@ -65,16 +66,17 @@ other-modules: Data.Geo.Jord.AngleSpec Data.Geo.Jord.EvalSpec - Data.Geo.Jord.GeoPosSpec Data.Geo.Jord.GreatCircleSpec + Data.Geo.Jord.LatLongSpec Data.Geo.Jord.LengthSpec + Data.Geo.Jord.PositionSpec Paths_jord hs-source-dirs: test ghc-options: -Wall build-depends: HUnit ==1.6.* - , base >=4.7 && <5 + , base >=4.9 && <5 , hspec ==2.* , jord default-language: Haskell2010 diff --git a/src/Data/Geo/Jord.hs b/src/Data/Geo/Jord.hs --- a/src/Data/Geo/Jord.hs +++ b/src/Data/Geo/Jord.hs @@ -16,22 +16,24 @@ module Data.Geo.Jord ( module Data.Geo.Jord.Angle , module Data.Geo.Jord.Eval - , module Data.Geo.Jord.GeoPos , module Data.Geo.Jord.GreatCircle + , module Data.Geo.Jord.LatLong , module Data.Geo.Jord.Length , module Data.Geo.Jord.NVector + , module Data.Geo.Jord.Position , module Data.Geo.Jord.Quantity , jordVersion ) where import Data.Geo.Jord.Angle import Data.Geo.Jord.Eval -import Data.Geo.Jord.GeoPos import Data.Geo.Jord.GreatCircle +import Data.Geo.Jord.LatLong import Data.Geo.Jord.Length import Data.Geo.Jord.NVector +import Data.Geo.Jord.Position import Data.Geo.Jord.Quantity -- | version. jordVersion :: String -jordVersion = "0.1.0.0" +jordVersion = "0.2.0.0" diff --git a/src/Data/Geo/Jord/Eval.hs b/src/Data/Geo/Jord/Eval.hs --- a/src/Data/Geo/Jord/Eval.hs +++ b/src/Data/Geo/Jord/Eval.hs @@ -1,6 +1,6 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# OPTIONS_GHC -fno-warn-orphans #-} - +{-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + -- | -- Module: Data.Geo.Jord.Eval -- Copyright: (c) 2018 Cedric Liegeois @@ -11,61 +11,62 @@ -- -- Types and functions for evaluating expressions in textual form. -- -module Data.Geo.Jord.Eval - ( Value(..) - , Vault - , Result - , emptyVault - , eval - , functions - , insert - , delete - , lookup - ) where - -import Control.Monad.Fail -import Data.Bifunctor -import Data.Geo.Jord.Angle -import Data.Geo.Jord.GeoPos -import Data.Geo.Jord.GreatCircle -import Data.Geo.Jord.Length -import Data.Geo.Jord.NVector -import Data.List hiding (delete, insert, lookup) -import Data.Maybe -import Prelude hiding (fail, lookup) -import Text.ParserCombinators.ReadP -import Text.Read (readMaybe) - +module Data.Geo.Jord.Eval + ( Value(..) + , Vault + , Result + , emptyVault + , eval + , functions + , insert + , delete + , lookup + ) where + +import Control.Monad.Fail +import Data.Bifunctor +import Data.Geo.Jord.Angle +import Data.Geo.Jord.GreatCircle +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Position +import Data.List hiding (delete, insert, lookup) +import Data.Maybe +import Prelude hiding (fail, lookup) +import Text.ParserCombinators.ReadP +import Text.Read (readMaybe) + -- | A value accepted and returned by 'eval'. -data Value +data Value = Ang Angle -- ^ 'Angle' | AngDec Double -- ^ 'Angle' in decimal degrees | Bool Bool -- ^ boolean | Double Double -- ^ double | Len Length -- ^ 'Length' | Gc GreatCircle -- ^ 'GreatCircle' - | Geo GeoPos -- ^ 'GeoPos' - | Geos [GeoPos] -- ^ list of 'GeoPos' - | GeoDec (Double, Double) -- ^ latitude and longitude in decimal degrees - | GeosDec [(Double, Double)] -- ^ list of latitude and longitude in decimal degrees + | Ll LatLong -- ^ 'LatLong' + | Lls [LatLong] -- ^ list of 'LatLong' + | LlDec (Double, Double) -- ^ latitude and longitude in decimal degrees + | LlsDec [(Double, Double)] -- ^ list of latitude and longitude in decimal degrees | Vec NVector -- ^ 'NVector' | Vecs [NVector] -- ^ list of 'NVector's - deriving (Eq, Show) - + deriving (Eq, Show) + -- | 'Either' an error or a 'Value'. -type Result = Either String Value - +type Result = Either String Value + -- | A location for 'Value's to be shared by successive evalations. -newtype Vault = - Vault [(String, Value)] - +newtype Vault = + Vault [(String, Value)] + -- | An empty 'Vault'. -emptyVault :: Vault -emptyVault = Vault [] - -instance MonadFail (Either String) where - fail = Left - +emptyVault :: Vault +emptyVault = Vault [] + +instance MonadFail (Either String) where + fail = Left + -- | Evaluates @s@, an expression of the form @"(f x y ..)"@. -- -- >>> eval "finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E" @@ -73,7 +74,7 @@ -- -- @f@ must be one of the supported 'functions' and each parameter @x@, @y@, .. , is either another function call -- or a 'String' parameter. Parameters are either resolved by name using the 'Resolve' --- function @r@ or if it returns 'Nothing', 'read' to an 'Angle', a 'Length' or a 'GeoPos'. +-- function @r@ or if it returns 'Nothing', 'read' to an 'Angle', a 'Length' or a 'LatLong'. -- -- If the evaluation is successful, returns the resulting 'Value' ('Right') otherwise -- a description of the error ('Left'). @@ -88,11 +89,11 @@ -- a2 = eval "(finalBearing a1 54S154W)" vault -- @ -- --- By default, all returned positions are 'Geo' ('GeoPos'), to get back a 'Vec' ('NVector'), the +-- All returned positions are 'LatLong' by default, to get back a 'NVector' the -- expression must be wrapped by 'toNVector'. -- -- @ --- dest = eval "destination 54°N,154°E 54° 1000m" -- Right Geo +-- dest = eval "destination 54°N,154°E 54° 1000m" -- Right Ll -- dest = eval "toNVector (destination 54°N,154°E 54° 1000m)" -- Right Vec -- @ -- @@ -105,20 +106,20 @@ -- length = eval "distance antipode 54N154E 54S154W" -- Left String -- @ -- -eval :: String -> Vault -> Result -eval s r = - case expr s of - Left err -> Left err - Right (rvec, ex) -> convert (evalExpr ex r) rvec - -convert :: Result -> Bool -> Result -convert r True = r -convert r False = - case r of - Right (Vec v) -> Right (Geo (fromNVector v)) - Right (Vecs vs) -> Right (Geos (map fromNVector vs)) - oth -> oth - +eval :: String -> Vault -> Result +eval s r = + case expr s of + Left err -> Left err + Right (rvec, ex) -> convert (evalExpr ex r) rvec + +convert :: Result -> Bool -> Result +convert r True = r +convert r False = + case r of + Right (Vec v) -> Right (Ll (fromNVector v)) + Right (Vecs vs) -> Right (Lls (map fromNVector vs)) + oth -> oth + -- | All supported functions: -- -- * 'antipode' @@ -149,7 +150,7 @@ -- -- * 'latLongDecimal' -- --- * 'readGeoPos' +-- * 'readLatLong' -- -- * 'toDecimalDegrees' -- @@ -161,360 +162,360 @@ -- -- * 'toNVector' -- -functions :: [String] -functions = - [ "antipode" - , "crossTrackDistance" - , "destination" - , "decimalDegrees" - , "distance" - , "finalBearing" - , "greatCircle" - , "initialBearing" - , "interpolate" - , "intersections" - , "isInside" - , "latLong" - , "latLongDecimal" - , "mean" - , "readGeoPos" - , "toDecimalDegrees" - , "toKilometres" - , "toMetres" - , "toNauticalMiles" - , "toNVector" - ] - +functions :: [String] +functions = + [ "antipode" + , "crossTrackDistance" + , "destination" + , "decimalDegrees" + , "distance" + , "finalBearing" + , "greatCircle" + , "initialBearing" + , "interpolate" + , "intersections" + , "isInside" + , "latLong" + , "latLongDecimal" + , "mean" + , "readLatLong" + , "toDecimalDegrees" + , "toKilometres" + , "toMetres" + , "toNauticalMiles" + , "toNVector" + ] + -- | @insert k v vault@ inserts value @v@ for key @k@. Overwrites any previous value. -insert :: String -> Value -> Vault -> Vault -insert k v vault = Vault (e ++ [(k, v)]) - where - Vault e = delete k vault - +insert :: String -> Value -> Vault -> Vault +insert k v vault = Vault (e ++ [(k, v)]) + where + Vault e = delete k vault + -- | @lookup k vault@ looks up the value of key @k@ in the vault. -lookup :: String -> Vault -> Maybe Value -lookup k (Vault es) = fmap snd (find (\e -> fst e == k) es) - +lookup :: String -> Vault -> Maybe Value +lookup k (Vault es) = fmap snd (find (\e -> fst e == k) es) + -- | @delete k vault@ deletes key @k@ from the vault. -delete :: String -> Vault -> Vault -delete k (Vault es) = Vault (filter (\e -> fst e /= k) es) - -expr :: (MonadFail m) => String -> m (Bool, Expr) -expr s = do - ts <- tokenise s - ast <- parse ts - fmap (\a -> (expectVec ts, a)) (transform ast) - -expectVec :: [Token] -> Bool -expectVec (_:Func "toNVector":_) = True -expectVec _ = False - -evalExpr :: Expr -> Vault -> Result -evalExpr (Param p) vault = - case lookup p vault of - Just (Geo g) -> Right (Vec (toNVector g)) - Just v -> Right v - Nothing -> tryRead p -evalExpr (Antipode a) vault = - case evalExpr a vault of - (Right (Vec p)) -> Right (Vec (antipode p)) - r -> Left ("Call error: antipode " ++ showErr [r]) -evalExpr (CrossTrackDistance a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p), Right (Gc gc)] -> Right (Len (crossTrackDistance p gc)) - r -> Left ("Call error: crossTrackDistance " ++ showErr r) -evalExpr (DecimalDegrees d) _ = Right (Ang (decimalDegrees d)) -evalExpr (Destination a b c) vault = - case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of - [Right (Vec p), Right (Ang a'), Right (Len l)] -> Right (Vec (destination p a' l)) - r -> Left ("Call error: destination " ++ showErr r) -evalExpr (Distance a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p1), Right (Vec p2)] -> Right (Len (distance p1 p2)) - r -> Left ("Call error: distance " ++ showErr r) -evalExpr (FinalBearing a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p1), Right (Vec p2)] -> Right (Ang (finalBearing p1 p2)) - r -> Left ("Call error: finalBearing " ++ showErr r) -evalExpr (GreatCircleSC a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p1), Right (Vec p2)] -> bimap id Gc (greatCircleE p1 p2) - [Right (Vec p), Right (Ang a')] -> Right (Gc (greatCircleBearing p a')) - r -> Left ("Call error: greatCircle " ++ showErr r) -evalExpr (InitialBearing a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p1), Right (Vec p2)] -> Right (Ang (initialBearing p1 p2)) - r -> Left ("Call error: initialBearing " ++ showErr r) -evalExpr (Interpolate a b c) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Vec p1), Right (Vec p2)] -> Right (Vec (interpolate p1 p2 c)) - r -> Left ("Call error: interpolate " ++ showErr r) -evalExpr (Intersections a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Gc gc1), Right (Gc gc2)] -> - maybe - (Right (Vecs [])) - (\is -> Right (Vecs [fst is, snd is])) - (intersections gc1 gc2 :: Maybe (NVector, NVector)) - r -> Left ("Call error: intersections " ++ showErr r) -evalExpr (IsInside as) vault = - let m = map (`evalExpr` vault) as - ps = [p | Right (Vec p) <- m] - in if length m == length ps && length ps > 3 - then Right (Bool (isInside (head ps) (tail ps))) - else Left ("Call error: isInside " ++ showErr m) -evalExpr (Mean as) vault = - let m = map (`evalExpr` vault) as - ps = [p | Right (Vec p) <- m] - in if length m == length ps - then maybe (Left ("Call error: mean " ++ showErr m)) (Right . Vec) (mean ps) - else Left ("Call error: mean " ++ showErr m) -evalExpr (LatLong a b) vault = - case [evalExpr a vault, evalExpr b vault] of - [Right (Ang lat), Right (Ang lon)] -> - bimap (\e -> "Call error: latLong : " ++ e) (Vec . toNVector) (latLongE lat lon) - r -> Left ("Call error: latLong " ++ showErr r) -evalExpr (LatLongDecimal a b) _ = - bimap (\e -> "Call error: LatLongDecimal : " ++ e) (Vec . toNVector) (latLongDecimalE a b) -evalExpr (ReadGeoPos s) _ = - bimap (\e -> "Call error: readGeoPos : " ++ e) (Vec . toNVector) (readGeoPosE s) -evalExpr (ToDecimalDegrees e) vault = - case evalExpr e vault of - (Right (Ang a)) -> Right (AngDec (toDecimalDegrees a)) - (Right (Geo p)) -> Right (GeoDec (toDecimalDegrees' p)) - (Right (Geos ps)) -> Right (GeosDec (map toDecimalDegrees' ps)) - (Right (Vec p)) -> Right (GeoDec ((toDecimalDegrees' . fromNVector) p)) - (Right (Vecs ps)) -> Right (GeosDec (map (toDecimalDegrees' . fromNVector) ps)) - r -> Left ("Call error: toDecimalDegrees" ++ showErr [r]) -evalExpr (ToKilometres e) vault = - case evalExpr e vault of - (Right (Len l)) -> Right (Double (toKilometres l)) - r -> Left ("Call error: toKilometres" ++ showErr [r]) -evalExpr (ToMetres e) vault = - case evalExpr e vault of - (Right (Len l)) -> Right (Double (toMetres l)) - r -> Left ("Call error: toMetres" ++ showErr [r]) -evalExpr (ToNauticalMiles e) vault = - case evalExpr e vault of - (Right (Len l)) -> Right (Double (toNauticalMiles l)) - r -> Left ("Call error: toNauticalMiles" ++ showErr [r]) -evalExpr (ToNVector a) vault = - case evalExpr a vault of - r@(Right (Vec _)) -> r - r -> Left ("Call error: toNVector " ++ showErr [r]) - -showErr :: [Result] -> String -showErr rs = " > " ++ intercalate " & " (map (either id show) rs) - -tryRead :: String -> Result -tryRead s = - case r of - [a@(Right (Ang _)), _, _] -> a - [_, l@(Right (Len _)), _] -> l - [_, _, Right (Geo g)] -> Right (Vec (toNVector g)) - _ -> Left ("couldn't read " ++ s) - where - r = map ($ s) [readE readAngleE Ang, readE readLengthE Len, readE readGeoPosE Geo] - -readE :: (String -> Either String a) -> (a -> Value) -> String -> Either String Value -readE p v s = bimap id v (p s) - +delete :: String -> Vault -> Vault +delete k (Vault es) = Vault (filter (\e -> fst e /= k) es) + +expr :: (MonadFail m) => String -> m (Bool, Expr) +expr s = do + ts <- tokenise s + ast <- parse ts + fmap (\a -> (expectVec ts, a)) (transform ast) + +expectVec :: [Token] -> Bool +expectVec (_:Func "toNVector":_) = True +expectVec _ = False + +evalExpr :: Expr -> Vault -> Result +evalExpr (Param p) vault = + case lookup p vault of + Just (Ll ll) -> Right (Vec (toNVector ll)) + Just v -> Right v + Nothing -> tryRead p +evalExpr (Antipode a) vault = + case evalExpr a vault of + (Right (Vec p)) -> Right (Vec (antipode p)) + r -> Left ("Call error: antipode " ++ showErr [r]) +evalExpr (CrossTrackDistance a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p), Right (Gc gc)] -> Right (Len (crossTrackDistance p gc)) + r -> Left ("Call error: crossTrackDistance " ++ showErr r) +evalExpr (DecimalDegrees d) _ = Right (Ang (decimalDegrees d)) +evalExpr (Destination a b c) vault = + case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of + [Right (Vec p), Right (Ang a'), Right (Len l)] -> Right (Vec (destination p a' l)) + r -> Left ("Call error: destination " ++ showErr r) +evalExpr (Distance a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p1), Right (Vec p2)] -> Right (Len (distance p1 p2)) + r -> Left ("Call error: distance " ++ showErr r) +evalExpr (FinalBearing a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p1), Right (Vec p2)] -> Right (Ang (finalBearing p1 p2)) + r -> Left ("Call error: finalBearing " ++ showErr r) +evalExpr (GreatCircleSC a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p1), Right (Vec p2)] -> bimap id Gc (greatCircleE p1 p2) + [Right (Vec p), Right (Ang a')] -> Right (Gc (greatCircleBearing p a')) + r -> Left ("Call error: greatCircle " ++ showErr r) +evalExpr (InitialBearing a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p1), Right (Vec p2)] -> Right (Ang (initialBearing p1 p2)) + r -> Left ("Call error: initialBearing " ++ showErr r) +evalExpr (Interpolate a b c) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Vec p1), Right (Vec p2)] -> Right (Vec (interpolate p1 p2 c)) + r -> Left ("Call error: interpolate " ++ showErr r) +evalExpr (Intersections a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Gc gc1), Right (Gc gc2)] -> + maybe + (Right (Vecs [])) + (\is -> Right (Vecs [fst is, snd is])) + (intersections gc1 gc2 :: Maybe (NVector, NVector)) + r -> Left ("Call error: intersections " ++ showErr r) +evalExpr (IsInside as) vault = + let m = map (`evalExpr` vault) as + ps = [p | Right (Vec p) <- m] + in if length m == length ps && length ps > 3 + then Right (Bool (isInside (head ps) (tail ps))) + else Left ("Call error: isInside " ++ showErr m) +evalExpr (Mean as) vault = + let m = map (`evalExpr` vault) as + ps = [p | Right (Vec p) <- m] + in if length m == length ps + then maybe (Left ("Call error: mean " ++ showErr m)) (Right . Vec) (mean ps) + else Left ("Call error: mean " ++ showErr m) +evalExpr (LatLong a b) vault = + case [evalExpr a vault, evalExpr b vault] of + [Right (Ang lat), Right (Ang lon)] -> + bimap (\e -> "Call error: latLong : " ++ e) (Vec . toNVector) (latLongE lat lon) + r -> Left ("Call error: latLong " ++ showErr r) +evalExpr (LatLongDecimal a b) _ = + bimap (\e -> "Call error: LatLongDecimal : " ++ e) (Vec . toNVector) (latLongDecimalE a b) +evalExpr (ReadLatLong s) _ = + bimap (\e -> "Call error: readLatLong : " ++ e) (Vec . toNVector) (readLatLongE s) +evalExpr (ToDecimalDegrees e) vault = + case evalExpr e vault of + (Right (Ang a)) -> Right (AngDec (toDecimalDegrees a)) + (Right (Ll p)) -> Right (LlDec (toDecimalDegrees' p)) + (Right (Lls ps)) -> Right (LlsDec (map toDecimalDegrees' ps)) + (Right (Vec p)) -> Right (LlDec ((toDecimalDegrees' . fromNVector) p)) + (Right (Vecs ps)) -> Right (LlsDec (map (toDecimalDegrees' . fromNVector) ps)) + r -> Left ("Call error: toDecimalDegrees" ++ showErr [r]) +evalExpr (ToKilometres e) vault = + case evalExpr e vault of + (Right (Len l)) -> Right (Double (toKilometres l)) + r -> Left ("Call error: toKilometres" ++ showErr [r]) +evalExpr (ToMetres e) vault = + case evalExpr e vault of + (Right (Len l)) -> Right (Double (toMetres l)) + r -> Left ("Call error: toMetres" ++ showErr [r]) +evalExpr (ToNauticalMiles e) vault = + case evalExpr e vault of + (Right (Len l)) -> Right (Double (toNauticalMiles l)) + r -> Left ("Call error: toNauticalMiles" ++ showErr [r]) +evalExpr (ToNVector a) vault = + case evalExpr a vault of + r@(Right (Vec _)) -> r + r -> Left ("Call error: toNVector " ++ showErr [r]) + +showErr :: [Result] -> String +showErr rs = " > " ++ intercalate " & " (map (either id show) rs) + +tryRead :: String -> Result +tryRead s = + case r of + [a@(Right (Ang _)), _, _] -> a + [_, l@(Right (Len _)), _] -> l + [_, _, Right (Ll ll)] -> Right (Vec (toNVector ll)) + _ -> Left ("couldn't read " ++ s) + where + r = map ($ s) [readE readAngleE Ang, readE readLengthE Len, readE readLatLongE Ll] + +readE :: (String -> Either String a) -> (a -> Value) -> String -> Either String Value +readE p v s = bimap id v (p s) + ------------------------------------------ -- Lexical Analysis: String -> [Token] -- ------------------------------------------ -data Token - = Paren Char - | Func String - | Str String - deriving (Show) - -tokenise :: (MonadFail m) => String -> m [Token] -tokenise s - | null r = fail ("Lexical error: " ++ s) - | (e, "") <- last r = return (wrap e) - | otherwise = fail ("Lexical error: " ++ snd (last r)) - where - r = readP_to_S tokens s - +data Token + = Paren Char + | Func String + | Str String + deriving (Show) + +tokenise :: (MonadFail m) => String -> m [Token] +tokenise s + | null r = fail ("Lexical error: " ++ s) + | (e, "") <- last r = return (wrap e) + | otherwise = fail ("Lexical error: " ++ snd (last r)) + where + r = readP_to_S tokens s + -- | wraps top level expression between () if needed. -wrap :: [Token] -> [Token] -wrap ts - | null ts = ts - | (Paren '(') <- head ts = ts - | otherwise = Paren '(' : ts ++ [Paren ')'] - -tokens :: ReadP [Token] -tokens = many1 token - -token :: ReadP Token -token = (<++) ((<++) paren func) str - -paren :: ReadP Token -paren = (<++) parenO parenC - -parenO :: ReadP Token -parenO = do - optional (char ' ') - c <- char '(' - return (Paren c) - -parenC :: ReadP Token -parenC = do - c <- char ')' - optional (char ' ') - return (Paren c) - -func :: ReadP Token -func = do - n <- choice (map string functions) - _ <- char ' ' - return (Func n) - -str :: ReadP Token -str = do - optional (char ' ') - v <- munch1 (\c -> c /= '(' && c /= ')' && c /= ' ') - if v `elem` functions - then pfail - else return (Str v) - +wrap :: [Token] -> [Token] +wrap ts + | null ts = ts + | (Paren '(') <- head ts = ts + | otherwise = Paren '(' : ts ++ [Paren ')'] + +tokens :: ReadP [Token] +tokens = many1 token + +token :: ReadP Token +token = (<++) ((<++) paren func) str + +paren :: ReadP Token +paren = (<++) parenO parenC + +parenO :: ReadP Token +parenO = do + optional (char ' ') + c <- char '(' + return (Paren c) + +parenC :: ReadP Token +parenC = do + c <- char ')' + optional (char ' ') + return (Paren c) + +func :: ReadP Token +func = do + n <- choice (map string functions) + _ <- char ' ' + return (Func n) + +str :: ReadP Token +str = do + optional (char ' ') + v <- munch1 (\c -> c /= '(' && c /= ')' && c /= ' ') + if v `elem` functions + then pfail + else return (Str v) + ----------------------------------------- -- Syntactic Analysis: [Token] -> Ast -- ----------------------------------------- -data Ast - = Call String - [Ast] - | Lit String - deriving (Show) - +data Ast + = Call String + [Ast] + | Lit String + deriving (Show) + -- | syntax is (f x y) where x and y can be function themselves. -parse :: (MonadFail m) => [Token] -> m Ast -parse ts = fmap fst (walk ts) - -walk :: (MonadFail m) => [Token] -> m (Ast, [Token]) -walk [] = fail "Syntax error: empty" -walk (h:t) - | (Str s) <- h = return (Lit s, t) - | (Paren '(') <- h = walkFunc t - | otherwise = fail ("Syntax error: expected String or '(' but got " ++ show h) - -walkFunc :: (MonadFail m) => [Token] -> m (Ast, [Token]) -walkFunc [] = fail "Syntax error: '(' unexpected" -walkFunc (h:t) - | (Func n) <- h = walkParams n t [] - | otherwise = fail ("Syntax error: expected Function but got " ++ show h) - -walkParams :: (MonadFail m) => String -> [Token] -> [Ast] -> m (Ast, [Token]) -walkParams _ [] _ = fail "Syntax error: ')' not found" -walkParams n ts@(h:t) acc - | (Paren ')') <- h = return (Call n (reverse acc), t) - | otherwise = do - (el, t') <- walk ts - walkParams n t' (el : acc) - +parse :: (MonadFail m) => [Token] -> m Ast +parse ts = fmap fst (walk ts) + +walk :: (MonadFail m) => [Token] -> m (Ast, [Token]) +walk [] = fail "Syntax error: empty" +walk (h:t) + | (Str s) <- h = return (Lit s, t) + | (Paren '(') <- h = walkFunc t + | otherwise = fail ("Syntax error: expected String or '(' but got " ++ show h) + +walkFunc :: (MonadFail m) => [Token] -> m (Ast, [Token]) +walkFunc [] = fail "Syntax error: '(' unexpected" +walkFunc (h:t) + | (Func n) <- h = walkParams n t [] + | otherwise = fail ("Syntax error: expected Function but got " ++ show h) + +walkParams :: (MonadFail m) => String -> [Token] -> [Ast] -> m (Ast, [Token]) +walkParams _ [] _ = fail "Syntax error: ')' not found" +walkParams n ts@(h:t) acc + | (Paren ')') <- h = return (Call n (reverse acc), t) + | otherwise = do + (el, t') <- walk ts + walkParams n t' (el : acc) + ------------------------------------- -- Semantic Analysis: Ast -> Expr -- ------------------------------------- -data Expr - = Param String - | Antipode Expr - | CrossTrackDistance Expr - Expr - | DecimalDegrees Double - | Destination Expr - Expr - Expr - | Distance Expr - Expr - | FinalBearing Expr - Expr - | GreatCircleSC Expr - Expr - | InitialBearing Expr - Expr - | Interpolate Expr - Expr - Double - | Intersections Expr - Expr - | IsInside [Expr] - | Mean [Expr] - | LatLong Expr - Expr - | LatLongDecimal Double - Double - | ReadGeoPos String - | ToDecimalDegrees Expr - | ToKilometres Expr - | ToMetres Expr - | ToNauticalMiles Expr - | ToNVector Expr - deriving (Show) - -transform :: (MonadFail m) => Ast -> m Expr -transform (Call "antipode" [e]) = fmap Antipode (transform e) -transform (Call "crossTrackDistance" [e1, e2]) = do - p <- transform e1 - gc <- transform e2 - return (CrossTrackDistance p gc) -transform (Call "decimalDegrees" [Lit s]) = fmap DecimalDegrees (readDouble s) -transform (Call "destination" [e1, e2, e3]) = do - p1 <- transform e1 - p2 <- transform e2 - p3 <- transform e3 - return (Destination p1 p2 p3) -transform (Call "distance" [e1, e2]) = do - p1 <- transform e1 - p2 <- transform e2 - return (Distance p1 p2) -transform (Call "finalBearing" [e1, e2]) = do - p1 <- transform e1 - p2 <- transform e2 - return (FinalBearing p1 p2) -transform (Call "greatCircle" [e1, e2]) = do - p1 <- transform e1 - p2 <- transform e2 - return (GreatCircleSC p1 p2) -transform (Call "initialBearing" [e1, e2]) = do - p1 <- transform e1 - p2 <- transform e2 - return (InitialBearing p1 p2) -transform (Call "interpolate" [e1, e2, Lit s]) = do - p1 <- transform e1 - p2 <- transform e2 - d <- readDouble s - if d >= 0.0 && d <= 1.0 - then return (Interpolate p1 p2 d) - else fail "Semantic error: interpolate expects [0..1] as last argument" -transform (Call "intersections" [e1, e2]) = do - gc1 <- transform e1 - gc2 <- transform e2 - return (Intersections gc1 gc2) -transform (Call "isInside" e) = do - ps <- mapM transform e - return (IsInside ps) -transform (Call "latLong" [e1, e2]) = do - gc1 <- transform e1 - gc2 <- transform e2 - return (LatLong gc1 gc2) -transform (Call "latLongDecimal" [Lit s1, Lit s2]) = do - d1 <- readDouble s1 - d2 <- readDouble s2 - return (LatLongDecimal d1 d2) -transform (Call "mean" e) = do - ps <- mapM transform e - return (Mean ps) -transform (Call "readGeoPos" [Lit s]) = return (ReadGeoPos s) -transform (Call "toDecimalDegrees" [e]) = fmap ToDecimalDegrees (transform e) -transform (Call "toKilometres" [e]) = fmap ToKilometres (transform e) -transform (Call "toMetres" [e]) = fmap ToMetres (transform e) -transform (Call "toNauticalMiles" [e]) = fmap ToNauticalMiles (transform e) -transform (Call "toNVector" [e]) = fmap ToNVector (transform e) -transform (Call f e) = fail ("Semantic error: " ++ f ++ " does not accept " ++ show e) -transform (Lit s) = return (Param s) - -readDouble :: (MonadFail m) => String -> m Double -readDouble s = - case readMaybe s of - Just d -> return d - Nothing -> fail ("Unparsable double: " ++ s) +data Expr + = Param String + | Antipode Expr + | CrossTrackDistance Expr + Expr + | DecimalDegrees Double + | Destination Expr + Expr + Expr + | Distance Expr + Expr + | FinalBearing Expr + Expr + | GreatCircleSC Expr + Expr + | InitialBearing Expr + Expr + | Interpolate Expr + Expr + Double + | Intersections Expr + Expr + | IsInside [Expr] + | Mean [Expr] + | LatLong Expr + Expr + | LatLongDecimal Double + Double + | ReadLatLong String + | ToDecimalDegrees Expr + | ToKilometres Expr + | ToMetres Expr + | ToNauticalMiles Expr + | ToNVector Expr + deriving (Show) + +transform :: (MonadFail m) => Ast -> m Expr +transform (Call "antipode" [e]) = fmap Antipode (transform e) +transform (Call "crossTrackDistance" [e1, e2]) = do + p <- transform e1 + gc <- transform e2 + return (CrossTrackDistance p gc) +transform (Call "decimalDegrees" [Lit s]) = fmap DecimalDegrees (readDouble s) +transform (Call "destination" [e1, e2, e3]) = do + p1 <- transform e1 + p2 <- transform e2 + p3 <- transform e3 + return (Destination p1 p2 p3) +transform (Call "distance" [e1, e2]) = do + p1 <- transform e1 + p2 <- transform e2 + return (Distance p1 p2) +transform (Call "finalBearing" [e1, e2]) = do + p1 <- transform e1 + p2 <- transform e2 + return (FinalBearing p1 p2) +transform (Call "greatCircle" [e1, e2]) = do + p1 <- transform e1 + p2 <- transform e2 + return (GreatCircleSC p1 p2) +transform (Call "initialBearing" [e1, e2]) = do + p1 <- transform e1 + p2 <- transform e2 + return (InitialBearing p1 p2) +transform (Call "interpolate" [e1, e2, Lit s]) = do + p1 <- transform e1 + p2 <- transform e2 + d <- readDouble s + if d >= 0.0 && d <= 1.0 + then return (Interpolate p1 p2 d) + else fail "Semantic error: interpolate expects [0..1] as last argument" +transform (Call "intersections" [e1, e2]) = do + gc1 <- transform e1 + gc2 <- transform e2 + return (Intersections gc1 gc2) +transform (Call "isInside" e) = do + ps <- mapM transform e + return (IsInside ps) +transform (Call "latLong" [e1, e2]) = do + gc1 <- transform e1 + gc2 <- transform e2 + return (LatLong gc1 gc2) +transform (Call "latLongDecimal" [Lit s1, Lit s2]) = do + d1 <- readDouble s1 + d2 <- readDouble s2 + return (LatLongDecimal d1 d2) +transform (Call "mean" e) = do + ps <- mapM transform e + return (Mean ps) +transform (Call "readLatLong" [Lit s]) = return (ReadLatLong s) +transform (Call "toDecimalDegrees" [e]) = fmap ToDecimalDegrees (transform e) +transform (Call "toKilometres" [e]) = fmap ToKilometres (transform e) +transform (Call "toMetres" [e]) = fmap ToMetres (transform e) +transform (Call "toNauticalMiles" [e]) = fmap ToNauticalMiles (transform e) +transform (Call "toNVector" [e]) = fmap ToNVector (transform e) +transform (Call f e) = fail ("Semantic error: " ++ f ++ " does not accept " ++ show e) +transform (Lit s) = return (Param s) + +readDouble :: (MonadFail m) => String -> m Double +readDouble s = + case readMaybe s of + Just d -> return d + Nothing -> fail ("Unparsable double: " ++ s) diff --git a/src/Data/Geo/Jord/GeoPos.hs b/src/Data/Geo/Jord/GeoPos.hs deleted file mode 100644 --- a/src/Data/Geo/Jord/GeoPos.hs +++ /dev/null @@ -1,221 +0,0 @@ --- | --- Module: Data.Geo.Jord.GeoPos --- Copyright: (c) 2018 Cedric Liegeois --- License: BSD3 --- Maintainer: Cedric Liegeois --- Stability: experimental --- Portability: portable --- --- Types to represent a geographic position by its latitude and longitude. --- -module Data.Geo.Jord.GeoPos - ( - -- * The 'GeoPos' type - GeoPos(latitude, longitude) - -- * Smart constructors - , latLong - , latLongE - , latLongF - , latLongDecimal - , latLongDecimalE - , latLongDecimalF - -- * read - , readGeoPos - , readGeoPosE - , readGeoPosF - -- * Misc. - , toDecimalDegrees' - ) where - -import Control.Applicative hiding (many) -import Control.Monad.Fail -import Data.Char -import Data.Geo.Jord.Angle -import Data.Geo.Jord.Parse -import Data.Maybe -import Prelude hiding (fail) -import Text.ParserCombinators.ReadP -import Text.Read hiding (pfail) - --- | A geographic position (latitude and longitude). -data GeoPos = GeoPos - { latitude :: Angle - , longitude :: Angle - } deriving (Eq) - --- | See 'readGeoPos'. -instance Read GeoPos where - readsPrec _ = readP_to_S geo - --- | Produced string format: d°(m')(s'')[N|S],d°(m')(s'')[E|W] - e.g. 55°36'21''N,13°0'2''E. -instance Show GeoPos where - show (GeoPos lat lon) = showLat lat ++ "," ++ showLon lon - --- | 'GeoPos' from given latitude and longitude. --- 'error's if given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLong :: Angle -> Angle -> GeoPos -latLong lat lon = - fromMaybe - (error ("Invalid latitude=" ++ show lat ++ " or longitude=" ++ show lon)) - (latLongF lat lon) - --- | 'GeoPos' from given latitude and longitude. --- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLongE :: Angle -> Angle -> Either String GeoPos -latLongE lat lon - | not (isWithin lat (decimalDegrees (-90)) (decimalDegrees 90)) = - Left ("Invalid latitude=" ++ show lat) - | not (isWithin lon (decimalDegrees (-180)) (decimalDegrees 180)) = - Left ("Invalid longitude=" ++ show lon) - | otherwise = Right (GeoPos lat lon) - --- | 'GeoPos' from given latitude and longitude. --- 'fail's if given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLongF :: (MonadFail m) => Angle -> Angle -> m GeoPos -latLongF lat lon = - case e of - Left err -> fail err - Right g -> return g - where - e = latLongE lat lon - --- | 'GeoPos' from given latitude and longitude in decimal degrees. --- 'error's if given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLongDecimal :: Double -> Double -> GeoPos -latLongDecimal lat lon = latLong (decimalDegrees lat) (decimalDegrees lon) - --- | 'GeoPos' from given latitude and longitude in decimal degrees. --- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLongDecimalE :: Double -> Double -> Either String GeoPos -latLongDecimalE lat lon = latLongE (decimalDegrees lat) (decimalDegrees lon) - --- | 'GeoPos' from given latitude and longitude in decimal degrees. --- 'fail's if given latitude is outisde [-90..90]° and/or --- given longitude is outisde [-180..180]°. -latLongDecimalF :: (MonadFail m) => Double -> Double -> m GeoPos -latLongDecimalF lat lon = latLongF (decimalDegrees lat) (decimalDegrees lon) - --- | Obtains a 'GeoPos' from the given string formatted as either: --- --- * DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - e.g. 553621N0130002E or 0116S03649E or 47N122W --- --- * 'Angle'[N|S] 'Angle'[E|W] - e.g. 55°36'21''N 13°0'02''E or 11°16'S 36°49'E or 47°N 122°W --- --- This simply calls @read s :: GeoPos@ so 'error' should be handled at the call site. --- -readGeoPos :: String -> GeoPos -readGeoPos s = read s :: GeoPos - --- | Same as 'readGeoPos' but returns a 'Either'. -readGeoPosE :: String -> Either String GeoPos -readGeoPosE s = - case readMaybe s of - Nothing -> Left ("couldn't read geo pos " ++ s) - Just g -> Right g - --- | Same as 'readGeoPos' but returns a 'MonadFail'. -readGeoPosF :: (MonadFail m) => String -> m GeoPos -readGeoPosF s = - let pg = readGeoPosE s - in case pg of - Left e -> fail e - Right g -> return g - --- | Converts the given 'GeoPos' to tuple of latitude and longitude in decimal degrees. -toDecimalDegrees' :: GeoPos -> (Double, Double) -toDecimalDegrees' g = (toDecimalDegrees (latitude g), toDecimalDegrees (longitude g)) - --- | Parses and returns a 'GeoPos'. -geo :: ReadP GeoPos -geo = block <|> human - --- | Parses and returns a 'GeoPos' - DD(D)MMSS. -block :: ReadP GeoPos -block = do - lat <- blat - lon <- blon - latLongF lat lon - --- | Parses and returns a latitude, DDMMSS expected. -blat :: ReadP Angle -blat = do - d' <- digits 2 - (m', s') <- option (0, 0) (ms <|> m) - h <- hemisphere - if h == 'N' - then dmsF d' m' s' 0 - else dmsF (-d') m' s' 0 - --- | Parses and returns a longitude, DDDMMSS expected. -blon :: ReadP Angle -blon = do - d' <- digits 3 - (m', s') <- option (0, 0) (ms <|> m) - m'' <- meridian - if m'' == 'E' - then dmsF d' m' s' 0 - else dmsF (-d') m' s' 0 - --- | Parses N or S char. -hemisphere :: ReadP Char -hemisphere = char 'N' <|> char 'S' - --- | Parses E or W char. -meridian :: ReadP Char -meridian = char 'E' <|> char 'W' - --- | Parses minutes and seconds. -ms :: ReadP (Int, Int) -ms = do - m' <- digits 2 - s' <- digits 2 - return (m', s') - --- | Parses minutes. -m :: ReadP (Int, Int) -m = do - m' <- digits 2 - return (m', 0) - --- | Parses and returns a 'GeoPos' from a human friendly text - see 'Angle'. -human :: ReadP GeoPos -human = do - lat <- hlat - _ <- char ' ' <|> char ',' - lon <- hlon - latLongF lat lon - --- | Parses and returns a latitude, 'Angle'N|S expected. -hlat :: ReadP Angle -hlat = do - lat <- angle - h <- hemisphere - if h == 'N' - then return lat - else return (negate' lat) - --- | Parses and returns a longitude, 'Angle'E|W expected. -hlon :: ReadP Angle -hlon = do - lon <- angle - m' <- meridian - if m' == 'E' - then return lon - else return (negate' lon) - --- | Latitude to string. -showLat :: Angle -> String -showLat lat - | isNegative lat = show (negate' lat) ++ "S" - | otherwise = show lat ++ "N" - --- | Longitude to string. -showLon :: Angle -> String -showLon lon - | isNegative lon = show (negate' lon) ++ "W" - | otherwise = show lon ++ "E" diff --git a/src/Data/Geo/Jord/GreatCircle.hs b/src/Data/Geo/Jord/GreatCircle.hs --- a/src/Data/Geo/Jord/GreatCircle.hs +++ b/src/Data/Geo/Jord/GreatCircle.hs @@ -13,59 +13,32 @@ -- -- This module assumes a spherical earth. -- --- TODO: --- --- * alongTrackDistance :: Position -> GreatArc -> Length --- --- * intersection :: GreatArc -> GreatArc -> Maybe Position --- --- * nearestPoint :: Position -> GreatArc -> Position --- --- * area :: [Position] -> Surface --- --- * closestApproach --- -module Data.Geo.Jord.GreatCircle - ( - -- * The 'GreatCircle' type - GreatCircle - -- * The 'Position' type - , Position(..) - -- * Smart constructors - , greatCircle - , greatCircleE - , greatCircleF - , greatCircleBearing - -- * Geodesic calculations - , antipode - , crossTrackDistance - , crossTrackDistance' - , destination - , destination' - , distance - , distance' - , finalBearing - , initialBearing - , interpolate - , intersections - , isInside - , mean - -- * Misc. - , meanEarthRadius - , northPole - , southPole - ) where - -import Control.Monad.Fail -import Data.Geo.Jord.Angle -import Data.Geo.Jord.GeoPos -import Data.Geo.Jord.Length -import Data.Geo.Jord.NVector -import Data.Geo.Jord.Quantity -import Data.List (subsequences) -import Data.Maybe (fromMaybe) -import Prelude hiding (fail) - +module Data.Geo.Jord.GreatCircle + ( + -- * The 'GreatCircle' type + GreatCircle + -- * Smart constructors + , greatCircle + , greatCircleE + , greatCircleF + , greatCircleBearing + -- * Geodesic calculations + , crossTrackDistance + , crossTrackDistance' + , intersections + , isInside + ) where + +import Control.Monad.Fail +import Data.Geo.Jord.Angle +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Position +import Data.Geo.Jord.Quantity +import Data.Maybe (fromMaybe) +import Prelude hiding (fail) + -- | A circle on the surface of the Earth which lies in a plane passing through -- the Earth's centre. Every two distinct and non-antipodal points on the surface -- of the Earth define a Great Circle. @@ -75,267 +48,86 @@ -- -- See 'greatCircle', 'greatCircleE', 'greatCircleF' or 'greatCircleBearing' constructors. -- -data GreatCircle = GreatCircle - { normal :: NVector - , dscr :: String - } deriving (Eq) - -instance Show GreatCircle where - show = dscr - --- | The 'Position' class defines 2 functions to convert a position to and from a 'NVector'. --- All functions in this module first convert 'Position' to 'NVector' and any resulting 'NVector' back --- to a 'Position'. This allows the call site to pass either 'NVector' or 'GeoPos' and to get back --- the same class instance. -class Position a where - -- | Converts a 'NVector' into 'Position' instance. - fromNVector :: NVector -> a - -- | Converts the 'Position' instance into a 'NVector'. - toNVector :: a -> NVector - --- | 'GeoPos' to/from 'NVector'. -instance Position GeoPos where - fromNVector v = latLong lat lon - where - lat = atan2' (z v) (sqrt (x v * x v + y v * y v)) - lon = atan2' (y v) (x v) - toNVector g = nvector x' y' z' - where - lat = latitude g - lon = longitude g - cl = cos' lat - x' = cl * cos' lon - y' = cl * sin' lon - z' = sin' lat - --- | Identity. -instance Position NVector where - fromNVector v = v - toNVector v = v - --- | 'GreateCircle' passing by both given 'Position's. 'error's if given positions are --- equal or antipodal. -greatCircle :: (Eq a, Position a, Show a) => a -> a -> GreatCircle -greatCircle p1 p2 = - fromMaybe - (error (show p1 ++ " and " ++ show p2 ++ " do not define a unique Great Circle")) - (greatCircleF p1 p2) - --- | 'GreateCircle' passing by both given 'Position's. A 'Left' indicates that given positions are --- equal or antipodal. -greatCircleE :: (Eq a, Position a) => a -> a -> Either String GreatCircle -greatCircleE p1 p2 - | p1 == p2 = Left "Invalid Great Circle: positions are equal" - | p1 == antipode p2 = Left "Invalid Great Circle: positions are antipodal" - | otherwise = - Right - (GreatCircle - (cross v1 v2) - ("passing by " ++ - show (fromNVector v1 :: GeoPos) ++ " & " ++ show (fromNVector v2 :: GeoPos))) - where - v1 = toNVector p1 - v2 = toNVector p2 - --- | 'GreateCircle' passing by both given 'Position's. 'fail's if given positions are --- equal or antipodal. -greatCircleF :: (Eq a, MonadFail m, Position a) => a -> a -> m GreatCircle -greatCircleF p1 p2 = - case e of - Left err -> fail err - Right gc -> return gc - where - e = greatCircleE p1 p2 - --- | 'GreatCircle' passing by the given 'Position' and heading on given bearing. -greatCircleBearing :: (Position a) => a -> Angle -> GreatCircle -greatCircleBearing p b = - GreatCircle - (sub n' e') - ("passing by " ++ show (fromNVector v :: GeoPos) ++ " heading on " ++ show b) - where - v = toNVector p - e = cross northPole v -- easting - n = cross v e -- northing - e' = scale e (cos' b / norm e) - n' = scale n (sin' b / norm n) - --- | Returns the antipodal 'Position' of the given 'Position' - i.e. the position on the surface --- of the Earth which is diametrically opposite to the given position. -antipode :: (Position a) => a -> a -antipode p = fromNVector (scale (toNVector p) (-1.0)) - --- | 'crossTrackDistance'' assuming a radius of 'meanEarthRadius'. -crossTrackDistance :: (Position a) => a -> GreatCircle -> Length -crossTrackDistance p gc = crossTrackDistance' p gc meanEarthRadius - --- | Signed distance from given 'Position' to given 'GreatCircle'. --- Returns a negative 'Length' if position if left of great circle, --- positive 'Length' if position if right of great circle; the orientation of the --- great circle is therefore important: --- --- @ --- let gc1 = greatCircle (latLongDecimal 51 0) (latLongDecimal 52 1) --- let gc2 = greatCircle (latLongDecimal 52 1) (latLongDecimal 51 0) --- crossTrackDistance p gc1 == (- crossTrackDistance p gc2) --- @ -crossTrackDistance' :: (Position a) => a -> GreatCircle -> Length -> Length -crossTrackDistance' p gc = - arcLength (sub (angleBetween (normal gc) (toNVector p) Nothing) (decimalDegrees 90)) - --- | 'destination'' assuming a radius of 'meanEarthRadius'. -destination :: (Position a) => a -> Angle -> Length -> a -destination p b d = destination' p b d meanEarthRadius - --- | Computes the destination 'Position' from the given 'Position' having travelled the given distance on the --- given initial bearing (bearing will normally vary before destination is reached) and using the given earth radius. --- --- This is known as the direct geodetic problem. -destination' :: (Position a) => a -> Angle -> Length -> Length -> a -destination' p b d r - | isZero d = p - | otherwise = fromNVector (add (scale v (cos' ta)) (scale de (sin' ta))) - where - v = toNVector p - ed = unit (cross northPole v) -- east direction vector at v - nd = cross v ed -- north direction vector at v - ta = central d r -- central angle - de = add (scale nd (cos' b)) (scale ed (sin' b)) -- unit vector in the direction of the azimuth - --- | 'distance'' assuming a radius of 'meanEarthRadius'. -distance :: (Position a) => a -> a -> Length -distance p1 p2 = distance' p1 p2 meanEarthRadius - --- | Computes the surface distance (length of geodesic) in 'Meters' assuming a --- spherical Earth between the two given 'Position's and using the given earth radius. -distance' :: (Position a) => a -> a -> Length -> Length -distance' p1 p2 = arcLength (angleBetween v1 v2 Nothing) - where - v1 = toNVector p1 - v2 = toNVector p2 - --- | Computes the final bearing arriving at given destination @p2@ 'Position' from given 'Position' @p1@. --- the final bearing will differ from the 'initialBearing' by varying degrees according to distance and latitude. --- Returns 180 if both position are equals. -finalBearing :: (Position a) => a -> a -> Angle -finalBearing p1 p2 = normalise (initialBearing p2 p1) (decimalDegrees 180) - --- | Computes the initial bearing from given @p1@ 'Position' to given @p2@ 'Position', in compass degrees. --- Returns 0 if both position are equals. -initialBearing :: (Position a) => a -> a -> Angle -initialBearing p1 p2 = normalise (angleBetween gc1 gc2 (Just v1)) (decimalDegrees 360) - where - v1 = toNVector p1 - v2 = toNVector p2 - gc1 = cross v1 v2 -- great circle through p1 & p2 - gc2 = cross v1 northPole -- great circle through p1 & north pole - --- | Computes the 'Position' at given fraction @f@ between the two given 'Position's @p0@ and @p1@. --- --- Special conditions: +data GreatCircle = GreatCircle + { normal :: NVector + , dscr :: String + } deriving (Eq) + +instance Show GreatCircle where + show = dscr + +-- | 'GreateCircle' passing by both given 'Position's. 'error's if given positions are +-- equal or antipodal. +greatCircle :: (Eq a, Position a, Show a) => a -> a -> GreatCircle +greatCircle p1 p2 = + fromMaybe + (error (show p1 ++ " and " ++ show p2 ++ " do not define a unique Great Circle")) + (greatCircleF p1 p2) + +-- | 'GreateCircle' passing by both given 'Position's. A 'Left' indicates that given positions are +-- equal or antipodal. +greatCircleE :: (Eq a, Position a) => a -> a -> Either String GreatCircle +greatCircleE p1 p2 + | p1 == p2 = Left "Invalid Great Circle: positions are equal" + | p1 == antipode p2 = Left "Invalid Great Circle: positions are antipodal" + | otherwise = + Right + (GreatCircle + (cross v1 v2) + ("passing by " ++ + show (fromNVector v1 :: LatLong) ++ " & " ++ show (fromNVector v2 :: LatLong))) + where + v1 = toNVector p1 + v2 = toNVector p2 + +-- | 'GreateCircle' passing by both given 'Position's. 'fail's if given positions are +-- equal or antipodal. +greatCircleF :: (Eq a, MonadFail m, Position a) => a -> a -> m GreatCircle +greatCircleF p1 p2 = + case e of + Left err -> fail err + Right gc -> return gc + where + e = greatCircleE p1 p2 + +-- | 'GreatCircle' passing by the given 'Position' and heading on given bearing. +greatCircleBearing :: (Position a) => a -> Angle -> GreatCircle +greatCircleBearing p b = + GreatCircle + (sub n' e') + ("passing by " ++ show (fromNVector v :: LatLong) ++ " heading on " ++ show b) + where + v = toNVector p + e = cross northPole v -- easting + n = cross v e -- northing + e' = scale e (cos' b / norm e) + n' = scale n (sin' b / norm n) + +-- | 'crossTrackDistance'' assuming a radius of 'meanEarthRadius'. +crossTrackDistance :: (Position a) => a -> GreatCircle -> Length +crossTrackDistance p gc = crossTrackDistance' p gc meanEarthRadius + +-- | Signed distance from given 'Position' to given 'GreatCircle'. +-- Returns a negative 'Length' if position if left of great circle, +-- positive 'Length' if position if right of great circle; the orientation of the +-- great circle is therefore important: -- -- @ --- interpolate p0 p1 0.0 => p0 --- interpolate p0 p1 1.0 => p1 +-- let gc1 = greatCircle (latLongDecimal 51 0) (latLongDecimal 52 1) +-- let gc2 = greatCircle (latLongDecimal 52 1) (latLongDecimal 51 0) +-- crossTrackDistance p gc1 == (- crossTrackDistance p gc2) -- @ --- --- 'error's if @f < 0 || f > 1.0@ --- -interpolate :: (Position a) => a -> a -> Double -> a -interpolate p0 p1 f - | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f) - | f == 0 = p0 - | f == 1 = p1 - | otherwise = fromNVector (unit (add v0 (scale (sub v1 v0) f))) - where - v0 = toNVector p0 - v1 = toNVector p1 - +crossTrackDistance' :: (Position a) => a -> GreatCircle -> Length -> Length +crossTrackDistance' p gc = + arcLength (sub (angularDistance (normal gc) (toNVector p) Nothing) (decimalDegrees 90)) + -- | Computes the intersections between the two given 'GreatCircle's. -- Two 'GreatCircle's intersect exactly twice unless there are equal (regardless of orientation), -- in which case 'Nothing' is returned. -intersections :: (Position a) => GreatCircle -> GreatCircle -> Maybe (a, a) -intersections gc1 gc2 - | norm i == 0.0 = Nothing - | otherwise - , let ni = unit i = Just (fromNVector ni, fromNVector (antipode ni)) - where - i = cross (normal gc1) (normal gc2) - --- | Determines whether the given 'Position' is inside the polygon defined by the given list of 'Position's. --- The polygon is closed if needed (i.e. if @head ps /= last ps@). --- --- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles --- will sum to less than 360°, and exterior point angles will be small but non-zero. --- --- Always returns 'False' if positions does not at least defines a triangle. --- -isInside :: (Eq a, Position a) => a -> [a] -> Bool -isInside p ps - | null ps = False - | head ps == last ps = isInside p (init ps) - | length ps < 3 = False - | otherwise = - let aSum = foldl (\a v' -> add a (uncurry angleBetween v' (Just v))) (decimalDegrees 0) es - in abs (toDecimalDegrees aSum) > 180.0 - where - v = toNVector p - es = egdes (map (sub v . toNVector) ps) - --- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] -egdes :: [NVector] -> [(NVector, NVector)] -egdes ps = zip ps ps' - where - ps' = tail ps ++ [head ps] - --- | Computes the geographic mean 'Position' of the given 'Position's if it is defined. --- --- The geographic mean is not defined for the antipodals positions (since they --- cancel each other). --- --- Special conditions: --- --- @ --- mean [] == Nothing --- mean [p] == Just p --- mean [p1, p2, p3] == Just circumcentre --- mean [p1, .., antipode p1] == Nothing --- @ --- -mean :: (Position a) => [a] -> Maybe a -mean [] = Nothing -mean [p] = Just p -mean ps = - if null antipodals - then Just (fromNVector (unit (foldl add zero vs))) - else Nothing - where - vs = map toNVector ps - ts = filter (\l -> length l == 2) (subsequences vs) - antipodals = - filter - (\t -> (fromNVector (antipode (head t)) :: GeoPos) == (fromNVector (last t) :: GeoPos)) - ts - --- | a, b,c => a b, a, c, b, c --- | Mean Earth radius: 6,371,008.8 metres. -meanEarthRadius :: Length -meanEarthRadius = metres 6371008.8 - --- | 'Position' of the North Pole. -northPole :: (Position a) => a -northPole = fromNVector (nvector 0.0 0.0 1.0) - --- | 'Position' of the South Pole. -southPole :: (Position a) => a -southPole = fromNVector (nvector 0.0 0.0 (-1.0)) - --- | Angle bteween the tow given 'NVector's. --- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], --- signed + if @v1@ is clockwise looking along @n@, - in opposite direction. -angleBetween :: NVector -> NVector -> Maybe NVector -> Angle -angleBetween v1 v2 n = atan2' sinO cosO - where - sign = maybe 1 (signum . dot (cross v1 v2)) n - sinO = sign * norm (cross v1 v2) - cosO = dot v1 v2 +intersections :: (Position a) => GreatCircle -> GreatCircle -> Maybe (a, a) +intersections gc1 gc2 + | norm i == 0.0 = Nothing + | otherwise + , let ni = unit i = Just (fromNVector ni, fromNVector (antipode ni)) + where + i = cross (normal gc1) (normal gc2) diff --git a/src/Data/Geo/Jord/LatLong.hs b/src/Data/Geo/Jord/LatLong.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geo/Jord/LatLong.hs @@ -0,0 +1,221 @@ +-- | +-- Module: Data.Geo.Jord.GeoPos +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois +-- Stability: experimental +-- Portability: portable +-- +-- Types to represent a geographic position by its latitude and longitude. +-- +module Data.Geo.Jord.LatLong + ( + -- * The 'LatLong' type + LatLong(latitude, longitude) + -- * Smart constructors + , latLong + , latLongE + , latLongF + , latLongDecimal + , latLongDecimalE + , latLongDecimalF + -- * read + , readLatLong + , readLatLongE + , readLatLongF + -- * Misc. + , toDecimalDegrees' + ) where + +import Control.Applicative hiding (many) +import Control.Monad.Fail +import Data.Char +import Data.Geo.Jord.Angle +import Data.Geo.Jord.Parse +import Data.Maybe +import Prelude hiding (fail) +import Text.ParserCombinators.ReadP +import Text.Read hiding (pfail) + +-- | A geographic position (latitude and longitude). +data LatLong = LatLong + { latitude :: Angle + , longitude :: Angle + } deriving (Eq) + +-- | See 'readLatLong'. +instance Read LatLong where + readsPrec _ = readP_to_S ll + +-- | Produced string format: d°(m')(s'')[N|S],d°(m')(s'')[E|W] - e.g. 55°36'21''N,13°0'2''E. +instance Show LatLong where + show (LatLong lat lon) = showLat lat ++ "," ++ showLon lon + +-- | 'LatLong' from given latitude and longitude. +-- 'error's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLong :: Angle -> Angle -> LatLong +latLong lat lon = + fromMaybe + (error ("Invalid latitude=" ++ show lat ++ " or longitude=" ++ show lon)) + (latLongF lat lon) + +-- | 'LatLong' from given latitude and longitude. +-- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLongE :: Angle -> Angle -> Either String LatLong +latLongE lat lon + | not (isWithin lat (decimalDegrees (-90)) (decimalDegrees 90)) = + Left ("Invalid latitude=" ++ show lat) + | not (isWithin lon (decimalDegrees (-180)) (decimalDegrees 180)) = + Left ("Invalid longitude=" ++ show lon) + | otherwise = Right (LatLong lat lon) + +-- | 'LatLong' from given latitude and longitude. +-- 'fail's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLongF :: (MonadFail m) => Angle -> Angle -> m LatLong +latLongF lat lon = + case e of + Left err -> fail err + Right g -> return g + where + e = latLongE lat lon + +-- | 'LatLong' from given latitude and longitude in decimal degrees. +-- 'error's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLongDecimal :: Double -> Double -> LatLong +latLongDecimal lat lon = latLong (decimalDegrees lat) (decimalDegrees lon) + +-- | 'LatLong' from given latitude and longitude in decimal degrees. +-- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLongDecimalE :: Double -> Double -> Either String LatLong +latLongDecimalE lat lon = latLongE (decimalDegrees lat) (decimalDegrees lon) + +-- | 'LatLong' from given latitude and longitude in decimal degrees. +-- 'fail's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +latLongDecimalF :: (MonadFail m) => Double -> Double -> m LatLong +latLongDecimalF lat lon = latLongF (decimalDegrees lat) (decimalDegrees lon) + +-- | Obtains a 'LatLong' from the given string formatted as either: +-- +-- * DD(MM)(SS)[N|S]DDD(MM)(SS)[E|W] - e.g. 553621N0130002E or 0116S03649E or 47N122W +-- +-- * 'Angle'[N|S] 'Angle'[E|W] - e.g. 55°36'21''N 13°0'02''E or 11°16'S 36°49'E or 47°N 122°W +-- +-- This simply calls @read s :: GeoPos@ so 'error' should be handled at the call site. +-- +readLatLong :: String -> LatLong +readLatLong s = read s :: LatLong + +-- | Same as 'readLatLong' but returns a 'Either'. +readLatLongE :: String -> Either String LatLong +readLatLongE s = + case readMaybe s of + Nothing -> Left ("couldn't read geo pos " ++ s) + Just g -> Right g + +-- | Same as 'readLatLong' but returns a 'MonadFail'. +readLatLongF :: (MonadFail m) => String -> m LatLong +readLatLongF s = + let pg = readLatLongE s + in case pg of + Left e -> fail e + Right g -> return g + +-- | Converts the given 'LatLong' to tuple of latitude and longitude in decimal degrees. +toDecimalDegrees' :: LatLong -> (Double, Double) +toDecimalDegrees' g = (toDecimalDegrees (latitude g), toDecimalDegrees (longitude g)) + +-- | Parses and returns a 'LatLong'. +ll :: ReadP LatLong +ll = block <|> human + +-- | Parses and returns a 'LatLong' - DD(D)MMSS. +block :: ReadP LatLong +block = do + lat <- blat + lon <- blon + latLongF lat lon + +-- | Parses and returns a latitude, DDMMSS expected. +blat :: ReadP Angle +blat = do + d' <- digits 2 + (m', s') <- option (0, 0) (ms <|> m) + h <- hemisphere + if h == 'N' + then dmsF d' m' s' 0 + else dmsF (-d') m' s' 0 + +-- | Parses and returns a longitude, DDDMMSS expected. +blon :: ReadP Angle +blon = do + d' <- digits 3 + (m', s') <- option (0, 0) (ms <|> m) + m'' <- meridian + if m'' == 'E' + then dmsF d' m' s' 0 + else dmsF (-d') m' s' 0 + +-- | Parses N or S char. +hemisphere :: ReadP Char +hemisphere = char 'N' <|> char 'S' + +-- | Parses E or W char. +meridian :: ReadP Char +meridian = char 'E' <|> char 'W' + +-- | Parses minutes and seconds. +ms :: ReadP (Int, Int) +ms = do + m' <- digits 2 + s' <- digits 2 + return (m', s') + +-- | Parses minutes. +m :: ReadP (Int, Int) +m = do + m' <- digits 2 + return (m', 0) + +-- | Parses and returns a 'LatLong' from a human friendly text - see 'Angle'. +human :: ReadP LatLong +human = do + lat <- hlat + _ <- char ' ' <|> char ',' + lon <- hlon + latLongF lat lon + +-- | Parses and returns a latitude, 'Angle'N|S expected. +hlat :: ReadP Angle +hlat = do + lat <- angle + h <- hemisphere + if h == 'N' + then return lat + else return (negate' lat) + +-- | Parses and returns a longitude, 'Angle'E|W expected. +hlon :: ReadP Angle +hlon = do + lon <- angle + m' <- meridian + if m' == 'E' + then return lon + else return (negate' lon) + +-- | Latitude to string. +showLat :: Angle -> String +showLat lat + | isNegative lat = show (negate' lat) ++ "S" + | otherwise = show lat ++ "N" + +-- | Longitude to string. +showLon :: Angle -> String +showLon lon + | isNegative lon = show (negate' lon) ++ "W" + | otherwise = show lon ++ "E" diff --git a/src/Data/Geo/Jord/Position.hs b/src/Data/Geo/Jord/Position.hs new file mode 100644 --- /dev/null +++ b/src/Data/Geo/Jord/Position.hs @@ -0,0 +1,224 @@ +-- | +-- Module: Data.Geo.Jord.GreatCircle +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with positions. +-- +-- All functions are implemented using the vector-based approached described in +-- +-- +-- This module assumes a spherical earth. +-- +module Data.Geo.Jord.Position + ( + -- * The 'Position' type + Position(..) + -- * Geodetic calculations + , angularDistance + , antipode + , destination + , destination' + , distance + , distance' + , finalBearing + , initialBearing + , interpolate + , isInside + , mean + -- * Misc. + , meanEarthRadius + , northPole + , southPole + ) where + +import Data.Geo.Jord.Angle +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Quantity +import Data.List (subsequences) +import Prelude hiding (fail) + +-- | The 'Position' class defines 2 functions to convert a position to and from a 'NVector'. +-- All functions in this module first convert 'Position' to 'NVector' and any resulting 'NVector' back +-- to a 'Position'. This allows the call site to pass either 'NVector' or 'LatLong' and to get back +-- the same class instance. +class Position a where + -- | Converts a 'NVector' into 'Position' instance. + fromNVector :: NVector -> a + -- | Converts the 'Position' instance into a 'NVector'. + toNVector :: a -> NVector + +-- | 'LatLong' to/from 'NVector'. +instance Position LatLong where + fromNVector v = latLong lat lon + where + lat = atan2' (z v) (sqrt (x v * x v + y v * y v)) + lon = atan2' (y v) (x v) + toNVector g = nvector x' y' z' + where + lat = latitude g + lon = longitude g + cl = cos' lat + x' = cl * cos' lon + y' = cl * sin' lon + z' = sin' lat + +-- | Identity. +instance Position NVector where + fromNVector v = v + toNVector v = v + +-- | Angle between the two given 'NVector's. +-- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], +-- signed + if @v1@ is clockwise looking along @n@, - in opposite direction. +angularDistance :: NVector -> NVector -> Maybe NVector -> Angle +angularDistance v1 v2 n = atan2' sinO cosO + where + sign = maybe 1 (signum . dot (cross v1 v2)) n + sinO = sign * norm (cross v1 v2) + cosO = dot v1 v2 + +-- | Returns the antipodal 'Position' of the given 'Position' - i.e. the position on the surface +-- of the Earth which is diametrically opposite to the given position. +antipode :: (Position a) => a -> a +antipode p = fromNVector (scale (toNVector p) (-1.0)) + +-- | 'destination'' assuming a radius of 'meanEarthRadius'. +destination :: (Position a) => a -> Angle -> Length -> a +destination p b d = destination' p b d meanEarthRadius + +-- | Computes the destination 'Position' from the given 'Position' having travelled the given distance on the +-- given initial bearing (bearing will normally vary before destination is reached) and using the given earth radius. +-- +-- This is known as the direct geodetic problem. +destination' :: (Position a) => a -> Angle -> Length -> Length -> a +destination' p b d r + | isZero d = p + | otherwise = fromNVector (add (scale v (cos' ta)) (scale de (sin' ta))) + where + v = toNVector p + ed = unit (cross northPole v) -- east direction vector at v + nd = cross v ed -- north direction vector at v + ta = central d r -- central angle + de = add (scale nd (cos' b)) (scale ed (sin' b)) -- unit vector in the direction of the azimuth + +-- | 'distance'' assuming a radius of 'meanEarthRadius'. +distance :: (Position a) => a -> a -> Length +distance p1 p2 = distance' p1 p2 meanEarthRadius + +-- | Computes the surface distance (length of geodesic) in 'Meters' assuming a +-- spherical Earth between the two given 'Position's and using the given earth radius. +distance' :: (Position a) => a -> a -> Length -> Length +distance' p1 p2 = arcLength (angularDistance v1 v2 Nothing) + where + v1 = toNVector p1 + v2 = toNVector p2 + +-- | Computes the final bearing arriving at given destination @p2@ 'Position' from given 'Position' @p1@. +-- the final bearing will differ from the 'initialBearing' by varying degrees according to distance and latitude. +-- Returns 180 if both position are equals. +finalBearing :: (Position a) => a -> a -> Angle +finalBearing p1 p2 = normalise (initialBearing p2 p1) (decimalDegrees 180) + +-- | Computes the initial bearing from given @p1@ 'Position' to given @p2@ 'Position', in compass degrees. +-- Returns 0 if both position are equals. +initialBearing :: (Position a) => a -> a -> Angle +initialBearing p1 p2 = normalise (angularDistance gc1 gc2 (Just v1)) (decimalDegrees 360) + where + v1 = toNVector p1 + v2 = toNVector p2 + gc1 = cross v1 v2 -- great circle through p1 & p2 + gc2 = cross v1 northPole -- great circle through p1 & north pole + +-- | Computes the 'Position' at given fraction @f@ between the two given 'Position's @p0@ and @p1@. +-- +-- Special conditions: +-- +-- @ +-- interpolate p0 p1 0.0 => p0 +-- interpolate p0 p1 1.0 => p1 +-- @ +-- +-- 'error's if @f < 0 || f > 1.0@ +-- +interpolate :: (Position a) => a -> a -> Double -> a +interpolate p0 p1 f + | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f) + | f == 0 = p0 + | f == 1 = p1 + | otherwise = fromNVector (unit (add v0 (scale (sub v1 v0) f))) + where + v0 = toNVector p0 + v1 = toNVector p1 + +-- | Determines whether the given 'Position' is inside the polygon defined by the given list of 'Position's. +-- The polygon is closed if needed (i.e. if @head ps /= last ps@). +-- +-- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles +-- will sum to less than 360°, and exterior point angles will be small but non-zero. +-- +-- Always returns 'False' if positions does not at least defines a triangle. +-- +isInside :: (Eq a, Position a) => a -> [a] -> Bool +isInside p ps + | null ps = False + | head ps == last ps = isInside p (init ps) + | length ps < 3 = False + | otherwise = + let aSum = foldl (\a v' -> add a (uncurry angularDistance v' (Just v))) (decimalDegrees 0) es + in abs (toDecimalDegrees aSum) > 180.0 + where + v = toNVector p + es = egdes (map (sub v . toNVector) ps) + +-- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] +egdes :: [NVector] -> [(NVector, NVector)] +egdes ps = zip ps ps' + where + ps' = tail ps ++ [head ps] + +-- | Computes the geographic mean 'Position' of the given 'Position's if it is defined. +-- +-- The geographic mean is not defined for the antipodals positions (since they +-- cancel each other). +-- +-- Special conditions: +-- +-- @ +-- mean [] == Nothing +-- mean [p] == Just p +-- mean [p1, p2, p3] == Just circumcentre +-- mean [p1, .., antipode p1] == Nothing +-- @ +-- +mean :: (Position a) => [a] -> Maybe a +mean [] = Nothing +mean [p] = Just p +mean ps = + if null antipodals + then Just (fromNVector (unit (foldl add zero vs))) + else Nothing + where + vs = map toNVector ps + ts = filter (\l -> length l == 2) (subsequences vs) + antipodals = + filter + (\t -> (fromNVector (antipode (head t)) :: LatLong) == (fromNVector (last t) :: LatLong)) + ts + +-- | Mean Earth radius: 6,371,008.8 metres. +meanEarthRadius :: Length +meanEarthRadius = metres 6371008.8 + +-- | 'Position' of the North Pole. +northPole :: (Position a) => a +northPole = fromNVector (nvector 0.0 0.0 1.0) + +-- | 'Position' of the South Pole. +southPole :: (Position a) => a +southPole = fromNVector (nvector 0.0 0.0 (-1.0)) diff --git a/test/Data/Geo/Jord/AngleSpec.hs b/test/Data/Geo/Jord/AngleSpec.hs --- a/test/Data/Geo/Jord/AngleSpec.hs +++ b/test/Data/Geo/Jord/AngleSpec.hs @@ -1,84 +1,89 @@ -module Data.Geo.Jord.AngleSpec - ( spec - ) where - -import Data.Geo.Jord -import Test.Hspec - -spec :: Spec -spec = do - describe "Reading valid angles" $ do - it "reads 55°36'21\"" $ readAngle "55°36'21\"" `shouldBe` decimalDegrees 55.6058333 - it "reads 55°36'21''" $ readAngle "55°36'21''" `shouldBe` decimalDegrees 55.6058333 - it "reads 55d36m21.0s" $ readAngle "55d36m21.0s" `shouldBe` decimalDegrees 55.6058333 - it "reads 55.6058333°" $ readAngle "55.6058333°" `shouldBe` decimalDegrees 55.6058333 - it "reads -55.6058333°" $ readAngle "-55.6058333°" `shouldBe` decimalDegrees (-55.6058333) - it "reads 96°01′18″" $ readAngle "96°01′18″" `shouldBe` decimalDegrees 96.02166666 - describe "Adding/Subtracting angles" $ do - it "adds angles" $ - add (decimalDegrees 55.6058333) (decimalDegrees 5.0) `shouldBe` - decimalDegrees 60.6058333 - it "subtracts angles" $ - sub (decimalDegrees 5.0) (decimalDegrees 55.6058333) `shouldBe` - decimalDegrees (-50.6058333) - describe "Angle normalisation" $ do - it "370 degrees normalised to [0..360] = 10" $ - normalise (decimalDegrees 370) (decimalDegrees 360) `shouldBe` decimalDegrees 10 - it "350 degrees normalised to [0..360] = 350" $ - normalise (decimalDegrees 350) (decimalDegrees 360) `shouldBe` decimalDegrees 350 - describe "Angle equality" $ do - it "considers 59.9999999° == 60.0°" $ decimalDegrees 59.9999999 `shouldBe` decimalDegrees 60 - it "considers 59.9999998° /= 60.0°" $ - decimalDegrees 59.9999998 `shouldNotBe` decimalDegrees 60 - describe "Showing angles" $ do - it "shows 59.99999999999999 as 60°0'0.0\"" $ - show (decimalDegrees 59.99999999999999) `shouldBe` "60°0'0.0\"" - it "shows 154.915 as 154°54'54.0\"" $ - show (decimalDegrees 154.915) `shouldBe` "154°54'54.0\"" - it "shows -154.915 as -154°54'54.0\"" $ - show (decimalDegrees (-154.915)) `shouldBe` "-154°54'54.0\"" - describe "Angle from decimal degrees" $ do - it "returns 1 millisecond when called with 1 / 3600000" $ do - let actual = decimalDegrees (1 / 3600000) - getDegrees actual `shouldBe` 0 - getMinutes actual `shouldBe` 0 - getSeconds actual `shouldBe` 0 - getMilliseconds actual `shouldBe` 1 - it "returns 1 second when called with 1000 / 3600000" $ do - let actual = decimalDegrees (1000 / 3600000) - getDegrees actual `shouldBe` 0 - getMinutes actual `shouldBe` 0 - getSeconds actual `shouldBe` 1 - getMilliseconds actual `shouldBe` 0 - it "returns 1 minute when called with 60000 / 3600000" $ do - let actual = decimalDegrees (60000 / 3600000) - getDegrees actual `shouldBe` 0 - getMinutes actual `shouldBe` 1 - getSeconds actual `shouldBe` 0 - getMilliseconds actual `shouldBe` 0 - it "returns 1 degree when called with 1" $ do - let actual = decimalDegrees 1 - getDegrees actual `shouldBe` 1 - getMinutes actual `shouldBe` 0 - getSeconds actual `shouldBe` 0 - getMilliseconds actual `shouldBe` 0 - it "accepts positve values" $ do - let actual = decimalDegrees 154.9150300 - getDegrees actual `shouldBe` 154 - getMinutes actual `shouldBe` 54 - getSeconds actual `shouldBe` 54 - getMilliseconds actual `shouldBe` 108 - it "accepts negative values" $ do - let actual = decimalDegrees (-154.915) - getDegrees actual `shouldBe` (-154) - getMinutes actual `shouldBe` 54 - getSeconds actual `shouldBe` 54 - getMilliseconds actual `shouldBe` 0 - describe "Arc length" $ do - it "computes the length of an arc with a central angle of 1 milliseconds" $ - arcLength (decimalDegrees (1.0 / 3600000.0)) meanEarthRadius `shouldBe` metres 0.031 - it - "arc length with central angle of 0.6 milliseconds == arc length with central angle of 1 milliseconds" $ - arcLength (decimalDegrees (0.6 / 3600000.0)) meanEarthRadius `shouldBe` metres 0.031 - it "arc length with central angle of 0.5 milliseconds == 0" $ - arcLength (decimalDegrees (0.4 / 3600000.0)) meanEarthRadius `shouldBe` metres 0 +module Data.Geo.Jord.AngleSpec + ( spec + ) where + +import Data.Geo.Jord +import System.IO +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid angles" $ do + it "reads 55°36'21\"" $ readAngle "55°36'21\"" `shouldBe` decimalDegrees 55.6058333 + it "reads 55°36'21''" $ readAngle "55°36'21''" `shouldBe` decimalDegrees 55.6058333 + it "reads 55d36m21.0s" $ readAngle "55d36m21.0s" `shouldBe` decimalDegrees 55.6058333 + it "reads 55.6058333°" $ readAngle "55.6058333°" `shouldBe` decimalDegrees 55.6058333 + it "reads -55.6058333°" $ readAngle "-55.6058333°" `shouldBe` decimalDegrees (-55.6058333) + it "reads 96°01′18″" $ do + hSetEncoding stdin utf8 + hSetEncoding stdout utf8 + hSetEncoding stderr utf8 + readAngle "96°01′18″" `shouldBe` decimalDegrees 96.02166666 + describe "Adding/Subtracting angles" $ do + it "adds angles" $ + add (decimalDegrees 55.6058333) (decimalDegrees 5.0) `shouldBe` + decimalDegrees 60.6058333 + it "subtracts angles" $ + sub (decimalDegrees 5.0) (decimalDegrees 55.6058333) `shouldBe` + decimalDegrees (-50.6058333) + describe "Angle normalisation" $ do + it "370 degrees normalised to [0..360] = 10" $ + normalise (decimalDegrees 370) (decimalDegrees 360) `shouldBe` decimalDegrees 10 + it "350 degrees normalised to [0..360] = 350" $ + normalise (decimalDegrees 350) (decimalDegrees 360) `shouldBe` decimalDegrees 350 + describe "Angle equality" $ do + it "considers 59.9999999° == 60.0°" $ decimalDegrees 59.9999999 `shouldBe` decimalDegrees 60 + it "considers 59.9999998° /= 60.0°" $ + decimalDegrees 59.9999998 `shouldNotBe` decimalDegrees 60 + describe "Showing angles" $ do + it "shows 59.99999999999999 as 60°0'0.0\"" $ + show (decimalDegrees 59.99999999999999) `shouldBe` "60°0'0.0\"" + it "shows 154.915 as 154°54'54.0\"" $ + show (decimalDegrees 154.915) `shouldBe` "154°54'54.0\"" + it "shows -154.915 as -154°54'54.0\"" $ + show (decimalDegrees (-154.915)) `shouldBe` "-154°54'54.0\"" + describe "Angle from decimal degrees" $ do + it "returns 1 millisecond when called with 1 / 3600000" $ do + let actual = decimalDegrees (1 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 1 + it "returns 1 second when called with 1000 / 3600000" $ do + let actual = decimalDegrees (1000 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 1 + getMilliseconds actual `shouldBe` 0 + it "returns 1 minute when called with 60000 / 3600000" $ do + let actual = decimalDegrees (60000 / 3600000) + getDegrees actual `shouldBe` 0 + getMinutes actual `shouldBe` 1 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 0 + it "returns 1 degree when called with 1" $ do + let actual = decimalDegrees 1 + getDegrees actual `shouldBe` 1 + getMinutes actual `shouldBe` 0 + getSeconds actual `shouldBe` 0 + getMilliseconds actual `shouldBe` 0 + it "accepts positve values" $ do + let actual = decimalDegrees 154.9150300 + getDegrees actual `shouldBe` 154 + getMinutes actual `shouldBe` 54 + getSeconds actual `shouldBe` 54 + getMilliseconds actual `shouldBe` 108 + it "accepts negative values" $ do + let actual = decimalDegrees (-154.915) + getDegrees actual `shouldBe` (-154) + getMinutes actual `shouldBe` 54 + getSeconds actual `shouldBe` 54 + getMilliseconds actual `shouldBe` 0 + describe "Arc length" $ do + it "computes the length of an arc with a central angle of 1 milliseconds" $ + arcLength (decimalDegrees (1.0 / 3600000.0)) meanEarthRadius `shouldBe` metres 0.031 + it + "arc length with central angle of 0.6 milliseconds == arc length with central angle of 1 milliseconds" $ + arcLength (decimalDegrees (0.6 / 3600000.0)) meanEarthRadius `shouldBe` metres 0.031 + it "arc length with central angle of 0.5 milliseconds == 0" $ + arcLength (decimalDegrees (0.4 / 3600000.0)) meanEarthRadius `shouldBe` metres 0 diff --git a/test/Data/Geo/Jord/EvalSpec.hs b/test/Data/Geo/Jord/EvalSpec.hs --- a/test/Data/Geo/Jord/EvalSpec.hs +++ b/test/Data/Geo/Jord/EvalSpec.hs @@ -1,41 +1,41 @@ -module Data.Geo.Jord.EvalSpec - ( spec - ) where - -import Data.Geo.Jord -import Test.Hspec - -spec :: Spec -spec = do - describe "Expression evaluation" $ do - it "evaluates simple expression" $ - case eval "antipode 54N154E" emptyVault of - (Right (Geo g)) -> g `shouldBe` latLongDecimal (-54.0) (-26.0) - r -> fail (show r) - it "evaluates an expression with one function call" $ - case eval "distance 54N154E (antipode 54N154E)" emptyVault of - (Right (Len l)) -> l `shouldBe` kilometres 20015.114442000002 - r -> fail (show r) - it "evaluates an expression with one function call" $ - case eval "distance (antipode 54N154E) 54N154E" emptyVault of - (Right (Len l)) -> l `shouldBe` kilometres 20015.114442000002 - r -> fail (show r) - it "evaluates expression with nested function calls" $ - case eval - "finalBearing (destination (antipode 54°N,154°E) 54° 1000m) (readGeoPos 54°N,154°E)" - emptyVault of - (Right (Ang a)) -> a `shouldBe` decimalDegrees 126 - r -> fail (show r) - it "resolves variables" $ do - let vault = insert "a" (Geo (latLongDecimal 54.0 154.0)) emptyVault - case eval "antipode a" vault of - (Right (Geo g)) -> g `shouldBe` latLongDecimal (-54.0) (-26.0) - r -> fail (show r) - it "rejects expression with lexical error" $ - case eval "finalBearing(destination" emptyVault of - (Left e) -> e `shouldBe` "Lexical error: finalBearing(destination" - r -> fail (show r) - it "rejects expression with syntaxic error" $ - case eval "finalBearing (destination a" emptyVault of - (Left e) -> e `shouldBe` "Syntax error: ')' not found" - r -> fail (show r) +module Data.Geo.Jord.EvalSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Expression evaluation" $ do + it "evaluates simple expression" $ + case eval "antipode 54N154E" emptyVault of + (Right (Ll ll)) -> ll `shouldBe` latLongDecimal (-54.0) (-26.0) + r -> fail (show r) + it "evaluates an expression with one function call" $ + case eval "distance 54N154E (antipode 54N154E)" emptyVault of + (Right (Len l)) -> l `shouldBe` kilometres 20015.114442000002 + r -> fail (show r) + it "evaluates an expression with one function call" $ + case eval "distance (antipode 54N154E) 54N154E" emptyVault of + (Right (Len l)) -> l `shouldBe` kilometres 20015.114442000002 + r -> fail (show r) + it "evaluates expression with nested function calls" $ + case eval + "finalBearing (destination (antipode 54°N,154°E) 54° 1000m) (readLatLong 54°N,154°E)" + emptyVault of + (Right (Ang a)) -> a `shouldBe` decimalDegrees 126 + r -> fail (show r) + it "resolves variables" $ do + let vault = insert "a" (Ll (latLongDecimal 54.0 154.0)) emptyVault + case eval "antipode a" vault of + (Right (Ll ll)) -> ll `shouldBe` latLongDecimal (-54.0) (-26.0) + r -> fail (show r) + it "rejects expression with lexical error" $ + case eval "finalBearing(destination" emptyVault of + (Left e) -> e `shouldBe` "Lexical error: finalBearing(destination" + r -> fail (show r) + it "rejects expression with syntaxic error" $ + case eval "finalBearing (destination a" emptyVault of + (Left e) -> e `shouldBe` "Syntax error: ')' not found" + r -> fail (show r) diff --git a/test/Data/Geo/Jord/GeoPosSpec.hs b/test/Data/Geo/Jord/GeoPosSpec.hs deleted file mode 100644 --- a/test/Data/Geo/Jord/GeoPosSpec.hs +++ /dev/null @@ -1,63 +0,0 @@ -module Data.Geo.Jord.GeoPosSpec - ( spec - ) where - -import Data.Geo.Jord -import Test.Hspec - -spec :: Spec -spec = do - describe "Reading valid DMS text" $ do - it "reads 553621N0130002E" $ - readGeoPos "553621N0130002E" `shouldBe` latLongDecimal 55.6058333 13.0005555 - it "reads 55°36'21''N 013°00'02''E" $ - readGeoPos "55°36'21''N 013°00'02''E" `shouldBe` latLongDecimal 55.6058333 13.0005555 - it "reads 5536N01300E" $ readGeoPos "5536N01300E" `shouldBe` latLongDecimal 55.6 13.0 - it "reads 55N013E" $ readGeoPos "55N013E" `shouldBe` latLongDecimal 55.0 13.0 - it "reads 011659S0364900E" $ - readGeoPos "011659S0364900E" `shouldBe` latLongDecimal (-1.2830555) 36.8166666 - it "reads 0116S03649E" $ - readGeoPos "0116S03649E" `shouldBe` latLongDecimal (-1.2666666) 36.8166666 - it "reads 1°16'S,36°49'E" $ - readGeoPos "1°16'S,36°49'E" `shouldBe` latLongDecimal (-1.2666666) 36.8166666 - it "reads 01S036E" $ readGeoPos "01S036E" `shouldBe` latLongDecimal (-1.0) 36.0 - it "reads 473622N1221955W" $ - readGeoPos "473622N1221955W" `shouldBe` latLongDecimal 47.6061111 (-122.3319444) - it "reads 4736N12219W" $ - readGeoPos "4736N12219W" `shouldBe` latLongDecimal 47.6 (-122.3166666) - it "reads 47N122W" $ readGeoPos "47N122W" `shouldBe` latLongDecimal 47.0 (-122.0) - it "reads 47°N 122°W" $ readGeoPos "47°N 122°W" `shouldBe` latLongDecimal 47.0 (-122.0) - it "reads 544807S0681811W" $ - readGeoPos "544807S0681811W" `shouldBe` latLongDecimal (-54.8019444) (-68.3030555) - it "reads 5448S06818W" $ readGeoPos "5448S06818W" `shouldBe` latLongDecimal (-54.8) (-68.3) - it "reads 54S068W" $ readGeoPos "54S068W" `shouldBe` latLongDecimal (-54.0) (-68.0) - describe "Reading invalid DMS text" $ do - it "fails to read 553621K0130002E" $ - readGeoPosE "553621K0130002E" `shouldBe` Left "couldn't read geo pos 553621K0130002E" - it "fails to read 011659S0364900Z" $ - readGeoPosE "011659S0364900Z" `shouldBe` Left "couldn't read geo pos 011659S0364900Z" - it "fails to read 4736221221955W" $ - readGeoPosE "4736221221955W" `shouldBe` Left "couldn't read geo pos 4736221221955W" - it "fails to read 54480S0681811W" $ - readGeoPosE "54480S0681811W" `shouldBe` Left "couldn't read geo pos 54480S0681811W" - it "fails to read 553621N013000E" $ - readGeoPosE "553621N013000E" `shouldBe` Left "couldn't read geo pos 553621N013000E" - it "fails to read 914807S0681811W" $ - readGeoPosE "914807S0681811W" `shouldBe` Left "couldn't read geo pos 914807S0681811W" - it "fails to read 544807S1811811W" $ - readGeoPosE "544807S1811811W" `shouldBe` Left "couldn't read geo pos 544807S1811811W" - it "fails to read 546007S1801811W" $ - readGeoPosE "546007S1801811W" `shouldBe` Left "couldn't read geo pos 546007S1801811W" - it "fails to read 545907S1801860W" $ - readGeoPosE "545907S1801860W" `shouldBe` Left "couldn't read geo pos 545907S1801860W" - describe "Showing geographic positions" $ do - it "shows the N/E position formatted in DMS with symbols" $ - show (latLongDecimal 55.60583333 13.00055556) `shouldBe` "55°36'21.0\"N,13°0'2.0\"E" - it "shows the S/E position formatted in DMS with symbols" $ - show (latLongDecimal (-1.28305556) 36.81666) `shouldBe` "1°16'59.0\"S,36°48'59.976\"E" - it "shows the N/W position formatted in DMS with symbols" $ - show (latLongDecimal 47.60611 (-122.33194)) `shouldBe` - "47°36'21.996\"N,122°19'54.984\"W" - it "shows the S/W position formatted in DMS with symbols" $ - show (latLongDecimal (-54.80194) (-68.30305)) `shouldBe` - "54°48'6.984\"S,68°18'10.980\"W" diff --git a/test/Data/Geo/Jord/GreatCircleSpec.hs b/test/Data/Geo/Jord/GreatCircleSpec.hs --- a/test/Data/Geo/Jord/GreatCircleSpec.hs +++ b/test/Data/Geo/Jord/GreatCircleSpec.hs @@ -1,192 +1,48 @@ -module Data.Geo.Jord.GreatCircleSpec - ( spec - ) where - -import Control.Exception.Base -import Data.Geo.Jord -import Data.Maybe -import Test.Hspec - -spec :: Spec -spec = do - describe "Antipode" $ do - it "returns the antipodal point" $ - antipode (readGeoPos "484137N0061105E") `shouldBe` - latLongDecimal (-48.6936111) (-173.8152777) - it "returns the south pole when called with the north pole" $ - antipode (northPole :: GeoPos) `shouldBe` latLongDecimal (-90.0) (-180.0) - it "returns the north pole when called with the south pole" $ - antipode (southPole :: GeoPos) `shouldBe` latLongDecimal 90.0 (-180.0) - describe "Cross Track Distance" $ do - it "returns a negative length when position is left of great circle (bearing)" $ do - let p = latLongDecimal 53.2611 (-0.7972) - let gc = greatCircleBearing (latLongDecimal 53.3206 (-1.7297)) (decimalDegrees 96.0) - crossTrackDistance p gc `shouldBe` metres (-305.663) - it "returns a negative length when position is left of great circle" $ do - let p = latLongDecimal 53.2611 (-0.7972) - let gc = greatCircle (latLongDecimal 53.3206 (-1.7297)) (latLongDecimal 53.1887 0.1334) - crossTrackDistance p gc `shouldBe` metres (-307.547) - it "returns a positve length when position is right of great circle (bearing)" $ do - let p = readGeoPos "531540N0014750W" - let gc = greatCircleBearing (readGeoPos "531914N0014347W") (readAngle "96d01m18s") - crossTrackDistance p gc `shouldBe` metres 7042.324 - it "returns a positive length when position is left of great circle" $ do - let p = antipode (latLongDecimal 53.2611 (-0.7972)) - let gc = greatCircle (latLongDecimal 53.3206 (-1.7297)) (latLongDecimal 53.1887 0.1334) - crossTrackDistance p gc `shouldBe` metres 307.547 - describe "Distance" $ do - it "returns 0 if both points are equal" $ - distance (readGeoPos "500359N1795959W") (readGeoPos "500359N1795959W") `shouldBe` - metres 0.0 - it "returns the distance between 2 points" $ - distance (readGeoPos "500359N0054253W") (readGeoPos "583838N0030412W") `shouldBe` - metres 968854.873 - it "handles singularity at the pole" $ - distance (northPole :: GeoPos) (southPole :: GeoPos) `shouldBe` - metres 2.00151144420359e7 - it "handles the discontinuity at the Date Line" $ - distance (readGeoPos "500359N1795959W") (readGeoPos "500359N1795959E") `shouldBe` - metres 39.66 - describe "Destination" $ do - it "return the given point if distance is 0 meter" $ - destination (readGeoPos "531914N0014347W") (decimalDegrees 96.0217) (metres 0) `shouldBe` - readGeoPos "531914N0014347W" - it "return the destination point along great-circle at distance and bearing" $ - destination (readGeoPos "531914N0014347W") (decimalDegrees 96.0217) (metres 124800) `shouldBe` - latLongDecimal 53.1882691 0.1332744 - describe "Initial bearing" $ do - it "returns the 0 if both point are the same" $ - initialBearing (readGeoPos "500359N0054253W") (readGeoPos "500359N0054253W") `shouldBe` - decimalDegrees 0 - it "returns the initial bearing in compass degrees" $ - initialBearing (readGeoPos "500359N0054253W") (readGeoPos "583838N0030412W") `shouldBe` - decimalDegrees 9.1198181 - it "returns the initial bearing in compass degrees" $ - initialBearing (readGeoPos "583838N0030412W") (readGeoPos "500359N0054253W") `shouldBe` - decimalDegrees 191.2752013 - describe "Interpolate" $ do - it "fails if f < 0.0" $ - evaluate (interpolate (readGeoPos "44N044E") (readGeoPos "46N046E") (-0.5)) `shouldThrow` - errorCall "fraction must be in range [0..1], was -0.5" - it "fails if f > 1.0" $ - evaluate (interpolate (readGeoPos "44N044E") (readGeoPos "46N046E") 1.1) `shouldThrow` - errorCall "fraction must be in range [0..1], was 1.1" - it "returns p0 if f == 0" $ - interpolate (readGeoPos "44N044E") (readGeoPos "46N046E") 0.0 `shouldBe` - readGeoPos "44N044E" - it "returns p1 if f == 1" $ - interpolate (readGeoPos "44N044E") (readGeoPos "46N046E") 1.0 `shouldBe` - readGeoPos "46N046E" - it "returns the interpolated position" $ - interpolate - (readGeoPos "53°28'46''N 2°14'43''W") - (readGeoPos "55°36'21''N 13°02'09''E") - 0.5 `shouldBe` - latLongDecimal 54.7835574 5.1949856 - describe "Intersections" $ do - it "returns nothing if both great circle are equals" $ do - let gc = greatCircleBearing (latLongDecimal 51.885 0.235) (decimalDegrees 108.63) - (intersections gc gc :: Maybe (GeoPos, GeoPos)) `shouldBe` Nothing - it "returns nothing if both great circle are equals (opposite orientation)" $ do - let gc1 = greatCircle (latLongDecimal 51.885 0.235) (latLongDecimal 52.885 1.235) - let gc2 = greatCircle (latLongDecimal 52.885 1.235) (latLongDecimal 51.885 0.235) - (intersections gc1 gc2 :: Maybe (GeoPos, GeoPos)) `shouldBe` Nothing - it "returns the two points where the two great circles intersects" $ do - let gc1 = greatCircleBearing (latLongDecimal 51.885 0.235) (decimalDegrees 108.63) - let gc2 = greatCircleBearing (latLongDecimal 49.008 2.549) (decimalDegrees 32.72) - let (i1, i2) = fromJust (intersections gc1 gc2) - i1 `shouldBe` latLongDecimal 50.9017226 4.4942782 - i2 `shouldBe` antipode i1 - describe "isInside" $ do - it "return False if polygon is empty" $ isInside (latLongDecimal 45 1) [] `shouldBe` False - it "return False if polygon does not define at least a triangle" $ - isInside (latLongDecimal 45 1) [latLongDecimal 45 1, latLongDecimal 45 2] `shouldBe` - False - it "returns True if point is inside polygon" $ do - let polygon = - [ latLongDecimal 45 1 - , latLongDecimal 45 2 - , latLongDecimal 46 2 - , latLongDecimal 46 1 - ] - let p = latLongDecimal 45.1 1.1 - isInside p polygon `shouldBe` True - it "returns False if point is inside polygon" $ do - let polygon = - [ latLongDecimal 45 1 - , latLongDecimal 45 2 - , latLongDecimal 46 2 - , latLongDecimal 46 1 - ] - let p = antipode (latLongDecimal 45.1 1.1) - isInside p polygon `shouldBe` False - it "returns False if point is a vertex of the polygon" $ do - let polygon = - [ latLongDecimal 45 1 - , latLongDecimal 45 2 - , latLongDecimal 46 2 - , latLongDecimal 46 1 - ] - let p = latLongDecimal 45 1 - isInside p polygon `shouldBe` False - it "handles closed polygons" $ do - let polygon = - [ latLongDecimal 45 1 - , latLongDecimal 45 2 - , latLongDecimal 46 2 - , latLongDecimal 46 1 - , latLongDecimal 45 1 - ] - let p = latLongDecimal 45.1 1.1 - isInside p polygon `shouldBe` True - it "handles concave polygons" $ do - let malmo = latLongDecimal 55.6050 13.0038 - let ystad = latLongDecimal 55.4295 13.82 - let lund = latLongDecimal 55.7047 13.1910 - let helsingborg = latLongDecimal 56.0465 12.6945 - let kristianstad = latLongDecimal 56.0294 14.1567 - let polygon = [malmo, ystad, kristianstad, helsingborg, lund] - let hoor = latLongDecimal 55.9295 13.5297 - let hassleholm = latLongDecimal 56.1589 13.7668 - isInside hoor polygon `shouldBe` True - isInside hassleholm polygon `shouldBe` False - describe "Final bearing" $ do - it "returns the 180.0 if both point are the same" $ - finalBearing (readGeoPos "500359N0054253W") (readGeoPos "500359N0054253W") `shouldBe` - decimalDegrees 180 - it "returns the final bearing in compass degrees" $ - finalBearing (readGeoPos "500359N0054253W") (readGeoPos "583838N0030412W") `shouldBe` - decimalDegrees 11.2752013 - it "returns the final bearing in compass degrees" $ - finalBearing (readGeoPos "583838N0030412W") (readGeoPos "500359N0054253W") `shouldBe` - decimalDegrees 189.1198181 - it "returns the final bearing in compass degrees" $ - finalBearing (readGeoPos "535941S0255915W") (readGeoPos "54N154E") `shouldBe` - decimalDegrees 125.6839436 - describe "Great Circle Smart constructors" $ do - it "fails if both positions are equal" $ - greatCircleE (latLongDecimal 3 154) (latLongDecimal 3 154) `shouldBe` - Left "Invalid Great Circle: positions are equal" - it "fails if both positions are antipodal" $ - greatCircleE (latLongDecimal 3 154) (antipode (latLongDecimal 3 154)) `shouldBe` - Left "Invalid Great Circle: positions are antipodal" - describe "Mean" $ do - it "returns Nothing if no point is given" $ (mean [] :: Maybe GeoPos) `shouldBe` Nothing - it "returns the unique given point" $ - mean [readGeoPos "500359N0054253W"] `shouldBe` Just (readGeoPos "500359N0054253W") - it "returns the geographical mean" $ - mean [readGeoPos "500359N0054253W", readGeoPos "583838N0030412W"] `shouldBe` - Just (latLongDecimal 54.3622869 (-4.5306725)) - it "returns Nothing if list contains antipodal points" $ do - let points = - [ latLongDecimal 45 1 - , latLongDecimal 45 2 - , latLongDecimal 46 2 - , latLongDecimal 46 1 - , antipode (latLongDecimal 45 2) - ] - mean points `shouldBe` Nothing - describe "North pole" $ - it "returns (90, 0)" $ (northPole :: GeoPos) `shouldBe` latLongDecimal 90.0 0.0 - describe "South pole" $ - it "returns (-90, 0)" $ (southPole :: GeoPos) `shouldBe` latLongDecimal (-90.0) 0.0 +module Data.Geo.Jord.GreatCircleSpec + ( spec + ) where + +import Data.Geo.Jord +import Data.Maybe +import Test.Hspec + +spec :: Spec +spec = do + describe "Cross Track Distance" $ do + it "returns a negative length when position is left of great circle (bearing)" $ do + let p = latLongDecimal 53.2611 (-0.7972) + let gc = greatCircleBearing (latLongDecimal 53.3206 (-1.7297)) (decimalDegrees 96.0) + crossTrackDistance p gc `shouldBe` metres (-305.663) + it "returns a negative length when position is left of great circle" $ do + let p = latLongDecimal 53.2611 (-0.7972) + let gc = greatCircle (latLongDecimal 53.3206 (-1.7297)) (latLongDecimal 53.1887 0.1334) + crossTrackDistance p gc `shouldBe` metres (-307.547) + it "returns a positve length when position is right of great circle (bearing)" $ do + let p = readLatLong "531540N0014750W" + let gc = greatCircleBearing (readLatLong "531914N0014347W") (readAngle "96d01m18s") + crossTrackDistance p gc `shouldBe` metres 7042.324 + it "returns a positive length when position is left of great circle" $ do + let p = antipode (latLongDecimal 53.2611 (-0.7972)) + let gc = greatCircle (latLongDecimal 53.3206 (-1.7297)) (latLongDecimal 53.1887 0.1334) + crossTrackDistance p gc `shouldBe` metres 307.547 + describe "Intersections" $ do + it "returns nothing if both great circle are equals" $ do + let gc = greatCircleBearing (latLongDecimal 51.885 0.235) (decimalDegrees 108.63) + (intersections gc gc :: Maybe (LatLong, LatLong)) `shouldBe` Nothing + it "returns nothing if both great circle are equals (opposite orientation)" $ do + let gc1 = greatCircle (latLongDecimal 51.885 0.235) (latLongDecimal 52.885 1.235) + let gc2 = greatCircle (latLongDecimal 52.885 1.235) (latLongDecimal 51.885 0.235) + (intersections gc1 gc2 :: Maybe (LatLong, LatLong)) `shouldBe` Nothing + it "returns the two points where the two great circles intersects" $ do + let gc1 = greatCircleBearing (latLongDecimal 51.885 0.235) (decimalDegrees 108.63) + let gc2 = greatCircleBearing (latLongDecimal 49.008 2.549) (decimalDegrees 32.72) + let (i1, i2) = fromJust (intersections gc1 gc2) + i1 `shouldBe` latLongDecimal 50.9017226 4.4942782 + i2 `shouldBe` antipode i1 + describe "Great Circle Smart constructors" $ do + it "fails if both positions are equal" $ + greatCircleE (latLongDecimal 3 154) (latLongDecimal 3 154) `shouldBe` + Left "Invalid Great Circle: positions are equal" + it "fails if both positions are antipodal" $ + greatCircleE (latLongDecimal 3 154) (antipode (latLongDecimal 3 154)) `shouldBe` + Left "Invalid Great Circle: positions are antipodal" diff --git a/test/Data/Geo/Jord/LatLongSpec.hs b/test/Data/Geo/Jord/LatLongSpec.hs new file mode 100644 --- /dev/null +++ b/test/Data/Geo/Jord/LatLongSpec.hs @@ -0,0 +1,63 @@ +module Data.Geo.Jord.LatLongSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Reading valid DMS text" $ do + it "reads 553621N0130002E" $ + readLatLong "553621N0130002E" `shouldBe` latLongDecimal 55.6058333 13.0005555 + it "reads 55°36'21''N 013°00'02''E" $ + readLatLong "55°36'21''N 013°00'02''E" `shouldBe` latLongDecimal 55.6058333 13.0005555 + it "reads 5536N01300E" $ readLatLong "5536N01300E" `shouldBe` latLongDecimal 55.6 13.0 + it "reads 55N013E" $ readLatLong "55N013E" `shouldBe` latLongDecimal 55.0 13.0 + it "reads 011659S0364900E" $ + readLatLong "011659S0364900E" `shouldBe` latLongDecimal (-1.2830555) 36.8166666 + it "reads 0116S03649E" $ + readLatLong "0116S03649E" `shouldBe` latLongDecimal (-1.2666666) 36.8166666 + it "reads 1°16'S,36°49'E" $ + readLatLong "1°16'S,36°49'E" `shouldBe` latLongDecimal (-1.2666666) 36.8166666 + it "reads 01S036E" $ readLatLong "01S036E" `shouldBe` latLongDecimal (-1.0) 36.0 + it "reads 473622N1221955W" $ + readLatLong "473622N1221955W" `shouldBe` latLongDecimal 47.6061111 (-122.3319444) + it "reads 4736N12219W" $ + readLatLong "4736N12219W" `shouldBe` latLongDecimal 47.6 (-122.3166666) + it "reads 47N122W" $ readLatLong "47N122W" `shouldBe` latLongDecimal 47.0 (-122.0) + it "reads 47°N 122°W" $ readLatLong "47°N 122°W" `shouldBe` latLongDecimal 47.0 (-122.0) + it "reads 544807S0681811W" $ + readLatLong "544807S0681811W" `shouldBe` latLongDecimal (-54.8019444) (-68.3030555) + it "reads 5448S06818W" $ readLatLong "5448S06818W" `shouldBe` latLongDecimal (-54.8) (-68.3) + it "reads 54S068W" $ readLatLong "54S068W" `shouldBe` latLongDecimal (-54.0) (-68.0) + describe "Reading invalid DMS text" $ do + it "fails to read 553621K0130002E" $ + readLatLongE "553621K0130002E" `shouldBe` Left "couldn't read geo pos 553621K0130002E" + it "fails to read 011659S0364900Z" $ + readLatLongE "011659S0364900Z" `shouldBe` Left "couldn't read geo pos 011659S0364900Z" + it "fails to read 4736221221955W" $ + readLatLongE "4736221221955W" `shouldBe` Left "couldn't read geo pos 4736221221955W" + it "fails to read 54480S0681811W" $ + readLatLongE "54480S0681811W" `shouldBe` Left "couldn't read geo pos 54480S0681811W" + it "fails to read 553621N013000E" $ + readLatLongE "553621N013000E" `shouldBe` Left "couldn't read geo pos 553621N013000E" + it "fails to read 914807S0681811W" $ + readLatLongE "914807S0681811W" `shouldBe` Left "couldn't read geo pos 914807S0681811W" + it "fails to read 544807S1811811W" $ + readLatLongE "544807S1811811W" `shouldBe` Left "couldn't read geo pos 544807S1811811W" + it "fails to read 546007S1801811W" $ + readLatLongE "546007S1801811W" `shouldBe` Left "couldn't read geo pos 546007S1801811W" + it "fails to read 545907S1801860W" $ + readLatLongE "545907S1801860W" `shouldBe` Left "couldn't read geo pos 545907S1801860W" + describe "Showing geographic positions" $ do + it "shows the N/E position formatted in DMS with symbols" $ + show (latLongDecimal 55.60583333 13.00055556) `shouldBe` "55°36'21.0\"N,13°0'2.0\"E" + it "shows the S/E position formatted in DMS with symbols" $ + show (latLongDecimal (-1.28305556) 36.81666) `shouldBe` "1°16'59.0\"S,36°48'59.976\"E" + it "shows the N/W position formatted in DMS with symbols" $ + show (latLongDecimal 47.60611 (-122.33194)) `shouldBe` + "47°36'21.996\"N,122°19'54.984\"W" + it "shows the S/W position formatted in DMS with symbols" $ + show (latLongDecimal (-54.80194) (-68.30305)) `shouldBe` + "54°48'6.984\"S,68°18'10.980\"W" diff --git a/test/Data/Geo/Jord/PositionSpec.hs b/test/Data/Geo/Jord/PositionSpec.hs new file mode 100644 --- /dev/null +++ b/test/Data/Geo/Jord/PositionSpec.hs @@ -0,0 +1,153 @@ +module Data.Geo.Jord.PositionSpec + ( spec + ) where + +import Control.Exception.Base +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Antipode" $ do + it "returns the antipodal point" $ + antipode (readLatLong "484137N0061105E") `shouldBe` + latLongDecimal (-48.6936111) (-173.8152777) + it "returns the south pole when called with the north pole" $ + antipode (northPole :: LatLong) `shouldBe` latLongDecimal (-90.0) (-180.0) + it "returns the north pole when called with the south pole" $ + antipode (southPole :: LatLong) `shouldBe` latLongDecimal 90.0 (-180.0) + describe "Distance" $ do + it "returns 0 if both points are equal" $ + distance (readLatLong "500359N1795959W") (readLatLong "500359N1795959W") `shouldBe` + metres 0.0 + it "returns the distance between 2 points" $ + distance (readLatLong "500359N0054253W") (readLatLong "583838N0030412W") `shouldBe` + metres 968854.873 + it "handles singularity at the pole" $ + distance (northPole :: LatLong) (southPole :: LatLong) `shouldBe` + metres 2.00151144420359e7 + it "handles the discontinuity at the Date Line" $ + distance (readLatLong "500359N1795959W") (readLatLong "500359N1795959E") `shouldBe` + metres 39.66 + describe "Destination" $ do + it "return the given point if distance is 0 meter" $ + destination (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (metres 0) `shouldBe` + readLatLong "531914N0014347W" + it "return the destination point along great-circle at distance and bearing" $ + destination (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (metres 124800) `shouldBe` + latLongDecimal 53.1882691 0.1332744 + describe "Initial bearing" $ do + it "returns the 0 if both point are the same" $ + initialBearing (readLatLong "500359N0054253W") (readLatLong "500359N0054253W") `shouldBe` + decimalDegrees 0 + it "returns the initial bearing in compass degrees" $ + initialBearing (readLatLong "500359N0054253W") (readLatLong "583838N0030412W") `shouldBe` + decimalDegrees 9.1198181 + it "returns the initial bearing in compass degrees" $ + initialBearing (readLatLong "583838N0030412W") (readLatLong "500359N0054253W") `shouldBe` + decimalDegrees 191.2752013 + describe "Interpolate" $ do + it "fails if f < 0.0" $ + evaluate (interpolate (readLatLong "44N044E") (readLatLong "46N046E") (-0.5)) `shouldThrow` + errorCall "fraction must be in range [0..1], was -0.5" + it "fails if f > 1.0" $ + evaluate (interpolate (readLatLong "44N044E") (readLatLong "46N046E") 1.1) `shouldThrow` + errorCall "fraction must be in range [0..1], was 1.1" + it "returns p0 if f == 0" $ + interpolate (readLatLong "44N044E") (readLatLong "46N046E") 0.0 `shouldBe` + readLatLong "44N044E" + it "returns p1 if f == 1" $ + interpolate (readLatLong "44N044E") (readLatLong "46N046E") 1.0 `shouldBe` + readLatLong "46N046E" + it "returns the interpolated position" $ + interpolate + (readLatLong "53°28'46''N 2°14'43''W") + (readLatLong "55°36'21''N 13°02'09''E") + 0.5 `shouldBe` + latLongDecimal 54.7835574 5.1949856 + describe "isInside" $ do + it "return False if polygon is empty" $ isInside (latLongDecimal 45 1) [] `shouldBe` False + it "return False if polygon does not define at least a triangle" $ + isInside (latLongDecimal 45 1) [latLongDecimal 45 1, latLongDecimal 45 2] `shouldBe` + False + it "returns True if point is inside polygon" $ do + let polygon = + [ latLongDecimal 45 1 + , latLongDecimal 45 2 + , latLongDecimal 46 2 + , latLongDecimal 46 1 + ] + let p = latLongDecimal 45.1 1.1 + isInside p polygon `shouldBe` True + it "returns False if point is inside polygon" $ do + let polygon = + [ latLongDecimal 45 1 + , latLongDecimal 45 2 + , latLongDecimal 46 2 + , latLongDecimal 46 1 + ] + let p = antipode (latLongDecimal 45.1 1.1) + isInside p polygon `shouldBe` False + it "returns False if point is a vertex of the polygon" $ do + let polygon = + [ latLongDecimal 45 1 + , latLongDecimal 45 2 + , latLongDecimal 46 2 + , latLongDecimal 46 1 + ] + let p = latLongDecimal 45 1 + isInside p polygon `shouldBe` False + it "handles closed polygons" $ do + let polygon = + [ latLongDecimal 45 1 + , latLongDecimal 45 2 + , latLongDecimal 46 2 + , latLongDecimal 46 1 + , latLongDecimal 45 1 + ] + let p = latLongDecimal 45.1 1.1 + isInside p polygon `shouldBe` True + it "handles concave polygons" $ do + let malmo = latLongDecimal 55.6050 13.0038 + let ystad = latLongDecimal 55.4295 13.82 + let lund = latLongDecimal 55.7047 13.1910 + let helsingborg = latLongDecimal 56.0465 12.6945 + let kristianstad = latLongDecimal 56.0294 14.1567 + let polygon = [malmo, ystad, kristianstad, helsingborg, lund] + let hoor = latLongDecimal 55.9295 13.5297 + let hassleholm = latLongDecimal 56.1589 13.7668 + isInside hoor polygon `shouldBe` True + isInside hassleholm polygon `shouldBe` False + describe "Final bearing" $ do + it "returns the 180.0 if both point are the same" $ + finalBearing (readLatLong "500359N0054253W") (readLatLong "500359N0054253W") `shouldBe` + decimalDegrees 180 + it "returns the final bearing in compass degrees" $ + finalBearing (readLatLong "500359N0054253W") (readLatLong "583838N0030412W") `shouldBe` + decimalDegrees 11.2752013 + it "returns the final bearing in compass degrees" $ + finalBearing (readLatLong "583838N0030412W") (readLatLong "500359N0054253W") `shouldBe` + decimalDegrees 189.1198181 + it "returns the final bearing in compass degrees" $ + finalBearing (readLatLong "535941S0255915W") (readLatLong "54N154E") `shouldBe` + decimalDegrees 125.6839436 + describe "Mean" $ do + it "returns Nothing if no point is given" $ (mean [] :: Maybe LatLong) `shouldBe` Nothing + it "returns the unique given point" $ + mean [readLatLong "500359N0054253W"] `shouldBe` Just (readLatLong "500359N0054253W") + it "returns the geographical mean" $ + mean [readLatLong "500359N0054253W", readLatLong "583838N0030412W"] `shouldBe` + Just (latLongDecimal 54.3622869 (-4.5306725)) + it "returns Nothing if list contains antipodal points" $ do + let points = + [ latLongDecimal 45 1 + , latLongDecimal 45 2 + , latLongDecimal 46 2 + , latLongDecimal 46 1 + , antipode (latLongDecimal 45 2) + ] + mean points `shouldBe` Nothing + describe "North pole" $ + it "returns (90, 0)" $ (northPole :: LatLong) `shouldBe` latLongDecimal 90.0 0.0 + describe "South pole" $ + it "returns (-90, 0)" $ (southPole :: LatLong) `shouldBe` latLongDecimal (-90.0) 0.0