{-# 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 ) where import Data.Ord (comparing) import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy) import Data.Geo.GPX.PtType import Data.Geo.GPX.LongitudeType import Data.Geo.GPX.LatitudeType import Data.Geo.GPX.Accessor.Time import Data.Geo.GPX.Accessor.Lon import Data.Geo.GPX.Accessor.Lat 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