gps 0.7 → 0.8
raw patch · 5 files changed
+638/−269 lines, 5 filesdep +statisticsdep +vectordep ~timesetup-changed
Dependencies added: statistics, vector
Dependency ranges changed: time
Files
- Data/GPS.hs +5/−266
- Data/GPS/Core.hs +183/−0
- Data/GPS/Trail.hs +441/−0
- Setup.hs +0/−0
- gps.cabal +9/−3
Data/GPS.hs view
@@ -1,273 +1,12 @@-{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls, BangPatterns #-}+{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls, BangPatterns, TupleSections #-} -- |A basic GPS library with calculations for distance and speed along -- with helper functions for filtering/smoothing trails. All distances -- are in meters and time is in seconds. Speed is thus meters/second module Data.GPS- ( -- * Types- Distance- , Heading- , Speed- , Vector- , Trail- -- * Constants- , north- , south- , east- , west- , radiusOfEarth- -- * Coordinate Functions- , heading- , distance- , speed- , addVector- , getRadianPair- , getDMSPair- , divideArea- -- * Trail Functions- , totalDistance- , restLocations- , closestDistance- , filterByMaxSpeed- , convexHull- -- * Other helpers- , readGPX- , module Data.Geo.GPX+ ( module Data.GPS.Core+ , module Data.GPS.Trail ) where -import Data.Function (on)-import Data.Ord (comparing)-import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy)-import Data.Geo.GPX hiding (none, cmt)--import Text.XML.HXT.Arrow-import Text.XML.XSD.DateTime(DateTime, toUTCTime)--import Data.Time-import Data.Maybe (listToMaybe)-import Data.Fixed-import Control.Monad---- |Distances are expressed in meters-type Distance = Double---- |Angles are expressed in radians from North.--- 0 == North--- pi/2 == West--- pi == South--- (3/2)pi == East == - (pi / 2)-type Heading = Double---- |Speed is hard coded as meters per second-type Speed = Double--type Vector = (Distance, Heading)-type Trail a = [a]--getUTCTime :: (Lat a, Lon a, Time a) => a -> Maybe UTCTime-getUTCTime = fmap toUTCTime . time--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)- x' = if x > 1 then 1 else x- in radiusOfEarth * acos x'- where- (lat1, lon1) = getRadianPairD a- (lat2, lon2) = getRadianPairD b---- | Find the total distance traveled-totalDistance :: (Lat a, Lon a) => [a] -> Distance-totalDistance as = sum $ zipWith distance as (drop 1 as)---- | 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 -- ^ 0 = North, pi/2 = West...-heading a b =- 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--getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector-getVector a b = (distance a b, heading a b)---- | 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-speed a b = - case (getUTCTime b, getUTCTime a) of- (Just x, Just y) -> - let timeDiff = realToFrac (diffUTCTime x y)- in if timeDiff == 0 then Nothing else Just $ (distance a b) / timeDiff- _ -> Nothing---- |Filter out all points that result in a speed greater than a given--- value (the second point is dropped)-filterByMaxSpeed :: (Lat loc, Lon loc, Time loc) => Speed -> Trail loc -> Trail loc-filterByMaxSpeed mx xs =- let ss = zipWith speed xs (drop 1 xs)- fs = filter ((< Just mx) . fst) (zip ss $ drop 1 xs)- in take 1 xs ++ map snd fs--data TempTrail a = T (Trail a) a---- |radius of the earth in meters-radiusOfEarth :: Double-radiusOfEarth = 6378700---- |North is 0 radians-north :: Heading-north = 0---- |South, being 180 degrees from North, is pi.-south :: Heading-south = pi---- |East is 270 degrees from North-east :: Heading-east = (3 / 2) * pi---- |West is 90 degrees (pi/2)-west :: Heading-west = pi / 2--toDecimal = (*) (180 / pi)----- |Given a vector and coordinate, computes a new coordinate.--- Within some epsilon it should hold that if------ @dest = addVector (dist,heading) start@------ then------ @heading == dmsHeading start dest@--- --- @dist == distance start dest@-addVector :: (Lat c, Lon c) => Vector -> c -> c-addVector (d,h) p = setLon (longitudeType lon2) . setLat (latitudeType 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))--getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double)-getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair--getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)-getDMSPair c = (lat c, lon c)---- |Provides a lat/lon pair of doubles in radians-getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)-getRadianPair p = (toRadians (lat p), toRadians (lon p))--toRadians :: Floating f => f -> f-toRadians = (*) (pi / 180)---- |Filters out any points that go backward in time (thus must not be valid if this is a trail)-linearTime :: (Lat a, Lon a, Time a) => Trail a -> Trail a-linearTime [] = []-linearTime (p:ps) = go (getUTCTime p) ps- where- go _ [] = []- go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps---- |Creates a list of trails all of which are within the given distance of each--- other spanning atleast the given amount of time.------ For example @restLocations 50 600@--- would return lists of all points that are within 50 meters of each other and--- span at least 10 minutes (600 seconds).------ Note this gives points within fifty meters of the earliest point - wandering--- in a rest area with a 50 meter radius could result in several rest points--- ([a,b..]) or even none if the distance between individual points exceeds 50m.-restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a]-restLocations d s xs = go xs- where- go [] = []- go (a:as) =- case takeWhileLast ((<=) d . distance a) as of- (Just l, close, far) ->- case (getUTCTime a, getUTCTime l) of- (Just t1, Just t2) ->- let diff = diffUTCTime t2 t1- in if diff >= s then (a:close) : go far else go as- _ -> go as- _ -> go as--takeWhileLast :: (a -> Bool) -> [a] -> (Maybe a, [a], [a])-takeWhileLast p [] = (Nothing, [], [])-takeWhileLast p (x:xs)- | not (p x) = (Nothing, [], x:xs)- | otherwise = go x xs- where- go !a [] = (Just a, [], [])- go !a (b:bs)- | p b = let (c,d,f) = go b bs in (c, a:d, f)- | otherwise = (Just a, [], b:bs)---- |Returns the closest distance between two trails (or Nothing if a trail is empty)--- O( (n * m) * log (n * m) )-closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance-closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs]----- |@divideArea vDist hDist nw se@ divides an area into a grid of equally--- spaced coordinates within the box drawn by the northwest point (nw) and--- southeast point (se). Because this uses floating point there might be a--- different number of points in some rows (the last might be too far east based--- on a heading from the se point).-divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]-divideArea vDist hDist nw se =- let (top,left) = (lat nw, lon nw)- (btm,right) = (lat se, lon se)- columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw- buildRow = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))- in map buildRow columnOne---- |Reads a GPX file (using the GPX library) by simply concatenating all the--- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.-readGPX :: FilePath -> IO (Trail WptType)-readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile---- | Uses Grahams scan to compute the convex hull of the given points.--- This operation requires sorting of the points, so don't try it unless--- you have notably more memory than the list of points will consume.-convexHull :: (Eq c, Lat c, Lon c) => [c] -> [c]-convexHull xs =- let first = southMost xs- in case first of- Nothing -> []- Just f ->- let sorted = sortBy (comparing (eastZeroHeading f)) (filter (/= f) xs)- in case sorted of- (a:b:cs) -> grahamScan (b:a:f:[]) cs- cs -> f : cs- where- grahamScan [] _ = []- grahamScan ps [] = ps- grahamScan (x:[]) _ = [x]- grahamScan (p2:p1:ps) (x:xs) =- case turn p1 p2 x of- LeftTurn -> grahamScan (x:p2:p1:ps) xs- Straight -> grahamScan (x:p2:p1:ps) xs- _ -> grahamScan (p1:ps) (x:xs)--eastZeroHeading :: (Lat c, Lon c) => c -> c -> Heading-eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s--data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)--turn :: (Lat c, Lon c) => c -> c -> c -> Turn-turn a b c =- let h1 = eastZeroHeading a b- h2 = eastZeroHeading b c- d = h2 - h1- in if d >= 0 && d < pi then LeftTurn else RightTurn---- | Find the southmost point-southMost :: (Lat c) => [c] -> Maybe c-southMost [] = Nothing-southMost cs = Just . minimumBy (comparing lat) $ cs+import Data.GPS.Core+import Data.GPS.Trail
+ Data/GPS/Core.hs view
@@ -0,0 +1,183 @@+module Data.GPS.Core+ ( -- * Types+ Distance+ , Heading+ , Speed+ , Vector+ , Trail+ -- * Constants+ , north+ , south+ , east+ , west+ , radiusOfEarth+ -- * Coordinate Functions+ , heading+ , distance+ , speed+ , getVector+ , addVector+ , getRadianPair+ , getDMSPair+ , divideArea+ , interpolate+ -- * IO helpers+ , writeGPX+ , readGPX+ , readGPXSegments+ -- * Utility+ , getUTCTime+ , module Data.Geo.GPX+ ) where++import Data.Time+import Data.Maybe+import Control.Monad+import Text.XML.HXT.Arrow+import Text.XML.XSD.DateTime(DateTime,toUTCTime)+import Data.Geo.GPX++-- |Distances are expressed in meters+type Distance = Double++-- |Angles are expressed in radians from North.+-- 0 == North+-- pi/2 == West+-- pi == South+-- (3/2)pi == East == - (pi / 2)+type Heading = Double++-- |Speed is hard coded as meters per second+type Speed = Double+type Vector = (Distance, Heading)++type Trail a = [a]++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++-- | 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)) + (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))+ where+ (lat1, lon1) = getRadianPairD a+ (lat2, lon2) = getRadianPairD b+ diffLon = lon1 - lon2++getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector+getVector a b = (distance a b, heading a b)++-- |Given a vector and coordinate, computes a new coordinate.+-- Within some epsilon it should hold that if+--+-- @dest = addVector (dist,heading) start@+--+-- then+--+-- @heading == heading start dest@+-- +-- @dist == distance start dest@+addVector :: (Lat c, Lon c) => Vector -> c -> c+addVector (d,h) p = setLon (longitudeType $ toDegrees lon2) + . 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))++-- | 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+speed a b = + case (getUTCTime b, getUTCTime a) of+ (Just x, Just y) -> + let timeDiff = realToFrac (diffUTCTime x y)+ in if timeDiff == 0 then Nothing else Just $ (distance a b) / timeDiff+ _ -> Nothing++-- |radius of the earth in meters+radiusOfEarth :: Double+radiusOfEarth = 6378700++-- |North is 0 radians+north :: Heading+north = 0++-- |South, being 180 degrees from North, is pi.+south :: Heading+south = pi++-- |East is 270 degrees (3 pi / 2)+east :: Heading+east = (3 / 2) * pi++-- |West is 90 degrees (pi/2)+west :: Heading+west = pi / 2++toDegrees = (*) (180 / pi)++getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double)+getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair++getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)+getDMSPair c = (lat c, lon c)++-- |Provides a lat/lon pair of doubles in radians+getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)+getRadianPair p = (toRadians (lat p), toRadians (lon p))++toRadians :: Floating f => f -> f+toRadians = (*) (pi / 180)++-- | @interpolate c1 c2 w@ where 0 <= w <= 1 Gives a point on the line+-- between c1 and c2 equal to @c1 when @w == 0@ (weighted linearly+-- toward c2).+interpolate :: (Lat a, Lon a) => a -> a -> Double -> a+interpolate c1 c2 w+ | w < 0 || w > 1 = error "Interpolate only works with a weight between zero and one"+ | otherwise = + let (h,d) = (heading c1 c2, distance c1 c2)+ v = (d * w, h)+ in addVector v c1+ +-- |@divideArea vDist hDist nw se@ divides an area into a grid of equally+-- spaced coordinates within the box drawn by the northwest point (nw) and+-- southeast point (se). Because this uses floating point there might be a+-- different number of points in some rows (the last might be too far east based+-- on a heading from the se point).+divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]+divideArea vDist hDist nw se =+ let (top,left) = (lat nw, lon nw)+ (btm,right) = (lat se, lon se)+ columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw+ buildRow = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))+ in map buildRow columnOne++-- |Reads a GPX file (using the GPX library) by simply concatenating all the+-- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.+readGPX :: FilePath -> IO (Trail WptType)+readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile++writeGPX :: FilePath -> Trail WptType -> IO ()+writeGPX fp ps = writeGpxFile fp $ gpx $ gpxType "1.0" "Haskell GPS Package (via the GPX package)" Nothing [] [] [trkType Nothing Nothing Nothing Nothing [] Nothing Nothing Nothing [trksegType ps Nothing]] Nothing++-- writeGpxFile should go in the GPX package+writeGpxFile :: FilePath -> Gpx -> IO ()+writeGpxFile fp gpx = runX_ (constA gpx >>> xpickleDocument (xpickle :: PU Gpx) [] fp)++runX_ t = runX t >> return ()++readGPXSegments :: FilePath -> IO [Trail WptType]+readGPXSegments = liftM (map (concatMap trkpts) . map trksegs . concatMap trks) . readGpxFile
+ Data/GPS/Trail.hs view
@@ -0,0 +1,441 @@+{-# LANGUAGE TupleSections #-}+module Data.GPS.Trail+ ( -- * Types + AvgMethod(..)+ , Selected(..)+ -- * Utility Functions+ , isSelected+ , isNotSelected+ , onSelected+ , selLength+ -- * Trail Functions+ -- ** Queries+ , totalDistance+ , avgSpeeds+ , slidingAverageSpeed+ , closestDistance+ , convexHull+ -- ** Transformations+-- , bezierCurveAt+-- , bezierCurve+ , linearTime+ , filterPoints+ -- ** Grouping Methods+ , betweenSpeeds+ , restLocations+ , spansTime+ , everyNPoints+ -- ** Group Transformations + , intersectionOf+ , invertSelection+ , firstGrouping+ , lastGrouping+ , unionOf+ , refineGrouping+ , (/\), (\/)+ -- -- ** Composite Operations (Higher Level)+-- , smoothRests+-- , smoothSegments+-- , smoothPath+ ) where++import Text.Show.Functions ()+import Data.GPS.Core hiding (fix)++import Control.Arrow (first, second)+import Control.Monad+import Data.Fixed (mod')+import Data.Function (on,fix)+import Data.List as L+import Data.Maybe+import Data.Ord+import Data.Time++import Statistics.Function as F+import Statistics.Sample+import qualified Data.Vector.Unboxed as V++takeWhileEnd :: (a -> Bool) -> [a] -> (Maybe a, [a],[a])+takeWhileEnd p xs = go xs Nothing+ where+ go [] e = (e, [], [])+ go (a:as) e+ | p a = let (e',xs,zs) = go as (Just a) in (e',a:xs,zs)+ | otherwise = (e,[], a:as)++data AvgMethod c+ = AvgMean -- ^ Obtain the 'mean' of the considered points+ | AvgHarmonicMean -- ^ Obtain the 'harmonicMean'+ | AvgGeometricMean -- ^ Obtain the 'geometricMean'+ | AvgMedian -- ^ Obtain the median of the considered points+ | AvgEndPoints -- ^ Compute the speed considering only the given endpoints+ | AvgMinOf [AvgMethod c] -- ^ Take the minimum of the speeds from the given methods+ | AvgWith ([c] -> Speed)+ +-- | @avgSpeeds n points@+-- Average speed using a window of up to @n@ seconds and averaging by taking the+-- Median ('AvgMedian').+avgSpeeds :: (Lat a, Lon a, Time a) => NominalDiffTime -> Trail a -> [(UTCTime, Speed)]+avgSpeeds = slidingAverageSpeed AvgHarmonicMean++-- | @slidingAverageSpeed m n@ Average speed using a moving window of up to @n@ seconds+-- and an 'AvgMethod' of @m@.+slidingAverageSpeed :: (Lat a, Lon a, Time a) => + AvgMethod a -> NominalDiffTime -> Trail a -> [(UTCTime, Speed)]+slidingAverageSpeed _ _ [] = []+slidingAverageSpeed m n (x:xs) =+ let avg = getAvg (x:xs') m+ avgTime = getAvgTime x (fromMaybe x e)+ in case avgTime of+ Nothing -> []+ Just t -> (t,avg) : slidingAverageSpeed m n xs+ where+ (e,xs',rest) = takeWhileEnd (\c -> getTimeDiff c x <= Just n) xs+ getTimeDiff a b = on (liftM2 diffUTCTime) getUTCTime a b+ + -- getAvg :: [] -> AvgMethod -> Speed+ getAvg cs AvgMean =+ let ss = getSpeedsV cs+ in case m of+ AvgMean -> mean ss+ AvgHarmonicMean -> harmonicMean ss+ AvgGeometricMean -> geometricMean ss+ AvgMedian ->+ let ss' = F.sort $ getSpeedsV cs+ len = V.length ss'+ mid = len `div` 2+ in if V.length ss' < 3+ then mean ss'+ else if odd len then ss' V.! mid else mean (V.slice mid 2 ss')+ AvgEndPoints -> fromMaybe 0 . join . fmap (speed x) $ e+ AvgMinOf as -> minimum $ map (getAvg cs) as+ AvgWith f -> f cs+ getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (getUTCTime a)+ getSpeedsV = V.fromList . getSpeeds+ getSpeeds zs = concatMap (maybeToList . uncurry speed) $ zip zs (drop 1 zs)++type TrailTransformation c = [Selected (Trail c)] -> Trail c++-- | A PointGrouping is a function that selects segments of a trail.+-- +-- Grouping point _does not_ result in deleted points. It is always true that:+--+-- forall g : PointGrouping c -->+-- concatMap unSelect (g ts) == ts+--+-- The purpose of grouping is usually for later processing. Any desire to drop+-- points that didn't meet a particular grouping criteria can be filled with+-- a composition with 'filter' (or directly via 'filterPoints'.+type PointGrouping c = Trail c -> [Selected (Trail c)]++type TransformGrouping c = [Selected (Trail c)] -> [Selected (Trail c)]++-- | When grouping points, lists of points are either marked as 'Select' or 'NotSelect'.+data Selected a = Select {unSelect :: a} | NotSelect {unSelect :: a}+ deriving (Eq, Ord, Show)++isSelected :: Selected a -> Bool+isSelected (Select _) = True+isSelected _ = False++isNotSelected :: Selected a -> Bool+isNotSelected = not . isSelected++selLength :: Selected [a] -> Int+selLength = length . unSelect++onSelected :: (a -> b) -> (a -> b) -> Selected a -> b+onSelected f _ (Select a) = f a+onSelected _ g (NotSelect a) = g a++instance Functor Selected where+ fmap f (Select x) = Select $ f x+ fmap f (NotSelect x) = NotSelect $ f x++dropExact :: Int -> [Selected [a]] -> [Selected [a]]+dropExact i [] = []+dropExact i (x:xs) =+ case compare (selLength x) i of+ EQ -> xs+ LT -> dropExact (i - selLength x) xs+ GT -> fmap (drop i) x : xs++-- | Groups trail segments into contiguous points within the speed+-- and all others outside of the speed. The "speed" from point p(i)+-- to p(i+1) is associated with p(i) (execpt for the first speed+-- value, which is associated with both the first and second point)+betweenSpeeds :: (Lat a, Lon a, Time a) => Double -> Double -> PointGrouping a+betweenSpeeds low hi ps =+ let spds = concatMap maybeToList $ zipWith speed ps (drop 1 ps)+ psSpds = [(p,s) | p <- ps, s <- maybeToList (listToMaybe spds) ++ spds]+ inRange x = x >= low && x <= hi+ chunk [] = []+ chunk xs@(x:_) =+ let op p = if inRange (snd x) then first Select . span p else first NotSelect . break p+ (r,rest) = op (inRange . snd) xs+ in r : chunk xs+ in map (fmap (map fst)) $ chunk psSpds++-- | A "rest point" means the coordinates remain within a given distance+-- for at least a particular amount of time.+restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> PointGrouping a+restLocations d s ps =+ let consToFirst x [] = [NotSelect [x]]+ consToFirst x (a:as) = (fmap (x:) a) : as+ go [] [] = []+ go [] nonRests = [NotSelect $ reverse nonRests]+ go (a:as) nonRests =+ case takeWhileEnd ((<=) d . distance a) as of+ (Just l, close, far) ->+ case (getUTCTime a, getUTCTime l) of+ (Just t1, Just t2) ->+ let diff = diffUTCTime t2 t1+ in if diff >= s then NotSelect (reverse nonRests) : Select (a:close) : go far [] else go as (a:nonRests)+ _ -> consToFirst a $ go as nonRests+ _ -> consToFirst a $ go as nonRests+ in go ps []+ +-- |chunking points into groups spanning at most the given time+-- interval.+spansTime :: (Lat a, Lon a, Time a) => NominalDiffTime -> PointGrouping a+spansTime n ps =+ let times = mkTimePair ps+ chunk [] = []+ chunk (x:xs) =+ let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs in good : chunk rest+ in map (Select . map fst) $ chunk times++-- | intersects the given groupings+intersectionOf :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a+intersectionOf gs ps =+ let groupings = map ($ ps) gs+ -- chunk :: [[Selected [pnts]]] -> pnts -> [pnts]+ chunk _ [] = []+ chunk ggs xs = + let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs -- FIXME this is all manner of broken+ sel = if all isSelected (concatMap (take 1) ggs) then Select else NotSelect+ (c,rest) = splitAt minLen xs+ in sel c : chunk (filter (not . null) $ map (dropExact minLen) ggs) rest+ in chunk groupings ps++-- | Union all the groupings+unionOf :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a+unionOf gs ps =+ let groupings = map ($ ps) gs+ chunk _ [] = []+ chunk ggs xs =+ let getSegs = concatMap (take 1)+ segs = getSegs ggs+ len =+ if any isSelected segs+ then max 1 . maximum . getSegs . map (map selLength) . map (filter isSelected) $ ggs+ else max 1 . minimum . getSegs . map (map selLength) $ ggs+ sel = if any isSelected segs then Select else NotSelect+ (c,rest) = splitAt len xs+ in sel c : chunk (filter (not . null) $ map (dropExact len) ggs) rest+ in chunk groupings ps+ +-- | Intersection binary operator+(/\) :: [Selected (Trail a)] -> TransformGrouping a+(/\) xs [] = []+(/\) [] ys = []+(/\) xsL@(Select x:xs) ysL@(Select y:ys) =+ let z = if length x < length y then x else y+ xs' = selListDrop (length z) xsL+ ys' = selListDrop (length z) ysL+ in Select z : (xs' /\ ys')+(/\) xs (NotSelect y:ys) = NotSelect y : (selListDrop (length y) xs /\ ys)+(/\) (NotSelect x:xs) ys = NotSelect x : (xs /\ selListDrop (length x) ys)++-- | Union binary operator+(\/) :: [Selected (Trail a)] -> TransformGrouping a+(\/) xs [] = xs+(\/) [] ys = ys+(\/) (Select x:xs) (Select y : ys) =+ let xLen = length x+ yLen = length y+ in if xLen < yLen+ then (Select y :) (selListDrop (yLen - xLen) xs \/ ys)+ else (Select x :) (xs \/ selListDrop (xLen - yLen) ys)+(\/) (Select x:xs) ys = Select x : selListDrop (length x) ys+(\/) xs (Select y:ys) = Select y : selListDrop (length y) xs+(\/) xsL@(NotSelect x:xs) ysL@(NotSelect y:ys) =+ let xLen = length x+ yLen = length y+ in if xLen < yLen+ then (NotSelect x:) (xs \/ selListDrop xLen ysL)+ else (NotSelect y:) (selListDrop yLen xsL \/ ys)+ +selListTake :: Int -> [Selected [a]] -> [Selected [a]]+selListTake 0 _ = []+selListTake n [] = []+selListTake n (x:xs) =+ let x' = take n (unSelect x)+ in fmap (const x') x : selListTake (n - length x') xs++selListDrop :: Int -> [Selected [a]] -> [Selected [a]]+selListDrop 0 xs = xs+selListDrop n [] = []+selListDrop n (x:xs) =+ let x' = drop n (unSelect x)+ in fmap (const x') x : selListDrop (n - (selLength x - length x')) xs++-- |Inverts the selected/nonselected segments+invertSelection :: TransformGrouping a+invertSelection = map (onSelected NotSelect Select)++-- |@firstGrouping f ps@ only the first segment remains 'Select'ed, and only+-- if it was already selected by @f@.+firstGrouping :: TransformGrouping a+firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps)++-- | Only the last segment, if any, is selected (note: the current+-- implementation is inefficient, using 'reverse')+lastGrouping :: TransformGrouping a+lastGrouping ps = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps')++-- | chunk the trail into groups of N points+everyNPoints :: Int -> PointGrouping a+everyNPoints n ps+ | n <= 0 = [NotSelect ps]+ | otherwise = go ps+ where+ go [] = []+ go xs = let (h,t) = splitAt n xs in Select h : go t+ +-- |For every selected group, refine the selection using the second+-- grouping method. This differs from 'IntersectionOf' by restarting+-- the second grouping algorithm at the beginning each group selected+-- by the first algorithm.+refineGrouping :: PointGrouping a -> TransformGrouping a+refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x]))++-- |Remove all points that remain 'NotSelect'ed by the given grouping algorithm.+filterPoints :: PointGrouping a -> Trail a -> Trail a+filterPoints g = concatMap unSelect . filter isSelected . g++-- Extract the time from each coordinate. If no time is available then+-- the coordinate is dropped!+mkTimePair :: (Time a) => Trail a -> [(a,UTCTime)]+mkTimePair xs =+ let timesM = map (\x-> fmap (x,) $ getUTCTime x) xs+ in concatMap maybeToList timesM++-- |Construct a bezier curve using the provided trail. Construct a+-- new trail by sampling the given bezier curve at the given times.+-- The current implementation assumes the times of the input+-- coordinates are available and all equal (Ex: all points are 5+-- seconds apart), the results will be poor if this is not the case!+bezierCurveAt :: (Lat a, Lon a, Time a) => [UTCTime] -> Trail a -> Trail a+bezierCurveAt selectedTimes xs = + let timesDef = mkTimePair xs+ end = last timesDef+ 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++bezierPoint :: (Lat a, Lon a) => [a] -> Double -> a+bezierPoint pnts t = go pnts+ where+ go [] = error "GPS Package: Can not create a bezier point from an empty list"+ go [p] = p+ go ps = interpolate (go (init ps)) (go (tail ps)) t++-- |Interpolate selected points onto a bezier curve. Note this gets+-- exponentially more expensive with the length of the segement being+-- transformed - it is not advisable to perform this operation on+-- trail segements with more than ten points!+bezierCurve :: (Lat a, Lon a, Time a) => [Selected (Trail a)] -> Trail a+bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id)++-- |Filters out any points that go backward in time (thus must not be+-- valid if this is a trail)+linearTime :: (Lon a, Lat a, Time a) => [a] -> [a]+linearTime [] = []+linearTime (p:ps) = go (getUTCTime p) ps+ where+ go _ [] = []+ go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps++-- |Returns the closest distance between two trails (or Nothing if a+-- trail is empty). Inefficient implementation:+-- O( (n * m) * log (n * m) )+closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance+closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]++-- | Find the total distance traveled+totalDistance :: (Lat a, Lon a) => [a] -> Distance+totalDistance as = sum $ zipWith distance as (drop 1 as)++-- | Uses Grahams scan to compute the convex hull of the given points.+-- This operation requires sorting of the points, so don't try it unless+-- you have notably more memory than the list of points will consume.+convexHull :: (Eq c, Lat c, Lon c) => [c] -> [c]+convexHull xs =+ let first = southMost xs+ in case first of+ Nothing -> []+ Just f ->+ let sorted = sortBy (comparing (eastZeroHeading f)) (filter (/= f) xs)+ in case sorted of+ (a:b:cs) -> grahamScan (b:a:f:[]) cs+ cs -> f : cs+ where+ grahamScan [] _ = []+ grahamScan ps [] = ps+ grahamScan (x:[]) _ = [x]+ grahamScan (p2:p1:ps) (x:xs) =+ case turn p1 p2 x of+ LeftTurn -> grahamScan (x:p2:p1:ps) xs+ Straight -> grahamScan (x:p2:p1:ps) xs+ _ -> grahamScan (p1:ps) (x:xs)++eastZeroHeading :: (Lat c, Lon c) => c -> c -> Heading+eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s++data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)++turn :: (Lat c, Lon c) => c -> c -> c -> Turn+turn a b c =+ let h1 = eastZeroHeading a b+ h2 = eastZeroHeading b c+ d = h2 - h1+ in if d >= 0 && d < pi then LeftTurn else RightTurn++-- | Find the southmost point+southMost :: (Lat c) => [c] -> Maybe c+southMost [] = Nothing+southMost cs = Just . minimumBy (comparing lat) $ cs++---------- COMPOSIT OPERATIONS ---------------+-- These operations are simply implemented using the previously+-- defined functions. They can serve either for concise use for novice+-- users or as instructional examples.+------------------------------------------++smoothRests :: (Lat a, Lon a, Time a) => Trail a -> Trail a+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''++smoothPath :: (Lat a, Lon a, Time a) => Trail a -> Trail a+smoothPath ps = undefined+ +slidingWindow :: Int -> Int -> ([a] -> [a]) -> [a] -> [a]+slidingWindow width step f trail = go trail+ where+ go [] = []+ go xs =+ let (window,e) = splitAt width xs+ (hT,eT) = splitAt step (f window)+ in hT ++ go (eT ++ e)
Setup.hs view
gps.cabal view
@@ -1,10 +1,10 @@ name: gps-version: 0.7+version: 0.8 license: BSD3 license-file: LICENSE author: Thomas DuBuisson <thomas.dubuisson@gmail.com> maintainer: Thomas DuBuisson-description: Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations+description: Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations. NOTE: Version range 0.8.* won't strictly follow PVP - I will be adding additional functions in minor releases 0.8.x. synopsis: For manipulating GPS coordinates and trails. category: Data stability: stable@@ -16,7 +16,13 @@ Library Build-Depends: base >= 3 && < 6, pretty >= 1.0 , prettyclass >= 1.0,- time >= 1.1 && < 1.3, GPX == 0.4.*, hxt >= 8.5, xsd == 0.3.*+ time >= 1.1, GPX == 0.4.*, hxt >= 8.5, 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:++source-repository head+ type: git+ location: git@github.com:TomMD/gps.git