gps-0.4.0: Data/GPS.hs

{-# 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
	  Coordinate(..)
	, Location(..)
	, DMSCoordinate(..)
	, DMS(..)
	, Latitude
	, Longitude
	, Distance
	, Heading
	, Speed
	, Vector
	, Trail
	-- * Constants
	, north
	, south
	, east
	, west
	, radiusOfEarth
	-- * Helper Functions
	, smoothTrails
	, longestSmoothTrail
	, restLocations
	, closestDistance
	, normalizeDMS
	, dmsToDegreePair
	, dmsToRadianPair
	, degreePairToDMS
	, addVector
	, divideArea
	) 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)

-- |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]

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

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

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

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

-- |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.
-- Within some epsilon it should hold that if
--
-- 	@dest = addVector (dist,heading) start@
--
-- then
--
-- 	@heading == dmsHeading start dest@
-- 	
-- 	@dist    == distance start dest@
addVector :: Vector -> DMSCoordinate -> DMSCoordinate
addVector (d,h) dms = degreePairToDMS (toDecimal  lat2, toDecimal lon2)
  where
	(lat,lon) = dmsToRadianPair dms
	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)

-- |@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)

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

-- |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 [] = []
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

-- |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 :: 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]


-- |@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 :: Distance -> Distance -> DMSCoordinate -> DMSCoordinate -> [[DMSCoordinate]]
divideArea vDist hDist nw se =
	let (top,left)  = dmsToDegreePair nw
	    (btm,right) = dmsToDegreePair 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))

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)