gps-0.2.4: Data/GPS.hs

{-# 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)