module Geo.Computations.Basic
(
Distance
, Heading
, Speed
, Vector
, Trail
, Circle
, Arc
, Point(..)
, pt
, north
, south
, east
, west
, radiusOfEarth
, circumferenceOfEarth
, heading
, distance
, speed
, getVector
, addVector
, interpolate
, circleIntersectionPoints
, intersectionArcsOf
, maximumDistanceOfArc
, readGPXFile
) where
import Data.Time
import Data.Maybe
import Geo.GPX.Conduit
type Distance = Double
type Heading = Double
type Speed = Double
type Vector = (Distance, Heading)
type Circle a = (a, Distance)
type Arc a = (Circle a, Heading, Heading)
type Trail a = [a]
distance :: Point -> Point -> Distance
distance x y =
let (lat1,lon1) = getRadianPair x
(lat2,lon2) = getRadianPair y
deltaLat = lat2 lat1
deltaLon = lon2 lon1
a = (sin (deltaLat / 2))^(2::Int) + cos lat1 * cos lat2 * (sin (deltaLon / 2))^(2::Int)
c = 2 * atan2 (a**0.5) ((1a)**0.5)
in radiusOfEarth * c
heading :: Point -> Point -> Heading
heading a b =
atan2 (sin (diffLon) * cos (lat2))
(cos(lat1) * sin (lat2) sin(lat1) * cos lat2 * cos (diffLon))
where
(lat1, lon1) = getRadianPair a
(lat2, lon2) = getRadianPair b
diffLon = lon2 lon1
getVector :: Point -> Point -> Vector
getVector a b = (distance a b, heading a b)
addVector :: Vector -> Point -> Point
addVector (d,h) p =
p { pntLon = toDegrees lon2
, pntLat = toDegrees lat2
}
where
(lat,lon) = getRadianPair p
lat2 = asin (sin lat * cos (d / radiusOfEarth) + cos lat
* sin(d/radiusOfEarth) * cos h)
lon2 = lon atan2 (sin h * sin (d / radiusOfEarth) * cos lat)
(cos (d/radiusOfEarth) sin lat * sin lat2)
speed :: Point -> Point -> Maybe Speed
speed a b =
case (pntTime b, pntTime a) of
(Just x, Just y) ->
let timeDiff = realToFrac (diffUTCTime x y)
in if timeDiff == 0 then Nothing else Just $ (distance a b) / timeDiff
_ -> Nothing
radiusOfEarth :: Double
radiusOfEarth = 6378700
circumferenceOfEarth :: Double
circumferenceOfEarth = radiusOfEarth * 2 * pi
north :: Heading
north = 0
south :: Heading
south = pi
east :: Heading
east = (3 / 2) * pi
west :: Heading
west = pi / 2
toDegrees :: Double -> Double
toDegrees = (*) (180 / pi)
getRadianPair :: Point -> (Double,Double)
getRadianPair p = (toRadians (pntLat p), toRadians (pntLon p))
toRadians :: Floating f => f -> f
toRadians = (*) (pi / 180)
interpolate :: Point -> Point -> Double -> Point
interpolate c1 c2 w
| w < 0 = c1
| w > 1 = c2
| otherwise =
let (h,d) = (heading c1 c2, distance c1 c2)
v = (d * w, h)
in addVector v c1
circleIntersectionPoints :: (Point, Distance) -> (Point, Distance) -> Maybe (Point,Point)
circleIntersectionPoints (a,r1) (b,r2)
| pntLat a == pntLat b && pntLon a == pntLon b && r1 == r2 = Nothing
| r1 + r2 < ab = Nothing
| any isNaN (map pntLat pts) || any isNaN (map pntLon pts) = Nothing
| otherwise = Just (p1, p2)
where
ab = distance a b
angABX = acos ( (r1^(2::Int) + ab^(2::Int) r2^(2::Int)) / (2 * r1 * ab) )
ang1 = heading a b + angABX
ang2 = heading a b angABX
p1 = addVector (r1, ang1) a
p2 = addVector (r1, ang2) a
pts = [p1,p2]
intersectionArcsOf :: [Circle Point] -> [Arc Point]
intersectionArcsOf cs =
let isArcWithinCircle circ arc = maximumDistanceOfArc (fst circ) arc <= (snd circ)
isArcWithinAllCircles arc = all ($ arc) (map isArcWithinCircle cs)
getArcs c1 c2 = concatMap (buildArcsFromPoints c1 c2) . maybeToList $ circleIntersectionPoints c1 c2
buildArcsFromPoints c1 c2 (p1,p2) =
let c1h1 = heading (fst c1) p1
c1h2 = heading (fst c1) p2
c2h1 = heading (fst c2) p1
c2h2 = heading (fst c2) p2
in [(c1,c1h1,c1h2), (c1,c1h2, c1h1), (c2,c2h1,c2h2), (c2,c2h2,c2h1)]
in filter isArcWithinAllCircles . concatMap (uncurry getArcs) . choose2 $ cs
maximumDistanceOfArc :: Point -> Arc Point -> Distance
maximumDistanceOfArc pnt ((c,r), h1, h2) =
let pcHeading = heading pnt c
in if ((pcHeading < h1 || pcHeading > h2) && h1 < h2) || ((pcHeading > h2 && pcHeading < h1) && h1 > h2)
then max (distance pnt (addVector (r,h1) c)) (distance pnt (addVector (r,h2) c))
else distance pnt c + r
choose2 :: [a] -> [(a,a)]
choose2 [] = []
choose2 (x:xs) = map (x,) xs ++ choose2 xs