jord (empty) → 0.1.0.0
raw patch · 21 files changed
+2413/−0 lines, 21 filesdep +HUnitdep +basedep +haskelinesetup-changed
Dependencies added: HUnit, base, haskeline, hspec, jord
Files
- ChangeLog.md +3/−0
- LICENSE +30/−0
- README.md +41/−0
- Setup.hs +2/−0
- app/Main.hs +169/−0
- jord.cabal +80/−0
- src/Data/Geo/Jord.hs +37/−0
- src/Data/Geo/Jord/Angle.hs +286/−0
- src/Data/Geo/Jord/Eval.hs +520/−0
- src/Data/Geo/Jord/GeoPos.hs +221/−0
- src/Data/Geo/Jord/GreatCircle.hs +341/−0
- src/Data/Geo/Jord/Length.hs +121/−0
- src/Data/Geo/Jord/NVector.hs +80/−0
- src/Data/Geo/Jord/Parse.hs +52/−0
- src/Data/Geo/Jord/Quantity.hs +17/−0
- test/Data/Geo/Jord/AngleSpec.hs +84/−0
- test/Data/Geo/Jord/EvalSpec.hs +41/−0
- test/Data/Geo/Jord/GeoPosSpec.hs +63/−0
- test/Data/Geo/Jord/GreatCircleSpec.hs +192/−0
- test/Data/Geo/Jord/LengthSpec.hs +32/−0
- test/Spec.hs +1/−0
+ ChangeLog.md view
@@ -0,0 +1,3 @@+### 0.1.0.0++- Initial version
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Cedric Liegeois (c) 2018++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Author name here nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,41 @@+# Jord + +[](https://travis-ci.org/ofmooseandmen/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). + +## How do I build it? + +If you have [Stack](https://docs.haskellstack.org/en/stable/README/), +then: +```sh +$ stack build --test +``` + +## How do I use it? + +```haskell +import Data.Geo.Jord + +-- destination position from 531914N0014347W having travelled 500Nm on a heading of 96.0217° +destination (readGeoPos "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) + +-- distance between 54°N,154°E and its antipodal position +let p = latLongDecimal 54 154 +distance p (antipode p) +``` + +Jord comes with a REPL (built with [haskeline](https://github.com/judah/haskeline)): + +```sh +$ jord-exe +jord> finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E +jord> angle: 126°0'0.0" +```
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/Main.hs view
@@ -0,0 +1,169 @@+-- | +-- Module: Main +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- 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+ 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
+ jord.cabal view
@@ -0,0 +1,80 @@+-- This file has been generated from package.yaml by hpack version 0.28.2.+--+-- see: https://github.com/sol/hpack+--+-- hash: c02feb8438ae6672aeef69101f3c30532d6cf7132fece36f1fb9b67b207c9332++name: jord+version: 0.1.0.0+synopsis: Geographic position calculations on Great Circles+description: Please see the README on GitHub at <https://github.com/ofmooseandmen/jord#readme>+category: Geography+stability: experimental+homepage: https://github.com/ofmooseandmen/jord+bug-reports: https://github.com/ofmooseandmen/jord/issues+author: Cedric Liegeois+maintainer: Cedric Liegeois <ofmooseandmen@yahoo.com>+copyright: 2018 Cedric Liegeois+license: BSD3+license-file: LICENSE+build-type: Simple+cabal-version: >= 1.10+extra-source-files:+ ChangeLog.md+ README.md++source-repository head+ type: git+ location: https://github.com/ofmooseandmen/jord++library+ exposed-modules:+ Data.Geo.Jord+ Data.Geo.Jord.Angle+ Data.Geo.Jord.Eval+ Data.Geo.Jord.Length+ Data.Geo.Jord.GeoPos+ Data.Geo.Jord.GreatCircle+ Data.Geo.Jord.Quantity+ Data.Geo.Jord.NVector+ other-modules:+ Data.Geo.Jord.Parse+ hs-source-dirs:+ src+ ghc-options: -Wall+ build-depends:+ base >=4.7 && <5+ default-language: Haskell2010++executable jord-exe+ main-is: Main.hs+ other-modules:+ Paths_jord+ hs-source-dirs:+ app+ ghc-options: -Wall+ build-depends:+ base >=4.7 && <5+ , haskeline >=0.7 && <0.8+ , jord+ default-language: Haskell2010++test-suite jord-test+ type: exitcode-stdio-1.0+ main-is: Spec.hs+ other-modules:+ Data.Geo.Jord.AngleSpec+ Data.Geo.Jord.EvalSpec+ Data.Geo.Jord.GeoPosSpec+ Data.Geo.Jord.GreatCircleSpec+ Data.Geo.Jord.LengthSpec+ Paths_jord+ hs-source-dirs:+ test+ ghc-options: -Wall+ build-depends:+ HUnit ==1.6.*+ , base >=4.7 && <5+ , hspec ==2.*+ , jord+ default-language: Haskell2010
+ src/Data/Geo/Jord.hs view
@@ -0,0 +1,37 @@+-- | +-- Module: Data.Geo.Jord +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Geographic position calculations (distance, bearing, intersection, etc...) on great circles using +-- the algorithms described in <http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf Gade, K. (2010). A Non-singular Horizontal Position Representation>. +-- +-- See <http://www.navlab.net/nvector Position calculations - simple and exact solutions> +-- +-- See <http://www.movable-type.co.uk/scripts/latlong-vectors.html Vector-based geodesy> +-- +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.Length + , module Data.Geo.Jord.NVector + , 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.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Quantity + +-- | version. +jordVersion :: String +jordVersion = "0.1.0.0"
+ src/Data/Geo/Jord/Angle.hs view
@@ -0,0 +1,286 @@+-- | +-- Module: Data.Geo.Jord.Angle +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with angles representing latitudes, longitude and bearings. +-- +module Data.Geo.Jord.Angle + (+ -- * The 'Angle' type+ Angle+ -- * Smart constructors + , decimalDegrees+ , dms+ , dmsE+ , dmsF+ -- * Calculations + , arcLength+ , central+ , isNegative+ , isWithin+ , negate'+ , normalise+ -- * Trigonometric functions + , atan2'+ , cos'+ , sin'+ -- * Accessors + , getDegrees+ , getMinutes+ , getSeconds+ , getMilliseconds+ , toDecimalDegrees+ -- * Read + , angle+ , readAngle+ , readAngleE+ , readAngleF+ ) where++import Control.Applicative+import Control.Monad.Fail+import Data.Fixed+import Data.Geo.Jord.Length+import Data.Geo.Jord.Parse+import Data.Geo.Jord.Quantity+import Data.Maybe+import Prelude hiding (fail, length)+import Text.ParserCombinators.ReadP+import Text.Read hiding (get, look, pfail)++-- | An angle with a resolution of a milliseconds of a degree. +-- When used as a latitude/longitude this roughly translate to a precision +-- of 30 millimetres at the equator. +newtype Angle = Angle+ { milliseconds :: Int+ } deriving (Eq)++-- | See 'readAngle'. +instance Read Angle where+ readsPrec _ = readP_to_S angle++-- | Angle is shown degrees, minutes, seconds and milliseconds - e.g. 154°25'43.5". +instance Show Angle where+ show a =+ show (getDegrees a) +++ "°" +++ show (getMinutes a) ++ "'" ++ show (getSeconds a) ++ "." ++ show (getMilliseconds a) ++ "\""++-- | Add/Subtract 'Angle'. +instance Quantity Angle where+ add (Angle millis1) (Angle millis2) = Angle (millis1 + millis2)+ sub (Angle millis1) (Angle millis2) = Angle (millis1 - millis2)++-- | 'Angle' from given decimal degrees. Any 'Double' is accepted: it must be +-- validated by the call site when used to represent a latitude or longitude. +decimalDegrees :: Double -> Angle+decimalDegrees dec = Angle (round (dec * 3600000.0))++-- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. +-- 'error's if given minutes, seconds and/or milliseconds are invalid. +-- Degrees are not validated and can be any 'Int': they must be validated by the call site +-- when used to represent a latitude or longitude. +dms :: Int -> Int -> Int -> Int -> Angle+dms degs mins secs millis =+ fromMaybe+ (error+ ("Invalid minutes=" +++ show mins ++ " or seconds=" ++ show secs ++ " or milliseconds=" ++ show millis))+ (dmsF degs mins secs millis)++-- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. +-- A 'Left' indicates that given minutes, seconds and/or milliseconds are invalid. +-- Degrees are not validated and can be any 'Int': they must be validated by the call site +-- when used to represent a latitude or longitude. +dmsE :: Int -> Int -> Int -> Int -> Either String Angle+dmsE degs mins secs millis+ | mins < 0 || mins > 59 = Left ("Invalid minutes: " ++ show mins)+ | secs < 0 || secs >= 60 = Left ("Invalid seconds: " ++ show secs)+ | millis < 0 || millis >= 1000 = Left ("Invalid milliseconds: " ++ show millis)+ | otherwise = Right (decimalDegrees ms)+ where+ ms =+ signed+ (fromIntegral (abs degs) + (fromIntegral mins / 60.0 :: Double) ++ (fromIntegral secs / 3600.0 :: Double) ++ (fromIntegral millis / 3600000.0 :: Double))+ (signum degs)++-- | 'Angle' from the given given degrees, minutes, seconds and milliseconds. +-- 'fail's if given minutes, seconds and/or milliseconds are invalid. +-- Degrees are not validated and can be any 'Int': they must be validated by the call site +-- when used to represent a latitude or longitude. +dmsF :: (MonadFail m) => Int -> Int -> Int -> Int -> m Angle+dmsF degs mins secs millis =+ case e of+ Left err -> fail err+ Right a -> return a+ where+ e = dmsE degs mins secs millis++-- | @arcLength a r@ computes the 'Length' of the arc that subtends the angle @a@ for radius @r@. +arcLength :: Angle -> Length -> Length+arcLength a r = metres (toMetres r * toRadians a)++-- | @central l r@ computes the central 'Angle' from the arc length @l@ and radius @r@. +central :: Length -> Length -> Angle+central s r = fromRadians (toMetres s / toMetres r)++-- | Returns the given 'Angle' negated. +negate' :: Angle -> Angle+negate' (Angle millis) = Angle (-millis)++-- | @normalise a n@ normalises @a@ to [0, @n@]. +normalise :: Angle -> Angle -> Angle+normalise a n = decimalDegrees dec+ where+ dec = mod' (toDecimalDegrees a + toDecimalDegrees n) 360.0++-- | Is given 'Angle' < 0? +isNegative :: Angle -> Bool+isNegative (Angle millis) = millis < 0++-- | Is given 'Angle' within range [@low@..@high@] inclusive? +isWithin :: Angle -> Angle -> Angle -> Bool+isWithin (Angle millis) (Angle low) (Angle high) = millis >= low && millis <= high++-- | @atan2' y x@ computes the 'Angle' (from the positive x-axis) of the vector from the origin to the point (x,y). +atan2' :: Double -> Double -> Angle+atan2' y x = fromRadians (atan2 y x)++-- | @cos' a@ returns the cosinus of @a@. +cos' :: Angle -> Double+cos' a = cos (toRadians a)++-- | @sin' a@ returns the sinus of @a@. +sin' :: Angle -> Double+sin' a = sin (toRadians a)++-- | radians to degrees. +fromRadians :: Double -> Angle+fromRadians r = decimalDegrees (r / pi * 180.0)++-- | degrees to radians. +toRadians :: Angle -> Double+toRadians a = toDecimalDegrees a * pi / 180.0++-- | Converts the given 'Angle' to decimal degrees. +toDecimalDegrees :: Angle -> Double+toDecimalDegrees (Angle millis) = fromIntegral millis / 3600000.0++-- | @getDegrees a@ returns the degree component of @a@. +getDegrees :: Angle -> Int+getDegrees a = signed (field a 3600000.0 360.0) (signum (milliseconds a))++-- | @getMinutes a@ returns the minute component of @a@. +getMinutes :: Angle -> Int+getMinutes a = field a 60000.0 60.0++-- | @getSeconds a@ returns the second component of @a@. +getSeconds :: Angle -> Int+getSeconds a = field a 1000.0 60.0++-- | @getMilliseconds a@ returns the milliseconds component of @a@. +getMilliseconds :: Angle -> Int+getMilliseconds (Angle millis) = mod millis 1000++field :: Angle -> Double -> Double -> Int+field (Angle millis) divisor modulo =+ truncate (mod' (fromIntegral (abs millis) / divisor) modulo) :: Int++signed :: (Num a, Num b, Ord b) => a -> b -> a+signed n s+ | s < 0 = -n+ | otherwise = n++-- | Parses and returns an 'Angle'. +angle :: ReadP Angle+angle = degsMinsSecs <|> decimal++-- | Obtains a 'Angle' from the given string formatted as either: +-- +-- * d°m′s.ms″ - e.g. 55°36'21.3", where minutes, seconds and milliseconds are optional. +-- +-- * decimal° - e.g. 55.6050° or -133° +-- +-- Symbols: +-- +-- * degree: ° or d +-- +-- * minute: ', ′ or m +-- +-- * second: ", ″, '' or s +-- +-- This simply calls @read s :: Angle@ so 'error' should be handled at the call site. +-- +readAngle :: String -> Angle+readAngle s = read s :: Angle++-- | Same as 'readAngle' but returns an 'Either'. +readAngleE :: String -> Either String Angle+readAngleE s =+ case readMaybe s of+ Nothing -> Left ("couldn't read angle " ++ s)+ Just a -> Right a++-- | Same as 'readAngle' but returns a 'MonadFail'. +readAngleF :: (MonadFail m) => String -> m Angle+readAngleF s =+ let p = readAngleE s+ in case p of+ Left e -> fail e+ Right l -> return l++-- | Parses DMS.MS and returns an 'Angle'. +degsMinsSecs :: ReadP Angle+degsMinsSecs = do+ d' <- fmap fromIntegral integer+ degSymbol+ (m', s', ms') <- option (0, 0, 0) (minsSecs <|> minsOnly)+ dmsF d' m' s' ms'++-- | Parses minutes, seconds with optionally milliseconds. +minsSecs :: ReadP (Int, Int, Int)+minsSecs = do+ m' <- natural+ minSymbol+ s' <- natural+ ms' <- option 0 (char '.' >> natural)+ secSymbol+ return (m', s', ms')++-- | Parses minutes. +minsOnly :: ReadP (Int, Int, Int)+minsOnly = do+ m' <- natural+ minSymbol+ return (m', 0, 0)++-- | Parses decimal degrees. +decimal :: ReadP Angle+decimal = do+ d <- double+ degSymbol+ return (decimalDegrees d)++-- | skips degree symbol. +degSymbol :: ReadP ()+degSymbol = do+ _ <- char '°' <|> char 'd'+ return ()++-- | skips minute symbol. +minSymbol :: ReadP ()+minSymbol = do+ _ <- char '\'' <|> char '′' <|> char 'm'+ return ()++-- | skips second symbol. +secSymbol :: ReadP ()+secSymbol = do+ _ <- string "\"" <|> string "''" <|> string "″" <|> string "s"+ return ()
+ src/Data/Geo/Jord/Eval.hs view
@@ -0,0 +1,520 @@+{-# LANGUAGE FlexibleInstances #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- | +-- Module: Data.Geo.Jord.Eval +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- 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)++-- | A value accepted and returned by 'eval'. +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 + | Vec NVector -- ^ 'NVector' + | Vecs [NVector] -- ^ list of 'NVector's + deriving (Eq, Show)++-- | 'Either' an error or a 'Value'. +type Result = Either String Value++-- | A location for 'Value's to be shared by successive evalations. +newtype Vault =+ Vault [(String, Value)]++-- | An empty 'Vault'. +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" +-- 126° +-- +-- @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'. +-- +-- If the evaluation is successful, returns the resulting 'Value' ('Right') otherwise +-- a description of the error ('Left'). +-- +-- @ +-- vault = emptyVault +-- angle = eval "finalBearing 54N154E 54S154W" vault -- Right Ang +-- length = eval "distance (antipode 54N154E) 54S154W" vault -- Right Len +-- -- parameter resolution from vault +-- a1 = eval "finalBearing 54N154E 54S154W" vault +-- vault = insert "a1" vault +-- a2 = eval "(finalBearing a1 54S154W)" vault +-- @ +-- +-- By default, all returned positions are 'Geo' ('GeoPos'), to get back a 'Vec' ('NVector'), the +-- expression must be wrapped by 'toNVector'. +-- +-- @ +-- dest = eval "destination 54°N,154°E 54° 1000m" -- Right Geo +-- dest = eval "toNVector (destination 54°N,154°E 54° 1000m)" -- Right Vec +-- @ +-- +-- Every function call must be wrapped between parentheses, however they can be ommitted for the top level call. +-- +-- @ +-- angle = eval "finalBearing 54N154E 54S154W" -- Right Ang +-- angle = eval "(finalBearing 54N154E 54S154W)" -- Right Ang +-- length = eval "distance (antipode 54N154E) 54S154W" -- Right Len +-- 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++-- | All supported functions: +-- +-- * 'antipode' +-- +-- * 'crossTrackDistance' +-- +-- * 'decimalDegrees' +-- +-- * 'destination' +-- +-- * 'distance' +-- +-- * 'finalBearing' +-- +-- * 'greatCircle' +-- +-- * 'initialBearing' +-- +-- * 'interpolate' +-- +-- * 'intersections' +-- +-- * 'isInside' +-- +-- * 'mean' +-- +-- * 'latLong' +-- +-- * 'latLongDecimal' +-- +-- * 'readGeoPos' +-- +-- * 'toDecimalDegrees' +-- +-- * 'toKilometres' +-- +-- * 'toMetres' +-- +-- * 'toNauticalMiles' +-- +-- * 'toNVector' +-- +functions :: [String]+functions =+ [ "antipode"+ , "crossTrackDistance"+ , "destination"+ , "decimalDegrees"+ , "distance"+ , "finalBearing"+ , "greatCircle"+ , "initialBearing"+ , "interpolate"+ , "intersections"+ , "isInside"+ , "latLong"+ , "latLongDecimal"+ , "mean"+ , "readGeoPos"+ , "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++-- | @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)++-- | @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)++------------------------------------------ +-- 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++-- | 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)++----------------------------------------- +-- Syntactic Analysis: [Token] -> Ast -- +----------------------------------------- +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)++------------------------------------- +-- 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)
+ src/Data/Geo/Jord/GeoPos.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.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
@@ -0,0 +1,341 @@+-- | +-- 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 <https://en.wikipedia.org/wiki/Great_circle Great Circle>. +-- +-- 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. +-- +-- 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)++-- | 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. +-- +-- It is internally represented as its normal vector - i.e. the normal vector +-- to the plane containing the great circle. +-- +-- 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: +-- +-- @ +-- 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++-- | 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
+ src/Data/Geo/Jord/Length.hs view
@@ -0,0 +1,121 @@+-- | +-- Module: Data.Geo.Jord.Length +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with (signed) lengths in metres, kilometres or nautical miles. +-- +module Data.Geo.Jord.Length + ( + -- * The 'Length' type + Length(millimetres) + -- * Smart constructors + , kilometres + , metres + , nauticalMiles + -- * Read + , readLength + , readLengthE + , readLengthF + -- * Conversions + , toKilometres + , toMetres + , toNauticalMiles + -- * Misc. + , isZero + ) where + +import Control.Applicative +import Control.Monad.Fail +import Data.Geo.Jord.Parse +import Data.Geo.Jord.Quantity +import Prelude hiding (fail, length) +import Text.ParserCombinators.ReadP +import Text.Read hiding (pfail) + +-- | A length with a resolution of 1 millimetre. +newtype Length = Length + { millimetres :: Int + } deriving (Eq) + +-- | See 'readLength'. +instance Read Length where + readsPrec _ = readP_to_S length + +-- | Length is shown in metres when <= 10,000 m and in kilometres otherwise. +instance Show Length where + show l + | m <= 10000.0 = show m ++ "m" + | otherwise = show (m / 1000.0) ++ "km" + where + m = toMetres l + +-- | Add/Subtract Length. +instance Quantity Length where + add a b = Length (millimetres a + millimetres b) + sub a b = Length (millimetres a - millimetres b) + +-- | 'Length' from given amount of nautical miles. +nauticalMiles :: Double -> Length +nauticalMiles nm = metres (nm * 1852.0) + +-- | 'Length' from given amount of metres. +metres :: Double -> Length +metres m = Length (round (m * 1000.0)) + +-- | 'Length' from given amount of kilometres. +kilometres :: Double -> Length +kilometres km = metres (km * 1000.0) + +-- | Obtains a 'Length' from the given string formatted as (-)float[m|km|nm] - e.g. 3000m, 2.5km or -154nm. +-- +-- This simply calls @read s :: Length@ so 'error' should be handled at the call site. +-- +readLength :: String -> Length +readLength s = read s :: Length + +-- | Same as 'readLength' but returns a 'Either'. +readLengthE :: String -> Either String Length +readLengthE s = + case readMaybe s of + Nothing -> Left ("couldn't read length " ++ s) + Just l -> Right l + +-- | Same as 'readLength' but returns a 'MonadFail'. +readLengthF :: (MonadFail m) => String -> m Length +readLengthF s = + let p = readEither s + in case p of + Left e -> fail e + Right l -> return l + +-- | @toKilometres l@ converts @l@ to kilometres. +toKilometres :: Length -> Double +toKilometres l = toMetres l / 1000.0 + +-- | @toMetres l@ converts @l@ to metres. +toMetres :: Length -> Double +toMetres (Length mm) = fromIntegral mm / 1000.0 + +-- | @toNauticalMiles l@ converts @l@ to nautical miles. +toNauticalMiles :: Length -> Double +toNauticalMiles l = toMetres l / 1852.0 + +-- | Is given 'Length' == 0? +isZero :: Length -> Bool +isZero (Length mm) = mm == 0 + +-- | Parses and returns a 'Length'. +length :: ReadP Length +length = do + v <- number + skipSpaces + u <- string "m" <|> string "km" <|> string "Nm" + case u of + "m" -> return (metres v) + "km" -> return (kilometres v) + "Nm" -> return (nauticalMiles v) + _ -> pfail
+ src/Data/Geo/Jord/NVector.hs view
@@ -0,0 +1,80 @@+-- | +-- Module: Data.Geo.Jord.NVector +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Types and functions for working with n-vectors. +-- +module Data.Geo.Jord.NVector + ( NVector(x, y, z) + , nvector + , cross + , dot + , norm + , scale + , unit + , zero + ) where + +import Data.Geo.Jord.Quantity + +-- | Represents a position as the normal vector to the sphere. +data NVector = NVector + { x :: Double + , y :: Double + , z :: Double + } deriving (Eq, Show) + +-- | Add and subtract 'NVector's. +instance Quantity NVector where + add a b = NVector x' y' z' + where + x' = x a + x b + y' = y a + y b + z' = z a + z b + sub a b = NVector x' y' z' + where + x' = x a - x b + y' = y a - y b + z' = z a - z b + +-- | Smart 'NVector' constructor. The returned 'NVector' is a 'unit' vector. +nvector :: Double -> Double -> Double -> NVector +nvector x' y' z' = unit (NVector x' y' z') + +-- | Computes the cross product of the two given 'NVector's. +cross :: NVector -> NVector -> NVector +cross a b = NVector x' y' z' + where + x' = y a * z b - z a * y b + y' = z a * x b - x a * z b + z' = x a * y b - y a * x b + +-- | Computes the dot product of the two given 'NVector's. +dot :: NVector -> NVector -> Double +dot a b = x a * x b + y a * y b + z a * z b + +-- | Computes the norm of the given 'NVector'. +norm :: NVector -> Double +norm a = sqrt ((x a * x a) + (y a * y a) + (z a * z a)) + +-- | Multiplies each component of the given 'NVector' by the given value. +scale :: NVector -> Double -> NVector +scale a s = NVector x' y' z' + where + x' = x a * s + y' = y a * s + z' = z a * s + +-- | Normalises the given 'NVector'. +unit :: NVector -> NVector +unit a = scale a s + where + s = 1.0 / norm a + +-- | [0, 0, 0] - not a valid 'NVector', but can be used as the identity value during reduction. +zero :: NVector +zero = NVector 0.0 0.0 0.0
+ src/Data/Geo/Jord/Parse.hs view
@@ -0,0 +1,52 @@+-- |+-- Module: Data.Geo.Jord.Parse+-- Copyright: (c) 2018 Cedric Liegeois+-- License: BSD3+-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr>+-- Stability: experimental+-- Portability: portable+--+-- internal 'ReadP' parsers used by "Jord".+--+module Data.Geo.Jord.Parse+ ( digits+ , double+ , integer+ , natural+ , number+ ) where++import Control.Applicative+import Data.Char+import Text.ParserCombinators.ReadP++-- | Parses the given number of digits and returns the read 'Int'.+digits :: Int -> ReadP Int+digits n = fmap read (count n digit)++-- | Parses optionally a @-@ followed by a 'positive'.'positive' and returns the read 'Double'.+double :: ReadP Double+double = do+ s <- option 1.0 (fmap (\_ -> -1.0) (char '-'))+ i <- natural+ f <- char '.' >> natural+ return (s * (read (show i ++ "." ++ show f) :: Double))++-- | Parses optionally a @-@ followed by a 'positive' and returns the read 'Int'.+integer :: ReadP Int+integer = do+ s <- option 1 (fmap (\_ -> -1) (char '-'))+ p <- natural+ return (s * p)++-- | Parses 1 or more 'digit's and returns the read 'Int'.+natural :: ReadP Int+natural = fmap read (munch1 isDigit)++-- | Parses an 'integer' or 'double' and returns the read 'Double'.+number :: ReadP Double+number = double <|> fmap fromIntegral integer++-- | Parses and returns a digit.+digit :: ReadP Char+digit = satisfy isDigit
+ src/Data/Geo/Jord/Quantity.hs view
@@ -0,0 +1,17 @@+-- |+-- Module: Data.Geo.Jord.Quantity+-- Copyright: (c) 2018 Cedric Liegeois+-- License: BSD3+-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr>+-- Stability: experimental+-- Portability: portable+--+-- Defines the class 'Quantity' for something that can be added or subtracted.+--+module Data.Geo.Jord.Quantity+ ( Quantity(..)+ ) where++-- | Something that can be added or subtracted.+class Quantity a where+ add, sub :: a -> a -> a
+ test/Data/Geo/Jord/AngleSpec.hs view
@@ -0,0 +1,84 @@+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
+ test/Data/Geo/Jord/EvalSpec.hs view
@@ -0,0 +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)
+ test/Data/Geo/Jord/GeoPosSpec.hs view
@@ -0,0 +1,63 @@+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
@@ -0,0 +1,192 @@+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
+ test/Data/Geo/Jord/LengthSpec.hs view
@@ -0,0 +1,32 @@+module Data.Geo.Jord.LengthSpec+ ( spec+ ) where++import Data.Geo.Jord+import Test.Hspec++spec :: Spec+spec = do+ describe "Reading valid lengths" $ do+ it "reads -15.2m" $ readLength "-15.2m" `shouldBe` metres (-15.2)+ it "reads 154km" $ readLength "154km" `shouldBe` kilometres 154+ it "reads 1000Nm" $ readLength "1000Nm" `shouldBe` nauticalMiles 1000+ describe "Reading invalid lengths" $ do+ it "fails to read 5" $ readLengthE "5" `shouldBe` Left "couldn't read length 5"+ it "fails to read 5nmi" $ readLengthE "5nmi" `shouldBe` Left "couldn't read length 5nmi"+ describe "Showing lengths" $ do+ it "shows length in metres when <= 10000 m" $ show (metres 5) `shouldBe` "5.0m"+ it "shows length in kilometres when > 10000 m" $+ show (kilometres 1000) `shouldBe` "1000.0km"+ describe "Converting lengths" $ do+ it "converts metres to kilometres" $ toKilometres (metres 1000) `shouldBe` 1.0+ it "converts metres to nautical miles" $+ toNauticalMiles (metres 1000) `shouldBe` 0.5399568034557235+ it "converts kilometres to nautical miles" $+ toNauticalMiles (kilometres 1000) `shouldBe` 539.9568034557235+ it "converts nautical miles to metres" $ toMetres (nauticalMiles 10.5) `shouldBe` 19446+ it "converts nautical miles to kilometres" $+ toKilometres (nauticalMiles 10.5) `shouldBe` 19.446+ describe "Adding/Subtracting lengths" $ do+ it "adds lengths" $ add (kilometres 1000) (metres 1000) `shouldBe` metres 1001000+ it "subtracts lengths" $ sub (metres 1000) (nauticalMiles 10.5) `shouldBe` metres (-18446)
+ test/Spec.hs view
@@ -0,0 +1,1 @@+{-# OPTIONS_GHC -F -pgmF hspec-discover #-}