diff --git a/Data/GPS.hs b/Data/GPS.hs --- a/Data/GPS.hs +++ b/Data/GPS.hs @@ -17,13 +17,7 @@ module Data.GPS ( -- * Types and Classes - Coordinate(..) - , Location(..) - , DMSCoordinate(..) - , DMS(..) - , Latitude - , Longitude - , Distance + Distance , Heading , Speed , Vector @@ -35,30 +29,28 @@ , west , radiusOfEarth -- * Helper Functions - , smoothTrails - , longestSmoothTrail , restLocations , closestDistance - , normalizeDMS - , dmsToDegreePair - , dmsToRadianPair - , degreePairToDMS + , getRadianPair , addVector , divideArea + , convexHull ) where import Data.Ord (comparing) -import Data.List (sort, mapAccumL, maximumBy) -import Data.Binary -import Data.Binary.Put -import Data.Binary.Get +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.PrettyPrint.HughesPJClass -import Text.PrettyPrint -import Numeric (showFFloat) +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 @@ -76,66 +68,34 @@ type Vector = (Distance, Heading) type Trail a = [a] --- |Empty declaration for phantom type -data Latitude - --- |Empty declaration for phantom type -data Longitude - --- |DMSCoordinate is the typical degree minute second --- for latitude/longitude used by most civilian GPS devices. -data DMSCoordinate = DMSCoord - { latitude :: DMS Latitude - , longitude :: DMS Longitude - } deriving (Eq, Ord, Show, Read) - --- |DMS is the degrees, minutes, seconds used for both latitude and longitude -data DMS a = DMS - { degrees :: Double - , minutes :: Double - , seconds :: Double } - deriving (Eq, Ord, Show, Read) - --- |A coordinate is a place on earths surface. While twoDMS and fromDMS are the only --- required functions, coordinate systems may provide more accurate calculations --- of heading, distance, and vectors without the intermediate translation. -class Coordinate a where - toDMS :: a -> DMSCoordinate - fromDMS :: DMSCoordinate -> a - distance :: a -> a -> Distance - distance a b = - radiusOfEarth * acos( sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2 - lon1) ) - where - (lat1, lon1) = dmsToRadianPair (toDMS a) - (lat2, lon2) = dmsToRadianPair (toDMS b) - heading :: a -> a -> Heading -- 0 = North, pi/2 = West... - heading a b = dmsHeading (toDMS a) (toDMS b) - getVector :: a -> a -> Vector - getVector a b = (distance a b, heading a b) +getUTCTime :: (Lat a, Lon a, Time a) => a -> Maybe UTCTime +getUTCTime = fmap toUTCTime . time --- |A location is a coordinate at a specific time -class (Coordinate coord) => Location loc coord time | loc -> coord, loc -> time where - getCoordinate :: loc -> coord - getTime :: loc -> time - getUTC :: loc -> UTCTime - speed :: loc -> loc -> Speed -- meters per second - speed a b = distMeter / timeSec - where - timeSec = realToFrac $ diffUTCTime (getUTC b) (getUTC a) - distMeter = distance (getCoordinate a) (getCoordinate b) +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 -instance Coordinate DMSCoordinate where - toDMS = id - fromDMS = id +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 -instance (Coordinate c) => Location (c, UTCTime) c UTCTime where - getCoordinate = fst - getTime = snd - getUTC = snd +getVector :: (Lat a, Lon a) => a -> a -> Vector +getVector a b = (distance a b, heading a b) -instance (Coordinate c) => Coordinate (c,UTCTime) where - toDMS = toDMS . getCoordinate - fromDMS d = (fromDMS d , UTCTime (toEnum 0) 0) +-- | 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 @@ -161,17 +121,8 @@ toDecimal = (*) (180 / pi) --- |provides a heading in radians from North -dmsHeading :: DMSCoordinate -> DMSCoordinate -> Heading -dmsHeading a b = - atan2 (sin (diffLon) * cos (lat2)) - (cos(lat1) * sin (lat2) - sin(lat1) * cos lat2 * cos (diffLon)) - where - (lat1, lon1) = dmsToRadianPair a - (lat2, lon2) = dmsToRadianPair b - diffLon = lon1 - lon2 --- |Given a vector and coordinate, computes a new DMS coordinate. +-- |Given a vector and coordinate, computes a new coordinate. -- Within some epsilon it should hold that if -- -- @dest = addVector (dist,heading) start@ @@ -181,82 +132,33 @@ -- @heading == dmsHeading start dest@ -- -- @dist == distance start dest@ -addVector :: Vector -> DMSCoordinate -> DMSCoordinate -addVector (d,h) dms = degreePairToDMS (toDecimal lat2, toDecimal lon2) +addVector :: (Lat c, Lon c) => Vector -> c -> c +addVector (d,h) p = setLon (longitudeType lon2) . setLat (latitudeType lat2) $ p where - (lat,lon) = dmsToRadianPair dms + (lat,lon) = getRadianPairD p lat2 = lat + (cos h) * (d / radiusOfEarth) lon2 = lon + acos ( (cos (d/radiusOfEarth) - sin lat * sin lat2) / (cos lat * cos lat2)) --- |Provides a lat/lon pair of doubles in radians -dmsToRadianPair :: DMSCoordinate -> (Double, Double) -dmsToRadianPair (DMSCoord (DMS latD latM latS) (DMS lonD lonM lonS)) = (lat, lon) - where - lon = toRadians lonD lonM lonS - lat = toRadians latD latM latS - toRadians d m s = (d + m / 60 + s / 3600) * (pi / 180) - --- |Provides a lat/lon pair of doubles in degrees -dmsToDegreePair :: DMSCoordinate -> (Double,Double) -dmsToDegreePair (DMSCoord (DMS latD latM latS) (DMS lonD lonM lonS)) = (lat, lon) - where - lon = toDegrees lonD lonM lonS - lat = toDegrees latD latM latS - toDegrees d m s = (d + m / 60 + s / 3600) - --- Given a lat/long pair of doubles in degrees, produces a DMSCoordinate -degreePairToDMS :: (Double,Double) -> DMSCoordinate -degreePairToDMS (lat,lon) = DMSCoord (DMS lat 0 0) (DMS lon 0 0) - --- |Typically useful for printing, normalizes degrees, minutes, seconds --- into just degrees and decimal minutes: --- --- @DMSCoord (DMS 45 36.938455 0) (DMS ...)@ -normalizeDMS :: DMSCoordinate -> DMSCoordinate -normalizeDMS (DMSCoord lat lon) = - DMSCoord lat' lon' - where - lat' = normalize lat - lon' = normalize lon - normalize (DMS d m s) = - let all = d + (m / 60) + (s / 3600) - d' = fromIntegral (floor all) - m' = (all-d') * 60 - in DMS d' m' 0 - -fromTemp :: TempTrail a -> Trail a -fromTemp (T ps p) = reverse (p:ps) +getRadianPairD :: (Lat c, Lon c) => c -> (Double,Double) +getRadianPairD = (\(a,b) -> (realToFrac a, realToFrac b)) . getRadianPair --- |@smoothTrails speed trail@ should separate points that would mandate a speed --- in excess of the given rate into separate Trails. No point are deleted; --- if there is only one 'correct' trail expected then the largest resulting --- trail will likely be the one your looking for. -smoothTrails :: Location a b c => Speed -> Trail a -> [Trail a] -smoothTrails s ps = - let (trails, _) = mapAccumL go [] (linearTime ps) - in map fromTemp trails +-- |Provides a lat/lon pair of doubles in radians +getRadianPair :: (Lat p, Lon p) => p -> (LatitudeType, LongitudeType) +getRadianPair p = (t, g) where - go acc p = - let (i,trails) = mapAccumL (add p) 0 acc - in if i == 0 - then (T [] p : trails, ()) - else (trails, ()) - add p i (T ts l) = - if (speed l p) <= s - then (i + 1, T (l:ts) p) - else (i, T ts l) + t = toRadians (lat p) + g = toRadians (lon p) --- |Uses 'smoothTrails' but returns only the longest resulting trail -longestSmoothTrail :: Location a b c => Speed -> Trail a -> Trail a -longestSmoothTrail metersPerSec trail = maximumBy (comparing length) ([] : smoothTrails metersPerSec trail) +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 :: Location a b c => Trail a -> Trail a +linearTime :: (Lat a, Lon a, Time a) => Trail a -> Trail a linearTime [] = [] -linearTime (p:ps) = go (getUTC p) ps +linearTime (p:ps) = go (getUTCTime p) ps where go _ [] = [] - go t (p:ps) = if getUTC p < t then go t ps else p : go (getUTC p) ps + 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. @@ -268,19 +170,34 @@ -- 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 :: Location a b c => Distance -> NominalDiffTime -> Trail a -> [Trail a] -restLocations d s = filter timeSpan . span2 +restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a] +restLocations d s xs = go xs where - span2 :: Location a b c => [a] -> [[a]] - span2 [] = [] - span2 (a:as) = let (x,y) = span ((>) d . distance (getCoordinate a) . getCoordinate) as in (a : x) : (span2 y) - timeSpan :: Location a b c => Trail a -> Bool - timeSpan [] = False - timeSpan t = s <= diffUTCTime (getUTC (last t)) (getUTC (head t)) + 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 :: Coordinate a => Trail a -> Trail a -> Maybe Distance +closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs] @@ -289,47 +206,49 @@ -- 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 :: Distance -> Distance -> DMSCoordinate -> DMSCoordinate -> [[DMSCoordinate]] +divideArea :: (Lat c, Lon c) => Distance -> Distance -> c -> c -> [[c]] divideArea vDist hDist nw se = - let (top,left) = dmsToDegreePair nw - (btm,right) = dmsToDegreePair 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 -instance Pretty DMSCoordinate where - pPrint c = - let (DMSCoord (DMS d' m' _) (DMS d m _)) = normalizeDMS c - showF = showFFloat (Just 6) - showX = shows . floor - in text $ showX d' (' ' : showF m' ('\t' : showX d (' ' : showF m ""))) - -instance Binary DMSCoordinate where - put (DMSCoord (DMS tD tM tS) (DMS gD gM gS)) = do - put tD >> put tM >> put tS >> put gD >> put gM >> put gS - get = do - tD <- get - tM <- get - tS <- get - gD <- get - gM <- get - gS <- get - return (DMSCoord (DMS tD tM tS) (DMS gD gM gS)) +-- | 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) -putLocation :: Binary c => (c, UTCTime) -> Put -putLocation (c,t) = put c >> put t +data Turn = LeftTurn | RightTurn | Straight deriving (Eq, Ord, Show, Read, Enum) -getLocation :: Binary c => Get (c, UTCTime) -getLocation = do - c <- get - t <- get - return (c,t) +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 -instance Binary UTCTime where - get = do - day <- get - diff <- get - return (UTCTime (toEnum day) (fromRational diff)) - put (UTCTime day diff) = do - put (fromEnum day) - put (toRational diff) +-- | Find the southmost point +southMost :: (Lat c) => [c] -> Maybe c +southMost [] = Nothing +southMost cs = Just . minimumBy (comparing lat) $ cs diff --git a/Data/GPS/KML.hs b/Data/GPS/KML.hs deleted file mode 100644 --- a/Data/GPS/KML.hs +++ /dev/null @@ -1,83 +0,0 @@ -{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, EmptyDataDecls #-} --- |Author: Thomas DuBuisson --- Copyright: Thomas DuBuisson --- License: BSD3 -module Data.GPS.KML - ( -- * KML Operations (Open format used by GoogleEarth, among others) - KML - , trailToKML - , pointsToKML - , kmlToString - ) where - -import Text.XML.Light -import Data.GPS - --- |The KML type and operations might be moved into a Data.KML module in the future -type KML = Element - --- |Converts the KML elements to a string and prepends the proper XML header, --- thus making it the correct format for saving as a file and opening it --- with other programs such as GoogleEarth. -kmlToString :: KML -> String -kmlToString = (++) xml_header . ppElement - -kmlTop = Element (QName "kml" Nothing Nothing) [Attr (basicName "xmlns") "http://www.opengis.net/kml/2.2"] [] Nothing - -addElem :: Element -> Element -> Element -addElem top child = top { elContent = Elem child : elContent top } - -basicName :: String -> QName -basicName n = QName n Nothing Nothing - -basicElem :: String -> Element -basicElem n = Element (basicName n) [] [] Nothing - -filledElem :: String -> String -> Element -filledElem name val = Element (basicName name) [] [Text (CData CDataText val Nothing)] Nothing - --- |converts a given set of coordinates to a trail in KML format. --- Useful for saving as a file. -trailToKML :: Coordinate a => String -> Trail a -> KML -trailToKML name ps = addElem kmlTop doc - where - doc = foldl addElem (basicElem "Document") (nameElem : trailElem : []) - nameElem = filledElem "name" name - trailName = filledElem "name" (name ++ " trail") - trailElem = foldl addElem (basicElem "Placemark") [trailName, pointsElem] - pointsElem = addElem (basicElem "LineString") (filledElem "coordinates" points) - points :: String - points = concatMap packCoordinate ps - packCoordinate :: Coordinate a => a -> String - packCoordinate x = - let (lat, lon) = getDegreeLatLon x - in concat [lon, ",", lat, ",0 "] - --- |converts a given set of coordinates to points in KML format. --- Useful for saving as a file. -pointsToKML :: Coordinate a => String -> [a] -> [String] -> KML -pointsToKML name ps pointNames = addElem kmlTop doc - where - doc = foldl addElem (basicElem "Document") (nameElem : placemarkElems) - nameElem = filledElem "name" name - placemarkElems = map buildPoint (zip ps pointNames) - buildPoint (p,n) = - let (t,g) = getDegreeLatLon p - in foldl addElem (basicElem "Placemark") - [ filledElem "name" n - , foldl addElem (basicElem "LookAt") - [ filledElem "longitude" g - , filledElem "latitude" t - , filledElem "altitude" "0" - , filledElem "range" "500" - ] - , addElem (basicElem "Point") (filledElem "coordinates" (g ++ "," ++ t)) - ] - --- |Gets a pair of lat lon bytestrings (in degrees) -getDegreeLatLon :: Coordinate a => a -> (String,String) -getDegreeLatLon x = - let (lat,lon) = dmsToRadianPair . toDMS $ x - f = show . (*) (180 / pi) - in (f lat, f lon) - diff --git a/gps.cabal b/gps.cabal --- a/gps.cabal +++ b/gps.cabal @@ -1,5 +1,5 @@ name: gps -version: 0.4.0 +version: 0.5.0 license: BSD3 license-file: LICENSE author: Thomas DuBuisson @@ -9,7 +9,7 @@ category: Data stability: stable build-type: Simple -cabal-version: >= 1.2 +cabal-version: >= 1.6 tested-with: GHC == 6.10.3 extra-source-files: @@ -17,9 +17,9 @@ Description: Choose the split-up base package. Library - Build-Depends: base >= 3 && < 5, bytestring >= 0.9 && < 1.0, binary >= 0.4.0 && < 0.6.0, - pretty >= 1.0 && < 1.1 , prettyclass >= 1.0 && < 1.1, xml >= 1.3 && < 1.4, - time >= 1.1.2 && < 1.2 + Build-Depends: base >= 4 && < 5, bytestring >= 0.9 && < 1.0, + pretty >= 1.0 && < 1.1 , prettyclass >= 1.0 && < 1.1, + time >= 1.1 && < 1.2, GPX == 0.4.*, xsd == 0.3.* hs-source-dirs: - exposed-modules: Data.GPS, Data.GPS.KML + exposed-modules: Data.GPS ghc-options: