-- | -- Module: Data.Geo.Jord.GreatCircle -- Copyright: (c) 2020 Cedric Liegeois -- License: BSD3 -- Maintainer: Cedric Liegeois -- Stability: experimental -- Portability: portable -- -- Geographical Position calculations on great circles, i.e. using a __sphere__ to represent -- the celestial body that positions refer to. -- -- In order to use this module you should start with the following imports: -- -- @ -- import Data.Geo.Jord.GreatCircle -- import Data.Geo.Jord.Position -- @ -- -- If you wish to use both this module and the "Data.Geo.Jord.Geodesic" module you must qualify both imports. -- -- All functions are implemented using the vector-based approached described in -- -- module Data.Geo.Jord.GreatCircle ( -- * The 'GreatCircle' type GreatCircle , greatCircleThrough , greatCircleHeadingOn -- * The 'MinorArc' type , MinorArc , minorArc -- * Calculations , alongTrackDistance , alongTrackDistance' , angularDistance , crossTrackDistance , crossTrackDistance' , destination , finalBearing , initialBearing , interpolate , intersection , intersections , isBetween , isInsideSurface , mean , surfaceDistance ) where import Data.Fixed (Nano) import Data.List (subsequences) import Data.Geo.Jord.Internal import Data.Geo.Jord.Position -- | A circle on the __surface__ of a __sphere__ which lies in a plane -- passing through the sphere centre. Every two distinct and non-antipodal points -- define a unique Great Circle. -- -- It is internally represented as its normal vector - i.e. the normal vector -- to the plane containing the great circle. -- data GreatCircle a = GreatCircle !Vector3d !a String deriving (Eq) instance (Model a) => Show (GreatCircle a) where show (GreatCircle _ _ s) = s -- | @greatCircleThrough p1 p2@ returns the 'GreatCircle' passing by both positions @p1@ and @p2@. -- If positions are antipodal, any great circle passing through those positions will be returned. -- Returns 'Nothing' if given positions are equal. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p1 = latLongHeightPos 45.0 (-143.5) (metres 1500) S84 -- >>> let p2 = latLongHeightPos 46.0 14.5 (metres 3000) S84 -- >>> greatCircleThrough p1 p2 -- heights are ignored, great circle is always at surface. -- Just Great Circle { through 45°0'0.000"N,143°30'0.000"W 1500.0m (S84) & 46°0'0.000"N,14°30'0.000"E 3000.0m (S84) } -- greatCircleThrough :: (Spherical a) => Position a -> Position a -> Maybe (GreatCircle a) greatCircleThrough p1 p2 | llEq p1 p2 = Nothing | otherwise = Just (GreatCircle (normal' p1 p2) (model p1) dscr) where dscr = "Great Circle { through " ++ show p1 ++ " & " ++ show p2 ++ " }" -- | @greatCircleHeadingOn p b@ returns the 'GreatCircle' passing by position @p@ and -- heading on bearing @b@. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p = latLongPos 45.0 (-143.5) S84 -- >>> let b = decimalDegrees 33.0 -- >>> greatCircleHeadingOn p b -- Great Circle { by 45°0'0.000"N,143°30'0.000"W 0.0m (S84) & heading on 33°0'0.000" } -- greatCircleHeadingOn :: (Spherical a) => Position a -> Angle -> GreatCircle a greatCircleHeadingOn p b = GreatCircle (vsub n' e') (model p) dscr where v = nvec p e = vcross nvNorthPole v -- easting n = vcross v e -- northing e' = vscale e (cos' b / vnorm e) n' = vscale n (sin' b / vnorm n) dscr = "Great Circle { by " ++ show p ++ " & heading on " ++ show b ++ " }" -- | Oriented minor arc of a great circle between two positions: shortest path between -- positions on a great circle. data MinorArc a = MinorArc !Vector3d (Position a) (Position a) deriving (Eq) instance (Model a) => Show (MinorArc a) where show (MinorArc _ s e) = "Minor Arc { from: " ++ show s ++ ", to: " ++ show e ++ " }" -- | @minorArc p1 p2@ returns the 'MinorArc' from @p1@ to @p2@. -- Returns 'Nothing' if given positions are equal. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p1 = latLongHeightPos 45.0 (-143.5) (metres 1500) S84 -- >>> let p2 = latLongHeightPos 46.0 14.5 (metres 3000) S84 -- Just Minor Arc { from: 45°0'0.000"N,143°30'0.000"W 1500.0m (S84), to: 46°0'0.000"N,14°30'0.000"E 3000.0m (S84) } -- minorArc :: (Spherical a) => Position a -> Position a -> Maybe (MinorArc a) minorArc p1 p2 | llEq p1 p2 = Nothing | otherwise = Just (MinorArc (normal' p1 p2) p1 p2) -- | @alongTrackDistance p a@ computes how far Position @p@ is along a path described -- by the minor arc @a@: if a perpendicular is drawn from @p@ to the path, the -- along-track distance is the signed distance from the start point to where the -- perpendicular crosses the path. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p = s84Pos 53.2611 (-0.7972) zero -- >>> let g = minorArcBetween (s84Pos 53.3206 (-1.7297) zero) (s84Pos 53.1887 0.1334 zero) -- >>> fmap (alongTrackDistance p) a -- Right 62.3315757km -- alongTrackDistance :: (Spherical a) => Position a -> MinorArc a -> Length alongTrackDistance p (MinorArc n s _) = alongTrackDistance'' p s n -- | @alongTrackDistance' p s b@ computes how far Position @p@ is along a path starting -- at @s@ and heading on bearing @b@: if a perpendicular is drawn from @p@ to the path, the -- along-track distance is the signed distance from the start point to where the -- perpendicular crosses the path. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p = s84Pos 53.2611 (-0.7972) zero -- >>> let s = s84Pos 53.3206 (-1.7297) zero -- >>> let b = decimalDegrees 96.0017325 -- >>> alongTrackDistance' p s b -- 62.3315757km -- alongTrackDistance' :: (Spherical a) => Position a -> Position a -> Angle -> Length alongTrackDistance' p s b = alongTrackDistance'' p s n where (GreatCircle n _ _) = greatCircleHeadingOn s b -- | @angularDistance p1 p2 n@ computes the angle between the horizontal Points @p1@ and @p2@. -- If @n@ is 'Nothing', the angle is always in [0..180], otherwise it is in [-180, +180], -- signed + if @p1@ is clockwise looking along @n@, - in opposite direction. angularDistance :: (Spherical a) => Position a -> Position a -> Maybe (Position a) -> Angle angularDistance p1 p2 n = signedAngle v1 v2 vn where v1 = nvec p1 v2 = nvec p2 vn = fmap nvec n -- | @crossTrackDistance p gc@ computes the signed distance from horizontal Position @p@ to great circle @gc@. -- Returns a negative 'Length' if Position if left of great circle, -- positive 'Length' if Position if right of great circle; the orientation of the -- great circle is therefore important: -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let gc1 = greatCircleThrough (s84Pos 51 0 zero) (s84Pos 52 1 zero) -- >>> fmap (crossTrackDistance p) gc1 -- Right -176.7568725km -- >>> -- >>> let gc2 = greatCircleThrough (s84Pos 52 1 zero) (s84Pos 51 0 zero) -- >>> fmap (crossTrackDistance p) gc2 -- Right 176.7568725km -- >>> -- >>> let p = s84Pos 53.2611 (-0.7972) zero -- >>> let gc = greatCircleHeadingOn (s84Pos 53.3206 (-1.7297) zero) (decimalDegrees 96.0) -- >>> crossTrackDistance p gc -- -305.6629 metres -- crossTrackDistance :: (Spherical a) => Position a -> GreatCircle a -> Length crossTrackDistance p (GreatCircle n _ _) = arcLength (sub a (decimalDegrees 90)) (radius p) where a = radians (angleRadians n (nvec p)) -- | @crossTrackDistance' p s b@ computes the signed distance from horizontal Position @p@ to the -- great circle passing by @s@ and heading on bearing @b@. -- -- This is equivalent to: -- -- @ -- 'crossTrackDistance' p ('greatCircleHeadingOn' s b) -- @ -- crossTrackDistance' :: (Spherical a) => Position a -> Position a -> Angle -> Length crossTrackDistance' p s b = crossTrackDistance p (greatCircleHeadingOn s b) -- | @destination p b d@ computes the position along the great circle, reached from -- position @p@ having travelled the __surface__ distance @d@ on the initial bearing (compass angle) @b@ -- at __constant__ height. -- Note that the bearing will normally vary before destination is reached. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> destination (s84Pos 54 154 (metres 15000)) (decimalDegrees 33) (kilometres 1000) -- 61°10'44.188"N,164°10'19.254"E 15.0km (S84) -- destination :: (Spherical a) => Position a -> Angle -> Length -> Position a destination p b d | d == zero = p | otherwise = nvh nvd (height p) (model p) where nv = nvec p ed = vunit (vcross nvNorthPole nv) -- east direction vector at v nd = vcross nv ed -- north direction vector at v r = radius p ta = central d r -- central angle de = vadd (vscale nd (cos' b)) (vscale ed (sin' b)) -- vunit vector in the direction of the azimuth nvd = vadd (vscale nv (cos' ta)) (vscale de (sin' ta)) -- | @surfaceDistance p1 p2@ computes the surface distance on the great circle between the -- positions @p1@ and @p2@. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> surfaceDistance (northPole S84) (southPole S84) -- 20015.114352233km -- >>> -- >>> surfaceDistance (northPole S84) (northPole S84) -- 0.0m -- surfaceDistance :: (Spherical a) => Position a -> Position a -> Length surfaceDistance p1 p2 = arcLength a (radius p1) where a = radians (angleRadians (nvec p1) (nvec p2)) -- | @finalBearing p1 p2@ computes the final bearing arriving at @p2@ from @p1@ in compass angle. -- Compass angles are clockwise angles from true north: 0° = north, 90° = east, 180° = south, 270° = west. -- The final bearing will differ from the initial bearing by varying degrees according to distance and latitude. -- Returns 'Nothing' if both positions are equals. -- -- ==== __Examples__ -- -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p1 = s84Pos 0 1 (metres 12000) -- >>> let p2 = s84Pos 0 0 (metres 5000) -- >>> finalBearing p1 p2 -- Just 270°0'0.000" -- >>> -- >>> finalBearing p1 p1 -- Nothing -- finalBearing :: (Spherical a) => Position a -> Position a -> Maybe Angle finalBearing p1 p2 | llEq p1 p2 = Nothing | otherwise = Just (normalise (initialBearing' p2 p1) (decimalDegrees 180)) -- | @initialBearing p1 p2@ computes the initial bearing from @p1@ to @p2@ in compass angle. -- Compass angles are clockwise angles from true north: 0° = north, 90° = east, 180° = south, 270° = west. -- Returns 'Nothing' if both positions are equals. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p1 = s84Pos 58.643889 (-5.714722) (metres 12000) -- >>> let p2 = s84Pos 50.066389 (-5.714722) (metres 12000) -- >>> initialBearing p1 p2 -- Just 180°0'0.000" -- >>> -- >>> initialBearing p1 p1 -- Nothing -- initialBearing :: (Spherical a) => Position a -> Position a -> Maybe Angle initialBearing p1 p2 | llEq p1 p2 = Nothing | otherwise = Just (initialBearing' p1 p2) -- | @interpolate p0 p1 f# computes the position at fraction @f@ between the @p0@ and @p1@. -- -- Special conditions: -- -- @ -- interpolate p0 p1 0.0 = p0 -- interpolate p0 p1 1.0 = p1 -- @ -- -- 'error's if @f < 0 || f > 1@ -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let p1 = s84Pos 53.479444 (-2.245278) (metres 10000) -- >>> let p2 = s84Pos 55.605833 13.035833 (metres 20000) -- >>> interpolate p1 p2 0.5 -- 54°47'0.805"N,5°11'41.947"E 15.0km (S84) -- interpolate :: (Spherical a) => Position a -> Position a -> Double -> Position a interpolate p0 p1 f | f < 0 || f > 1 = error ("fraction must be in range [0..1], was " ++ show f) | f == 0 = p0 | f == 1 = p1 | otherwise = nvh iv ih (model p0) where nv0 = nvec p0 h0 = height p0 nv1 = nvec p1 h1 = height p1 iv = vunit (vadd nv0 (vscale (vsub nv1 nv0) f)) ih = lrph h0 h1 f -- | Computes the intersection between the two given minor arcs of great circle. -- -- see also 'intersections' -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let a1 = minorArcBetween (s84Pos 51.885 0.235 zero) (s84Pos 48.269 13.093 zero) -- >>> let a2 = minorArcBetween (s84Pos 49.008 2.549 zero) (s84Pos 56.283 11.304 zero) -- >>> join (intersection <$> a1 <*> a2) -- Just 50°54'6.260"N,4°29'39.052"E 0.0m (S84) -- intersection :: (Spherical a) => MinorArc a -> MinorArc a -> Maybe (Position a) intersection a1@(MinorArc n1 s1 _) a2@(MinorArc n2 _ _) = case intersections' n1 n2 (model s1) of Nothing -> Nothing (Just (i1, i2)) | isBetween i1 a1 && isBetween i1 a2 -> Just i1 | isBetween i2 a1 && isBetween i2 a2 -> Just i2 | otherwise -> Nothing -- | Computes the intersections between the two given 'GreatCircle's. -- Two great circles intersect exactly twice unless there are equal (regardless of orientation), -- in which case 'Nothing' is returned. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let gc1 = greatCircleHeadingOn (s84Pos 51.885 0.235 zero) (decimalDegrees 108.63) -- >>> let gc2 = greatCircleHeadingOn (s84Pos 49.008 2.549 zero) (decimalDegrees 32.72) -- >>> intersections gc1 gc2 -- Just (50°54'6.201"N,4°29'39.402"E 0.0m (S84),50°54'6.201"S,175°30'20.598"W 0.0m (S84)) -- >>> let i = intersections gc1 gc2 -- fmap fst i == fmap (antipode . snd) i -- >>> True -- intersections :: (Spherical a) => GreatCircle a -> GreatCircle a -> Maybe (Position a, Position a) intersections (GreatCircle n1 m _) (GreatCircle n2 _ _) = intersections' n1 n2 m -- | @isBetween p a@ determines whether position @p@ is within the minor arc -- of great circle @a@. -- -- If @p@ is not on the arc, returns whether @p@ is within the area bound -- by perpendiculars to the arc at each point (in the same hemisphere). -- isBetween :: (Spherical a) => Position a -> MinorArc a -> Bool isBetween p (MinorArc _ s e) = between && hemisphere where v0 = nvec p v1 = nvec s v2 = nvec e v10 = vsub v0 v1 v12 = vsub v2 v1 v20 = vsub v0 v2 v21 = vsub v1 v2 e1 = vdot v10 v12 -- p is on e side of s e2 = vdot v20 v21 -- p is on s side of e between = e1 >= 0 && e2 >= 0 hemisphere = vdot v0 v1 >= 0 && vdot v0 v2 >= 0 -- | @isInsideSurface p ps@ determines whether position @p@ is inside the __surface__ polygon defined by -- positions @ps@ (i.e. ignoring the height of the positions). -- The polygon can be opened or closed (i.e. if @head ps /= last ps@). -- -- Uses the angle summation test: on a sphere, due to spherical excess, enclosed point angles -- will sum to less than 360°, and exterior point angles will be small but non-zero. -- -- Always returns 'False' if @ps@ does not at least defines a triangle. -- -- ==== __Examples__ -- -- >>> import Data.Geo.Jord.GreatCircle -- >>> import Data.Geo.Jord.Position -- >>> -- >>> let malmo = s84Pos 55.6050 13.0038 zero -- >>> let ystad = s84Pos 55.4295 13.82 zero -- >>> let lund = s84Pos 55.7047 13.1910 zero -- >>> let helsingborg = s84Pos 56.0465 12.6945 zero -- >>> let kristianstad = s84Pos 56.0294 14.1567 zero -- >>> let polygon = [malmo, ystad, kristianstad, helsingborg, lund] -- >>> let hoor = s84Pos 55.9295 13.5297 zero -- >>> let hassleholm = s84Pos 56.1589 13.7668 zero -- >>> isInsideSurface hoor polygon -- True -- >>> isInsideSurface hassleholm polygon -- False -- isInsideSurface :: (Spherical a) => Position a -> [Position a] -> Bool isInsideSurface p ps | null ps = False | llEq (head ps) (last ps) = isInsideSurface p (init ps) | length ps < 3 = False | otherwise = let aSum = foldl (\a v' -> add a (uncurry signedAngle v' (Just v))) (decimalDegrees 0) (egdes (map (vsub v) vs)) in abs (toDecimalDegrees aSum) > 180.0 where v = nvec p vs = fmap nvec ps -- | @mean ps@ computes the geographic mean surface position of @ps@, if it is defined. -- -- The geographic mean is not defined for antipodals positions (since they -- cancel each other). -- -- Special conditions: -- -- @ -- mean [] = Nothing -- mean [p] = Just p -- mean [p1, p2, p3] = Just circumcentre -- mean [p1, .., antipode p1] = Nothing -- @ mean :: (Spherical a) => [Position a] -> Maybe (Position a) mean [] = Nothing mean [p] = Just p mean ps = if null antipodals then Just (nvh nv zero (model . head $ ps)) else Nothing where vs = fmap nvec ps ts = filter (\l -> length l == 2) (subsequences vs) antipodals = filter (\t -> (realToFrac (vnorm (vadd (head t) (last t)) :: Double) :: Nano) == 0) ts nv = vunit $ foldl vadd vzero vs -- private alongTrackDistance'' :: (Spherical a) => Position a -> Position a -> Vector3d -> Length alongTrackDistance'' p s n = arcLength a (radius s) where a = signedAngle (nvec s) (vcross (vcross n (nvec p)) n) (Just n) -- | [p1, p2, p3, p4] to [(p1, p2), (p2, p3), (p3, p4), (p4, p1)] egdes :: [Vector3d] -> [(Vector3d, Vector3d)] egdes ps = zip ps (tail ps ++ [head ps]) lrph :: Length -> Length -> Double -> Length lrph h0 h1 f = metres h where h0' = toMetres h0 h1' = toMetres h1 h = h0' + (h1' - h0') * f intersections' :: (Spherical a) => Vector3d -> Vector3d -> a -> Maybe (Position a, Position a) intersections' n1 n2 s | (vnorm i :: Double) == 0.0 = Nothing | otherwise , let ni = nvh (vunit i) zero s = Just (ni, antipode ni) where i = vcross n1 n2 initialBearing' :: Position a -> Position a -> Angle initialBearing' p1 p2 = normalise a (decimalDegrees 360) where v1 = nvec p1 v2 = nvec p2 gc1 = vcross v1 v2 -- great circle through p1 & p2 gc2 = vcross v1 nvNorthPole -- great circle through p1 & north pole a = radians (signedAngleRadians gc1 gc2 (Just v1)) -- | reference sphere radius. radius :: (Spherical a) => Position a -> Length radius = equatorialRadius . surface . model normal' :: (Spherical a) => Position a -> Position a -> Vector3d normal' p1 p2 = vcross (nvec p1) (nvec p2) signedAngle :: Vector3d -> Vector3d -> Maybe Vector3d -> Angle signedAngle v1 v2 n = radians (signedAngleRadians v1 v2 n)