jord 0.2.0.0 → 0.3.0.0
raw patch · 33 files changed
+2782/−1644 lines, 33 files
Files
- ChangeLog.md +5/−0
- README.md +18/−5
- app/Eval.hs +487/−0
- app/Main.hs +170/−157
- jord.cabal +18/−10
- src/Data/Geo/Jord.hs +19/−7
- src/Data/Geo/Jord/Angle.hs +19/−4
- src/Data/Geo/Jord/AngularPosition.hs +59/−0
- src/Data/Geo/Jord/Earth.hs +117/−0
- src/Data/Geo/Jord/EcefPosition.hs +59/−0
- src/Data/Geo/Jord/Eval.hs +0/−521
- src/Data/Geo/Jord/Frames.hs +290/−0
- src/Data/Geo/Jord/Geodetics.hs +317/−0
- src/Data/Geo/Jord/GreatCircle.hs +0/−133
- src/Data/Geo/Jord/LatLong.hs +15/−15
- src/Data/Geo/Jord/Length.hs +99/−93
- src/Data/Geo/Jord/NVector.hs +30/−69
- src/Data/Geo/Jord/Position.hs +0/−224
- src/Data/Geo/Jord/Quantity.hs +13/−12
- src/Data/Geo/Jord/Rotation.hs +138/−0
- src/Data/Geo/Jord/Transform.hs +192/−0
- src/Data/Geo/Jord/Vector3d.hs +124/−0
- test/Data/Geo/Jord/AngleSpec.hs +92/−89
- test/Data/Geo/Jord/EarthSpec.hs +27/−0
- test/Data/Geo/Jord/EvalSpec.hs +0/−41
- test/Data/Geo/Jord/FramesSpec.hs +77/−0
- test/Data/Geo/Jord/GeodeticsSpec.hs +213/−0
- test/Data/Geo/Jord/GreatCircleSpec.hs +0/−48
- test/Data/Geo/Jord/LatLongSpec.hs +63/−63
- test/Data/Geo/Jord/LengthSpec.hs +3/−0
- test/Data/Geo/Jord/PositionSpec.hs +0/−153
- test/Data/Geo/Jord/RotationSpec.hs +49/−0
- test/Data/Geo/Jord/TransformSpec.hs +69/−0
ChangeLog.md view
@@ -1,3 +1,8 @@+### 0.3.0.0 + +- Added ellipsoid, ECEF positions +- Added NEDVector + ### 0.2.0.0 - GeoPos -> LatLong
README.md view
@@ -10,6 +10,11 @@ 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). +- Transformation between ECEF (earth-centred, earth-fixed), Latitude/Longitude and N-Vector positions for spherical and ellipsoidal earth model +- Transformation between Latitude/Longitude and N-Vector positions +- Local, Body and North, East, Down Frames: delta between positions, target position from reference position and delta +- surface distance, initial & final bearing, interpolated position, great circle intersections, cross track distance, ... + ## How do I build it? If you have [Stack](https://docs.haskellstack.org/en/stable/README/), @@ -25,12 +30,20 @@ ```haskell import Data.Geo.Jord +-- Delta between positions in frameL +let p1 = decimalLatLongHeight 1 2 (metres (-3)) +let p2 = decimalLatLongHeight 4 5 (metres (-6)) +let w = decimalDegrees 5 -- wander azimuth +deltaBetween p1 p2 (frameL w) wgs84 -- = deltaMetres 359490.579 302818.523 17404.272 + -- destination position from 531914N0014347W having travelled 500Nm on a heading of 96.0217° -destination (readGeoPos "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) +-- using mean earth radius derived from the WG84 ellipsoid +destination (readLatLong "531914N0014347W") (decimalDegrees 96.0217) (nauticalMiles 500) r84 --- distance between 54°N,154°E and its antipodal position -let p = latLongDecimal 54 154 -distance p (antipode p) +-- surface distance between 54°N,154°E and its antipodal position +-- using mean earth radius derived from the WG84 ellipsoid +let p = decimalLatLong 54 154 +surfaceDistance p (antipode p) r84 ``` Jord comes with a REPL (built with [haskeline](https://github.com/judah/haskeline)): @@ -38,5 +51,5 @@ ```sh $ jord-exe jord> finalBearing (destination (antipode 54°N,154°E) 54° 1000m) 54°N,154°E -jord> angle: 126°0'0.0" +jord> angle: 126°0'0.0" (126.0) ```
+ app/Eval.hs view
@@ -0,0 +1,487 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE TupleSections #-}+{-# OPTIONS_GHC -fno-warn-orphans #-}++-- | +-- Module: 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 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.AngularPosition+import Data.Geo.Jord.Geodetics+import Data.Geo.Jord.LatLong+import Data.Geo.Jord.Length+import Data.Geo.Jord.NVector+import Data.Geo.Jord.Quantity+import Data.Geo.Jord.Transform+import Data.List hiding (delete, insert, lookup)+import Data.Maybe+import Prelude hiding (fail, lookup)+import Text.ParserCombinators.ReadP+import Text.Read (readEither, readMaybe)++-- | A value accepted and returned by 'eval'. +data Value+ = Ang Angle -- ^ angle + | Bool Bool -- ^ boolean + | Double Double -- ^ double + | Len Length -- ^ length + | Gc GreatCircle -- ^ great circle + | Geo (AngularPosition LatLong) -- ^ latitude, longitude and height + | NVec (AngularPosition NVector) -- ^ n-vector and height + | Vals [Value] -- array of values + 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 'LatLong'. +-- +-- 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 "surfaceDistance (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 +-- @ +-- +-- All returned positions are 'LatLong' by default, to get back a n-vector the +-- expression must be wrapped by 'toNVector'. +-- +-- @ +-- dest = eval "destination 54°N,154°E 54° 1000m" -- Right Ll +-- dest = eval "toNVector (destination 54°N,154°E 54° 1000m)" -- Right NVec +-- @ +-- +-- 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, expr') -> convert (evalExpr expr' r) rvec++convert :: Result -> Bool -> Result+convert r True = r+convert r False =+ case r of+ Right v@(NVec _) -> Right (toGeo v)+ Right (Vals vs) -> Right (Vals (map toGeo vs))+ oth -> oth++toGeo :: Value -> Value+toGeo (NVec v) = Geo (fromNVector v)+toGeo val = val++-- | All supported functions. +functions :: [String]+functions =+ [ "antipode"+ , "crossTrackDistance"+ , "destination"+ , "finalBearing"+ , "geoPos"+ , "greatCircle"+ , "initialBearing"+ , "interpolate"+ , "intersections"+ , "insideSurface"+ , "mean"+ , "surfaceDistance"+ , "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 (expectVec ts,) (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 geo) -> Right (NVec (toNVector geo))+ Just v -> Right v+ Nothing -> tryRead p+evalExpr (Antipode a) vault =+ case evalExpr a vault of+ (Right (NVec p)) -> Right (NVec (antipode p))+ r -> Left ("Call error: antipode " ++ showErr [r])+evalExpr (CrossTrackDistance a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (NVec p), Right (Gc gc)] -> Right (Len (crossTrackDistance84 p gc))+ r -> Left ("Call error: crossTrackDistance " ++ showErr r)+evalExpr (Destination a b c) vault =+ case [evalExpr a vault, evalExpr b vault, evalExpr c vault] of+ [Right (NVec p), Right (Ang a'), Right (Len l)] -> Right (NVec (destination84 p a' l))+ [Right (NVec p), Right (Double a'), Right (Len l)] ->+ Right (NVec (destination84 p (decimalDegrees a') l))+ r -> Left ("Call error: destination " ++ showErr r)+evalExpr (FinalBearing a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (NVec p1), Right (NVec p2)] ->+ maybe+ (Left "Call error: finalBearing identical points")+ (Right . Ang)+ (finalBearing p1 p2)+ r -> Left ("Call error: finalBearing " ++ showErr r)+evalExpr (GeoPos as) vault =+ case vs of+ [Right p@(NVec _)] -> Right p+ [Right (NVec v), Right (Len h)] -> Right (NVec (AngularPosition (pos v) h))+ [Right (Double lat), Right (Double lon)] ->+ bimap+ (\e -> "Call error: geoPos : " ++ e)+ (NVec . toNVector)+ (decimalLatLongHeightE lat lon zero)+ [Right (Double lat), Right (Double lon), Right (Len h)] ->+ bimap+ (\e -> "Call error: geoPos : " ++ e)+ (NVec . toNVector)+ (decimalLatLongHeightE lat lon h)+ [Right (Double lat), Right (Double lon), Right (Double h)] ->+ bimap+ (\e -> "Call error: geoPos : " ++ e)+ (NVec . toNVector)+ (decimalLatLongHeightE lat lon (metres h))+ r -> Left ("Call error: geoPos " ++ showErr r)+ where+ vs = map (`evalExpr` vault) as+evalExpr (GreatCircleSC a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (NVec p1), Right (NVec p2)] -> bimap id Gc (greatCircleE p1 p2)+ [Right (NVec 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 (NVec p1), Right (NVec p2)] ->+ maybe+ (Left "Call error: initialBearing identical points")+ (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 (NVec p1), Right (NVec p2)] -> Right (NVec (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 (Vals []))+ (\is -> Right (Vals [NVec (fst is), NVec (snd is)]))+ (intersections gc1 gc2 :: Maybe (AngularPosition NVector, AngularPosition NVector))+ r -> Left ("Call error: intersections " ++ showErr r)+evalExpr (InsideSurface as) vault =+ let m = map (`evalExpr` vault) as+ ps = [p | Right (NVec p) <- m]+ in if length m == length ps && length ps > 3+ then Right (Bool (insideSurface (head ps) (tail ps)))+ else Left ("Call error: insideSurface " ++ showErr m)+evalExpr (Mean as) vault =+ let m = map (`evalExpr` vault) as+ ps = [p | Right (NVec p) <- m]+ in if length m == length ps+ then maybe (Left ("Call error: mean " ++ showErr m)) (Right . NVec) (mean ps)+ else Left ("Call error: mean " ++ showErr m)+evalExpr (SurfaceDistance a b) vault =+ case [evalExpr a vault, evalExpr b vault] of+ [Right (NVec p1), Right (NVec p2)] -> Right (Len (surfaceDistance84 p1 p2))+ r -> Left ("Call error: surfaceDistance " ++ 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 (NVec _)) -> 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 geo), _] -> Right (NVec (toNVector geo))+ [_, _, _, d@(Right (Double _))] -> d+ _ -> Left ("couldn't read " ++ s)+ where+ r =+ map+ ($ s)+ [ readE readAngleE Ang+ , readE readLengthE Len+ , readE readLatLongE (\ll -> Geo (AngularPosition ll zero))+ , readE readEither Double+ ]++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+ | Destination Expr+ Expr+ Expr+ | FinalBearing Expr+ Expr+ | GeoPos [Expr]+ | GreatCircleSC Expr+ Expr+ | InitialBearing Expr+ Expr+ | Interpolate Expr+ Expr+ Double+ | Intersections Expr+ Expr+ | InsideSurface [Expr]+ | Mean [Expr]+ | SurfaceDistance Expr+ 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 "destination" [e1, e2, e3]) = do+ p1 <- transform e1+ p2 <- transform e2+ p3 <- transform e3+ return (Destination p1 p2 p3)+transform (Call "finalBearing" [e1, e2]) = do+ p1 <- transform e1+ p2 <- transform e2+ return (FinalBearing p1 p2)+transform (Call "geoPos" e) = do+ ps <- mapM transform e+ return (GeoPos ps)+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 "insideSurface" e) = do+ ps <- mapM transform e+ return (InsideSurface ps)+transform (Call "mean" e) = do+ ps <- mapM transform e+ return (Mean ps)+transform (Call "surfaceDistance" [e1, e2]) = do+ p1 <- transform e1+ p2 <- transform e2+ return (SurfaceDistance p1 p2)+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)
app/Main.hs view
@@ -8,162 +8,175 @@ -- -- 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 Eval+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" ++ - " 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 + 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) 54°N,154°E)\n" +++ "\n" +++ " Top level () can be ommitted: antipode 54N028E\n" +++ "\n Position calculations (Spherical Earth):\n\n" +++ " The following calculations assume a spherical earth model with a radius\n" +++ " derived from the WGS84 ellipsoid: " +++ show r84 +++ "\n" +++ "\n antipode pos antipodal point of pos\n" +++ " crossTrackDistance pos gc signed distance from pos to great circle gc\n" +++ " destination pos ang len destination position from pos having travelled len\n" +++ " on initial bearing ang (either in text form or decimal degrees)\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" +++ " insideSurface pos [pos] is p inside surface polygon?\n" +++ " mean [pos] geographical mean surface position of [pos]\n" +++ " surfaceDistance pos1 pos2 surface distance between pos1 and pos2\n" +++ "\n Constructors and conversions:\n\n" +++ " geoPos latlong surface geographic position from latlong\n" +++ " geoPos latlong height geographic position from latlong and height\n" +++ " geoPos lat long height geographic position from decimal latitude, longitude and height\n" +++ " geoPos lat long metres geographic position from decimal latitude, longitude and metres\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" +++ " 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 ++ " (" ++ show (toDecimalDegrees a) ++ ")"+showV (Bool b) = show b+showV (Double d) = show d+showV (Gc gc) = "great circle: " ++ show gc+showV (Len l) = "length: " ++ show l+showV (Geo g) =+ "latitude, longitude: " +++ show ll +++ " (" +++ show (toDecimalDegrees (latitude ll)) +++ ", " ++ show (toDecimalDegrees (longitude ll)) ++ ") - height: " ++ show h+ where+ ll = pos g+ h = height g+showV (NVec nv) =+ "n-vector: (" ++ show x ++ ", " ++ show y ++ ", " ++ show z ++ ") - height: " ++ show h+ where+ v = vec (pos nv)+ x = vx v+ y = vy v+ z = vz v+ h = height nv+showV (Vals []) = "empty"+showV (Vals vs) = "\n " ++ intercalate "\n " (map showV vs)++showVar :: String -> Value -> String+showVar n v = n ++ "=" ++ showV v
jord.cabal view
@@ -2,10 +2,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: c1197ea6dc351c977d808ecb2b0eef21c2993f1a1c14e40471ff7a0b18fcd501+-- hash: 0f7d62b4ffb4fdc3ff1a411f577d06587329b2b015bbfece734345148e0ebed8 name: jord-version: 0.2.0.0+version: 0.3.0.0 synopsis: Geographical Position Calculations description: Please see the README on GitHub at <https://github.com/ofmooseandmen/jord#readme> category: Geography@@ -31,13 +31,18 @@ exposed-modules: Data.Geo.Jord Data.Geo.Jord.Angle- Data.Geo.Jord.Eval- Data.Geo.Jord.Length+ Data.Geo.Jord.AngularPosition+ Data.Geo.Jord.Earth+ Data.Geo.Jord.EcefPosition+ Data.Geo.Jord.Frames+ Data.Geo.Jord.Geodetics Data.Geo.Jord.LatLong- Data.Geo.Jord.GreatCircle- Data.Geo.Jord.Position- Data.Geo.Jord.Quantity+ Data.Geo.Jord.Length Data.Geo.Jord.NVector+ Data.Geo.Jord.Quantity+ Data.Geo.Jord.Rotation+ Data.Geo.Jord.Transform+ Data.Geo.Jord.Vector3d other-modules: Data.Geo.Jord.Parse hs-source-dirs:@@ -50,6 +55,7 @@ executable jord-exe main-is: Main.hs other-modules:+ Eval Paths_jord hs-source-dirs: app@@ -65,11 +71,13 @@ main-is: Spec.hs other-modules: Data.Geo.Jord.AngleSpec- Data.Geo.Jord.EvalSpec- Data.Geo.Jord.GreatCircleSpec+ Data.Geo.Jord.EarthSpec+ Data.Geo.Jord.FramesSpec+ Data.Geo.Jord.GeodeticsSpec Data.Geo.Jord.LatLongSpec Data.Geo.Jord.LengthSpec- Data.Geo.Jord.PositionSpec+ Data.Geo.Jord.RotationSpec+ Data.Geo.Jord.TransformSpec Paths_jord hs-source-dirs: test
src/Data/Geo/Jord.hs view
@@ -13,27 +13,39 @@ -- -- See <http://www.movable-type.co.uk/scripts/latlong-vectors.html Vector-based geodesy> -- +-- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> +-- module Data.Geo.Jord ( module Data.Geo.Jord.Angle - , module Data.Geo.Jord.Eval - , module Data.Geo.Jord.GreatCircle + , module Data.Geo.Jord.AngularPosition + , module Data.Geo.Jord.Earth + , module Data.Geo.Jord.EcefPosition + , module Data.Geo.Jord.Frames + , module Data.Geo.Jord.Geodetics , 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 + , module Data.Geo.Jord.Rotation + , module Data.Geo.Jord.Transform + , module Data.Geo.Jord.Vector3d , jordVersion ) where import Data.Geo.Jord.Angle -import Data.Geo.Jord.Eval -import Data.Geo.Jord.GreatCircle +import Data.Geo.Jord.AngularPosition +import Data.Geo.Jord.Earth +import Data.Geo.Jord.EcefPosition +import Data.Geo.Jord.Frames +import Data.Geo.Jord.Geodetics 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.Geo.Jord.Rotation +import Data.Geo.Jord.Transform +import Data.Geo.Jord.Vector3d -- | version. jordVersion :: String -jordVersion = "0.2.0.0" +jordVersion = "0.3.0.0"
src/Data/Geo/Jord/Angle.hs view
@@ -8,7 +8,7 @@ -- -- Types and functions for working with angles representing latitudes, longitude and bearings. -- -module Data.Geo.Jord.Angle +module Data.Geo.Jord.Angle ( -- * The 'Angle' type Angle@@ -25,6 +25,7 @@ , negate' , normalise -- * Trigonometric functions + , asin' , atan2' , cos' , sin'@@ -50,6 +51,7 @@ import Data.Maybe import Prelude hiding (fail, length) import Text.ParserCombinators.ReadP+import Text.Printf import Text.Read hiding (get, look, pfail) -- | An angle with a resolution of a milliseconds of a degree. @@ -66,14 +68,23 @@ -- | Angle is shown degrees, minutes, seconds and milliseconds - e.g. 154°25'43.5". instance Show Angle where show a =- show (getDegrees a) +++ s +++ show d ++ "°" ++- show (getMinutes a) ++ "'" ++ show (getSeconds a) ++ "." ++ show (getMilliseconds a) ++ "\""+ show (getMinutes a) +++ "'" ++ show (getSeconds a) ++ "." ++ printf "%03d" (getMilliseconds a) ++ "\""+ where+ d = getDegrees a+ s =+ if d == 0 && milliseconds a < 0+ then "-"+ else "" -- | Add/Subtract 'Angle'. instance Quantity Angle where add (Angle millis1) (Angle millis2) = Angle (millis1 + millis2) sub (Angle millis1) (Angle millis2) = Angle (millis1 - millis2)+ zero = Angle 0 -- | '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. @@ -152,6 +163,10 @@ atan2' :: Double -> Double -> Angle atan2' y x = fromRadians (atan2 y x) +-- | @asin' a@ computes arc sine of @a@. +asin' :: Double -> Angle+asin' a = fromRadians (asin a)+ -- | @cos' a@ returns the cosinus of @a@. cos' :: Angle -> Double cos' a = cos (toRadians a)@@ -186,7 +201,7 @@ -- | @getMilliseconds a@ returns the milliseconds component of @a@. getMilliseconds :: Angle -> Int-getMilliseconds (Angle millis) = mod millis 1000+getMilliseconds (Angle millis) = mod (abs millis) 1000 field :: Angle -> Double -> Double -> Int field (Angle millis) divisor modulo =
+ src/Data/Geo/Jord/AngularPosition.hs view
@@ -0,0 +1,59 @@+-- | +-- Module: Data.Geo.Jord.AngularPosition +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Angular positions. +-- +-- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> +-- +module Data.Geo.Jord.AngularPosition+ ( AngularPosition(..)+ , latLongHeight+ , decimalLatLongHeight+ , decimalLatLongHeightE+ , decimalLatLongHeightF+ , nvectorHeight+ ) where++import Control.Monad.Fail+import Data.Geo.Jord.LatLong+import Data.Geo.Jord.Length+import Data.Geo.Jord.NVector++-- | An earth position defined by an horizontal position and height. +-- +-- horizontal position can be either a 'LatLong' or a 'NVector'. +data AngularPosition a = AngularPosition+ { pos :: a+ , height :: Length+ } deriving (Eq, Show)++-- | 'AngularPosition' from a 'LatLong' and height. +latLongHeight :: LatLong -> Length -> AngularPosition LatLong+latLongHeight = AngularPosition++-- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. +-- 'error's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +decimalLatLongHeight :: Double -> Double -> Length -> AngularPosition LatLong+decimalLatLongHeight lat lon = latLongHeight (decimalLatLong lat lon)++-- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. +-- A 'Left' indicates that the given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +decimalLatLongHeightE :: Double -> Double -> Length -> Either String (AngularPosition LatLong)+decimalLatLongHeightE lat lon h = fmap (`latLongHeight` h) (decimalLatLongE lat lon)++-- | 'AngularPosition' from given latitude and longitude in __decimal degrees__ and height. +-- 'fail's if given latitude is outisde [-90..90]° and/or +-- given longitude is outisde [-180..180]°. +decimalLatLongHeightF :: (MonadFail m) => Double -> Double -> Length -> m (AngularPosition LatLong)+decimalLatLongHeightF lat lon h = fmap (`latLongHeight` h) (decimalLatLongF lat lon)++-- | 'AngularPosition' from a 'NVector' and height. +nvectorHeight :: NVector -> Length -> AngularPosition NVector+nvectorHeight = AngularPosition
+ src/Data/Geo/Jord/Earth.hs view
@@ -0,0 +1,117 @@+-- | +-- Module: Data.Geo.Jord.Earth +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Ellipsoidal and derived spherical earth models. +-- +module Data.Geo.Jord.Earth+ ( Earth(..)+ , Ellipsoid(..)+ , eccentricity+ , meanRadius+ , polarRadius+ , spherical+ -- * Reference ellipsoids.+ , wgs84+ , grs80+ , wgs72+ -- * Spherical models dervived from reference ellipsoids.+ , s84+ , s80+ , s72+ , r84+ , r80+ , r72+ ) where++import Data.Geo.Jord.Length++-- | Earth model: ellipsoidal or spherical.+data Earth+ = Ellipsoidal Ellipsoid+ | Spherical Length+ deriving (Eq, Show)++-- | Primary ellipsoid parameters. +data Ellipsoid = Ellipsoid+ { equatorialRadius :: Length -- ^ equatorial radius or semi-major axis (a). + , inverseFlattening :: Double -- ^ inverse flattening. + } deriving (Eq, Show)++-- | Computes the eccentricity of the given 'Earth' model. +eccentricity :: Earth -> Double+eccentricity (Ellipsoidal e) = sqrt (1.0 - (b * b) / (a * a))+ where+ a = semiMajorAxis e+ b = semiMinorAxis a (inverseFlattening e)+eccentricity (Spherical _) = 0++-- | Computes the mean radius of the given 'Earth' model. +-- +-- This radius can be used for geodetic calculations assuming a spherical earth model. +meanRadius :: Earth -> Length+meanRadius (Ellipsoidal e) = metres ((2.0 * a + b) / 3.0)+ where+ a = semiMajorAxis e+ b = semiMinorAxis a (inverseFlattening e)+meanRadius (Spherical r) = r++-- | Computes the polar radius or semi-minor axis (b) of the given 'Earth' model. +polarRadius :: Earth -> Length+polarRadius (Ellipsoidal e) = metres (semiMinorAxis a f)+ where+ a = semiMajorAxis e+ f = inverseFlattening e+polarRadius (Spherical r) = r++-- | Spherical model derived from given model.+spherical :: Earth -> Earth+spherical e = Spherical (meanRadius e)++-- | World Geodetic System WGS84 ellipsoid. +wgs84 :: Earth+wgs84 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257223563))++-- | Geodetic Reference System 1980 ellipsoid. +grs80 :: Earth+grs80 = Ellipsoidal (Ellipsoid (metres 6378137.0) (1.0 / 298.257222101))++-- | World Geodetic System WGS72 ellipsoid. +wgs72 :: Earth+wgs72 = Ellipsoidal (Ellipsoid (metres 6378135.0) (1.0 / 298.26))++-- | Spherical earth model derived from 'wgs84'.+s84 :: Earth+s84 = spherical wgs84++-- | Spherical earth model derived from 'grs80'. +s80 :: Earth+s80 = spherical grs80++-- | Spherical earth model derived from 'wgs72'. +s72 :: Earth+s72 = spherical wgs72++-- | Mean earth radius derived from the 'wgs84' ellipsoid. +r84 :: Length+r84 = meanRadius s84++-- | Mean earth radius derived from the 'grs80' ellipsoid. +r80 :: Length+r80 = meanRadius s80++-- | Mean earth radius derived from the 'wgs72' ellipsoid. +r72 :: Length+r72 = meanRadius s72++-- | semi-major axis (a) in metres. +semiMajorAxis :: Ellipsoid -> Double+semiMajorAxis = toMetres . equatorialRadius++-- | Computes the polar semi-minor axis (b) from @a@ anf @f@. +semiMinorAxis :: Double -> Double -> Double+semiMinorAxis a f = a * (1.0 - f)
+ src/Data/Geo/Jord/EcefPosition.hs view
@@ -0,0 +1,59 @@+-- | +-- Module: Data.Geo.Jord.EcefPosition +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Earth Centred, Earth Fixed (ECEF) position. +-- +-- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> +-- +module Data.Geo.Jord.EcefPosition+ ( EcefPosition+ , ecef+ , ecefMetres+ , ex+ , ey+ , ez+ ) where++import Data.Geo.Jord.Length+import Data.Geo.Jord.Vector3d++-- | An earth position expressed in the Earth Centred, Earth Fixed (ECEF) coordinates system. +-- +-- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. +-- +-- Note: on a spherical model earth, an n-vector is equivalent to a normalised version of an (ECEF) cartesian coordinate. +newtype EcefPosition =+ EcefPosition Vector3d+ deriving (Eq, Show)++instance IsVector3d EcefPosition where+ vec (EcefPosition v) = v++-- | 'EcefPosition' from given x, y and z length. +-- +-- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. +ecef :: Length -> Length -> Length -> EcefPosition+ecef x y z = EcefPosition (Vector3d (toMetres x) (toMetres y) (toMetres z))++-- | 'EcefPosition' from given x, y and z length in _metres_. +-- +-- @ex-ey@ plane is the equatorial plane, @ey@ is on the prime meridian, and @ez@ on the polar axis. +ecefMetres :: Double -> Double -> Double -> EcefPosition+ecefMetres x y z = ecef (metres x) (metres y) (metres z)++-- | x coordinate of the given '$tc'EcefPosition'. +ex :: EcefPosition -> Length+ex (EcefPosition v) = metres (vx v)++-- | y coordinate of the given '$tc'EcefPosition'. +ey :: EcefPosition -> Length+ey (EcefPosition v) = metres (vy v)++-- | z coordinate of the given '$tc'EcefPosition'. +ez :: EcefPosition -> Length+ez (EcefPosition v) = metres (vz v)
− src/Data/Geo/Jord/Eval.hs
@@ -1,521 +0,0 @@-{-# 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.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 - = Ang Angle -- ^ 'Angle' - | AngDec Double -- ^ 'Angle' in decimal degrees - | Bool Bool -- ^ boolean - | Double Double -- ^ double - | Len Length -- ^ 'Length' - | Gc GreatCircle -- ^ 'GreatCircle' - | 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) - --- | '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 'LatLong'. --- --- 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 --- @ --- --- 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 Ll --- 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 (Ll (fromNVector v)) - Right (Vecs vs) -> Right (Lls (map fromNVector vs)) - oth -> oth - --- | All supported functions: --- --- * 'antipode' --- --- * 'crossTrackDistance' --- --- * 'decimalDegrees' --- --- * 'destination' --- --- * 'distance' --- --- * 'finalBearing' --- --- * 'greatCircle' --- --- * 'initialBearing' --- --- * 'interpolate' --- --- * 'intersections' --- --- * 'isInside' --- --- * 'mean' --- --- * 'latLong' --- --- * 'latLongDecimal' --- --- * 'readLatLong' --- --- * '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 - --- | @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 (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 - --- | 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 - | 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/Frames.hs view
@@ -0,0 +1,290 @@+{-# LANGUAGE FlexibleInstances #-} + +-- | +-- Module: Data.Geo.Jord.Frames +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Type and functions for working with delta vectors in different reference frames. +-- +-- 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> +-- +module Data.Geo.Jord.Frames + ( + -- * Reference Frames + Frame(..) + -- ** Body frame + , FrameB + , yaw + , pitch + , roll + , frameB + -- ** Local frame + , FrameL + , wanderAzimuth + , frameL + -- ** North-East-Down frame + , FrameN + , frameN + -- * Deltas + , Delta + , delta + , deltaMetres + -- * Delta in the north, east, down frame + , Ned + , ned + , nedMetres + , north + , east + , down + , bearing + , elevation + , norm + -- * Calculations + , deltaBetween + , nedBetween + , target + , targetN + ) where + +import Data.Geo.Jord.Angle +import Data.Geo.Jord.AngularPosition +import Data.Geo.Jord.Earth +import Data.Geo.Jord.EcefPosition +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Rotation +import Data.Geo.Jord.Transform +import Data.Geo.Jord.Vector3d + +-- | class for reference frames. +-- +-- Supported frames: +-- +-- * 'FrameB': 'rEF' returns R_EB +-- +-- * 'FrameL': 'rEF' returns R_EL +-- +-- * 'FrameN': 'rEF' returns R_EN +-- +class Frame a where + rEF :: a -> [Vector3d] -- ^ rotation matrix to transform vectors decomposed in frame @a@ to vectors decomposed Earth-Fixed frame. + +-- | Body frame (typically of a vehicle). +-- +-- * Position: The origin is in the vehicle’s reference point. +-- +-- * Orientation: The x-axis points forward, the y-axis to the right (starboard) and the z-axis +-- in the vehicle’s down direction. +-- +-- * Comments: The frame is fixed to the vehicle. +-- +data FrameB = + FrameB Angle + Angle + Angle + Vector3d + deriving (Eq, Show) + +-- | body yaw angle (vertical axis). +yaw :: FrameB -> Angle +yaw (FrameB a _ _ _) = a + +-- | body pitch angle (transverse axis). +pitch :: FrameB -> Angle +pitch (FrameB _ a _ _) = a + +-- | body roll angle (longitudinal axis). +roll :: FrameB -> Angle +roll (FrameB _ _ a _) = a + +-- | 'FrameB' from given yaw, pitch, roll, position (origin) and earth model. +frameB :: (ETransform a) => Angle -> Angle -> Angle -> a -> Earth -> FrameB +frameB yaw' pitch' roll' p e = FrameB yaw' pitch' roll' (nvec p e) + +-- | R_EB: frame B to Earth +instance Frame FrameB where + rEF (FrameB y p r o) = rm + where + rNB = zyx2r y p r + n = FrameN o + rEN = rEF n + rm = mdot rEN rNB -- closest frames cancel: N + +-- | Local level, Wander azimuth frame. +-- +-- * Position: The origin is directly beneath or above the vehicle (B), at Earth’s surface (surface +-- of ellipsoid model). +-- +-- * Orientation: The z-axis is pointing down. Initially, the x-axis points towards north, and the +-- y-axis points towards east, but as the vehicle moves they are not rotating about the z-axis +-- (their angular velocity relative to the Earth has zero component along the z-axis). +-- (Note: Any initial horizontal direction of the x- and y-axes is valid for L, but if the +-- initial position is outside the poles, north and east are usually chosen for convenience.) +-- +-- * Comments: The L-frame is equal to the N-frame except for the rotation about the z-axis, +-- which is always zero for this frame (relative to Earth). Hence, at a given time, the only +-- difference between the frames is an angle between the x-axis of L and the north direction; +-- this angle is called the wander azimuth angle. The L-frame is well suited for general +-- calculations, as it is non-singular. +-- +data FrameL = + FrameL Angle + LatLong + deriving (Eq, Show) + +-- | wander azimuth: angle between x-axis of the frame L and the north direction. +wanderAzimuth :: FrameL -> Angle +wanderAzimuth (FrameL a _) = a + +-- | R_EL: frame L to Earth +instance Frame FrameL where + rEF (FrameL w o) = rm + where + r = xyz2r (longitude o) (negate' (latitude o)) w + rEe' = [Vector3d 0 0 (-1), Vector3d 0 1 0, Vector3d 1 0 0] + rm = mdot rEe' r + +-- | 'FrameL' from given wander azimuth, position (origin) and earth model. +frameL :: (ETransform a) => Angle -> a -> Earth -> FrameL +frameL w p e = FrameL w ll + where + v = pos (ecefToNVector (toEcef p e) e) + ll = nvectorToLatLong v + +-- | North-East-Down (local level) frame. +-- +-- * Position: The origin is directly beneath or above the vehicle (B), at Earth’s surface (surface +-- of ellipsoid model). +-- +-- * Orientation: The x-axis points towards north, the y-axis points towards east (both are +-- horizontal), and the z-axis is pointing down. +-- +-- * Comments: When moving relative to the Earth, the frame rotates about its z-axis to allow the +-- x-axis to always point towards north. When getting close to the poles this rotation rate +-- will increase, being infinite at the poles. The poles are thus singularities and the direction of +-- the x- and y-axes are not defined here. Hence, this coordinate frame is not suitable for +-- general calculations. +-- +newtype FrameN = + FrameN Vector3d + deriving (Eq, Show) + +-- | R_EN: frame N to Earth +instance Frame FrameN where + rEF (FrameN o) = transpose rm + where + np = vec northPole + rd = vscale o (-1) -- down (pointing opposite to n-vector) + re = vunit (vcross np o) -- east (pointing perpendicular to the plane) + rn = vcross re rd -- north (by right hand rule) + rm = [rn, re, rd] + +-- | 'FrameN' from given position (origin) and earth model. +frameN :: (ETransform a) => a -> Earth -> FrameN +frameN p e = FrameN (nvec p e) + +-- | delta between position in one of the reference frames. +newtype Delta = + Delta Vector3d + deriving (Eq, Show) + +-- | 'Delta' from given x, y and z length. +delta :: Length -> Length -> Length -> Delta +delta x y z = Delta (Vector3d (toMetres x) (toMetres y) (toMetres z)) + +-- | 'Delta' from given x, y and z length in _metres_. +deltaMetres :: Double -> Double -> Double -> Delta +deltaMetres x y z = delta (metres x) (metres y) (metres z) + +-- | North, east and down delta (thus in frame 'FrameN'). +newtype Ned = + Ned Vector3d + deriving (Eq, Show) + +-- | 'Ned' from given north, east and down. +ned :: Length -> Length -> Length -> Ned +ned n e d = Ned (Vector3d (toMetres n) (toMetres e) (toMetres d)) + +-- | 'Ned' from given north, east and down in _metres_. +nedMetres :: Double -> Double -> Double -> Ned +nedMetres n e d = ned (metres n) (metres e) (metres d) + +-- | North component of given 'Ned'. +north :: Ned -> Length +north (Ned v) = metres (vx v) + +-- | East component of given 'Ned'. +east :: Ned -> Length +east (Ned v) = metres (vy v) + +-- | Down component of given 'Ned'. +down :: Ned -> Length +down (Ned v) = metres (vz v) + +-- | @bearing v@ computes the bearing in compass angle of the NED vector @v@ from north. +-- +-- Compass angles are clockwise angles from true north: 0 = north, 90 = east, 180 = south, 270 = west. +-- +bearing :: Ned -> Angle +bearing v = + let a = atan2' (toMetres (east v)) (toMetres (north v)) + in normalise a (decimalDegrees 360.0) + +-- | @elevation v@ computes the elevation of the NED vector @v@ from horizontal (ie tangent to ellipsoid surface). +elevation :: Ned -> Angle +elevation (Ned v) = negate' (asin' (vz v / vnorm v)) + +-- | @norm v@ computes the norm of the NED vector @v@. +norm :: Ned -> Length +norm (Ned v) = metres (vnorm v) + +-- | @deltaBetween p1 p2 f e@ computes the exact 'Delta' between the two positions @p1@ and @p2@ in frame @f@ +-- using earth model @e@. +deltaBetween :: (ETransform a, Frame c) => a -> a -> (a -> Earth -> c) -> Earth -> Delta +deltaBetween p1 p2 f e = deltaMetres (vx d) (vy d) (vz d) + where + e1 = ecefvec p1 e + e2 = ecefvec p2 e + de = vsub e2 e1 + -- rotation matrix to go from Earth Frame to Frame at p1 + rm = transpose (rEF (f p1 e)) + d = vrotate de rm + +-- | @nedBetween p1 p2 e@ computes the exact 'Ned' vector between the two positions @p1@ and @p2@, in north, east, and down +-- using earth model @e@. +-- +-- Produced 'Ned' delta is relative to @p1@: Due to the curvature of Earth and different directions to the North Pole, +-- the north, east, and down directions will change (relative to Earth) for different places. +-- +-- Position @p1@ must be outside the poles for the north and east directions to be defined. +nedBetween :: (ETransform a) => a -> a -> Earth -> Ned +nedBetween p1 p2 e = nedMetres (vx d) (vy d) (vz d) + where + (Delta d) = deltaBetween p1 p2 frameN e + +-- | @target p0 f d e@ computes the target position from position @p0@ and delta @d@ using in frame @f@ +-- and using earth model @e@. +target :: (ETransform a, Frame c) => a -> (a -> Earth -> c) -> Delta -> Earth -> a +target p0 f (Delta d) e = fromEcef (ecefMetres (vx e0 + vx c) (vy e0 + vy c) (vz e0 + vz c)) e + where + e0 = ecefvec p0 e + rm = rEF (f p0 e) + c = vrotate d rm + +-- | @targetN p0 d e@ computes the target position from position @p0@ and north, east, down @d@ using earth model @e@. +targetN :: (ETransform a) => a -> Ned -> Earth -> a +targetN p0 (Ned d) = target p0 frameN (Delta d) + +-- | ECEF position (as a 'Vector3d') from given position. +ecefvec :: (ETransform a) => a -> Earth -> Vector3d +ecefvec p m = vec (toEcef p m) + +-- | NVector (as a 'Vector3d') from given positon. +nvec :: (ETransform a) => a -> Earth -> Vector3d +nvec p e = vec (pos (ecefToNVector (toEcef p e) e))
+ src/Data/Geo/Jord/Geodetics.hs view
@@ -0,0 +1,317 @@+-- | +-- Module: Data.Geo.Jord.Geodetics +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Geodetic calculations assuming a __spherical__ earth model. +-- +-- 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> +-- +module Data.Geo.Jord.Geodetics + ( + -- * The 'GreatCircle' type + GreatCircle + -- * Smart constructors + , greatCircle + , greatCircleE + , greatCircleF + , greatCircleBearing + -- * Calculations + , angularDistance + , antipode + , crossTrackDistance + , crossTrackDistance84 + , destination + , destination84 + , finalBearing + , initialBearing + , interpolate + , intersections + , insideSurface + , mean + , surfaceDistance + , surfaceDistance84 + ) where + +import Control.Monad.Fail +import Data.Fixed +import Data.Geo.Jord.Angle +import Data.Geo.Jord.AngularPosition +import Data.Geo.Jord.Earth (r84) +import Data.Geo.Jord.LatLong +import Data.Geo.Jord.Length +import Data.Geo.Jord.NVector +import Data.Geo.Jord.Quantity +import Data.Geo.Jord.Transform +import Data.Geo.Jord.Vector3d +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 :: Vector3d + , dscr :: String + } deriving (Eq) + +instance Show GreatCircle where + show = dscr + +-- | 'GreatCircle' passing by both given positions. 'error's if given positions are +-- equal or antipodal. +greatCircle :: (NTransform 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) + +-- | 'GreatCircle' passing by both given positions. A 'Left' indicates that given positions are +-- equal or antipodal. +greatCircleE :: (NTransform a) => a -> a -> Either String GreatCircle +greatCircleE p1 p2 + | v1 == v2 = Left "Invalid Great Circle: positions are equal" + | (realToFrac (vnorm (vadd v1 v2)) :: Nano) == 0 = + Left "Invalid Great Circle: positions are antipodal" + | otherwise = + Right (GreatCircle (vcross v1 v2) ("passing by " ++ show (ll p1) ++ " & " ++ show (ll p2))) + where + v1 = vector3d p1 + v2 = vector3d p2 + +-- | 'GreatCircle' passing by both given positions. 'fail's if given positions are +-- equal or antipodal. +greatCircleF :: (NTransform a, MonadFail m) => 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 :: (NTransform a) => a -> Angle -> GreatCircle +greatCircleBearing p b = + GreatCircle (vsub n' e') ("passing by " ++ show (ll p) ++ " heading on " ++ show b) + where + v = vector3d p + e = vcross (vec northPole) v -- easting + n = vcross v e -- northing + e' = vscale e (cos' b / vnorm e) + n' = vscale n (sin' b / vnorm n) + +-- | @angularDistance p1 p2 n@ computes the angle between the horizontal positions @p1@ and @p2@. +-- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], +-- signed + if @p1@ is clockwise looking along @n@, - in opposite direction. +angularDistance :: (NTransform a) => a -> a -> Maybe a -> Angle +angularDistance p1 p2 n = angularDistance' v1 v2 vn + where + v1 = vector3d p1 + v2 = vector3d p2 + vn = fmap vector3d n + +-- | @antipode p@ computes the antipodal horizontal position of @p@: +-- the horizontal position on the surface of the Earth which is diametrically opposite to @p@. +antipode :: (NTransform a) => a -> a +antipode p = fromNVector (angular (vscale (vector3d nv) (-1.0)) h) + where + (AngularPosition nv h) = toNVector p + +-- | @crossTrackDistance p gc@ computes the signed distance horizontal position @p@ to great circle @gc@. +-- 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 (decimalLatLong 51 0) (decimalLatLong 52 1) +-- let gc2 = greatCircle (decimalLatLong 52 1) (decimalLatLong 51 0) +-- crossTrackDistance p gc1 == (- crossTrackDistance p gc2) +-- @ +crossTrackDistance :: (NTransform a) => a -> GreatCircle -> Length -> Length +crossTrackDistance p gc = + arcLength (sub (angularDistance' (normal gc) (vector3d p) Nothing) (decimalDegrees 90)) + +-- | 'crossTrackDistance' using the mean radius of the WGS84 reference ellipsoid. +crossTrackDistance84 :: (NTransform a) => a -> GreatCircle -> Length +crossTrackDistance84 p gc = crossTrackDistance p gc r84 + +-- | @destination p b d r@ computes the destination position from position @p@ having +-- travelled the distance @d@ on the initial bearing (compass angle) @b@ (bearing will normally vary +-- before destination is reached) and using the earth radius @r@. +destination :: (NTransform a) => a -> Angle -> Length -> Length -> a +destination p b d r + | toMetres d == 0.0 = p + | otherwise = fromNVector (angular vd h) + where + (AngularPosition nv h) = toNVector p + v = vec nv + ed = vunit (vcross (vec northPole) v) -- east direction vector at v + nd = vcross v ed -- north direction vector at v + ta = central d r -- central angle + de = vadd (vscale nd (cos' b)) (vscale ed (sin' b)) -- vunit vector in the direction of the azimuth + vd = vadd (vscale v (cos' ta)) (vscale de (sin' ta)) + +-- | 'destination' using the mean radius of the WGS84 reference ellipsoid. +destination84 :: (NTransform a) => a -> Angle -> Length -> a +destination84 p b d = destination p b d r84 + +-- | @finalBearing p1 p2@ computes the final bearing arriving at @p2@ from @p1@ in compass angle. +-- +-- Compass angles are clockwise angles from true north: 0 = north, 90 = east, 180 = south, 270 = west. +-- +-- The final bearing will differ from the 'initialBearing' by varying degrees according to distance and latitude. +-- +-- Returns 'Nothing' if both horizontal positions are equals. +finalBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle +finalBearing p1 p2 = fmap (\b -> normalise b (decimalDegrees 180)) (initialBearing p2 p1) + +-- | @initialBearing p1 p2@ computes the initial bearing from @p1@ to @p2@ in compass angle. +-- +-- Compass angles are clockwise angles from true north: 0 = north, 90 = east, 180 = south, 270 = west. +-- +-- Returns 'Nothing' if both horizontal positions are equals. +initialBearing :: (Eq a, NTransform a) => a -> a -> Maybe Angle +initialBearing p1 p2 + | p1 == p2 = Nothing + | otherwise = Just (normalise (angularDistance' gc1 gc2 (Just v1)) (decimalDegrees 360)) + where + v1 = vector3d p1 + v2 = vector3d p2 + gc1 = vcross v1 v2 -- great circle through p1 & p2 + gc2 = vcross v1 (vec northPole) -- great circle through p1 & north pole + +-- | @interpolate p0 p1 f# computes the horizontal position at fraction @f@ between the @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 :: (NTransform 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 (angular iv ih) + where + (AngularPosition nv0 h0) = toNVector p0 + (AngularPosition nv1 h1) = toNVector p1 + v0 = vec nv0 + v1 = vec nv1 + iv = vunit (vadd v0 (vscale (vsub v1 v0) f)) + ih = lrph h0 h1 f + +-- | 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 :: (NTransform a) => GreatCircle -> GreatCircle -> Maybe (a, a) +intersections gc1 gc2 + | (vnorm i :: Double) == 0.0 = Nothing + | otherwise + , let ni = fromNVector (angular (vunit i) zero) = Just (ni, antipode ni) + where + i = vcross (normal gc1) (normal gc2) + +-- | @insideSurface p ps@ determines whether the @p@ is inside the polygon defined by the list of positions @ps@. +-- 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 @ps@ does not at least defines a triangle. +-- +insideSurface :: (Eq a, NTransform a) => a -> [a] -> Bool +insideSurface p ps + | null ps = False + | head ps == last ps = insideSurface p (init ps) + | length ps < 3 = False + | otherwise = + let aSum = + foldl + (\a v' -> add a (uncurry angularDistance' v' (Just v))) + (decimalDegrees 0) + (egdes (map (vsub v) vs)) + in abs (toDecimalDegrees aSum) > 180.0 + where + v = vector3d p + vs = fmap vector3d ps + +-- | @mean ps@ computes the mean geographic horitzontal position of @ps@, if it is defined. +-- +-- The geographic mean is not defined for antipodals position (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 :: (NTransform a) => [a] -> Maybe a +mean [] = Nothing +mean [p] = Just p +mean ps = + if null antipodals + then Just (fromNVector (angular (vunit (foldl vadd vzero vs)) zero)) + else Nothing + where + vs = fmap vector3d ps + ts = filter (\l -> length l == 2) (subsequences vs) + antipodals = + filter (\t -> (realToFrac (vnorm (vadd (head t) (last t)) :: Double) :: Nano) == 0) ts + +-- | @surfaceDistance p1 p2@ computes the surface distance (length of geodesic) between the positions @p1@ and @p2@. +surfaceDistance :: (NTransform a) => a -> a -> Length -> Length +surfaceDistance p1 p2 = arcLength (angularDistance p1 p2 Nothing) + +-- | 'surfaceDistance' using the mean radius of the WGS84 reference ellipsoid. +surfaceDistance84 :: (NTransform a) => a -> a -> Length +surfaceDistance84 p1 p2 = surfaceDistance p1 p2 r84 + +-- | Angle between the two given n-vectors. +-- 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' :: Vector3d -> Vector3d -> Maybe Vector3d -> Angle +angularDistance' v1 v2 n = atan2' sinO cosO + where + sign = maybe 1 (signum . vdot (vcross v1 v2)) n + sinO = sign * vnorm (vcross v1 v2) + cosO = vdot v1 v2 + +-- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] +egdes :: [Vector3d] -> [(Vector3d, Vector3d)] +egdes ps = zip ps (tail ps ++ [head ps]) + +lrph :: Length -> Length -> Double -> Length +lrph h0 h1 f = metres h + where + h0' = toMetres h0 + h1' = toMetres h1 + h = h0' + (h1' - h0') * f + +vector3d :: (NTransform a) => a -> Vector3d +vector3d = vec . pos . toNVector + +angular :: Vector3d -> Length -> AngularPosition NVector +angular v = nvectorHeight (nvector (vx v) (vy v) (vz v)) + +ll :: (NTransform a) => a -> LatLong +ll = nvectorToLatLong . pos . toNVector
− src/Data/Geo/Jord/GreatCircle.hs
@@ -1,133 +0,0 @@--- | --- 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. --- -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. --- --- 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 - --- | '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: --- --- @ --- 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 (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)
src/Data/Geo/Jord/LatLong.hs view
@@ -1,12 +1,12 @@ -- | --- Module: Data.Geo.Jord.GeoPos +-- Module: Data.Geo.Jord.LatLong -- 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. +-- Geodetic latitude and longitude. -- module Data.Geo.Jord.LatLong ( @@ -16,9 +16,9 @@ , latLong , latLongE , latLongF - , latLongDecimal - , latLongDecimalE - , latLongDecimalF + , decimalLatLong + , decimalLatLongE + , decimalLatLongF -- * read , readLatLong , readLatLongE @@ -37,7 +37,7 @@ import Text.ParserCombinators.ReadP import Text.Read hiding (pfail) --- | A geographic position (latitude and longitude). +-- | Horizontal position defined by its geodetic latitude and longitude. data LatLong = LatLong { latitude :: Angle , longitude :: Angle @@ -82,23 +82,23 @@ where e = latLongE lat lon --- | 'LatLong' from given latitude and longitude in decimal degrees. +-- | '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) +decimalLatLong :: Double -> Double -> LatLong +decimalLatLong lat lon = latLong (decimalDegrees lat) (decimalDegrees lon) --- | 'LatLong' from given latitude and longitude in decimal degrees. +-- | '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) +decimalLatLongE :: Double -> Double -> Either String LatLong +decimalLatLongE lat lon = latLongE (decimalDegrees lat) (decimalDegrees lon) --- | 'LatLong' from given latitude and longitude in decimal degrees. +-- | '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) +decimalLatLongF :: (MonadFail m) => Double -> Double -> m LatLong +decimalLatLongF lat lon = latLongF (decimalDegrees lat) (decimalDegrees lon) -- | Obtains a 'LatLong' from the given string formatted as either: --
src/Data/Geo/Jord/Length.hs view
@@ -8,114 +8,120 @@ -- -- Types and functions for working with (signed) lengths in metres, kilometres or nautical miles. -- -module Data.Geo.Jord.Length - ( - -- * The 'Length' type - Length(millimetres) +module Data.Geo.Jord.Length+ (+ -- * The 'Length' type+ Length(millimetres) -- * Smart constructors - , kilometres - , metres - , nauticalMiles + , feet+ , kilometres+ , metres+ , nauticalMiles -- * Read - , readLength - , readLengthE - , readLengthF + , 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) - + , toFeet+ , toKilometres+ , toMetres+ , toNauticalMiles+ ) 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) - +newtype Length = Length+ { millimetres :: Int+ } deriving (Eq)+ -- | See 'readLength'. -instance Read Length where - readsPrec _ = readP_to_S length - +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 - +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) - +instance Quantity Length where+ add a b = Length (millimetres a + millimetres b)+ sub a b = Length (millimetres a - millimetres b)+ zero = Length 0++-- | 'Length' from given amount of feet.+feet :: Double -> Length+feet ft = metres (ft * 0.3048)++-- | 'Length' from given amount of kilometres.+kilometres :: Double -> Length+kilometres km = metres (km * 1000.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) - +metres :: Double -> Length+metres m = Length (round (m * 1000.0))++-- | 'Length' from given amount of nautical miles.+nauticalMiles :: Double -> Length+nauticalMiles nm = metres (nm * 1852.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 - +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 - +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 - +readLengthF :: (MonadFail m) => String -> m Length+readLengthF s =+ let p = readEither s+ in case p of+ Left e -> fail e+ Right l -> return l++-- | @toFeet l@ converts @l@ to feet.+toFeet :: Length -> Double+toFeet l = toMetres l / 0.3048+ -- | @toKilometres l@ converts @l@ to kilometres. -toKilometres :: Length -> Double -toKilometres l = toMetres l / 1000.0 - +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 - +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 - +toNauticalMiles :: Length -> Double+toNauticalMiles l = toMetres l / 1852.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 +length :: ReadP Length+length = do+ v <- number+ skipSpaces+ u <- string "m" <|> string "km" <|> string "Nm" <|> string "ft"+ case u of+ "m" -> return (metres v)+ "km" -> return (kilometres v)+ "Nm" -> return (nauticalMiles v)+ "ft" -> return (feet v)+ _ -> pfail
src/Data/Geo/Jord/NVector.hs view
@@ -8,73 +8,34 @@ -- -- 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 - +module Data.Geo.Jord.NVector+ ( NVector+ , nvector+ , northPole+ , southPole+ ) where++import Data.Geo.Jord.Vector3d+ -- | 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 +-- +-- Orientation: z-axis points to the North Pole along the Earth's rotation axis, +-- x-axis points towards the point where latitude = longitude = 0. +newtype NVector =+ NVector Vector3d+ deriving (Eq, Show)++instance IsVector3d NVector where+ vec (NVector v) = v++-- | Unit 'NVector' from given x, y and z. +nvector :: Double -> Double -> Double -> NVector+nvector x y z = NVector (vunit (Vector3d x y z))++-- | Horizontal position of the North Pole. +northPole :: NVector+northPole = NVector (Vector3d 0.0 0.0 1.0)++-- | Horizontal position of the North Pole. +southPole :: NVector+southPole = NVector (Vector3d 0.0 0.0 (-1.0))
− src/Data/Geo/Jord/Position.hs
@@ -1,224 +0,0 @@--- | --- 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))
src/Data/Geo/Jord/Quantity.hs view
@@ -1,17 +1,18 @@--- |--- 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 +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Classes for working with quantities. +-- module Data.Geo.Jord.Quantity ( Quantity(..) ) where --- | Something that can be added or subtracted.-class Quantity a where+-- | Something that can be added or subtracted. +class (Eq a) => Quantity a where add, sub :: a -> a -> a+ zero :: a
+ src/Data/Geo/Jord/Rotation.hs view
@@ -0,0 +1,138 @@+-- | +-- Module: Data.Geo.Jord.Rotation +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Rotation matrices from/to 3 angles about new axes +-- +-- 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> +-- +module Data.Geo.Jord.Rotation + ( r2xyz + , r2zyx + , xyz2r + , zyx2r + ) where + +import Data.Geo.Jord.Angle +import Data.Geo.Jord.Vector3d + +-- | Angles about new axes in the xyz-order from a rotation matrix. +-- +-- The produced list contains 3 'Angle's of rotation about new axes. +-- +-- The x, y, z angles are called Euler angles or Tait-Bryan angles and are +-- defined by the following procedure of successive rotations: +-- Given two arbitrary coordinate frames A and B. Consider a temporary frame +-- T that initially coincides with A. In order to make T align with B, we +-- first rotate T an angle x about its x-axis (common axis for both A and T). +-- Secondly, T is rotated an angle y about the NEW y-axis of T. Finally, T +-- is rotated an angle z about its NEWEST z-axis. The final orientation of +-- T now coincides with the orientation of B. +-- The signs of the angles are given by the directions of the axes and the +-- right hand rule. +r2xyz :: [Vector3d] -> [Angle] +r2xyz [v0, v1, v2] = [x, y, z] + where + v00 = vx v0 + v01 = vy v0 + v12 = vz v1 + v22 = vz v2 + z = atan2' (-v01) v00 + x = atan2' (-v12) v22 + sy = vz v0 + -- cos y is based on as many elements as possible, to average out + -- numerical errors. It is selected as the positive square root since + -- y: [-pi/2 pi/2] + cy = sqrt ((v00 * v00 + v01 * v01 + v12 * v12 + v22 * v22) / 2.0) + y = atan2' sy cy +r2xyz m = error ("Invalid rotation matrix " ++ show m) + +-- | Angles about new axes in the xyz-order from a rotation matrix. +-- +-- The produced list contains 3 'Angle's of rotation about new axes. +-- The z, x, y angles are called Euler angles or Tait-Bryan angles and are +-- defined by the following procedure of successive rotations: +-- Given two arbitrary coordinate frames A and B. Consider a temporary frame +-- T that initially coincides with A. In order to make T align with B, we +-- first rotate T an angle z about its z-axis (common axis for both A and T). +-- Secondly, T is rotated an angle y about the NEW y-axis of T. Finally, T +-- is rotated an angle x about its NEWEST x-axis. The final orientation of +-- T now coincides with the orientation of B. +-- The signs of the angles are given by the directions of the axes and the +-- right hand rule. +-- Note that if A is a north-east-down frame and B is a body frame, we +-- have that z=yaw, y=pitch and x=roll. +r2zyx :: [Vector3d] -> [Angle] +r2zyx m = [z, y, x] + where + [x, y, z] = fmap negate' (r2xyz (transpose m)) + +-- | Rotation matrix (direction cosine matrix) from 3 angles about new axes in the xyz-order. +-- +-- The produced (no unit) rotation matrix is such +-- that the relation between a vector v decomposed in A and B is given by: +-- @v_A = mdot R_AB v_B@ +-- +-- The rotation matrix R_AB is created based on 3 angles x,y,z about new axes +-- (intrinsic) in the order x-y-z. The angles are called Euler angles or +-- Tait-Bryan angles and are defined by the following procedure of successive +-- rotations: +-- Given two arbitrary coordinate frames A and B. Consider a temporary frame +-- T that initially coincides with A. In order to make T align with B, we +-- first rotate T an angle x about its x-axis (common axis for both A and T). +-- Secondly, T is rotated an angle y about the NEW y-axis of T. Finally, T +-- is rotated an angle z about its NEWEST z-axis. The final orientation of +-- T now coincides with the orientation of B. +-- The signs of the angles are given by the directions of the axes and the +-- right hand rule. +xyz2r :: Angle -> Angle -> Angle -> [Vector3d] +xyz2r x y z = [v1, v2, v3] + where + cx = cos' x + sx = sin' x + cy = cos' y + sy = sin' y + cz = cos' z + sz = sin' z + v1 = Vector3d (cy * cz) ((-cy) * sz) sy + v2 = Vector3d (sy * sx * cz + cx * sz) ((-sy) * sx * sz + cx * cz) ((-cy) * sx) + v3 = Vector3d ((-sy) * cx * cz + sx * sz) (sy * cx * sz + sx * cz) (cy * cx) + +-- | rotation matrix (direction cosine matrix) from 3 angles about new axes in the zyx-order. +-- +-- The produced (no unit) rotation matrix is such +-- that the relation between a vector v decomposed in A and B is given by: +-- @v_A = mdot R_AB v_B@ +-- +-- The rotation matrix R_AB is created based on 3 angles +-- z,y,x about new axes (intrinsic) in the order z-y-x. The angles are called +-- Euler angles or Tait-Bryan angles and are defined by the following +-- procedure of successive rotations: +-- Given two arbitrary coordinate frames A and B. Consider a temporary frame +-- T that initially coincides with A. In order to make T align with B, we +-- first rotate T an angle z about its z-axis (common axis for both A and T). +-- Secondly, T is rotated an angle y about the NEW y-axis of T. Finally, T +-- is rotated an angle x about its NEWEST x-axis. The final orientation of +-- T now coincides with the orientation of B. +-- The signs of the angles are given by the directions of the axes and the +-- right hand rule. +-- +-- Note that if A is a north-east-down frame and B is a body frame, we +-- have that z=yaw, y=pitch and x=roll. +zyx2r :: Angle -> Angle -> Angle -> [Vector3d] +zyx2r z y x = [v1, v2, v3] + where + cx = cos' x + sx = sin' x + cy = cos' y + sy = sin' y + cz = cos' z + sz = sin' z + v1 = Vector3d (cz * cy) ((-sz) * cx + cz * sy * sx) (sz * sx + cz * sy * cx) + v2 = Vector3d (sz * cy) (cz * cx + sz * sy * sx) ((-cz) * sx + sz * sy * cx) + v3 = Vector3d (-sy) (cy * sx) (cy * cx)
+ src/Data/Geo/Jord/Transform.hs view
@@ -0,0 +1,192 @@+{-# LANGUAGE FlexibleInstances #-}++-- | +-- Module: Data.Geo.Jord.Transform +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- Transformations between coordinates systems both in spherical and ellipsoidal form. +-- +-- 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> +-- +-- See <http://clynchg3c.com/Technote/geodesy/coorddef.pdf Earth Coordinates> +-- +module Data.Geo.Jord.Transform+ ( NTransform(..)+ , ETransform(..)+ , nvectorToLatLong+ , latLongToNVector+ , ecefToNVector+ , nvectorToEcef+ , geodeticHeight+ ) where++import Data.Geo.Jord.Angle+import Data.Geo.Jord.AngularPosition+import Data.Geo.Jord.Earth+import Data.Geo.Jord.EcefPosition+import Data.Geo.Jord.LatLong+import Data.Geo.Jord.Length+import Data.Geo.Jord.NVector+import Data.Geo.Jord.Quantity+import Data.Geo.Jord.Vector3d++-- | Transformation between positions and 'AngularPosition' of 'NVector'. +class NTransform a where+ toNVector :: a -> AngularPosition NVector -- ^ position to 'AngularPosition' of 'NVector'. + fromNVector :: AngularPosition NVector -> a -- ^ 'AngularPosition' of 'NVector' and height to position. ++-- | 'NVector' <-> 'AngularPosition' of 'NVector'.+instance NTransform NVector where+ toNVector nv = AngularPosition nv zero+ fromNVector = pos++-- | 'LatLong' <-> 'AngularPosition' of 'NVector'. +instance NTransform LatLong where+ toNVector ll = AngularPosition (latLongToNVector ll) zero+ fromNVector = nvectorToLatLong . pos++-- | 'NTransform' identity. +instance NTransform (AngularPosition NVector) where+ toNVector = id+ fromNVector = id++-- | 'AngularPosition' of 'LatLong' <-> 'AngularPosition' of 'NVector'. +instance NTransform (AngularPosition LatLong) where+ toNVector (AngularPosition ll h) = AngularPosition (latLongToNVector ll) h+ fromNVector (AngularPosition nv h) = AngularPosition (nvectorToLatLong nv) h++-- | Transformation between 'EcefPosition' and angular or n-vector positions. +class ETransform a where+ toEcef :: a -> Earth -> EcefPosition -- ^ position and earth model to to 'EcefPosition'. + fromEcef :: EcefPosition -> Earth -> a -- ^ 'EcefPosition' and earth model to position. ++-- | 'NVector' <-> 'EcefPosition'. +instance ETransform NVector where+ fromEcef p e = pos (ecefToNVector p e)+ toEcef v = nvectorToEcef (nvectorHeight v zero)++-- | 'LatLong' <-> 'EcefPosition'. +instance ETransform LatLong where+ fromEcef p e = fromNVector (nvectorHeight (fromEcef p e :: NVector) zero)+ toEcef = toEcef . toNVector++-- | 'AngularPosition' of 'NVector' <-> 'EcefPosition'. +instance ETransform (AngularPosition NVector) where+ fromEcef = ecefToNVector+ toEcef = nvectorToEcef++-- | 'AngularPosition' of 'LatLong' <-> 'EcefPosition'. +instance ETransform (AngularPosition LatLong) where+ fromEcef p e = fromNVector (ecefToNVector p e)+ toEcef = nvectorToEcef . toNVector++-- | 'ETransform' identity.+instance ETransform EcefPosition where+ fromEcef p _ = p+ toEcef p _ = p++-- | @nvectorToLatLong v@ transforms 'NVector' @v@ to an equivalent 'LatLong'. +-- +-- Same as 'toNVector'. +nvectorToLatLong :: NVector -> LatLong+nvectorToLatLong nv = latLong lat lon+ where+ v = vec nv+ lat = atan2' (vz v) (sqrt (vx v * vx v + vy v * vy v))+ lon = atan2' (vy v) (vx v)++-- | @latLongToNVector ll@ transforms 'LatLong' @ll@ to an equivalent 'NVector'. +-- +-- See also 'fromNVector'. +latLongToNVector :: LatLong -> NVector+latLongToNVector ll = nvector x' y' z'+ where+ lat = latitude ll+ lon = longitude ll+ cl = cos' lat+ x' = cl * cos' lon+ y' = cl * sin' lon+ z' = sin' lat++-- | @ecefToNVector p e@ transforms 'EcefPosition' @p@ to an equivalent 'NVector' and geodetic height +-- using earth model @e@. +-- +-- See also 'fromEcef' +ecefToNVector :: EcefPosition -> Earth -> AngularPosition NVector+-- Ellipsoidal+ecefToNVector ep e@(Ellipsoidal el) = nvectorHeight (nvecEllipsoidal d e2 k px py pz) (metres h)+ where+ ev = vec ep+ e' = eccentricity e+ e2 = e' * e'+ e4 = e2 * e2+ a = toMetres (equatorialRadius el)+ a2 = a * a+ px = vx ev+ py = vy ev+ pz = vz ev+ p = (px * px + py * py) / a2+ q = ((1 - e2) / a2) * (pz * pz)+ r = (p + q - e4) / 6.0+ s = (e4 * p * q) / (4.0 * r * r * r)+ t = (1.0 + s + sqrt (s * (2.0 + s))) ** (1 / 3)+ u = r * (1.0 + t + 1.0 / t)+ v = sqrt (u * u + q * e4)+ w = e2 * (u + v - q) / (2.0 * v)+ k = sqrt (u + v + w * w) - w+ d = k * sqrt (px * px + py * py) / (k + e2)+ h = ((k + e2 - 1.0) / k) * sqrt (d * d + pz * pz)+-- Spherical+ecefToNVector p (Spherical r) = nvectorHeight (nvector (vx nv) (vy nv) (vz nv)) h+ where+ ev = vec p+ nv = vunit ev+ h = sub (metres (vnorm ev)) r++nvecEllipsoidal :: Double -> Double -> Double -> Double -> Double -> Double -> NVector+nvecEllipsoidal d e2 k px py pz = nvector nx' ny' nz'+ where+ s = 1.0 / sqrt (d * d + pz * pz)+ a = k / (k + e2)+ nx' = s * a * px+ ny' = s * a * py+ nz' = s * pz++-- | @nvectorToEcef (n, h) e@ transforms 'NVector' @n@ and geodetic height @h@ +-- to an equivalent 'EcefPosition' using earth model @e@. +-- +-- See also 'toEcef' +nvectorToEcef :: AngularPosition NVector -> Earth -> EcefPosition+-- Ellipsoidal+nvectorToEcef (AngularPosition nv h) e@(Ellipsoidal el) = ecef ex' ey' ez'+ where+ v = vec nv+ uv = vunit v+ a = toMetres (equatorialRadius el)+ b = toMetres (polarRadius e)+ nx' = vx uv+ ny' = vy uv+ nz' = vz uv+ m = (a * a) / (b * b)+ n = b / sqrt ((nx' * nx' * m) + (ny' * ny' * m) + (nz' * nz'))+ h' = toMetres h+ ex' = metres (n * m * nx' + h' * nx')+ ey' = metres (n * m * ny' + h' * ny')+ ez' = metres (n * nz' + h' * nz')+-- Spherical+nvectorToEcef (AngularPosition nv h) (Spherical r) = ecefMetres (vx ev) (vy ev) (vz ev)+ where+ unv = vunit . vec $ nv+ n = add h r+ ev = vscale unv (toMetres n)++-- | @geodeticHeight p e@ computes the geodetic height of 'EcefPosition' @p@ using earth model @e@. +-- +-- The geodetic height (or ellipsoidal height) is __not__ the mean sea level (MSL) height. +geodeticHeight :: EcefPosition -> Earth -> Length+geodeticHeight p e = height (ecefToNVector p e)
+ src/Data/Geo/Jord/Vector3d.hs view
@@ -0,0 +1,124 @@+-- | +-- Module: Data.Geo.Jord.Vector3d +-- Copyright: (c) 2018 Cedric Liegeois +-- License: BSD3 +-- Maintainer: Cedric Liegeois <ofmooseandmen@yahoo.fr> +-- Stability: experimental +-- Portability: portable +-- +-- 3-element vectors. +-- +module Data.Geo.Jord.Vector3d + ( Vector3d(..) + , IsVector3d(..) + , vadd + , vsub + , vdot + , vnorm + , vcross + , vrotate + , vscale + , vunit + , vzero + , transpose + , mdot + ) where + +-- | 3-element vector. +data Vector3d = Vector3d + { vx :: Double + , vy :: Double + , vz :: Double + } deriving (Eq, Show) + +-- | class for data types assimilated to 'Vector3d'. +class IsVector3d a where + vec :: a -> Vector3d + +-- | Adds 2 vectors. +vadd :: Vector3d -> Vector3d -> Vector3d +vadd v1 v2 = Vector3d x y z + where + x = vx v1 + vx v2 + y = vy v1 + vy v2 + z = vz v1 + vz v2 + +-- | Subtracts 2 vectors. +vsub :: Vector3d -> Vector3d -> Vector3d +vsub v1 v2 = Vector3d x y z + where + x = vx v1 - vx v2 + y = vy v1 - vy v2 + z = vz v1 - vz v2 + +-- | Computes the cross product of 2 vectors: the vector perpendicular to given vectors. +vcross :: Vector3d -> Vector3d -> Vector3d +vcross v1 v2 = Vector3d x y z + where + x = vy v1 * vz v2 - vz v1 * vy v2 + y = vz v1 * vx v2 - vx v1 * vz v2 + z = vx v1 * vy v2 - vy v1 * vx v2 + +-- | Computes the dot product of 2 vectors. +vdot :: Vector3d -> Vector3d -> Double +vdot v1 v2 = vx v1 * vx v2 + vy v1 * vy v2 + vz v1 * vz v2 + +-- | Computes the norm of a vector. +vnorm :: Vector3d -> Double +vnorm v = sqrt (x * x + y * y + z * z) + where + x = vx v + y = vy v + z = vz v + +-- | @vrotate v rm@ applies rotation matrix @rm@ to @v@. +vrotate :: Vector3d -> [Vector3d] -> Vector3d +vrotate v rm + | length rm /= 3 = error ("Invalid rotation matrix" ++ show rm) + | otherwise = Vector3d x y z + where + [x, y, z] = map (vdot v) rm + +-- | @vscale v s@ multiplies each component of @v@ by @s@. +vscale :: Vector3d -> Double -> Vector3d +vscale v s = Vector3d x y z + where + x = vx v * s + y = vy v * s + z = vz v * s + +-- | Normalises a vector. The 'vnorm' of the produced vector is @1@. +vunit :: Vector3d -> Vector3d +vunit v + | s == 1.0 = v + | otherwise = vscale v s + where + s = 1.0 / vnorm v + +-- | vector of vnorm 0. +vzero :: Vector3d +vzero = Vector3d 0 0 0 + +-- | transpose __square__ matrix made of 'Vector3d'. +transpose :: [Vector3d] -> [Vector3d] +transpose m = fmap ds2v (transpose' xs) + where + xs = fmap v2ds m + +-- | transpose matrix. +transpose' :: [[Double]] -> [[Double]] +transpose' ([]:_) = [] +transpose' x = map head x : transpose' (map tail x) + +-- | multiplies 2 matrices of 'Vector3d'. +mdot :: [Vector3d] -> [Vector3d] -> [Vector3d] +mdot a b = fmap ds2v [[vdot ar bc | bc <- transpose b] | ar <- a] + +-- | 'Vector3d' to list of doubles. +v2ds :: Vector3d -> [Double] +v2ds (Vector3d x' y' z') = [x', y', z'] + +-- | list of doubles to 'Vector3d'. +ds2v :: [Double] -> Vector3d +ds2v [x', y', z'] = Vector3d x' y' z' +ds2v xs = error ("Invalid list: " ++ show xs)
test/Data/Geo/Jord/AngleSpec.hs view
@@ -1,89 +1,92 @@-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 +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.000\"" $+ show (decimalDegrees 59.99999999999999) `shouldBe` "60°0'0.000\""+ it "shows 154.915 as 154°54'54.000\"" $+ show (decimalDegrees 154.915) `shouldBe` "154°54'54.000\""+ it "shows -154.915 as -154°54'54.000\"" $+ show (decimalDegrees (-154.915)) `shouldBe` "-154°54'54.000\""+ it "show 0.5245 as 0°31'28.800\"" $ show (decimalDegrees 0.5245) `shouldBe` "0°31'28.200\""+ it "show -0.5245 as -0°31'28.800\"" $+ show (decimalDegrees (-0.5245)) `shouldBe` "-0°31'28.200\""+ 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)) (meanRadius wgs84) `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)) (meanRadius wgs84) `shouldBe` metres 0.031+ it "arc length with central angle of 0.5 milliseconds == 0" $+ arcLength (decimalDegrees (0.4 / 3600000.0)) (meanRadius wgs84) `shouldBe` metres 0
+ test/Data/Geo/Jord/EarthSpec.hs view
@@ -0,0 +1,27 @@+module Data.Geo.Jord.EarthSpec+ ( spec+ ) where++import Data.Geo.Jord+import Test.Hspec++spec :: Spec+spec = do+ describe "Eccentricity" $ do+ it "returns 0.08181919084262157 for the WGS84 ellipsoid" $+ eccentricity wgs84 `shouldBe` 0.08181919084262157+ it "returns 0.08181919104281514 for the GRS80 ellipsoid" $+ eccentricity grs80 `shouldBe` 0.08181919104281514+ it "returns 0.08181881066274845 for the WG72 ellipsoid" $+ eccentricity wgs72 `shouldBe` 0.08181881066274845+ describe "Polar radius" $ do+ it "returns 6356752.314 m for the WGS84 ellipsoid" $+ polarRadius wgs84 `shouldBe` metres 6356752.314+ it "returns 6356752.314 m for the GRS80 ellipsoid" $+ polarRadius grs80 `shouldBe` metres 6356752.314+ it "returns 6356750.52 m for the WG72 ellipsoid" $+ polarRadius wgs72 `shouldBe` metres 6356750.52+ describe "Mean radius" $ do+ it "returns 6371008.771 m for the WGS84 ellipsoid" $ r84 `shouldBe` metres 6371008.771+ it "returns 6371008.771 m for the GRS80 ellipsoid" $ r80 `shouldBe` metres 6371008.771+ it "returns 6371006.84 m for the WG72 ellipsoid" $ r72 `shouldBe` metres 6371006.84
− test/Data/Geo/Jord/EvalSpec.hs
@@ -1,41 +0,0 @@-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/FramesSpec.hs view
@@ -0,0 +1,77 @@+module Data.Geo.Jord.FramesSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "Ellipsoidal earth model" $ do + describe "target" $ do + it "return the given point if NED norm = 0" $ do + let p0 = readLatLong "531914N0014347W" + let d = ned zero zero zero + targetN p0 d wgs84 `shouldBe` p0 + it "computes the target point from p0 and NED" $ do + let p0 = decimalLatLong 49.66618 3.45063 + let d = nedMetres (-86126) (-78900) 1069 + targetN p0 d wgs84 `shouldBe` decimalLatLong 48.8866688 2.374721111 + it "computes the target point from p0 and vector in Frame B" $ do + let p0 = decimalLatLongHeight 49.66618 3.45063 zero + let y = decimalDegrees 10 -- yaw + let r = decimalDegrees 20 -- roll + let p = decimalDegrees 30 -- pitch + let d = deltaMetres 3000 2000 100 + target p0 (frameB y r p) d wgs84 `shouldBe` + decimalLatLongHeight 49.6918016 3.4812669 (metres 6.007) + describe "nedBetween" $ do + it "computes NED between LatLong positions" $ do + let p1 = decimalLatLong 49.66618 3.45063 + let p2 = decimalLatLong 48.88667 2.37472 + let d = nedBetween p1 p2 wgs84 + d `shouldBe` nedMetres (-86125.88049540376) (-78900.08718759022) 1069.1981930266265 + it "computes NED between angular positions" $ do + let p1 = decimalLatLongHeight 49.66618 3.45063 zero + let p2 = decimalLatLongHeight 48.88667 2.37472 zero + let d = nedBetween p1 p2 wgs84 + d `shouldBe` nedMetres (-86125.88049540376) (-78900.08718759022) 1069.1981930266265 + describe "deltaBetween" $ + it "computes delta between angular positions in frame L" $ do + let p1 = decimalLatLongHeight 1 2 (metres (-3)) + let p2 = decimalLatLongHeight 4 5 (metres (-6)) + let w = decimalDegrees 5 -- wander azimuth + let d = deltaBetween p1 p2 (frameL w) wgs84 + d `shouldBe` deltaMetres 359490.579 302818.523 17404.272 + describe "deltaBetween and target consistency" $ + it "computes targetN p1 (nedBetween p1 p2) = p2" $ do + let p1 = decimalLatLongHeight 49.66618 3.45063 zero + let p2 = decimalLatLongHeight 48.88667 2.37472 zero + targetN p1 (nedBetween p1 p2 wgs84) wgs84 `shouldBe` p2 + describe "rotation matrix to go from/to earth-fixed frame to/from frame" $ do + it "computes the rotation matrix to go from Frame N to earth-fixed frame" $ do + let p = decimalLatLong 49.66618 3.45063 + let f = frameN p wgs84 + rEF f `shouldBe` + [ Vector3d (-0.7609044147490918) (-6.018845508229421e-2) (-0.6460664218872152) + , Vector3d (-4.588084170935741e-2) 0.9981870315100305 (-3.8956366478847045e-2) + , Vector3d 0.6472398473358879 0.0 (-0.7622864160016345) + ] + it "computes the rotation matrix to go from Frame B to earth-fixed frame" $ do + let p = decimalLatLong 49.66618 3.45063 + let f = frameB (decimalDegrees 10) (decimalDegrees 20) (decimalDegrees 30) p wgs84 + rEF f `shouldBe` + [ Vector3d (-0.4930071357732754) (-0.3703899170519431) (-0.7872453705312508) + , Vector3d 0.13374504887910177 0.8618333991702691 (-0.48923966925725315) + , Vector3d 0.8596837941807199 (-0.34648881860873165) (-0.3753521980782417) + ] + describe "North, East, Down delta" $ do + describe "norm" $ + it "computes the norm of a NED vector" $ + norm (nedMetres (-86126) (-78900) 1069) `shouldBe` metres 116807.708 + describe "bearing" $ + it "computes the bearing of a NED vector" $ + bearing (nedMetres (-86126) (-78900) 1069) `shouldBe` decimalDegrees 222.4927888 + describe "elevation" $ + it "computes the elevation of a NED vector from horizontal" $ + elevation (nedMetres (-86126) (-78900) 1069) `shouldBe` decimalDegrees (-0.5243663)
+ test/Data/Geo/Jord/GeodeticsSpec.hs view
@@ -0,0 +1,213 @@+module Data.Geo.Jord.GeodeticsSpec+ ( spec+ ) where++import Control.Exception.Base+import Data.Geo.Jord+import Data.Maybe (fromJust)+import Test.Hspec++spec :: Spec+spec = do+ describe "antipode" $ do+ it "returns the antipodal point" $ do+ let p = latLongHeight (readLatLong "484137N0061105E") (metres 15000)+ let e = decimalLatLongHeight (-48.6936111) (-173.8152777) (metres 15000)+ antipode p `shouldBe` e+ it "returns the south pole when called with the north pole" $+ antipode northPole `shouldBe` southPole+ it "returns the north pole when called with the south pole" $+ antipode southPole `shouldBe` northPole+ describe "crossTrackDistance" $ do+ it "returns a negative length when position is left of great circle (bearing)" $ do+ let p = decimalLatLong 53.2611 (-0.7972)+ let gc = greatCircleBearing (decimalLatLong 53.3206 (-1.7297)) (decimalDegrees 96.0)+ crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres (-305.663)+ it "returns a negative length when position is left of great circle" $ do+ let p = decimalLatLong 53.2611 (-0.7972)+ let gc = greatCircle (decimalLatLong 53.3206 (-1.7297)) (decimalLatLong 53.1887 0.1334)+ crossTrackDistance p gc (meanRadius wgs84) `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 (meanRadius wgs84) `shouldBe` metres 7042.324+ it "returns a positive length when position is left of great circle" $ do+ let p = antipode (decimalLatLong 53.2611 (-0.7972))+ let gc = greatCircle (decimalLatLong 53.3206 (-1.7297)) (decimalLatLong 53.1887 0.1334)+ crossTrackDistance p gc (meanRadius wgs84) `shouldBe` metres 307.547+ describe "destination" $ do+ it "return the given point if distance is 0 meter" $ do+ let p0 = readLatLong "531914N0014347W"+ destination p0 (decimalDegrees 96.0217) zero r84 `shouldBe` p0+ it "return the angular position along great-circle at distance and bearing" $ do+ let p0 = latLongHeight (readLatLong "531914N0014347W") (metres 15000.0)+ let p1 = decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0)+ destination p0 (decimalDegrees 96.0217) (metres 124800) r84 `shouldBe` p1+ it "return the ECEF position along great-circle at distance and bearing" $ do+ let p0 = ecefToNVector (ecefMetres 3812864.094 (-115142.863) 5121515.161) s84+ let p1 = ecefMetres 3826406.4710518294 8900.536398998282 5112694.233184049+ let p = destination84 p0 (decimalDegrees 96.0217) (metres 124800)+ nvectorToEcef p s84 `shouldBe` p1+ describe "finalBearing" $ do+ it "returns the Nothing if both point are the same" $ do+ let p = readLatLong "500359N0054253W"+ finalBearing p p `shouldBe` Nothing+ it "returns 0° if both point have the same longitude (going north)" $ do+ let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)+ let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000)+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 0)+ it "returns 180° if both point have the same longitude (going south)" $ do+ let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000)+ let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 180)+ it "returns 90° at the equator going east" $ do+ let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000)+ let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000)+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 90)+ it "returns 270° at the equator going west" $ do+ let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000)+ let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000)+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 270)+ it "returns the final bearing in compass angle" $ do+ let p1 = readLatLong "500359N0054253W"+ let p2 = readLatLong "583838N0030412W"+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 11.2752013)+ it "returns the final bearing in compass angle" $ do+ let p1 = readLatLong "583838N0030412W"+ let p2 = readLatLong "500359N0054253W"+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 189.1198181)+ it "returns the final bearing in compass angle" $ do+ let p1 = readLatLong "535941S0255915W"+ let p2 = readLatLong "54N154E"+ finalBearing p1 p2 `shouldBe` Just (decimalDegrees 125.6839436)+ describe "greatCircle smart constructors" $ do+ it "fails if both positions are equal" $+ greatCircleE (decimalLatLong 3 154) (decimalLatLong 3 154) `shouldBe`+ Left "Invalid Great Circle: positions are equal"+ it "fails if both positions are antipodal" $+ greatCircleE (decimalLatLong 3 154) (antipode (decimalLatLong 3 154)) `shouldBe`+ Left "Invalid Great Circle: positions are antipodal"+ describe "initialBearing" $ do+ it "returns Nothing if both point are the same" $ do+ let p = readLatLong "500359N1795959W"+ initialBearing p p `shouldBe` Nothing+ it "returns 0° if both point have the same longitude (going north)" $ do+ let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)+ let p2 = latLongHeight (readLatLong "583838N0054253W") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 0)+ it "returns 180° if both point have the same longitude (going south)" $ do+ let p1 = latLongHeight (readLatLong "583838N0054253W") (metres 12000)+ let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 180)+ it "returns 90° at the equator going east" $ do+ let p1 = latLongHeight (readLatLong "000000N0000000E") (metres 12000)+ let p2 = latLongHeight (readLatLong "000000N0010000E") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 90)+ it "returns 270° at the equator going west" $ do+ let p1 = latLongHeight (readLatLong "000000N0010000E") (metres 12000)+ let p2 = latLongHeight (readLatLong "000000N0000000E") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 270)+ it "returns the initial bearing in compass angle" $ do+ let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 12000)+ let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 9.1198181)+ it "returns the initial bearing in compass angle" $ do+ let p1 = latLongHeight (readLatLong "583838N0030412W") (metres 12000)+ let p2 = latLongHeight (readLatLong "500359N0054253W") (metres 5000)+ initialBearing p1 p2 `shouldBe` Just (decimalDegrees 191.2752013)+ describe "interpolate" $ do+ let p1 = readLatLong "44N044E"+ let p2 = readLatLong "46N046E"+ it "fails if f < 0.0" $+ evaluate (interpolate p1 p2 (-0.5)) `shouldThrow`+ errorCall "fraction must be in range [0..1], was -0.5"+ it "fails if f > 1.0" $+ evaluate (interpolate p1 p2 1.1) `shouldThrow`+ errorCall "fraction must be in range [0..1], was 1.1"+ it "returns p0 if f == 0" $ interpolate p1 p2 0.0 `shouldBe` p1+ it "returns p1 if f == 1" $ interpolate p1 p2 1.0 `shouldBe` p2+ it "returns the interpolated position" $ do+ let p3 = latLongHeight (readLatLong "53°28'46''N 2°14'43''W") (metres 10000)+ let p4 = latLongHeight (readLatLong "55°36'21''N 13°02'09''E") (metres 20000)+ interpolate p3 p4 0.5 `shouldBe`+ decimalLatLongHeight 54.7835574 5.1949856 (metres 15000)+ describe "insideSurface" $ do+ let p1 = decimalLatLong 45 1+ let p2 = decimalLatLong 45 2+ let p3 = decimalLatLong 46 1+ let p4 = decimalLatLong 46 2+ let p5 = decimalLatLong 45.1 1.1+ it "return False if polygon is empty" $ insideSurface p1 [] `shouldBe` False+ it "return False if polygon does not define at least a triangle" $+ insideSurface p1 [p1, p2] `shouldBe` False+ it "returns True if point is inside polygon" $ do+ let polygon = [p1, p2, p4, p3]+ insideSurface p5 polygon `shouldBe` True+ it "returns False if point is inside polygon" $ do+ let polygon = [p1, p2, p4, p3]+ let p = antipode p5+ insideSurface p polygon `shouldBe` False+ it "returns False if point is a vertex of the polygon" $ do+ let polygon = [p1, p2, p4, p3]+ insideSurface p1 polygon `shouldBe` False+ it "handles closed polygons" $ do+ let polygon = [p1, p2, p4, p3, p1]+ insideSurface p5 polygon `shouldBe` True+ it "handles concave polygons" $ do+ let malmo = decimalLatLong 55.6050 13.0038+ let ystad = decimalLatLong 55.4295 13.82+ let lund = decimalLatLong 55.7047 13.1910+ let helsingborg = decimalLatLong 56.0465 12.6945+ let kristianstad = decimalLatLong 56.0294 14.1567+ let polygon = [malmo, ystad, kristianstad, helsingborg, lund]+ let hoor = decimalLatLong 55.9295 13.5297+ let hassleholm = decimalLatLong 56.1589 13.7668+ insideSurface hoor polygon `shouldBe` True+ insideSurface hassleholm polygon `shouldBe` False+ describe "intersections" $ do+ it "returns nothing if both great circle are equals" $ do+ let gc = greatCircleBearing (decimalLatLong 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 (decimalLatLong 51.885 0.235) (decimalLatLong 52.885 1.235)+ let gc2 = greatCircle (decimalLatLong 52.885 1.235) (decimalLatLong 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 (decimalLatLong 51.885 0.235) (decimalDegrees 108.63)+ let gc2 = greatCircleBearing (decimalLatLong 49.008 2.549) (decimalDegrees 32.72)+ let (i1, i2) = fromJust (intersections gc1 gc2)+ i1 `shouldBe` decimalLatLong 50.9017226 4.4942782+ i2 `shouldBe` antipode i1+ describe "mean" $ do+ it "returns Nothing if no point is given" $ (mean [] :: Maybe NVector) `shouldBe` Nothing+ it "returns the unique given point" $ do+ let p = readLatLong "500359N0054253W"+ mean [p] `shouldBe` Just p+ it "returns the geographical mean" $ do+ let p1 = latLongHeight (readLatLong "500359N0054253W") (metres 15000.0)+ let p2 = latLongHeight (readLatLong "583838N0030412W") (metres 25000.0)+ let e = decimalLatLongHeight 54.3622869 (-4.5306725) zero+ mean [p1, p2] `shouldBe` Just e+ it "returns Nothing if list contains antipodal points" $ do+ let points =+ [ decimalLatLong 45 1+ , decimalLatLong 45 2+ , decimalLatLong 46 2+ , decimalLatLong 46 1+ , antipode (decimalLatLong 45 2)+ ]+ mean points `shouldBe` Nothing+ describe "surfaceDistance" $ do+ it "returns 0 if both points are equal" $ do+ let p = readLatLong "500359N1795959W"+ surfaceDistance p p r84 `shouldBe` zero+ it "returns the distance between 2 points" $ do+ let p1 = readLatLong "500359N0054253W"+ let p2 = readLatLong "583838N0030412W"+ surfaceDistance84 p1 p2 `shouldBe` metres 968854.868+ it "handles singularity at the pole" $+ surfaceDistance northPole southPole r84 `shouldBe` kilometres 20015.114351+ it "handles the discontinuity at the Date Line" $ do+ let p1 = readLatLong "500359N1795959W"+ let p2 = readLatLong "500359N1795959E"+ surfaceDistance p1 p2 (meanRadius wgs84) `shouldBe` metres 39.66
− test/Data/Geo/Jord/GreatCircleSpec.hs
@@ -1,48 +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
@@ -1,63 +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" +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` decimalLatLong 55.6058333 13.0005555+ it "reads 55°36'21''N 013°00'02''E" $+ readLatLong "55°36'21''N 013°00'02''E" `shouldBe` decimalLatLong 55.6058333 13.0005555+ it "reads 5536N01300E" $ readLatLong "5536N01300E" `shouldBe` decimalLatLong 55.6 13.0+ it "reads 55N013E" $ readLatLong "55N013E" `shouldBe` decimalLatLong 55.0 13.0+ it "reads 011659S0364900E" $+ readLatLong "011659S0364900E" `shouldBe` decimalLatLong (-1.2830555) 36.8166666+ it "reads 0116S03649E" $+ readLatLong "0116S03649E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666+ it "reads 1°16'S,36°49'E" $+ readLatLong "1°16'S,36°49'E" `shouldBe` decimalLatLong (-1.2666666) 36.8166666+ it "reads 01S036E" $ readLatLong "01S036E" `shouldBe` decimalLatLong (-1.0) 36.0+ it "reads 473622N1221955W" $+ readLatLong "473622N1221955W" `shouldBe` decimalLatLong 47.6061111 (-122.3319444)+ it "reads 4736N12219W" $+ readLatLong "4736N12219W" `shouldBe` decimalLatLong 47.6 (-122.3166666)+ it "reads 47N122W" $ readLatLong "47N122W" `shouldBe` decimalLatLong 47.0 (-122.0)+ it "reads 47°N 122°W" $ readLatLong "47°N 122°W" `shouldBe` decimalLatLong 47.0 (-122.0)+ it "reads 544807S0681811W" $+ readLatLong "544807S0681811W" `shouldBe` decimalLatLong (-54.8019444) (-68.3030555)+ it "reads 5448S06818W" $ readLatLong "5448S06818W" `shouldBe` decimalLatLong (-54.8) (-68.3)+ it "reads 54S068W" $ readLatLong "54S068W" `shouldBe` decimalLatLong (-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 (decimalLatLong 55.60583333 13.00055556) `shouldBe` "55°36'21.000\"N,13°0'2.000\"E"+ it "shows the S/E position formatted in DMS with symbols" $+ show (decimalLatLong (-1.28305556) 36.81666) `shouldBe` "1°16'59.000\"S,36°48'59.976\"E"+ it "shows the N/W position formatted in DMS with symbols" $+ show (decimalLatLong 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 (decimalLatLong (-54.80194) (-68.30305)) `shouldBe`+ "54°48'6.984\"S,68°18'10.980\"W"
test/Data/Geo/Jord/LengthSpec.hs view
@@ -11,6 +11,7 @@ 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+ it "reads 25000ft" $ readLength "25000ft" `shouldBe` feet 25000 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"@@ -27,6 +28,8 @@ 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+ it "converts feet to metres" $ toMetres (feet 25000) `shouldBe` 7620+ it "converts metres to feet" $ toFeet (metres 7620) `shouldBe` 25000 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/Data/Geo/Jord/PositionSpec.hs
@@ -1,153 +0,0 @@-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
+ test/Data/Geo/Jord/RotationSpec.hs view
@@ -0,0 +1,49 @@+module Data.Geo.Jord.RotationSpec + ( spec + ) where + +import Data.Geo.Jord +import Test.Hspec + +spec :: Spec +spec = do + describe "r2xyz" $ + it "computes the 3 angles about new axes in the xyz-order from rotation matrix" $ do + let xyz = [decimalDegrees 45, decimalDegrees 45, decimalDegrees 5] + let rm = + [ Vector3d 0.7044160264027587 (-6.162841671621935e-2) 0.7071067811865475 + , Vector3d 0.559725765762092 0.6608381550289296 (-0.5) + , Vector3d (-0.43646893232965345) 0.7479938977765876 0.5000000000000001 + ] + r2xyz rm `shouldBe` xyz + describe "r2xyz" $ + it "computes the 3 angles about new axes in the zyx-order from rotation matrix" $ do + let zyx = [decimalDegrees 10, decimalDegrees 20, decimalDegrees 30] + let rm = + [ Vector3d 0.9254165783983234 1.802831123629725e-2 0.37852230636979245 + , Vector3d 0.16317591116653482 0.8825641192593856 (-0.44096961052988237) + , Vector3d (-0.3420201433256687) 0.46984631039295416 0.8137976813493738 + ] + r2zyx rm `shouldBe` zyx + describe "xyz2r" $ + it "computes the rotation matrix from 3 angles about new axes in the xyz-order" $ do + let x = decimalDegrees 45 + let y = decimalDegrees 45 + let z = decimalDegrees 5 + let rm = + [ Vector3d 0.7044160264027587 (-6.162841671621935e-2) 0.7071067811865475 + , Vector3d 0.559725765762092 0.6608381550289296 (-0.5) + , Vector3d (-0.43646893232965345) 0.7479938977765876 0.5000000000000001 + ] + xyz2r x y z `shouldBe` rm + describe "zyx2r" $ + it "computes the rotation matrix from 3 angles about new axes in the zyx-order" $ do + let x = decimalDegrees 10 + let y = decimalDegrees 20 + let z = decimalDegrees 30 + let rm = + [ Vector3d 0.9254165783983234 1.802831123629725e-2 0.37852230636979245 + , Vector3d 0.16317591116653482 0.8825641192593856 (-0.44096961052988237) + , Vector3d (-0.3420201433256687) 0.46984631039295416 0.8137976813493738 + ] + zyx2r x y z `shouldBe` rm
+ test/Data/Geo/Jord/TransformSpec.hs view
@@ -0,0 +1,69 @@+module Data.Geo.Jord.TransformSpec+ ( spec+ ) where++import Data.Geo.Jord+import Test.Hspec++spec :: Spec+spec = do+ describe "Ellipsoidal transformation between coordinates systems" $ do+ it "transforms NVector position to ECEF position" $ do+ let p = nvector 0.5 0.5 0.7071+ toEcef p wgs84 `shouldBe` ecefMetres 3194434.411 3194434.411 4487326.819+ it "transforms angular position to ECEF position" $ do+ let refAngular =+ [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834)+ , decimalLatLongHeight 45.0 45.0 zero+ , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)+ ]+ let refEcefs =+ [ ecefMetres 5733855.775 (-6370998.38) 7008137.511+ , ecefMetres 3194419.145 3194419.145 4487348.409+ , ecefMetres 4200996.77 172460.321 4780102.808+ ]+ mapM_ (\(a, e) -> toEcef a wgs84 `shouldBe` e) (zip refAngular refEcefs)+ it "transforms ECEF position to angular position" $ do+ let refAngular =+ [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834050887)+ , decimalLatLongHeight 45.0 45.0 zero+ , decimalLatLongHeight 48.8562 2.3508 (metres 67.36990469945641)+ ]+ let refEcefs =+ [ ecefMetres 5733855.774881717 (-6370998.380260889) 7008137.510624695+ , ecefMetres 3194419.145121972 3194419.145121971 4487348.408606014+ , ecefMetres 4200996.769831858 172460.320727757 4780102.807914356+ ]+ mapM_ (\(a, e) -> fromEcef e wgs84 `shouldBe` a) (zip refAngular refEcefs)+ describe "Spherical transformation between coordinates systems" $ do+ it "transforms NVector position to ECEF position" $ do+ let p = nvector 0.5 0.5 0.7071+ toEcef p s84 `shouldBe` ecefMetres 3185519.660103391 3185519.660103391 4504961.903318216+ it "transforms angular position to ECEF position" $ do+ let refAngular =+ [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834)+ , decimalLatLongHeight 45.0 45.0 zero+ , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)+ , latLongHeight (readLatLong "531914N0014347W") (metres 15000.0)+ , decimalLatLongHeight 53.1882691 0.1332744 (metres 15000.0)+ ]+ let refEcefs =+ [ ecefMetres 5725717.354041086 (-6361955.622990872) 7025277.913631903+ , ecefMetres 3185504.385500001 3185504.3855 4504983.504973072+ , ecefMetres 4188328.8912726147 171940.27595767862 4797806.669141033+ , ecefMetres 3812864.094233316 (-115142.863124558) 5121515.160893968+ , ecefMetres 3826406.4642097903 8900.535428827865 5112694.238306408+ ]+ mapM_ (\(a, e) -> toEcef a s84 `shouldBe` e) (zip refAngular refEcefs)+ it "transforms ECEF position to angular position" $ do+ let refAngular =+ [ decimalLatLongHeight 39.379 (-48.013) (metres 4702059.834217537)+ , decimalLatLongHeight 45.0 45.0 (metres 1e-3)+ , decimalLatLongHeight 48.8562 2.3508 (metres 67.36972232195099)+ ]+ let refEcefs =+ [ ecefMetres 5725717.354 (-6361955.623) 7025277.914+ , ecefMetres 3185504.386 3185504.386 4504983.505+ , ecefMetres 4188328.891 171940.276 4797806.669+ ]+ mapM_ (\(a, e) -> fromEcef e s84 `shouldBe` a) (zip refAngular refEcefs)