module Data.GPS
(
Distance
, Heading
, Speed
, Vector
, Trail
, north
, south
, east
, west
, radiusOfEarth
, heading
, distance
, speed
, addVector
, getRadianPair
, getDMSPair
, divideArea
, totalDistance
, restLocations
, closestDistance
, filterByMaxSpeed
, convexHull
, readGPX
, module Data.Geo.GPX
) where
import Data.Function (on)
import Data.Ord (comparing)
import Data.List (sort, mapAccumL, minimumBy, maximumBy, sortBy)
import Data.Geo.GPX hiding (none, cmt)
import Text.XML.HXT.Arrow
import Text.XML.XSD.DateTime(DateTime, toUTCTime)
import Data.Time
import Data.Maybe (listToMaybe)
import Data.Fixed
import Control.Monad
type Distance = Double
type Heading = Double
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
totalDistance :: (Lat a, Lon a) => [a] -> Distance
totalDistance as = sum $ zipWith distance as (drop 1 as)
heading :: (Lat a, Lon a) => a -> a -> Heading
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 :: (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
filterByMaxSpeed :: (Lat loc, Lon loc, Time loc) => Speed -> Trail loc -> Trail loc
filterByMaxSpeed mx xs =
let ss = zipWith speed xs (drop 1 xs)
fs = filter ((< Just mx) . fst) (zip ss $ drop 1 xs)
in take 1 xs ++ map snd fs
data TempTrail a = T (Trail a) a
radiusOfEarth :: Double
radiusOfEarth = 6378700
north :: Heading
north = 0
south :: Heading
south = pi
east :: Heading
east = (3 / 2) * pi
west :: Heading
west = pi / 2
toDecimal = (*) (180 / pi)
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
getDMSPair :: (Lat c, Lon c) => c -> (LatitudeType, LongitudeType)
getDMSPair c = (lat c, lon c)
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)
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
restLocations :: (Lat a, Lon a, Time a) => Distance -> NominalDiffTime -> Trail a -> [Trail a]
restLocations d s xs = go xs
where
go [] = []
go (a:as) =
case takeWhileLast ((<=) d . distance a) as of
(Just l, close, far) ->
case (getUTCTime a, getUTCTime l) of
(Just t1, Just t2) ->
let diff = diffUTCTime t2 t1
in if diff >= s then (a:close) : go far else go as
_ -> go as
_ -> 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, [], [])
go !a (b:bs)
| p b = let (c,d,f) = go b bs in (c, a:d, f)
| otherwise = (Just a, [], b:bs)
closestDistance :: (Lat a, Lon a) => Trail a -> Trail a -> Maybe Distance
closestDistance as bs = listToMaybe $ sort [distance a b | a <- as, b <- bs]
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
readGPX :: FilePath -> IO (Trail WptType)
readGPX = liftM (concatMap trkpts . concatMap trksegs . concatMap trks) . readGpxFile
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 (eastZeroHeading 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)
eastZeroHeading :: (Lat c, Lon c) => c -> c -> Heading
eastZeroHeading s = (`mod'` (2*pi)) . (+ pi/2) . heading s
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 = eastZeroHeading a b
h2 = eastZeroHeading b c
d = h2 h1
in if d >= 0 && d < pi then LeftTurn else RightTurn
southMost :: (Lat c) => [c] -> Maybe c
southMost [] = Nothing
southMost cs = Just . minimumBy (comparing lat) $ cs