gps 0.8.4 → 1.2
raw patch · 11 files changed
Files
- Data/GPS.hs +0/−12
- Data/GPS/Core.hs +0/−192
- Data/GPS/Trail.hs +0/−467
- Geo/Computations.hs +14/−0
- Geo/Computations/Basic.hs +203/−0
- Geo/Computations/Trail.hs +447/−0
- Geo/Types.hs +35/−0
- LICENSE +30/−0
- Setup.hs +0/−0
- Test/GpsTest.hs +19/−45
- gps.cabal +25/−22
− Data/GPS.hs
@@ -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
− Data/GPS/Core.hs
@@ -1,192 +0,0 @@-{-# LANGUAGE CPP #-}-module Data.GPS.Core- ( -- * Types- Distance- , Heading- , Speed- , Vector- , Trail- -- * Constants- , north- , south- , east- , west- , radiusOfEarth- , circumferenceOfEarth- -- * Coordinate Functions- , heading- , distance- , speed- , getVector- , addVector- , getRadianPair- , getDMSPair- , divideArea- , interpolate- -- * IO helpers- , writeGPX- , readGPX- , readGPXSegments- -- * Utility- , getUTCTime- , module Data.Geo.GPX- ) where--import Data.Time-import Data.Maybe-import Control.Monad-import Text.XML.HXT.Core-import Text.XML.XSD.DateTime(DateTime,toUTCTime)-import Data.Geo.GPX---- |Distances are expressed in meters-type Distance = Double---- |Angles are expressed in radians from North.--- 0 == North--- pi/2 == West--- pi == South--- (3/2)pi == East == - (pi / 2)-type Heading = Double--type Angle = Double--- |Speed is hard coded as meters per second-type Speed = Double-type Vector = (Distance, Heading)--type Trail a = [a]--getUTCTime :: (Time a) => a -> Maybe UTCTime-getUTCTime = fmap toUTCTime . time--distance :: (Lat a, Lon a, Lat b, Lon 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 :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Heading-heading a b =- atan2 (sin (diffLon) * cos (lat2))- (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon))- where- (lat1, lon1) = getRadianPairD a- (lat2, lon2) = getRadianPairD b- diffLon = lon2 - lon1--getVector :: (Lat a, Lon a, Lat b, Lon b) => a -> b -> Vector-getVector a b = (distance a b, heading a b)---- |Given a vector and coordinate, computes a new coordinate.--- Within some epsilon it should hold that if------ @dest = addVector (dist,heading) start@------ then------ @heading == heading start dest@--- --- @dist == distance start dest@-addVector :: (Lat c, Lon c) => Vector -> c -> c-addVector (d,h) p = setLon (longitudeType $ toDegrees lon2) - . setLat (latitudeType $ toDegrees lat2) $ p- where- (lat,lon) = getRadianPairD p- lat2 = asin (sin (lat) * cos (d / radiusOfEarth) + cos(lat) - * sin(d/radiusOfEarth) * cos h)- lon2 = lon + atan2 (sin h * sin (d / radiusOfEarth) * cos lat)- (cos (d/radiusOfEarth) - sin lat * sin lat2)---- | Speed in meters per second, only if a 'Time' was recorded for each waypoint.-speed :: (Lat loc, Lon loc, Time loc, Lat b, Lon b, Time b) => loc -> b -> Maybe Speed-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 earht (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 :: (Lat c, Lon c) => c -> (Double,Double)-getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair--getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)-getDMSPair c = (lat c, lon c)---- |Provides a lat/lon pair of doubles in radians-getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType)-getRadianPair p = (toRadians (lat p), toRadians (lon p))--toRadians :: Floating f => f -> f-toRadians = (*) (pi / 180)---- | @interpolate c1 c2 w@ where 0 <= w <= 1 Gives a point on the line--- between c1 and c2 equal to @c1 when @w == 0@ (weighted linearly--- toward c2).-interpolate :: (Lat a, Lon a) => a -> a -> Double -> a-interpolate c1 c2 w- | w < 0 || w > 1 = error "Interpolate only works with a weight between zero and one"- | otherwise = - let (h,d) = (heading c1 c2, distance c1 c2)- v = (d * w, h)- in addVector v c1- --- |@divideArea vDist hDist nw se@ divides an area into a grid of equally--- spaced coordinates within the box drawn by the northwest point (nw) and--- southeast point (se). Because this uses floating point there might be a--- different number of points in some rows (the last might be too far east based--- on a heading from the se point).-divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]]-divideArea vDist hDist nw se =- let (top,left) = (lat nw, lon nw)- (btm,right) = (lat se, lon se)- columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw- buildRow = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east))- in map buildRow columnOne---- |Reads a GPX file (using the GPX library) by simply concatenating all the--- tracks, segments, and points ('trkpts', 'trksegs', 'trks') into a single 'Trail'.-readGPX :: FilePath -> IO (Trail WptType)-readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile--writeGPX :: FilePath -> Trail WptType -> IO ()-writeGPX fp ps = writeGpxFile fp $ gpx $ gpxType "1.0" "Haskell GPS Package (via the GPX package)" Nothing [] [] [trkType Nothing Nothing Nothing Nothing [] Nothing Nothing Nothing [trksegType ps Nothing]] Nothing---- writeGpxFile should go in the GPX package-writeGpxFile :: FilePath -> Gpx -> IO ()-writeGpxFile fp gpx = runX_ (constA gpx >>> xpickleDocument (xpickle :: PU Gpx) [] fp)--runX_ t = runX t >> return ()--readGPXSegments :: FilePath -> IO [Trail WptType]-readGPXSegments = liftM (map (concatMap trkpts) . map trksegs . concatMap trks) . readGpxFile
− Data/GPS/Trail.hs
@@ -1,467 +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 Statistics.Function as F-import Statistics.Sample-import qualified Data.Vector.Unboxed as V--takeWhileEnd :: (a -> Bool) -> [a] -> (Maybe a, [a],[a])-takeWhileEnd p xs = go xs Nothing- where- go [] e = (e, [], [])- go (a:as) e- | p a = let (e',xs,zs) = go as (Just a) in (e',a:xs,zs)- | otherwise = (e,[], a:as)--data AvgMethod c- = AvgMean -- ^ Obtain the 'mean' of the considered points- | AvgHarmonicMean -- ^ Obtain the 'harmonicMean'- | AvgGeometricMean -- ^ Obtain the 'geometricMean'- | AvgMedian -- ^ Obtain the median of the considered points- | AvgEndPoints -- ^ Compute the speed considering only the given endpoints- | AvgMinOf [AvgMethod c] -- ^ Take the minimum of the speeds from the given methods- | AvgWith ([c] -> Speed)- --- | @avgSpeeds n points@--- Average speed using a window of up to @n@ seconds and averaging by taking the--- Median ('AvgMedian').-avgSpeeds :: (Lat a, Lon a, Time a) => NominalDiffTime -> Trail a -> [(UTCTime, Speed)]-avgSpeeds = slidingAverageSpeed AvgHarmonicMean---- | @slidingAverageSpeed m n@ Average speed using a moving window of up to @n@ seconds--- and an 'AvgMethod' of @m@.-slidingAverageSpeed :: (Lat a, Lon a, Time a) => - AvgMethod a -> NominalDiffTime -> Trail a -> [(UTCTime, Speed)]-slidingAverageSpeed _ _ [] = []-slidingAverageSpeed m 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 :: (Lat a, Lon a, Time a) => Double -> Double -> PointGrouping a-betweenSpeeds low hi ps =- let spds = concatMap maybeToList $ zipWith speed ps (drop 1 ps)- psSpds = [(p,s) | p <- ps, s <- maybeToList (listToMaybe spds) ++ spds]- inRange x = x >= low && x <= hi- chunk [] = []- chunk xs@(x:_) =- let op p = if inRange (snd x) then first Select . span p else first NotSelect . break p- (r,rest) = op (inRange . snd) xs- in r : chunk xs- in map (fmap (map fst)) $ chunk psSpds---- | A "rest point" means the coordinates remain within a given distance--- for at least a particular amount of time.-restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> PointGrouping a-restLocations d s ps =- let consToFirst x [] = [NotSelect [x]]- consToFirst x (a:as) = (fmap (x:) a) : as- go [] [] = []- go [] nonRests = [NotSelect $ reverse nonRests]- go (a:as) nonRests =- case takeWhileEnd ((<=) d . distance a) as of- (Just l, close, far) ->- case (getUTCTime a, getUTCTime l) of- (Just t1, Just t2) ->- let diff = diffUTCTime t2 t1- in if diff >= s then NotSelect (reverse nonRests) : Select (a:close) : go far [] else go as (a:nonRests)- _ -> consToFirst a $ go as nonRests- _ -> consToFirst a $ go as nonRests- in go ps []- --- |chunking points into groups spanning at most the given time--- interval.-spansTime :: (Lat a, Lon a, Time a) => NominalDiffTime -> PointGrouping a-spansTime n ps =- let times = mkTimePair ps- chunk [] = []- chunk 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 :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a-intersectionOf gs ps =- let groupings = map ($ ps) gs- -- chunk :: [[Selected [pnts]]] -> pnts -> [pnts]- chunk _ [] = []- chunk ggs xs = - let minLen = max 1 . minimum . concatMap (take 1) $ map (map selLength) ggs -- FIXME this is all manner of broken- sel = if all isSelected (concatMap (take 1) ggs) then Select else NotSelect- (c,rest) = splitAt minLen xs- in sel c : chunk (filter (not . null) $ map (dropExact minLen) ggs) rest- in chunk groupings ps---- | Union all the groupings-unionOf :: (Lat a, Lon a, Time a) => [PointGrouping a] -> PointGrouping a-unionOf gs ps =- let groupings = map ($ ps) gs- chunk _ [] = []- chunk ggs xs =- let getSegs = concatMap (take 1)- segs = getSegs ggs- len =- if any isSelected segs- then max 1 . maximum . getSegs . map (map selLength) . map (filter isSelected) $ ggs- else max 1 . minimum . getSegs . map (map selLength) $ ggs- sel = if any isSelected segs then Select else NotSelect- (c,rest) = splitAt len xs- in sel c : chunk (filter (not . null) $ map (dropExact len) ggs) rest- in chunk groupings ps- --- | Intersection binary operator-(/\) :: [Selected (Trail a)] -> TransformGrouping a-(/\) xs [] = []-(/\) [] ys = []-(/\) xsL@(Select x:xs) ysL@(Select y:ys) =- let z = if length x < length y then x else y- xs' = selListDrop (length z) xsL- ys' = selListDrop (length z) ysL- in Select z : (xs' /\ ys')-(/\) xs (NotSelect y:ys) = NotSelect y : (selListDrop (length y) xs /\ ys)-(/\) (NotSelect x:xs) ys = NotSelect x : (xs /\ selListDrop (length x) ys)---- | Union binary operator-(\/) :: [Selected (Trail a)] -> TransformGrouping a-(\/) xs [] = xs-(\/) [] ys = ys-(\/) (Select x:xs) (Select y : ys) =- let xLen = length x- yLen = length y- in if xLen < yLen- then (Select y :) (selListDrop (yLen - xLen) xs \/ ys)- else (Select x :) (xs \/ selListDrop (xLen - yLen) ys)-(\/) (Select x:xs) ys = Select x : selListDrop (length x) ys-(\/) xs (Select y:ys) = Select y : selListDrop (length y) xs-(\/) xsL@(NotSelect x:xs) ysL@(NotSelect y:ys) =- let xLen = length x- yLen = length y- in if xLen < yLen- then (NotSelect x:) (xs \/ selListDrop xLen ysL)- else (NotSelect y:) (selListDrop yLen xsL \/ ys)- -selListTake :: Int -> [Selected [a]] -> [Selected [a]]-selListTake 0 _ = []-selListTake n [] = []-selListTake n (x:xs) =- let x' = take n (unSelect x)- in fmap (const x') x : selListTake (n - length x') xs--selListDrop :: Int -> [Selected [a]] -> [Selected [a]]-selListDrop 0 xs = xs-selListDrop n [] = []-selListDrop n (x:xs) =- let x' = drop n (unSelect x)- in fmap (const x') x : selListDrop (n - (selLength x - length x')) xs---- |Inverts the selected/nonselected segments-invertSelection :: TransformGrouping a-invertSelection = map (onSelected NotSelect Select)---- |@firstGrouping f ps@ only the first segment remains 'Select'ed, and only--- if it was already selected by @f@.-firstGrouping :: TransformGrouping a-firstGrouping ps = take 1 ps ++ map (NotSelect . unSelect) (drop 1 ps)---- | Only the last segment, if any, is selected (note: the current--- implementation is inefficient, using 'reverse')-lastGrouping :: TransformGrouping a-lastGrouping ps = let ps' = reverse ps in reverse $ take 1 ps' ++ map (NotSelect . unSelect) (drop 1 ps')---- | chunk the trail into groups of N points-everyNPoints :: Int -> PointGrouping a-everyNPoints n ps- | n <= 0 = [NotSelect ps]- | otherwise = go ps- where- go [] = []- go xs = let (h,t) = splitAt n xs in Select h : go t- --- |For every selected group, refine the selection using the second--- grouping method. This differs from 'IntersectionOf' by restarting--- the second grouping algorithm at the beginning each group selected--- by the first algorithm.-refineGrouping :: PointGrouping a -> TransformGrouping a-refineGrouping b = concatMap (onSelected b (\x -> [NotSelect x]))---- |Remove all points that remain 'NotSelect'ed by the given grouping algorithm.-filterPoints :: PointGrouping a -> Trail a -> Trail a-filterPoints g = concatMap unSelect . filter isSelected . g---- Extract the time from each coordinate. If no time is available then--- the coordinate is dropped!-mkTimePair :: (Time a) => Trail a -> [(a,UTCTime)]-mkTimePair xs =- let timesM = map (\x-> fmap (x,) $ getUTCTime x) xs- in concatMap maybeToList timesM---- |Construct a bezier curve using the provided trail. Construct a--- new trail by sampling the given bezier curve at the given times.--- The current implementation assumes the times of the input--- coordinates are available and all equal (Ex: all points are 5--- seconds apart), the results will be poor if this is not the case!-bezierCurveAt :: (Lat a, Lon a, Time a) => [UTCTime] -> Trail a -> Trail a-bezierCurveAt _ [] = []-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 (setTime . Just . fromUTCTime) newTimes curvePoints--bezierPoint :: (Lat a, Lon a) => [a] -> Double -> a-bezierPoint pnts t = go pnts- where- go [] = error "GPS Package: Can not create a bezier point from an empty list"- go [p] = p- go ps = interpolate (go (init ps)) (go (tail ps)) t---- |Interpolate selected points onto a bezier curve. Note this gets--- exponentially more expensive with the length of the segement being--- transformed - it is not advisable to perform this operation on--- trail segements with more than ten points!-bezierCurve :: (Lat a, Lon a, Time a) => [Selected (Trail a)] -> Trail a-bezierCurve = concatMap (onSelected (bezierCurveAt []) Prelude.id)---- |Filters out any points that go backward in time (thus must not be--- valid if this is a trail)-linearTime :: (Lon a, Lat a, Time a) => [a] -> [a]-linearTime [] = []-linearTime (p:ps) = go (getUTCTime p) ps- where- go _ [] = []- go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps---- |Returns the closest distance between two trails (or Nothing if a--- trail is empty). Inefficient implementation:--- O( (n * m) * log (n * m) )-closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance-closestDistance as bs = listToMaybe $ L.sort [distance a b | a <- as, b <- bs]---- | Find the total distance traveled-totalDistance :: (Lat a, Lon a) => [a] -> Distance-totalDistance as = sum $ zipWith distance as (drop 1 as)--totalTime :: Time 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, Lat c, Lon 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 :: (Lat c, Lon c) => c -> c -> Heading-eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s--data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum)--turn :: (Lat c, Lon c) => c -> c -> c -> Turn-turn a b c =- let h1 = eastZeroHeading a b- h2 = eastZeroHeading b c- d = h2 - h1- in if d >= 0 && d < pi then LeftTurn else RightTurn---- | Find the southmost point-southMost :: (Lat c) => [c] -> Maybe c-southMost [] = Nothing-southMost cs = Just . minimumBy (comparing lat) $ cs------------ COMPOSIT OPERATIONS ------------------ These operations are simply implemented using the previously--- defined functions. They can serve either for concise use for novice--- users or as instructional examples.---------------------------------------------smoothRests :: (Lat a, Lon a, Time a) => Trail a -> Trail a-smoothRests = bezierCurve . refineGrouping (everyNPoints 8) . restLocations 30 60--smoothSome :: (Lat a, Lon a, Time 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 :: (Lat a, Lon a, Time 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)
+ Geo/Computations.hs view
@@ -0,0 +1,14 @@+{-# 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+ , module Geo.Types+ ) where++import Geo.Computations.Basic+import Geo.Computations.Trail+import Geo.Types
+ Geo/Computations/Basic.hs view
@@ -0,0 +1,203 @@+{-# 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+ ) where++import Data.Time+import Data.Maybe+import Geo.Types++-- |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
+ Geo/Computations/Trail.hs view
@@ -0,0 +1,447 @@+{-# 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.Types++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
+ Geo/Types.hs view
@@ -0,0 +1,35 @@+module Geo.Types where+import Data.Text+import Data.Time++data Track = Track+ { trkName :: Maybe Text+ , trkDescription :: Maybe Text+ , segments :: [Segment]+ }+ deriving (Eq, Ord, Show, Read)++-- |A GPX segments is just a bundle of points.+data Segment = Segment { points :: [Point] }+ deriving (Eq, Ord, Show, Read)++type Latitude = Double+type Longitude = Double++-- |Track point is a full-fledged representation of all the data+-- available in most GPS loggers. It is possible you don't want+-- all this data and can just made do with coordinates (via 'Pnt')+-- or a custom derivative.+data Point = Point+ { pntLat :: Latitude+ , pntLon :: Longitude+ , pntEle :: Maybe Double -- ^ In meters+ , pntTime :: Maybe UTCTime+ -- , pntSpeed :: Maybe Double -- ^ Non-standard. Usually in meters/second.+ }+ deriving (Eq, Ord, Show, Read)++pt :: Latitude -> Longitude -> Maybe Double -> Maybe UTCTime -> Point+pt t g e m = Point t g e m++zeroPoint = Point 0 0 Nothing Nothing
LICENSE view
@@ -0,0 +1,30 @@+Copyright (c)2012, Thomas M. DuBuisson++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Thomas M. DuBuisson nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
Setup.hs view
Test/GpsTest.hs view
@@ -1,4 +1,4 @@-import Data.GPS+import Geo.Computations import Data.Time import Data.List import Data.Ord@@ -6,19 +6,9 @@ import Test.QuickCheck import Test.Framework (Test, defaultMain, testGroup) import Test.Framework.Providers.QuickCheck2 (testProperty)-import Text.XML.XSD.DateTime import Control.Applicative import Control.Monad -instance Arbitrary LatitudeType where- arbitrary = liftM (latitudeType . flip mod' 180) arbitrary--instance Arbitrary LongitudeType where- arbitrary = liftM (longitudeType . flip mod' 180) arbitrary--instance Arbitrary DateTime where- arbitrary = liftM fromUTCTime arbitrary- instance Arbitrary UTCTime where arbitrary = UTCTime <$> arbitrary <*> liftM (secondsToDiffTime . abs) (arbitrary :: Gen Integer) @@ -27,53 +17,35 @@ instance Arbitrary NominalDiffTime where arbitrary = liftM fromIntegral (arbitrary :: Gen Int)-instance Arbitrary WptType where++instance Arbitrary Point where arbitrary =- wptType <$> arbitrary -- Lat- <*> arbitrary -- Lon+ pt <$> fmap (`mod'` 90) arbitrary -- Lat+ <*> fmap (`mod'` 180) arbitrary -- Lon <*> arbitrary -- Time <*> arbitrary -- elevation- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return [] -- LinkType- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing- <*> return Nothing -newtype Trl = Trl [WptType]+newtype Trl = Trl [Point] deriving (Show) instance Arbitrary Trl where arbitrary = do- b <- arbitrary :: Gen [Int]- u_ts_d <- mapM (\i -> (,,) <$> arbitrary <*> replicateM i arbitrary <*> arbitrary) b :: Gen [(UTCTime, [WptType],NominalDiffTime)]- let u_ts_d' = sortBy (comparing (\(a,_,_) -> a)) u_ts_d- xs = concat [zipWith (setTime' . fromUTCTime) (iterate (addUTCTime d) u) x | (u,x,d) <- u_ts_d']- return $ Trl xs+ b <- (`mod` 5) `fmap` arbitrary :: Gen Int+ pnts <- mapM (\_ -> arbitrary) [0..abs b]+ return $ Trl (sortBy (comparing pntTime) pnts) -approxEq :: WptType -> WptType -> Bool+approxEq :: Point -> Point -> Bool approxEq a b = distance a b <= 0.2 -- error of 13cm has been observed due to floating point issues when using add vector. -pSaneDistance :: WptType -> WptType -> Bool+pSaneDistance :: Point -> Point -> Bool pSaneDistance a b = distance a b <= circumferenceOfEarth / 2 -pTriangleTheorem :: WptType -> WptType -> WptType -> Bool+pTriangleTheorem :: Point -> Point -> Point -> Bool pTriangleTheorem a b c = distance a b + distance b c >= distance a c -- Traditional flat-surface geometry || distance a b + distance b c + distance c a == 2 * pi * radiusOfEarth -pAddVector_DistanceHeading_ident :: WptType -> WptType -> Bool+pAddVector_DistanceHeading_ident :: Point -> Point -> Bool pAddVector_DistanceHeading_ident a b = let v = (distance a b, heading a b) c = addVector v a@@ -82,7 +54,7 @@ pConvexHull_Has_Extreme_Points :: Trl -> Bool pConvexHull_Has_Extreme_Points (Trl ts) = let ch = convexHull ts- ts' = sortBy (comparing lat) ts+ ts' = sortBy (comparing pntLat) ts northMost = last ts' southMost = head ts' in length ts < 3 || (northMost `elem` ch && southMost `elem` ch)@@ -105,9 +77,11 @@ , testProperty "TriangleTheorem" pTriangleTheorem , testProperty "Vector identity" pAddVector_DistanceHeading_ident ]- , testGroup "Trail Computations"- [ testProperty "Hull has extreme points" pConvexHull_Has_Extreme_Points- , testProperty "HullContainsBezier" pConvexHull_Bezier_Const]+-- These might make some sense in a local scope, but in a global range what+-- is the "left most" point? On a sphere what is a convex hull?+-- , testGroup "Trail Computations"+-- [ testProperty "Hull has extreme points" pConvexHull_Has_Extreme_Points+-- , testProperty "HullContainsBezier" pConvexHull_Bezier_Const] ]
gps.cabal view
@@ -1,27 +1,30 @@-name: gps-version: 0.8.4-license: BSD3-license-file: LICENSE-author: Thomas DuBuisson <thomas.dubuisson@gmail.com>-maintainer: Thomas DuBuisson-description: Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations. NOTE: Version range 0.8.* won't strictly follow PVP - I will be adding additional functions in minor releases 0.8.x.-synopsis: For manipulating GPS coordinates and trails.-category: Data-stability: stable-build-type: Simple-cabal-version: >= 1.8-tested-with: GHC == 6.10.3+name: gps+version: 1.2+license: BSD3+license-file: LICENSE+author: Thomas DuBuisson <thomas.dubuisson@gmail.com>+maintainer: Thomas DuBuisson+description: Useful for manipulating GPS coordinages (in various forms), building paths, and performing basic computations.+synopsis: For manipulating GPS coordinates and trails.+category: Data+stability: stable+build-type: Simple+cabal-version: >= 1.8+tested-with: GHC == 6.10.3 extra-source-files: Library Build-Depends: base >= 3 && < 6,- pretty >= 1.0 , prettyclass >= 1.0,- time >= 1.1, GPX >= 0.5, hxt >= 9.1, xsd >= 0.3,- vector >= 0.7, statistics >= 0.9+ pretty >= 1.0,+ prettyclass >= 1.0,+ time >= 1.4,+ vector >= 0.7,+ statistics >= 0.9,+ text hs-source-dirs:- exposed-modules: Data.GPS- other-modules: Data.GPS.Core Data.GPS.Trail- ghc-options:+ exposed-modules: Geo.Computations, Geo.Types+ other-modules: Geo.Computations.Basic, Geo.Computations.Trail+ ghc-options: -Wall test-suite gps-tests hs-source-dirs: Test@@ -31,9 +34,9 @@ build-depends: base >= 4.0, QuickCheck >= 2.4.0.1,- test-framework >= 0.3.3 && < 0.5,- test-framework-quickcheck2 >= 0.2.9 && < 0.3,- time, GPX, hxt, xsd, vector, statistics, gps+ test-framework >= 0.3.3,+ test-framework-quickcheck2 >= 0.2.9,+ time, vector, statistics, gps, gpx-conduit ghc-options: -Wall source-repository head