{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls #-} -- |Author: Thomas DuBuisson -- Copyright: Thomas DuBuisson -- License: BSD3 -- -- 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 -- -- The current intent of this library is to -- 1) Provide a standard type class interface for coordinates -- 2) Include fleshed out support for relevant libraries, a task -- too often neglected in modern Hackage packages. -- -- Integration includes KML support via the xml package, pretty printing, -- anb binary. module Data.GPS ( -- * Types and Classes Distance , Heading , Speed , Vector , Trail -- * Constants , north , south , east , west , radiusOfEarth -- * Helper Functions , restLocations , closestDistance , getRadianPair , addVector , divideArea , convexHull , module Data.Geo.GPX ) where import Data.Ord (comparing) import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy) import Data.Geo.GPX hiding (none, cmt) import Text.XML.HXT.Arrow import Text.XML.XSD.DateTime(DateTime, toUTCTime) import Data.Time import Data.Maybe (listToMaybe) import Data.Fixed (mod') -- |Distances are expressed in meters type Distance = Double -- |Angles are expressed in radians from North. -- 0 == North -- pi/2 == West -- pi == South -- (3/2)pi == East == - (pi / 2) type Heading = Double -- |Speed is hard coded as meters per second type Speed = Double type Vector = (Distance, Heading) type Trail a = [a] getUTCTime :: (Lat a, Lon a, Time a) => a -> Maybe UTCTime getUTCTime = fmap toUTCTime . time distance :: (Lat a, Lon a) => a -> a -> Distance distance a b = radiusOfEarth * acos( sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2 - lon1) ) where (lat1, lon1) = getRadianPairD a (lat2, lon2) = getRadianPairD b heading :: (Lat a, Lon a) => a -> a -> Heading -- ^ 0 = North, pi/2 = West... heading a b = atan2 (sin (diffLon) * cos (lat2)) (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon)) where (lat1, lon1) = getRadianPairD a (lat2, lon2) = getRadianPairD b diffLon = lon1 - lon2 getVector :: (Lat a, Lon a) => a -> a -> Vector getVector a b = (distance a b, heading a b) -- | Speed in meters per second speed :: (Lat loc, Lon loc, Time loc) => loc -> loc -> Maybe Speed speed a b = case (getUTCTime b, getUTCTime a) of (Just x, Just y) -> Just $ realToFrac (diffUTCTime x y) / (distance a b) _ -> Nothing data TempTrail a = T (Trail a) a -- |radius of the earth in meters radiusOfEarth :: Double radiusOfEarth = 6378700 -- |North is 0 radians north :: Heading north = 0 -- |South, being 180 degrees from North, is pi. south :: Heading south = pi -- |East is 270 degrees from North east :: Heading east = (3 / 2) * pi -- |West is 90 degrees (pi/2) west :: Heading west = pi / 2 toDecimal = (*) (180 / pi) -- |Given a vector and coordinate, computes a new coordinate. -- Within some epsilon it should hold that if -- -- @dest = addVector (dist,heading) start@ -- -- then -- -- @heading == dmsHeading start dest@ -- -- @dist == distance start dest@ addVector :: (Lat c, Lon c) => Vector -> c -> c addVector (d,h) p = setLon (longitudeType lon2) . setLat (latitudeType lat2) $ p where (lat,lon) = getRadianPairD p lat2 = lat + (cos h) * (d / radiusOfEarth) lon2 = lon + acos ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2)) getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double) getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair -- |Provides a lat/lon pair of doubles in radians getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType) getRadianPair p = (t, g) where t = toRadians (lat p) g = toRadians (lon p) toRadians :: Floating f => f -> f toRadians = (*) (pi / 180) -- |Filters out any points that go backward in time (thus must not be valid if this is a trail) linearTime :: (Lat a, Lon a, Time a) => Trail a -> Trail a linearTime [] = [] linearTime (p:ps) = go (getUTCTime p) ps where go _ [] = [] go t (p:ps) = if getUTCTime p < t then go t ps else p : go (getUTCTime p) ps -- |Creates a list of trails all of which are within the given distance of each -- other spanning atleast the given amount of time. -- -- For example @restLocations 50 600@ -- would return lists of all points that are within 50 meters of each other and -- span at least 10 minutes (600 seconds). -- -- Note this gives points within fifty meters of the earliest point - wandering -- in a rest area with a 50 meter radius could result in several rest points -- ([a,b..]) or even none if the distance between individual points exceeds 50m. restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a] restLocations d s xs = go xs where go [] = [] go (a:as) = let (lst, close, far) = takeWhileLast ((<= d) . distance a) as in case lst of Just x -> case (getUTCTime a, getUTCTime x) of (Just a', Just x') -> let d = diffUTCTime a' x' in if d >= s then close : go far else go as _ -> go as Nothing -> go as takeWhileLast :: (a -> Bool) -> [a] -> (Maybe a, [a], [a]) takeWhileLast p [] = (Nothing, [], []) takeWhileLast p (x:xs) | not (p x) = (Nothing, [], x:xs) | otherwise = go x xs where go a [] = (Just a, [a], []) go a (b:bs) | p b = let (c,d,f) = go b bs in (c, a:d, f) | otherwise = (Just a, [a], b:bs) -- |Returns the closest distance between two trails (or Nothing if a trail is empty) -- O( (n * m) * log (n * m) ) closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs] -- |@divideArea vDist hDist nw se@ divides an area into a grid of equally -- spaced coordinates within the box drawn by the northwest point (nw) and -- southeast point (se). Because this uses floating point there might be a -- different number of points in some rows (the last might be too far east based -- on a heading from the se point). divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]] divideArea vDist hDist nw se = let (top,left) = (lat nw, lon nw) (btm,right) = (lat se, lon se) columnOne = takeWhile ( (<= west) . heading se) . iterate (addVector (vDist, south)) $ nw buildRow = takeWhile ((>= north) . heading se) . iterate (addVector (hDist, east)) in map buildRow columnOne -- | Uses Grahams scan to compute the convex hull of the given points. -- This operation requires sorting of the points, so don't try it unless -- you have notably more memory than the list of points will consume. convexHull :: (Eq c, Lat c, Lon c) => [c] -> [c] convexHull xs = let first = southMost xs in case first of Nothing -> [] Just f -> let sorted = sortBy (comparing ((`mod'` (2*pi)). (+ pi/2). heading 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) 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 = heading a b h2 = heading b c in case compare (h2 - h1) 0 of LT -> RightTurn GT -> LeftTurn EQ -> Straight -- | Find the southmost point southMost :: (Lat c) => [c] -> Maybe c southMost [] = Nothing southMost cs = Just . minimumBy (comparing lat) $ cs