diff --git a/Data/GPS/Core.hs b/Data/GPS/Core.hs --- a/Data/GPS/Core.hs +++ b/Data/GPS/Core.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} module Data.GPS.Core ( -- * Types Distance @@ -11,6 +12,7 @@ , east , west , radiusOfEarth + , circumferenceOfEarth -- * Coordinate Functions , heading , distance @@ -33,7 +35,7 @@ import Data.Time import Data.Maybe import Control.Monad -import Text.XML.HXT.Arrow +import Text.XML.HXT.Core import Text.XML.XSD.DateTime(DateTime,toUTCTime) import Data.Geo.GPX @@ -47,6 +49,7 @@ -- (3/2)pi == East == - (pi / 2) type Heading = Double +type Angle = Double -- |Speed is hard coded as meters per second type Speed = Double type Vector = (Distance, Heading) @@ -56,25 +59,25 @@ getUTCTime :: (Time a) => a -> Maybe UTCTime getUTCTime = fmap toUTCTime . time -acos' x = if x > 1 then acos 1 else if x < (-1) then acos (-1) else acos x - distance :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Distance -distance a b = - let x = sin lat1 * sin lat2 + cos lat1 * cos lat2 * cos (lon2 - lon1) - in radiusOfEarth * acos' x - where - (lat1, lon1) = getRadianPairD a - (lat2, lon2) = getRadianPairD b +distance x y = + let (lat1,lon1) = getRadianPairD x + (lat2,lon2) = getRadianPairD y + deltaLat = lat2 - lat1 + deltaLon = lon2 - lon1 + a = (sin (deltaLat / 2))^2 + cos lat1 * cos lat2 * (sin (deltaLon / 2))^2 + c = 2 * atan2 (a**0.5) ((1-a)**0.5) + in radiusOfEarth * c -- | Direction two points aim toward (0 = North, pi/2 = West, pi = South, 3pi/2 = East) heading :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Heading heading a b = - atan2 (sin (diffLon) * cos (lat2)) + atan2 (sin (diffLon) * cos (lat2)) (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon)) where (lat1, lon1) = getRadianPairD a (lat2, lon2) = getRadianPairD b - diffLon = lon1 - lon2 + diffLon = lon2 - lon1 getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector getVector a b = (distance a b, heading a b) @@ -94,8 +97,10 @@ . setLat (latitudeType $ toDegrees lat2) $ p where (lat,lon) = getRadianPairD p - lat2 = lat + (cos h) * (d / radiusOfEarth) - lon2 = lon - acos' ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2)) + lat2 = asin (sin (lat) * cos (d / radiusOfEarth) + cos(lat) + * sin(d/radiusOfEarth) * cos h) + lon2 = lon + atan2 (sin h * sin (d / radiusOfEarth) * cos lat) + (cos (d/radiusOfEarth) - sin lat * sin lat2) -- | Speed in meters per second, only if a 'Time' was recorded for each waypoint. speed :: (Lat loc, Lon loc, Time loc, Lat b, Lon b, Time b) => loc -> b -> Maybe Speed @@ -109,6 +114,10 @@ -- |radius of the earth in meters radiusOfEarth :: Double radiusOfEarth = 6378700 + +-- |Circumference of earht (meters) +circumferenceOfEarth :: Double +circumferenceOfEarth = radiusOfEarth * 2 * pi -- |North is 0 radians north :: Heading diff --git a/Data/GPS/Trail.hs b/Data/GPS/Trail.hs --- a/Data/GPS/Trail.hs +++ b/Data/GPS/Trail.hs @@ -17,8 +17,8 @@ , closestDistance , convexHull -- ** Transformations - -- , bezierCurveAt - -- , bezierCurve + , bezierCurveAt + , bezierCurve , linearTime , filterPoints -- ** Grouping Methods @@ -35,14 +35,17 @@ , refineGrouping , (/\), (\/) -- ** Composite Operations (Higher Level) - -- , smoothRests - -- , smoothSegments - -- , smoothPath + , smoothRests + , smoothSegments +-- , smoothPath + -- * Misc + , bezierPoint ) where import Text.Show.Functions () import Data.GPS.Core hiding (fix) +import Text.XML.XSD.DateTime (fromUTCTime) import Control.Arrow (first, second) import Control.Monad import Data.Fixed (mod') @@ -94,6 +97,7 @@ -- getAvg :: [] -> AvgMethod -> Speed getAvg _ [] = 0 + getAvg _ [x] = 0 getAvg m cs = let ss = getSpeedsV cs in case m of @@ -336,11 +340,16 @@ top = head timesDef totalTime = diffUTCTime (snd end) (snd top) times = if null selectedTimes then map snd timesDef else selectedTimes - queryTimes = [fromTo (diffUTCTime (snd end) t / totalTime) | t <- times] - fromTo = fromRational . toRational - in if null timesDef || totalTime == 0 || any (\x -> x < 0 || x > 1) queryTimes - then xs - else map (bezierPoint xs) queryTimes + diffTimes = [diffUTCTime t (snd top) / totalTime | t <- times] + queryTimes = map realToFrac diffTimes + in if any (\x -> x < 0 || x > 1) queryTimes + then error "bezierCurveAt has a out-of-bound time!" + else + if null timesDef || totalTime == 0 || any (\x -> x < 0 || x > 1) queryTimes + then xs + else let curvePoints = (map (bezierPoint xs) queryTimes) + newTimes = [addUTCTime t (snd top) | t <- diffTimes] + in zipWith (setTime . Just . fromUTCTime) newTimes curvePoints bezierPoint :: (Lat a, Lon a) => [a] -> Double -> a bezierPoint pnts t = go pnts @@ -429,12 +438,15 @@ smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60 smoothSegments :: (Lat a, Lon a, Time a) => Trail a -> Trail a -smoothSegments ps = - let ps' = bezierCurve . everyNPoints 5 $ ps - (h,t) = splitAt 2 ps' - ps'' = bezierCurve . everyNPoints 5 $ t - in h ++ ps'' - +smoothSegments ps = bezierCurve . everyNPoints 7 $ ps +{- + let op xs = + let xs' = bezierCurve . everyNPoints 7 $ xs + (h,t) = splitAt 3 xs' + xs'' = bezierCurve . everyNPoints 7 $ xs' + in h ++ xs'' + in iterate op ps !! 10 +-} smoothPath :: (Lat a, Lon a, Time a) => Trail a -> Trail a smoothPath ps = undefined diff --git a/Test/GpsTest.hs b/Test/GpsTest.hs new file mode 100644 --- /dev/null +++ b/Test/GpsTest.hs @@ -0,0 +1,118 @@ +import Data.GPS +import Data.Time +import Data.List +import Data.Ord +import Data.Fixed +import Test.QuickCheck +import Test.Framework (Test, defaultMain, testGroup) +import Test.Framework.Providers.QuickCheck2 (testProperty) +import Text.XML.XSD.DateTime +import Control.Applicative +import Control.Monad + +instance Arbitrary LatitudeType where + arbitrary = liftM (latitudeType . flip mod' 180) arbitrary + +instance Arbitrary LongitudeType where + arbitrary = liftM (longitudeType . flip mod' 180) arbitrary + +instance Arbitrary DateTime where + arbitrary = liftM fromUTCTime arbitrary + +instance Arbitrary UTCTime where + arbitrary = UTCTime <$> arbitrary <*> liftM (secondsToDiffTime . abs) (arbitrary :: Gen Integer) + +instance Arbitrary Day where + arbitrary = fromGregorian <$> liftM ((+ 1000) . (`mod` 1100) .abs) arbitrary <*> liftM abs arbitrary <*> liftM abs arbitrary + +instance Arbitrary NominalDiffTime where + arbitrary = liftM fromIntegral (arbitrary :: Gen Int) +instance Arbitrary WptType where + arbitrary = + wptType <$> arbitrary -- Lat + <*> arbitrary -- Lon + <*> arbitrary -- Time + <*> arbitrary -- elevation + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return [] -- LinkType + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + <*> return Nothing + +newtype Trl = Trl [WptType] + deriving (Show) + +instance Arbitrary Trl where + arbitrary = do + b <- arbitrary :: Gen [Int] + u_ts_d <- mapM (\i -> (,,) <$> arbitrary <*> replicateM i arbitrary <*> arbitrary) b :: Gen [(UTCTime, [WptType],NominalDiffTime)] + let u_ts_d' = sortBy (comparing (\(a,_,_) -> a)) u_ts_d + xs = concat [zipWith (setTime' . fromUTCTime) (iterate (addUTCTime d) u) x | (u,x,d) <- u_ts_d'] + return $ Trl xs + +approxEq :: WptType -> WptType -> Bool +approxEq a b = distance a b <= 0.2 -- error of 13cm has been observed due to floating point issues when using add vector. + +pSaneDistance :: WptType -> WptType -> Bool +pSaneDistance a b = distance a b <= circumferenceOfEarth / 2 + +pTriangleTheorem :: WptType -> WptType -> WptType -> Bool +pTriangleTheorem a b c = + distance a b + distance b c >= distance a c -- Traditional flat-surface geometry + || distance a b + distance b c + distance c a == 2 * pi * radiusOfEarth + +pAddVector_DistanceHeading_ident :: WptType -> WptType -> Bool +pAddVector_DistanceHeading_ident a b = + let v = (distance a b, heading a b) + c = addVector v a + in (distance c b) <= 0.01 * (distance a b) + +pConvexHull_Has_Extreme_Points :: Trl -> Bool +pConvexHull_Has_Extreme_Points (Trl ts) = + let ch = convexHull ts + ts' = sortBy (comparing lat) ts + northMost = last ts' + southMost = head ts' + in length ts < 3 || (northMost `elem` ch && southMost `elem` ch) + +pConvexHull_Bezier_Const :: Trl -> Double -> Bool +pConvexHull_Bezier_Const (Trl ts) n = + let ts' = take 10 ts + ch = convexHull ts' + n' = abs (n `mod'` 1) + bp = bezierPoint ts' n' + ch' = convexHull (bp:ts') + in length ts < 3 || ch == ch' + +tests :: [Test] +tests = + [ + testGroup "Coordinate Computations" + [ testProperty "approxEq_id" (\x -> approxEq x x) + , testProperty "saneDistance" pSaneDistance + , testProperty "TriangleTheorem" pTriangleTheorem + , testProperty "Vector identity" pAddVector_DistanceHeading_ident + ] + , testGroup "Trail Computations" + [ testProperty "Hull has extreme points" pConvexHull_Has_Extreme_Points + , testProperty "HullContainsBezier" pConvexHull_Bezier_Const] + ] + + +---------------------------------------------------------- +------------- * HARNESS * -------------------------------- +---------------------------------------------------------- +main :: IO () +main = defaultMain tests diff --git a/gps.cabal b/gps.cabal --- a/gps.cabal +++ b/gps.cabal @@ -1,5 +1,5 @@ name: gps -version: 0.8.1 +version: 0.8.2 license: BSD3 license-file: LICENSE author: Thomas DuBuisson @@ -9,20 +9,33 @@ category: Data stability: stable build-type: Simple -cabal-version: >= 1.6 +cabal-version: >= 1.8 tested-with: GHC == 6.10.3 extra-source-files: Library Build-Depends: base >= 3 && < 6, pretty >= 1.0 , prettyclass >= 1.0, - time >= 1.1, GPX == 0.4.*, hxt >= 8.5, xsd == 0.3.*, + time >= 1.1, GPX >= 0.5, hxt >= 9.1, xsd >= 0.3, vector >= 0.7, statistics >= 0.9 hs-source-dirs: exposed-modules: Data.GPS other-modules: Data.GPS.Core Data.GPS.Trail ghc-options: +test-suite gps-tests + hs-source-dirs: Test + main-is: GpsTest.hs + type: exitcode-stdio-1.0 + + build-depends: + base >= 4.0, + QuickCheck >= 2.4.0.1, + test-framework >= 0.3.3 && < 0.5, + test-framework-quickcheck2 >= 0.2.9 && < 0.3, + time, GPX, hxt, xsd, vector, statistics, gps + + ghc-options: -Wall source-repository head type: git location: git@github.com:TomMD/gps.git