diff --git a/Data/GPS.hs b/Data/GPS.hs new file mode 100644 --- /dev/null +++ b/Data/GPS.hs @@ -0,0 +1,334 @@ +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances #-} +-- |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 + Coordinate(..) + , Location(..) + , DMSCoordinate(..) + , DMS(..) + , Heading + , Speed + , Vector + , Trail + -- * Helper Functions + , smoothTrails + , longestSmoothTrail + , restLocations + , closestDistance + , normalizeDMS + -- * KML Operations (Open format used by GoogleEarth, among others) + , KML + , trailToKML + , pointsToKML + , kmlToString + ) where + +import Data.Ord (comparing) +import Data.List (sort, mapAccumL, maximumBy) +import Data.Binary +import Data.Binary.Put +import Data.Binary.Get + +import Text.PrettyPrint.HughesPJClass +import Text.PrettyPrint +import Numeric (showFFloat) + +import Data.Time +import Data.Maybe (listToMaybe) + +import Text.XML.Light + +-- |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) = dmsToLatLon . toDMS $ x + f = show . toDecimal + in (f lat, f lon) + +-- |radius of the earth in meters +radiusOfEarth :: Double +radiusOfEarth = 6378700 + +-- |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 = + 3963000 * acos( sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2 - lon1) ) + where + (lat1, lon1) = dmsToLatLon (toDMS a) + (lat2, lon2) = dmsToLatLon (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) + +-- |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 (getDMS a) (getDMS b) + +toDecimal = (*) (180 / pi) + +-- |Extract a lat/lon from a location. +getDMS :: (Location loc x y) => loc -> DMSCoordinate +getDMS = toDMS . getCoordinate + +-- |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) = dmsToLatLon a + (lat2, lon2) = dmsToLatLon b + diffLon = lon1 - lon2 + +-- provides a lat/lon pair in radians +dmsToLatLon :: DMSCoordinate -> (Double, Double) +dmsToLatLon (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) + +-- |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 + +-- |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] + +-- |DMSCoordinate is the typical degree minute second +-- for latitude/longitude used by most civilian GPS devices. +data DMSCoordinate = DMSCoord + { latitude :: DMS + , longitude :: DMS + } deriving (Eq, Ord, Show, Read) + +-- |DMS is the degree, minute, seconds used for both latitude and longitude +data DMS = DMS { degrees :: Double + , minutes :: Double + , seconds :: Double } + deriving (Eq, Ord, Show, Read) + +instance Coordinate DMSCoordinate where + toDMS = id + fromDMS = id + +instance (Coordinate c) => Location (c, UTCTime) c UTCTime where + getCoordinate = fst + getTime = snd + getUTC = snd + +instance (Coordinate c) => Coordinate (c,UTCTime) where + toDMS = toDMS . getCoordinate + fromDMS d = (fromDMS d , UTCTime (toEnum 0) 0) + +data TempTrail a = T (Trail a) a + +fromTemp :: TempTrail a -> Trail a +fromTemp (T ps p) = reverse (p:ps) + +-- |'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 + 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) + +-- |longestSmoothTrail is ust 'smoothTrails' returning only the longest smooth trail. +longestSmoothTrail :: Location a b c => Speed -> Trail a -> Trail a +longestSmoothTrail metersPerSec trail = maximumBy (comparing length) ([] : smoothTrails metersPerSec trail) + +linearTime :: Location a b c => Trail a -> Trail a +linearTime [] = [] +linearTime (p:ps) = go (getUTC p) ps + where + go _ [] = [] + go t (p:ps) = if getUTC p < t then go t ps else p : go (getUTC p) ps + +-- |'restLocation dist timeDiff' will create a list of trails, all of which are +-- within '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 :: Location a b c => Distance -> NominalDiffTime -> Trail a -> [Trail a] +restLocations d s = filter timeSpan . span2 + 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)) + +-- |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 as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs] + +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)) + +putLocation :: Binary c => (c, UTCTime) -> Put +putLocation (c,t) = put c >> put t + +getLocation :: Binary c => Get (c, UTCTime) +getLocation = do + c <- get + t <- get + return (c,t) + +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) diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright (c) Thomas 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: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + +2. 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. + +3. Neither the name of the author nor the names of his contributors + may be used to endorse or promote products derived from this software + without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE 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 AUTHORS 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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/gps.cabal b/gps.cabal new file mode 100644 --- /dev/null +++ b/gps.cabal @@ -0,0 +1,29 @@ +name: gps +version: 0.2.4 +license: BSD3 +license-file: LICENSE +author: Thomas DuBuisson +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.2 +tested-with: GHC == 6.10.3 +extra-source-files: + +Flag small_base + Description: Choose the split-up base package. + +Library + if flag(small_base) + 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 + else + 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 + hs-source-dirs: + exposed-modules: Data.GPS + ghc-options: