jord 0.1.0.0 → 0.2.0.0
raw patch · 16 files changed
+1536/−1494 lines, 16 filesdep ~base
Dependency ranges changed: base
Files
- ChangeLog.md +9/−3
- README.md +5/−4
- app/Main.hs +157/−157
- jord.cabal +10/−8
- src/Data/Geo/Jord.hs +5/−3
- src/Data/Geo/Jord/Eval.hs +408/−407
- src/Data/Geo/Jord/GeoPos.hs +0/−221
- src/Data/Geo/Jord/GreatCircle.hs +103/−311
- src/Data/Geo/Jord/LatLong.hs +221/−0
- src/Data/Geo/Jord/Position.hs +224/−0
- test/Data/Geo/Jord/AngleSpec.hs +89/−84
- test/Data/Geo/Jord/EvalSpec.hs +41/−41
- test/Data/Geo/Jord/GeoPosSpec.hs +0/−63
- test/Data/Geo/Jord/GreatCircleSpec.hs +48/−192
- test/Data/Geo/Jord/LatLongSpec.hs +63/−0
- test/Data/Geo/Jord/PositionSpec.hs +153/−0
ChangeLog.md view
@@ -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
README.md view
@@ -1,15 +1,14 @@-# Jord +# Jord - Geographical Position Calculations [](https://travis-ci.org/ofmooseandmen/jord) +[](http://hackage.haskell.org/package/jord) [](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
app/Main.hs view
@@ -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
jord.cabal view
@@ -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 <https://github.com/ofmooseandmen/jord#readme> 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
src/Data/Geo/Jord.hs view
@@ -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"
src/Data/Geo/Jord/Eval.hs view
@@ -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)
− src/Data/Geo/Jord/GeoPos.hs
@@ -1,221 +0,0 @@--- | --- Module: Data.Geo.Jord.GeoPos --- Copyright: (c) 2018 Cedric Liegeois --- License: BSD3 --- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> --- 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"
src/Data/Geo/Jord/GreatCircle.hs view
@@ -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)
+ src/Data/Geo/Jord/LatLong.hs view
@@ -0,0 +1,221 @@+-- | +-- Module: Data.Geo.Jord.GeoPos +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- 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"
+ src/Data/Geo/Jord/Position.hs view
@@ -0,0 +1,224 @@+-- | +-- Module: Data.Geo.Jord.GreatCircle +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with positions. +-- +-- All functions are implemented using the vector-based approached described in +-- <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation> +-- +-- 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))
test/Data/Geo/Jord/AngleSpec.hs view
@@ -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
test/Data/Geo/Jord/EvalSpec.hs view
@@ -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)
− test/Data/Geo/Jord/GeoPosSpec.hs
@@ -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"
test/Data/Geo/Jord/GreatCircleSpec.hs view
@@ -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"
+ test/Data/Geo/Jord/LatLongSpec.hs view
@@ -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"
+ test/Data/Geo/Jord/PositionSpec.hs view
@@ -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