gps 1.0.3 → 1.1

raw patch · 8 files changed

+673/−754 lines, 8 filesdep +gpx-conduitdep −data-lensdep ~GPXdep ~hxtdep ~timesetup-changedPVP ok

version bump matches the API change (PVP)

Dependencies added: gpx-conduit

Dependencies removed: data-lens

Dependency ranges changed: GPX, hxt, time, xsd

API changes (from Hackage documentation)

Files

• Data/GPS.hs +0/−12
• Data/GPS/Core.hs +0/−266
• Data/GPS/Trail.hs +0/−468
• Geo/Computations.hs +12/−0
• Geo/Computations/Basic.hs +205/−0
• Geo/Computations/Trail.hs +446/−0
• Setup.hs +0/−0
• gps.cabal +10/−8

@@ -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

@@ -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]

@@ -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)

@@ -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

@@ -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

@@ -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

@@ -1,5 +1,5 @@ name:		gps-version:	1.0.3+version:	1.1 license:	BSD3 license-file:	LICENSE author:		Thomas DuBuisson <thomas.dubuisson@gmail.com>@@ -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