diff --git a/Data/GPS.hs b/Data/GPS.hs deleted file mode 100644 --- a/Data/GPS.hs +++ /dev/null @@ -1,12 +0,0 @@ -{-# 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 - ( module Data.GPS.Core - , module Data.GPS.Trail - ) where - -import Data.GPS.Core -import Data.GPS.Trail diff --git a/Data/GPS/Core.hs b/Data/GPS/Core.hs deleted file mode 100644 --- a/Data/GPS/Core.hs +++ /dev/null @@ -1,266 +0,0 @@ -{-# LANGUAGE TupleSections #-} -module Data.GPS.Core - ( -- * Types - Distance - , Heading - , Speed - , Vector - , Trail - , Circle - , Arc - , Coordinate (..) - -- * Constants - , north - , south - , east - , west - , radiusOfEarth - , circumferenceOfEarth - -- * Coordinate Functions - , heading - , distance - , speed - , getVector - , addVector - , getRadianPair - , getDMSPair - , divideArea - , interpolate - , circleIntersectionPoints - , intersectionArcsOf - , maximumDistanceOfArc - -- * IO helpers - , writeGPX - , readGPX - , readGPXSegments - -- * Utility - , getUTCTime - , module Data.Geo.GPX - ) where - -import Data.Time -import Data.Maybe -import Data.List (sortBy) -import Data.Ord (comparing) -import Control.Monad -import Text.XML.HXT.Core -import Text.XML.XSD.DateTime(DateTime,toUTCTime) -import Data.Geo.GPX -import Data.Lens.Common - -class (LatL a, LonL a) => Coordinate a where - lat :: a -> Double - lat = runLatitude . (^. latL) - lon :: a -> Double - lon = runLongitude . (^. lonL) - -instance Coordinate Wpt -instance Coordinate Pt - --- |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) - --- | Genearlly a circle indicates a known area in which we are searching --- (so a center point and maximum possible distance from that point) -type Circle a = (a, Distance) - --- | An arc is represented as a circle, starting heading and ending heading -type Arc a = (Circle a, Heading, Heading) - -type Trail a = [a] - -getUTCTime :: (TimeL a) => a -> Maybe UTCTime -getUTCTime = fmap toUTCTime . (^. timeL) - -distance :: (Coordinate a, Coordinate b) => a -> b -> Distance -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 :: (Coordinate a, Coordinate 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 = lon2 - lon1 - -getVector :: (Coordinate a, Coordinate 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 :: (Coordinate c) => Vector -> c -> c -addVector (d,h) p = (lonL ^= longitude (toDegrees lon2)) - . (latL ^= latitude (toDegrees lat2)) - $ p - where - (lat,lon) = getRadianPairD p - 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 :: (Coordinate loc, TimeL loc, Coordinate b, TimeL 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 - --- |Circumference of earth (meters) -circumferenceOfEarth :: Double -circumferenceOfEarth = radiusOfEarth * 2 * pi - --- |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 :: (Coordinate c) => c -> (Double,Double) -getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair - -getDMSPair :: (Coordinate c) => c -> (Latitude, Longitude) -getDMSPair c = (c ^. latL, c ^. lonL) - --- |Provides a lat/lon pair of doubles in radians -getRadianPair :: (Coordinate p) => p -> (Latitude, Longitude) -getRadianPair p = (toRadians (p ^. latL), toRadians (p ^. lonL)) - -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 :: (Coordinate 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 - --- | Compute the points at which two circles intersect (assumes a flat plain). If --- the circles do not intersect or are identical then the result is @Nothing@. -circleIntersectionPoints :: (Coordinate a) => (a, Distance) -> (a, Distance) -> Maybe (a,a) -circleIntersectionPoints (a,r1) (b,r2) - | a ^. latL == b ^. latL && a ^. lonL == b ^. lonL && r1 == r2 = Nothing -- FIXME need approx eq - | r1 + r2 < ab = Nothing - | any isNaN (map (^. latL) pts) || any isNaN (map (^. lonL) pts) = Nothing - | otherwise = Just (p1, p2) - where - ab = distance a b - angABX = acos ( (r1^2 + ab^2 - r2^2) / (2 * r1 * ab) ) - ang1 = heading a b + angABX - ang2 = heading a b - angABX - p1 = addVector (r1, ang1) a - p2 = addVector (r1, ang2) a - pts = [p1,p2] - --- | Find the area in which all given circles intersect. The resulting --- area is described in terms of the bounding arcs. All cirlces must --- intersect at two points. -intersectionArcsOf :: (Coordinate a) => [Circle a] -> [Arc a] -intersectionArcsOf cs = - let isArcWithinCircle circ arc = maximumDistanceOfArc (fst circ) arc <= (snd circ) - isArcWithinAllCircles arc = all ($ arc) (map isArcWithinCircle cs) - -- getArcs :: Circle a -> Circle a -> [Arc a] - getArcs c1 c2 = concatMap (buildArcsFromPoints c1 c2) . maybeToList $ circleIntersectionPoints c1 c2 - -- buildArcsFromPoints :: (a, a) -> [Arc a] - buildArcsFromPoints c1 c2 (p1,p2) = - let c1h1 = heading (fst c1) p1 - c1h2 = heading (fst c1) p2 - c2h1 = heading (fst c2) p1 - c2h2 = heading (fst c2) p2 - in [(c1,c1h1,c1h2), (c1,c1h2, c1h1), (c2,c2h1,c2h2), (c2,c2h2,c2h1)] - in filter isArcWithinAllCircles . concatMap (uncurry getArcs) . choose2 $ cs - -maximumDistanceOfArc :: (Coordinate a) => a -> Arc a -> Distance -maximumDistanceOfArc pnt ((c,r), h1, h2) = - let pcHeading = heading pnt c - in if ((pcHeading < h1 || pcHeading > h2) && h1 < h2) || ((pcHeading > h2 && pcHeading < h1) && h1 > h2) - then max (distance pnt (addVector (r,h1) c)) (distance pnt (addVector (r,h2) c)) - else distance pnt c + r - -choose2 :: [a] -> [(a,a)] -choose2 [] = [] -choose2 (x:xs) = map (x,) xs ++ choose2 xs - --- |@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 :: (Coordinate c) => Distance -> Distance -> c -> c -> [[c]] -divideArea vDist hDist nw se = - let (top,left) = (nw ^. latL, nw ^. lonL) - (btm,right) = (se ^. latL, se ^. lonL) - 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 Wpt) -readGPX = liftM (concatMap (^. trkptsL). concatMap (^. trksegsL) . concatMap (^. trksL)) . readGpxFile - -writeGPX :: FilePath -> Trail Wpt -> IO () -writeGPX fp ps = writeGpxFile fp $ gpx "1.0" "Haskell GPS Package (via the GPX package)" Nothing [] [] [trk Nothing Nothing Nothing Nothing [] Nothing Nothing Nothing [trkseg 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 Wpt] -readGPXSegments = liftM (map (concatMap (^. trkptsL)) . map (^. trksegsL) . concatMap (^. trksL)) . readGpxFile - -readGpxFile :: FilePath -> IO [GPX] -readGpxFile = runX . xunpickleDocument (xpickle :: PU GPX) [withRemoveWS yes, withValidate no] diff --git a/Data/GPS/Trail.hs b/Data/GPS/Trail.hs deleted file mode 100644 --- a/Data/GPS/Trail.hs +++ /dev/null @@ -1,468 +0,0 @@ -{-# LANGUAGE TupleSections #-} -module Data.GPS.Trail - ( -- * Types - AvgMethod(..) - , Selected(..) - , PointGrouping - , TransformGrouping - -- * Utility Functions - , isSelected - , isNotSelected - , onSelected - , selLength - -- * Trail Functions - -- ** Queries - , totalDistance - , totalTime - , 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 - , smoothSome - , smoothMore - -- * 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') -import Data.Function (on,fix) -import Data.List as L -import Data.Maybe -import Data.Ord -import Data.Time -import Data.Lens.Common - -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 :: (Coordinate a, TimeL 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 :: (Coordinate a, TimeL a) => - AvgMethod a -> NominalDiffTime -> Trail a -> [(UTCTime, Speed)] -slidingAverageSpeed _ _ [] = [] -slidingAverageSpeed m minTime xs = - let pts = map unSelect (spansTime minTime xs) - spds = map (getAvg m) pts - times = map getAvgTimes pts - in concatMap maybeToList $ zipWith (\t s -> fmap (,s) t) times spds - where - getTimeDiff a b = on (liftM2 diffUTCTime) getUTCTime a b - - -- getAvg :: [] -> AvgMethod -> Speed - getAvg _ [] = 0 - getAvg _ [x] = 0 - getAvg m cs = - 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 $ speed (head cs) (last cs) - AvgMinOf as -> minimum $ map (flip getAvg cs) as - AvgWith f -> f cs - getAvgTimes [] = Nothing - getAvgTimes [x] = getUTCTime x - getAvgTimes ps = getAvgTime (head ps) (last ps) - getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (getUTCTime a) - getSpeedsV = V.fromList . getSpeeds - getSpeeds zs = concatMap maybeToList $ zipWith speed zs (drop 1 zs) - --- | 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)] - --- | Given a selection of coordinates, transform the selected --- coordinates in some way (while leaving the non-selected --- coordinates unaffected). -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 :: (Coordinate a, TimeL 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 :: (Coordinate a, TimeL 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 :: (Coordinate a, TimeL a) => NominalDiffTime -> PointGrouping a -spansTime n ps = - let times = mkTimePair ps - chunk [] = [] - chunk xs@(x:_) = - let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs - in if null good then [xs] else good : chunk rest - in map (Select . map fst) $ chunk times - --- | intersects the given groupings -intersectionOf :: (Coordinate a, TimeL 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 :: (Coordinate a, TimeL 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 :: (TimeL 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 :: (Coordinate a, TimeL a) => [UTCTime] -> Trail a -> Trail a -bezierCurveAt _ [] = [] -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 - diffTimes = [diffUTCTime t (snd top) / totalTime | t <- times] - queryTimes = map realToFrac diffTimes - in if totalTime <= 0 || any (\x -> x < 0 || x > 1) queryTimes - then xs -- error "bezierCurveAt has a out-of-bound time!" - else - if null timesDef || 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 ((timeL ^=) . Just . fromUTCTime) newTimes curvePoints - -bezierPoint :: (Coordinate 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 :: (Coordinate a, TimeL 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 :: (Coordinate a, TimeL 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 :: (Coordinate 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 :: (Coordinate a) => [a] -> Distance -totalDistance as = sum $ zipWith distance as (drop 1 as) - -totalTime :: TimeL a => Trail a -> NominalDiffTime -totalTime [] = 0 -totalTime xs@(x:_) = fromMaybe 0 $ liftM2 diffUTCTime (getUTCTime x) (getUTCTime $ last xs) - --- | 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, Coordinate c) => [c] -> [c] -convexHull xs = - let first = southMost xs - in case first of - Nothing -> [] - Just f -> - let sorted = L.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 :: (Coordinate 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 :: (Coordinate 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 :: (LatL c) => [c] -> Maybe c -southMost [] = Nothing -southMost cs = Just . minimumBy (comparing (^. latL)) $ 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 :: (Coordinate a, TimeL a) => Trail a -> Trail a -smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60 - -smoothSome :: (Coordinate a, TimeL a) => Trail a -> Trail a -smoothSome = gSmoothSome 7 - -gSmoothSome n = bezierCurve . everyNPoints n - -gSmoothMore n k ps = - let ps' = gSmoothSome n ps - (h,t) = splitAt k ps' - in h ++ gSmoothSome n t - -smoothMore :: (Coordinate a, TimeL a) => Trail a -> Trail a -smoothMore - = gSmoothMore 7 3 . gSmoothMore 3 1 . gSmoothMore 5 2 - . gSmoothMore 3 1 . gSmoothMore 5 2 . gSmoothMore 7 3 - -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) diff --git a/Geo/Computations.hs b/Geo/Computations.hs new file mode 100644 --- /dev/null +++ b/Geo/Computations.hs @@ -0,0 +1,12 @@ +{-# 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 Geo.Computations + ( module Geo.Computations.Basic + , module Geo.Computations.Trail + ) where + +import Geo.Computations.Basic +import Geo.Computations.Trail diff --git a/Geo/Computations/Basic.hs b/Geo/Computations/Basic.hs new file mode 100644 --- /dev/null +++ b/Geo/Computations/Basic.hs @@ -0,0 +1,205 @@ +{-# LANGUAGE TupleSections #-} +module Geo.Computations.Basic + ( -- * Types + Distance + , Heading + , Speed + , Vector + , Trail + , Circle + , Arc + , Point(..) + , pt + -- * Constants + , north + , south + , east + , west + , radiusOfEarth + , circumferenceOfEarth + -- * Coordinate Functions + , heading + , distance + , speed + , getVector + , addVector + , interpolate + , circleIntersectionPoints + , intersectionArcsOf + , maximumDistanceOfArc + -- * IO helpers + , readGPXFile + ) where + +import Data.Time +import Data.Maybe +import Geo.GPX.Conduit + +-- |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) + +-- | Genearlly a circle indicates a known area in which we are searching +-- (so a center point and maximum possible distance from that point) +type Circle a = (a, Distance) + +-- | An arc is represented as a circle, starting heading and ending heading +type Arc a = (Circle a, Heading, Heading) + +type Trail a = [a] + +distance :: Point -> Point -> Distance +distance x y = + let (lat1,lon1) = getRadianPair x + (lat2,lon2) = getRadianPair y + deltaLat = lat2 - lat1 + deltaLon = lon2 - lon1 + a = (sin (deltaLat / 2))^(2::Int) + cos lat1 * cos lat2 * (sin (deltaLon / 2))^(2::Int) + 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 :: Point -> Point -> Heading +heading a b = + atan2 (sin (diffLon) * cos (lat2)) + (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon)) + where + (lat1, lon1) = getRadianPair a + (lat2, lon2) = getRadianPair b + diffLon = lon2 - lon1 + +getVector :: Point -> Point -> 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 :: Vector -> Point -> Point +addVector (d,h) p = + p { pntLon = toDegrees lon2 + , pntLat = toDegrees lat2 + } + where + (lat,lon) = getRadianPair p + 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 :: Point -> Point -> Maybe Speed +speed a b = + case (pntTime b, pntTime 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 + +-- |Circumference of earth (meters) +circumferenceOfEarth :: Double +circumferenceOfEarth = radiusOfEarth * 2 * pi + +-- |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 :: Double -> Double +toDegrees = (*) (180 / pi) + +-- get latitude and longituide in Radians as Double's +getRadianPair :: Point -> (Double,Double) +getRadianPair p = (toRadians (pntLat p), toRadians (pntLon 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 :: Point -> Point -> Double -> Point +interpolate c1 c2 w + | w < 0 = c1 + | w > 1 = c2 + | otherwise = + let (h,d) = (heading c1 c2, distance c1 c2) + v = (d * w, h) + in addVector v c1 + +-- | Compute the points at which two circles intersect (assumes a flat plain). If +-- the circles do not intersect or are identical then the result is @Nothing@. +circleIntersectionPoints :: (Point, Distance) -> (Point, Distance) -> Maybe (Point,Point) +circleIntersectionPoints (a,r1) (b,r2) + | pntLat a == pntLat b && pntLon a == pntLon b && r1 == r2 = Nothing -- FIXME need approx eq + | r1 + r2 < ab = Nothing + | any isNaN (map pntLat pts) || any isNaN (map pntLon pts) = Nothing + | otherwise = Just (p1, p2) + where + ab = distance a b + angABX = acos ( (r1^(2::Int) + ab^(2::Int) - r2^(2::Int)) / (2 * r1 * ab) ) + ang1 = heading a b + angABX + ang2 = heading a b - angABX + p1 = addVector (r1, ang1) a + p2 = addVector (r1, ang2) a + pts = [p1,p2] + +-- | Find the area in which all given circles intersect. The resulting +-- area is described in terms of the bounding arcs. All cirlces must +-- intersect at two points. +intersectionArcsOf :: [Circle Point] -> [Arc Point] +intersectionArcsOf cs = + let isArcWithinCircle circ arc = maximumDistanceOfArc (fst circ) arc <= (snd circ) + isArcWithinAllCircles arc = all ($ arc) (map isArcWithinCircle cs) + -- getArcs :: Circle a -> Circle a -> [Arc a] + getArcs c1 c2 = concatMap (buildArcsFromPoints c1 c2) . maybeToList $ circleIntersectionPoints c1 c2 + -- buildArcsFromPoints :: (a, a) -> [Arc a] + buildArcsFromPoints c1 c2 (p1,p2) = + let c1h1 = heading (fst c1) p1 + c1h2 = heading (fst c1) p2 + c2h1 = heading (fst c2) p1 + c2h2 = heading (fst c2) p2 + in [(c1,c1h1,c1h2), (c1,c1h2, c1h1), (c2,c2h1,c2h2), (c2,c2h2,c2h1)] + in filter isArcWithinAllCircles . concatMap (uncurry getArcs) . choose2 $ cs + +maximumDistanceOfArc :: Point -> Arc Point -> Distance +maximumDistanceOfArc pnt ((c,r), h1, h2) = + let pcHeading = heading pnt c + in if ((pcHeading < h1 || pcHeading > h2) && h1 < h2) || ((pcHeading > h2 && pcHeading < h1) && h1 > h2) + then max (distance pnt (addVector (r,h1) c)) (distance pnt (addVector (r,h2) c)) + else distance pnt c + r + +choose2 :: [a] -> [(a,a)] +choose2 [] = [] +choose2 (x:xs) = map (x,) xs ++ choose2 xs diff --git a/Geo/Computations/Trail.hs b/Geo/Computations/Trail.hs new file mode 100644 --- /dev/null +++ b/Geo/Computations/Trail.hs @@ -0,0 +1,446 @@ +{-# LANGUAGE TupleSections #-} +module Geo.Computations.Trail + ( -- * Types + AvgMethod(..) + , Selected(..) + , PointGrouping + , TransformGrouping + -- * Utility Functions + , isSelected + , isNotSelected + , onSelected + , selLength + -- * Trail Functions + -- ** Queries + , totalDistance + , totalTime + , 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 + , smoothTrail + -- * Misc + , bezierPoint + ) where + +import Text.Show.Functions () +import Geo.Computations.Basic +import Geo.GPX.Conduit () + +import Control.Arrow (first) +import Control.Monad +import Data.Fixed (mod') +import Data.Function (on) +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',xs2,zs) = go as (Just a) in (e',a:xs2,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 :: NominalDiffTime -> Trail Point -> [(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 :: AvgMethod Point -> NominalDiffTime -> Trail Point -> [(UTCTime, Speed)] +slidingAverageSpeed _ _ [] = [] +slidingAverageSpeed m minTime xs = + let pts = map unSelect (spansTime minTime xs) + spds = map (getAvg m) pts + times = map getAvgTimes pts + in concatMap maybeToList $ zipWith (\t s -> fmap (,s) t) times spds + where + getTimeDiff a b = on (liftM2 diffUTCTime) pntTime a b + + -- getAvg :: [] -> AvgMethod -> Speed + getAvg _ [] = 0 + getAvg _ [_] = 0 + getAvg m2 cs = + let ss = getSpeedsV cs + in case m2 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 $ speed (head cs) (last cs) + AvgMinOf as -> minimum $ map (flip getAvg cs) as + AvgWith f -> f cs + getAvgTimes [] = Nothing + getAvgTimes [x] = pntTime x + getAvgTimes ps = getAvgTime (head ps) (last ps) + getAvgTime a b = liftM2 addUTCTime (getTimeDiff b a) (pntTime a) + getSpeedsV = V.fromList . getSpeeds + getSpeeds zs = concatMap maybeToList $ zipWith speed zs (drop 1 zs) + +-- | 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)] + +-- | Given a selection of coordinates, transform the selected +-- coordinates in some way (while leaving the non-selected +-- coordinates unaffected). +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 _ [] = [] +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 :: Double -> Double -> PointGrouping Point +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 rest + 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 :: Distance -> NominalDiffTime -> PointGrouping Point +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 (pntTime a, pntTime 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 :: NominalDiffTime -> PointGrouping Point +spansTime n ps = + let times = mkTimePair ps + chunk [] = [] + chunk xs@(x:_) = + let (good,rest) = span ((<= addUTCTime n (snd x)) . snd) xs + in if null good then [xs] else good : chunk rest + in map (Select . map fst) $ chunk times + +-- | intersects the given groupings +intersectionOf :: [PointGrouping Point] -> PointGrouping Point +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 :: [PointGrouping Point] -> PointGrouping Point +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 +(/\) _ [] = [] +(/\) [] _ = [] +(/\) xsL@(Select x:_) ysL@(Select y:_) = + 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:_) ys = Select x : selListDrop (length x) ys +(\/) xs (Select y:_) = 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) + +selListDrop :: Int -> [Selected [a]] -> [Selected [a]] +selListDrop 0 xs = xs +selListDrop _ [] = [] +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 :: Trail Point -> [(Point,UTCTime)] +mkTimePair xs = + let timesM = map (\x-> fmap (x,) $ pntTime 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 :: [UTCTime] -> Trail Point -> Trail Point +bezierCurveAt _ [] = [] +bezierCurveAt selectedTimes xs = + let timesDef = mkTimePair xs + end = last timesDef + top = head timesDef + tTime = diffUTCTime (snd end) (snd top) + times = if null selectedTimes then map snd timesDef else selectedTimes + diffTimes = [diffUTCTime t (snd top) / tTime | t <- times] + queryTimes = map realToFrac diffTimes + in if tTime <= 0 || any (\x -> x < 0 || x > 1) queryTimes + then xs -- error "bezierCurveAt has a out-of-bound time!" + else + if null timesDef || 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 (\t p -> p { pntTime = Just t}) newTimes curvePoints + +bezierPoint :: [Point] -> Double -> Point +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 :: [Selected (Trail Point)] -> Trail Point +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 :: [Point] -> [Point] +linearTime [] = [] +linearTime (p:ps) = go (pntTime p) ps + where + go _ [] = [] + go t (x:xs) = if pntTime x < t then go t xs else x : go (pntTime x) xs + +-- |Returns the closest distance between two trails (or Nothing if a +-- trail is empty). Inefficient implementation: +-- O( (n * m) * log (n * m) ) +closestDistance :: Trail Point -> Trail Point -> Maybe Distance +closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs] + +-- | Find the total distance traveled +totalDistance :: [Point] -> Distance +totalDistance as = sum $ zipWith distance as (drop 1 as) + +totalTime :: Trail Point -> NominalDiffTime +totalTime [] = 0 +totalTime xs@(x:_) = fromMaybe 0 $ liftM2 diffUTCTime (pntTime x) (pntTime $ last xs) + +-- | 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 :: [Point] -> [Point] +convexHull lst = + let frst = southMost lst + in case frst of + Nothing -> [] + Just f -> + let sorted = L.sortBy (comparing (eastZeroHeading f)) (filter (/= f) lst) + 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 :: Point -> Point -> Heading +eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s + +data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum) + +turn :: Point -> Point -> Point -> 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 :: [Point] -> Maybe Point +southMost [] = Nothing +southMost cs = Just . minimumBy (comparing pntLat) $ 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. +------------------------------------------ + +-- | Smooth points with rest areas using a bezierCurve. +-- +-- Parameters: rest for 1 minute within 30 meters get smoothed +-- in a bezier curve over every 8 points. +smoothRests :: Trail Point -> Trail Point +smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60 + +-- |Smooth every 7 points using a bezier curve +smoothTrail :: Trail Point -> Trail Point +smoothTrail = gSmoothSome 7 + +-- |Smooth every n points using a bezier curve +gSmoothSome :: Int -> Trail Point -> Trail Point +gSmoothSome n = bezierCurve . everyNPoints n diff --git a/Setup.hs b/Setup.hs --- a/Setup.hs +++ b/Setup.hs diff --git a/gps.cabal b/gps.cabal --- a/gps.cabal +++ b/gps.cabal @@ -1,5 +1,5 @@ name: gps -version: 1.0.3 +version: 1.1 license: BSD3 license-file: LICENSE author: Thomas DuBuisson @@ -15,14 +15,16 @@ Library Build-Depends: base >= 3 && < 6, - pretty >= 1.0 , prettyclass >= 1.0, - time >= 1.1, GPX >= 0.7 && < 0.8 - , hxt >= 9.1 && < 9.3, xsd >= 0.3 && < 0.4, - vector >= 0.7, statistics >= 0.9, data-lens >= 2.0 && < 2.1 + pretty >= 1.0, + prettyclass >= 1.0, + time >= 1.4, + gpx-conduit >= 0.1, + vector >= 0.7, + statistics >= 0.9 hs-source-dirs: - exposed-modules: Data.GPS - other-modules: Data.GPS.Core Data.GPS.Trail - ghc-options: + exposed-modules: Geo.Computations + other-modules: Geo.Computations.Basic, Geo.Computations.Trail + ghc-options: -Wall test-suite gps-tests hs-source-dirs: Test