diff --git a/AddingProjections.txt b/AddingProjections.txt new file mode 100644 --- /dev/null +++ b/AddingProjections.txt @@ -0,0 +1,79 @@ +How to add projections to the Geodetics library +=============================================== + +The most useful thing you can do to help the Geodetics library is to +add new projections: + +New projection types +--------------------- + +A projection type is a generic mapping of an ellipsoid to a +plane. The library already contains the Stereographic and +Transverse Mercator projections. To add a new projection type: + +* Create a new module under Geodetics to contain your projection. + +* Define your new projection as a data type containing whatever + parameters your projection needs to define it. If your projection + includes a False Origin (most do) then use the GridOffset type. A + scaling factor should be a Dimensionless Double. + + If the projection has intermediate values derived purely from these + parameters then you can make it more efficient if you include them + in the data type rather than recalculating them every time. If you + do this then don't export the constructor, and have a separate + function to fill out the data structure correctly. See Stereographic + for an example. + +* Make your new projection an instance of GridClass and define the + "toGrid", "fromGrid" and "gridEllipsoid" functions. Generic grids + should have the ellipsoid as a parameter. + +National and other standard projections +--------------------------------------- + +A standard projection is a distinguished value of a generic +projection. For instance the UK National Grid is a value of the +TransverseMercator projection. + +A standard projection should be based on a generic projection. There +is little point in duplicating code for two standard projections if +they are based on the same underlying mathematics. Therefore the +design for a standard projection contains the following elements: + +* A representation of the projection as an untyped value of the + underlying generic projection. For instance the UK OS grid is + represented by "ukGrid :: GridTM OSGB36". + +* A unit type (that is, one constructor with no arguments) that is an + instance of the GridClass. Implement the "toGrid" and "fromGrid" + functions using the "unsafeGridCoerce" function and the untyped + value defined above. For instace the UkNationalGrid type has the + following implementations in its GridClass instance: + + toGrid _ = unsafeGridCoerce UkNationalGrid . toGrid ukGrid + fromGrid = fromGrid . unsafeGridCoerce ukGrid + +If your standard projection has a particular ellipsoid that is not +already defined (for instance, the UK national grid is based on +OSGB36) then this should be defined as a unit type in the same module +and made an instance of Ellipsoid. + +If there is a grid reference system associated with the projection +then this should be implemented as a pair of functions between strings +and gridpoints. When converting from a string make sure that some part +of the result indicates the precision of the reference. One way is to +return a gridpoint for the southwest corner of the grid square and an +offset to the centre. + +Testing +------- + +In the "test" folder create a new module for your projection and +import it into "Main.hs". All projections should have several test +points derived from a known oracle (such as a pre-existing +implementation or reference table). Cite the source in a comment to +your test data. In addition generic projections should have a "round +trip" test using QuickCheck to show that (fromGrid . toGrid) is an +identity to within some specified accuracy. + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,12 @@ +Copyright (c) 2014, Paul Johnson +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. + +3. Neither the name of Paul Johnson nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,21 @@ +geodetics +========= + +Haskell library of data types and calculations for positions on planet Earth + +This library provides "geodetic" positions. That is, latitude, longitude and altitude on a +specified Terrestrial Reference Frame (TRF). The basic TRF is the WGS84, which is the +one used by GPS and Google Earth. Others can be added by describing the underlying ellipsoid +and the difference in angle and centre with WGS84, and a position in one TRF can be +transformed into another. Given two points in the same TRF you can find the shortest distance between them +and the bearing from one to the other. + +Once you have a geodetic position defined you can project it onto a flat plane, or Grid. +At present Transverse Mercator and Oblique Stereographic grids are provided. More can be +added by defining new instances of the Grid typeclass: see "AddingProjections.txt" for +detais. + +The Paths module defines a path as a parametric function of distance that returns a +position and a bearing. Given two paths you can find their intersection using a fast +iterative algorithm. + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +import Distribution.Simple + +main = defaultMain diff --git a/ToDo.txt b/ToDo.txt new file mode 100644 --- /dev/null +++ b/ToDo.txt @@ -0,0 +1,6 @@ +To Do +===== + +* Add UTM and Universal Polar grids. + +* Add more national grid systems. diff --git a/geodetics.cabal b/geodetics.cabal new file mode 100644 --- /dev/null +++ b/geodetics.cabal @@ -0,0 +1,77 @@ +name: geodetics +version: 0.0.1 +cabal-version: >= 1.10 +build-type: Simple +author: Paul Johnson +data-files: + AddingProjections.txt, + LICENSE, + README.md, + ToDo.txt +license: BSD3 +copyright: Paul Johnson 2014. +synopsis: Terrestrial coordinate systems and associated calculations. +description: Precise geographical coordinates (latitude & longitude), with conversion between + different reference frames and projections. + . + Certain distinguished reference frames and grids are given distinct + types so that coordinates expressed within them cannot be confused with + from coordinates in other frames. +license-file: LICENSE +maintainer: Paul Johnson +homepage: https://github.com/PaulJohnson/geodetics +category: Geography +tested-with: GHC==7.6.3 + +source-repository head + type: git + location: https://github.com/PaulJohnson/geodetics + +library + hs-source-dirs: src + build-depends: + base >= 4.6 && < 4.7, + dimensional == 0.13.*, + array >= 0.4, + parsec >=3.1.3 && <3.2 + ghc-options: -Wall + exposed-modules: + Geodetics.Altitude, + Geodetics.Ellipsoids, + Geodetics.Geodetic, + Geodetics.Grid, + Geodetics.LatLongParser, + Geodetics.Path, + Geodetics.Stereographic, + Geodetics.TransverseMercator, + Geodetics.UK + Default-Language: Haskell2010 + +test-suite GeodeticTest + type: exitcode-stdio-1.0 + main-is: Main.hs + x-uses-tf: true + build-depends: + base >= 4, + HUnit >= 1.2 && < 2, + dimensional == 0.13.*, + QuickCheck >= 2.4, + test-framework >= 0.4.1, + test-framework-quickcheck2, + test-framework-hunit, + array >= 0.4 + hs-source-dirs: + src, + test + ghc-options: -Wall -rtsopts + other-modules: + ArbitraryInstances, + Main, + Geodetics.Geodetic, + Geodetics.Altitude, + Geodetics.Grid, + Geodetics.TransverseMercator, + Geodetics.Path, + Geodetics.Stereographic, + Geodetics.LatLongParser + Default-Language: Haskell2010 diff --git a/src/Geodetics/Altitude.hs b/src/Geodetics/Altitude.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Altitude.hs @@ -0,0 +1,17 @@ +module Geodetics.Altitude ( + HasAltitude (..) +) where + +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + +-- | All geographical coordinate systems need the concept of altitude above a reference point, usually associated with +-- local sea level. +-- +-- Minimum definition: altitude, setAltitude. +class HasAltitude a where + altitude :: a -> Length Double + setAltitude :: Length Double -> a -> a + -- | Set altitude to zero. + groundPosition :: a -> a + groundPosition = setAltitude _0 diff --git a/src/Geodetics/Ellipsoids.hs b/src/Geodetics/Ellipsoids.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Ellipsoids.hs @@ -0,0 +1,263 @@ +{-# LANGUAGE FlexibleContexts #-} + +{- | An Ellipsoid is a reasonable best fit for the surface of the +Earth over some defined area. WGS84 is the standard used for the whole +of the Earth. Other Ellipsoids are considered a best fit for some +specific area. +-} + +module Geodetics.Ellipsoids ( + -- ** Helmert transform between geodetic reference systems + Helmert (..), + inverseHelmert, + ECEF, + applyHelmert, + -- ** Ellipsoid models of the Geoid + Ellipsoid (..), + WGS84 (..), + LocalEllipsoid (..), + flattening, + minorRadius, + eccentricity2, + eccentricity'2, + -- ** Auxiliary latitudes and related Values + normal, + latitudeRadius, + meridianRadius, + primeVerticalRadius, + isometricLatitude, + -- ** Tiny linear algebra library for 3D vectors + Vec3, + Matrix3, + add3, + scale3, + negate3, + transform3, + invert3, + trans3, + dot3, + cross3 +) where + +import Data.Monoid +import Numeric.Units.Dimensional +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + + +-- | 3d vector as @(X,Y,Z)@. +type Vec3 a = (a,a,a) + +-- | 3x3 transform matrix for Vec3. +type Matrix3 a = Vec3 (Vec3 a) + + +-- | Multiply a vector by a scalar. +scale3 :: (Num a, Mul d d' d'') => + Vec3 (Dimensional DQuantity d a) -> Dimensional DQuantity d' a -> Vec3 (Dimensional DQuantity d'' a) +scale3 (x,y,z) s = (x*s, y*s, z*s) + + +-- | Negation of a vector. +negate3 :: (Num a) => Vec3 (Quantity d a) -> Vec3 (Quantity d a) +negate3 (x,y,z) = (negate x, negate y, negate z) + +-- | Add two vectors +add3 :: (Num a) => Vec3 (Quantity d a) -> Vec3 (Quantity d a) -> Vec3 (Quantity d a) +add3 (x1,y1,z1) (x2,y2,z2) = (x1+x2, y1+y2, z1+z2) + + +-- | Multiply a matrix by a vector in the Dimensional type system. +transform3 :: (Num a, Mul d d' d'') => + Matrix3 (Dimensional DQuantity d a) -> Vec3 (Dimensional DQuantity d' a) -> Vec3 (Dimensional DQuantity d'' a) +transform3 (tx,ty,tz) v = (t tx v, t ty v, t tz v) + where + t (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2 + + +-- | Inverse of a 3x3 matrix. +invert3 :: (Fractional a, Mul d d d2, Mul d2 d d3, Div d2 d3 d1') => + Matrix3 (Dimensional DQuantity d a) -> Matrix3 (Dimensional DQuantity d1' a) +invert3 ((x1,y1,z1), + (x2,y2,z2), + (x3,y3,z3)) = + ((det2 y2 z2 y3 z3 / det, det2 z1 y1 z3 y3 / det, det2 y1 z1 y2 z2 / det), + (det2 z2 x2 z3 x3 / det, det2 x1 z1 x3 z3 / det, det2 z1 x1 z2 x2 / det), + (det2 x2 y2 x3 y3 / det, det2 y1 x1 y3 x3 / det, det2 x1 y1 x2 y2 / det)) + where + det = (x1*y2*z3 + y1*z2*x3 + z1*x2*y3) - (z1*y2*x3 + y1*x2*z3 + x1*z2*y3) + det2 a b c d = a*d - b*c + +-- | Transpose of a 3x3 matrix. +trans3 :: Matrix3 a -> Matrix3 a +trans3 ((x1,y1,z1),(x2,y2,z2),(x3,y3,z3)) = ((x1,x2,x3),(y1,y2,y3),(z1,z2,z3)) + + +-- | Dot product of two vectors +dot3 :: (Num a, Mul d1 d2 d3) => + Vec3 (Dimensional DQuantity d1 a) -> Vec3 (Dimensional DQuantity d2 a) -> Dimensional DQuantity d3 a +dot3 (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2 + +-- | Cross product of two vectors +cross3 :: (Num a, Mul d1 d2 d3) => + Vec3 (Dimensional DQuantity d1 a) -> Vec3 (Dimensional DQuantity d2 a) -> Vec3 (Dimensional DQuantity d3 a) +cross3 (x1,y1,z1) (x2,y2,z2) = (y1*z2 - z1*y2, z1*x2 - x1*z2, x1*y2 - y1*x2) + + +-- | The 7 parameter Helmert transformation. The monoid instance allows composition. +data Helmert = Helmert { + cX, cY, cZ :: Length Double, + helmertScale :: Dimensionless Double, -- ^ Parts per million + rX, rY, rZ :: Dimensionless Double } deriving (Eq, Show) + +instance Monoid Helmert where + mempty = Helmert (0 *~ meter) (0 *~ meter) (0 *~ meter) _1 _0 _0 _0 + mappend h1 h2 = Helmert (cX h1 + cX h2) (cY h1 + cY h2) (cZ h1 + cZ h2) + (helmertScale h1 + helmertScale h2) + (rX h1 + rX h2) (rY h1 + rY h2) (rZ h1 + rZ h2) + + +-- | The inverse of a Helmert transformation. +inverseHelmert :: Helmert -> Helmert +inverseHelmert h = Helmert (negate $ cX h) (negate $ cY h) (negate $ cZ h) + (negate $ helmertScale h) + (negate $ rX h) (negate $ rY h) (negate $ rZ h) + + +-- | Earth-centred, Earth-fixed coordinates as a vector. The origin and axes are +-- not defined: use with caution. +type ECEF = Vec3 (Length Double) + +-- | Apply a Helmert transformation to earth-centered coordinates. +applyHelmert:: Helmert -> ECEF -> ECEF +applyHelmert h (x,y,z) = ( + cX h + s * ( x - rZ h * y + rY h * z), + cY h + s * ( rZ h * x + y - rX h * z), + cZ h + s * (negate (rY h) * x + rX h * y + z)) + where + s = _1 + helmertScale h * (1e-6 *~ one) + + +-- | An Ellipsoid is defined by the major radius and the inverse flattening (which define its shape), +-- and its Helmert transform relative to WGS84 (which defines its position and orientation). +-- +-- The inclusion of the Helmert parameters relative to WGS84 actually make this a Terrestrial +-- Reference Frame (TRF), but the term "Ellipsoid" will be used in this library for readability. +-- +-- Minimum definition: @majorRadius@, @flatR@ & @helmert@. +-- +-- Laws: +-- +-- > helmertToWGS84 = applyHelmert . helmert +-- > helmertFromWGS84 e . helmertToWGS84 e = id +class (Show a, Eq a) => Ellipsoid a where + majorRadius :: a -> Length Double + flatR :: a -> Dimensionless Double + -- ^ Inverse of the flattening. + helmert :: a -> Helmert + helmertToWSG84 :: a -> ECEF -> ECEF + -- ^ The Helmert transform that will convert a position wrt + -- this ellipsoid into a position wrt WGS84. + helmertToWSG84 e = applyHelmert (helmert e) + helmertFromWSG84 :: a -> ECEF -> ECEF + -- ^ And its inverse. + helmertFromWSG84 e = applyHelmert (inverseHelmert $ helmert e) + + +-- | The WGS84 geoid, major radius 6378137.0 meters, flattening = 1 / 298.257223563 +-- as defined in \"Technical Manual DMA TM 8358.1 - Datums, Ellipsoids, Grids, and +-- Grid Reference Systems\" at the National Geospatial-Intelligence Agency (NGA). +-- +-- The WGS84 has a special place in this library as the standard Ellipsoid against +-- which all others are defined. +data WGS84 = WGS84 + +instance Eq WGS84 where _ == _ = True + +instance Show WGS84 where + show _ = "WGS84" + +instance Ellipsoid WGS84 where + majorRadius _ = 6378137.0 *~ meter + flatR _ = 298.257223563 *~ one + helmert _ = mempty + helmertToWSG84 _ = id + helmertFromWSG84 _ = id + + +-- | Ellipsoids other than WGS84, used within a defined geographical area where +-- they are a better fit to the local geoid. Can also be used for historical ellipsoids. +-- +-- The @Show@ instance just returns the name. +-- Creating two different local ellipsoids with the same name is a Bad Thing. +data LocalEllipsoid = LocalEllipsoid { + nameLocal :: String, + majorRadiusLocal :: Length Double, + flatRLocal :: Dimensionless Double, + helmertLocal :: Helmert } deriving (Eq) + +instance Show LocalEllipsoid where + show = nameLocal + +instance Ellipsoid LocalEllipsoid where + majorRadius = majorRadiusLocal + flatR = flatRLocal + helmert = helmertLocal + + +-- | Flattening (f) of an ellipsoid. +flattening :: (Ellipsoid e) => e -> Dimensionless Double +flattening e = _1 / flatR e + +-- | The minor radius of an ellipsoid. +minorRadius :: (Ellipsoid e) => e -> Length Double +minorRadius e = majorRadius e * (_1 - flattening e) + + +-- | The eccentricity squared of an ellipsoid. +eccentricity2 :: (Ellipsoid e) => e -> Dimensionless Double +eccentricity2 e = _2 * f - (f * f) where f = flattening e + +-- | The second eccentricity squared of an ellipsoid. +eccentricity'2 :: (Ellipsoid e) => e -> Dimensionless Double +eccentricity'2 e = (f * (_2 - f)) / (_1 - f * f) where f = flattening e + + +-- | Distance from the surface at the specified latitude to the +-- axis of the Earth straight down. Also known as the radius of +-- curvature in the prime vertical, and often denoted @N@. +normal :: (Ellipsoid e) => e -> Angle Double -> Length Double +normal e lat = majorRadius e / sqrt (_1 - eccentricity2 e * sin lat ^ pos2) + + +-- | Radius of the circle of latitude: the distance from a point +-- at that latitude to the axis of the Earth. +latitudeRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double +latitudeRadius e lat = normal e lat * cos lat + + +-- | Radius of curvature in the meridian at the specified latitude. +-- Often denoted @M@. +meridianRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double +meridianRadius e lat = + majorRadius e * (_1 - eccentricity2 e) + / sqrt ((_1 - eccentricity2 e * sin lat ^ pos2) ^ pos3) + + +-- | Radius of curvature of the ellipsoid perpendicular to the meridian at the specified latitude. +primeVerticalRadius :: (Ellipsoid e) => e -> Angle Double -> Length Double +primeVerticalRadius e lat = + majorRadius e / sqrt (_1 - eccentricity2 e * sin lat ^ pos2) + + +-- | The isometric latitude. The isometric latitude is conventionally denoted by ψ +-- (not to be confused with the geocentric latitude): it is used in the development +-- of the ellipsoidal versions of the normal Mercator projection and the Transverse +-- Mercator projection. The name "isometric" arises from the fact that at any point +-- on the ellipsoid equal increments of ψ and longitude λ give rise to equal distance +-- displacements along the meridians and parallels respectively. +isometricLatitude :: (Ellipsoid e) => e -> Angle Double -> Angle Double +isometricLatitude ellipse lat = atanh sinLat - e * atanh (e * sinLat) + where + sinLat = sin lat + e = sqrt $ eccentricity2 ellipse diff --git a/src/Geodetics/Geodetic.hs b/src/Geodetics/Geodetic.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Geodetic.hs @@ -0,0 +1,298 @@ +module Geodetics.Geodetic ( + -- ** Geodetic Coordinates + Geodetic (..), + readGroundPosition, + toLocal, + toWGS84, + antipode, + geometricalDistance, + geometricalDistanceSq, + groundDistance, + properAngle, + showAngle, + -- ** Earth Centred Earth Fixed Coordinates + ECEF, + geoToEarth, + earthToGeo, + -- ** Re-exported for convenience + WGS84 (..) +) where + + +import Data.Char (chr) +import Data.Function +import Data.Maybe +import Data.Monoid +import Geodetics.Altitude +import Geodetics.Ellipsoids +import Geodetics.LatLongParser +import Numeric.Units.Dimensional.Prelude +import Text.ParserCombinators.ReadP +import qualified Prelude as P + +-- | Defines a three-D position on or around the Earth using latitude, +-- longitude and altitude with respect to a specified ellipsoid, with +-- positive directions being North and East. The default "show" +-- instance gives position in degrees, minutes and seconds to 5 decimal +-- places, which is a +-- resolution of about 1m on the Earth's surface. Internally latitude +-- and longitude are stored as double precision radians. Convert to +-- degrees using e.g. @latitude g /~ degree@. +-- +-- The functions here deal with altitude by assuming that the local +-- height datum is always co-incident with the ellipsoid in use, +-- even though the \"mean sea level\" (the usual height datum) can be tens +-- of meters above or below the ellipsoid, and two ellipsoids can +-- differ by similar amounts. This is because the altitude is +-- usually known with reference to a local datum regardless of the +-- ellipsoid in use, so it is simpler to preserve the altitude across +-- all operations. However if +-- you are working with ECEF coordinates from some other source then +-- this may give you the wrong results, depending on the altitude +-- correction your source has used. +-- +-- There is no "Eq" instance because comparing two arbitrary +-- co-ordinates on the Earth is a non-trivial exercise. Clearly if all +-- the parameters are equal on the same ellipsoid then they are indeed +-- in the same place. However if different ellipsoids are used then +-- two co-ordinates with different numbers can still refer to the same +-- physical location. If you want to find out if two co-ordinates are +-- the same to within a given tolerance then use "geometricDistance" +-- (or its squared variant to avoid an extra @sqrt@ operation). +data (Ellipsoid e) => Geodetic e = Geodetic { + latitude, longitude :: Angle Double, + geoAlt :: Length Double, + ellipsoid :: e +} + +instance (Ellipsoid e) => Show (Geodetic e) where + show g = concat [ + showAngle (abs $ latitude g), " ", letter "SN" (latitude g), ", ", + showAngle (abs $ longitude g), " ", letter "WE" (longitude g), ", ", + show (altitude g), " ", show (ellipsoid g)] + where letter s n = [s !! (if n < _0 then 0 else 1)] + + + +-- | Read the latitude and longitude of a ground position and +-- return a Geodetic position on the specified ellipsoid. +-- +-- The latitude and longitude may be in any of the following formats. +-- The comma between latitude and longitude is optional in all cases. +-- Latitude must always be first. +-- +-- * Signed decimal degrees: 34.52327, -46.23234 +-- +-- * Decimal degrees NSEW: 34.52327N, 46.23234W +-- +-- * Degrees and decimal minutes (units optional): 34° 31.43' N, 46° 13.92' +-- +-- * Degrees, minutes and seconds (units optional): 34° 31' 23.52\" N, 46° 13' 56.43\" W +-- +-- * DDDMMSS format with optional leading zeros: 343123.52N, 0461356.43W +readGroundPosition :: (Ellipsoid e) => e -> String -> Maybe (Geodetic e) +readGroundPosition e str = + case map fst $ filter (null . snd) $ readP_to_S latLong str of + [] -> Nothing + (lat,long) : _ -> Just $ groundPosition $ Geodetic (lat *~ degree) (long *~ degree) undefined e + + +-- | Show an angle as degrees, minutes and seconds to two decimal places. +showAngle :: Angle Double -> String +showAngle a + | isNaN a1 = "NaN" -- Not a Nangle + | isInfinite a1 = sgn ++ "Infinity" + | otherwise = concat [sgn, show d, [chr 0xB0, ' '], + show m, "' ", + show s, ".", dstr, "\"" ] + where + a1 = a /~ one + sgn = if a < _0 then "-" else "" + centisecs :: Integer + centisecs = P.abs $ P.round $ (a /~ degree) P.* 360000 -- hundredths of arcsec per degree. + (d, m1) = centisecs `P.divMod` 360000 + (m, s1) = m1 `P.divMod` 6000 -- hundredths of arcsec per arcmin + (s, ds) = s1 `P.divMod` 100 + dstr = reverse $ take 2 $ reverse (show ds) ++ "00" -- Decimal fraction with zero padding. + + +instance (Ellipsoid e) => HasAltitude (Geodetic e) where + altitude = geoAlt + setAltitude h g = g{geoAlt = h} + + + +-- | The point on the Earth diametrically opposite the argument, with +-- the same altitude. +antipode :: (Ellipsoid e) => Geodetic e -> Geodetic e +antipode g = Geodetic lat long (geoAlt g) (ellipsoid g) + where + lat = negate $ latitude g + long' = longitude g - 180 *~ degree + long | long' < _0 = long' + 360 *~ degree + | otherwise = long' + + + +-- | Convert a geodetic coordinate into earth centered, relative to the +-- ellipsoid in use. +geoToEarth :: (Ellipsoid e) => Geodetic e -> ECEF +geoToEarth geo = ( + (n + h) * coslat * coslong, + (n + h) * coslat * sinlong, + (n * (_1 - eccentricity2 e) + h) * sinlat) + where + n = normal e $ latitude geo + e = ellipsoid geo + coslat = cos $ latitude geo + coslong = cos $ longitude geo + sinlat = sin $ latitude geo + sinlong = sin $ longitude geo + h = altitude geo + + +-- | Convert an earth centred coordinate into a geodetic coordinate on +-- the specified geoid. +-- +-- Uses the closed form solution of H. Vermeille: Direct +-- transformation from geocentric coordinates to geodetic coordinates. +-- Journal of Geodesy Volume 76, Number 8 (2002), 451-454 +earthToGeo :: (Ellipsoid e) => e -> ECEF -> (Angle Double, Angle Double, Length Double) +earthToGeo e (x,y,z) = (phi, atan2 y x, sqrt (l ^ pos2 + p2) - norm) + where + -- Naming: numeric suffix inicates power. Hence x2 = x * x, x3 = x2 * x, etc. + p2 = x ^ pos2 + y ^ pos2 + a = majorRadius e + a2 = a ^ pos2 + e2 = eccentricity2 e + e4 = e2 ^ pos2 + zeta = (_1-e2) * (z ^ pos2 / a2) + rho = (p2 / a2 + zeta - e4) / _6 + rho2 = rho ^ pos2 + rho3 = rho * rho2 + s = e4 * zeta * p2 / (_4 * a2) + t = cbrt (s + rho3 + sqrt (s * (s + _2 * rho3))) + u = rho + t + rho2 / t + v = sqrt (u ^ pos2 + e4 * zeta) + w = e2 * (u + v - zeta) / (_2 * v) + kappa = _1 + e2 * (sqrt (u + v + w ^ pos2) + w) / (u + v) + phi = atan (kappa * z / sqrt p2) + norm = normal e phi + l = z + e2 * norm * sin phi + + +-- | Convert a position from any geodetic to another one, assuming local altitude stays constant. +toLocal :: (Ellipsoid e1, Ellipsoid e2) => e2 -> Geodetic e1 -> Geodetic e2 +toLocal e2 g = Geodetic lat lon alt e2 + where + alt = altitude g + (lat, lon, _) = earthToGeo e2 $ applyHelmert h $ geoToEarth g + h = helmert (ellipsoid g) `mappend` inverseHelmert (helmert e2) + +-- | Convert a position from any geodetic to WGS84, assuming local +-- altitude stays constant. +toWGS84 :: (Ellipsoid e) => Geodetic e -> Geodetic WGS84 +toWGS84 g = Geodetic lat lon alt WGS84 + where + alt = altitude g + (lat, lon, _) = earthToGeo WGS84 $ applyHelmert h $ geoToEarth g + h = helmert (ellipsoid g) + + +-- | The absolute distance in a straight line between two geodetic +-- points. They must be on the same ellipsoid. +-- Note that this is not the geodetic distance taken by following +-- the curvature of the earth. +geometricalDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Length Double +geometricalDistance g1 g2 = sqrt $ geometricalDistanceSq g1 g2 + +-- | The square of the absolute distance. Comes out as "Area" type of course. +geometricalDistanceSq :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Area Double +geometricalDistanceSq g1 g2 = (x1-x2) ^ pos2 + (y1-y2) ^ pos2 + (z1-z2) ^ pos2 + where + (x1,y1,z1) = geoToEarth g1 + (x2,y2,z2) = geoToEarth g2 + + +-- | The shortest ellipsoidal distance between two points on the +-- ground with reference to the same ellipsoid. Altitude is ignored. +-- +-- The results are the distance between the points, the bearing of +-- the second point from the first, and (180 degrees - the bearing +-- of the first point from the second). +-- +-- The algorithm can fail to converge where the arguments are near to +-- antipodal. In this case it returns @Nothing@. +-- +-- Uses Vincenty's formula. \"Direct and inverse solutions of +-- geodesics on the ellipsoid with application of nested +-- equations\". T. Vincenty. Survey Review XXII 176, April +-- 1975. http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf +groundDistance :: (Ellipsoid e) => Geodetic e -> Geodetic e -> + Maybe (Length Double, Dimensionless Double, Dimensionless Double) +groundDistance p1 p2 = do + (_, (lambda, (cos2Alpha, delta, sinDelta, cosDelta, cos2DeltaM))) <- + listToMaybe $ dropWhile converging $ take 100 $ zip lambdas $ tail lambdas + let + uSq = cos2Alpha * (a^pos2 - b^pos2) / b^pos2 + bigA = _1 + uSq/(16384*~one) * ((4096*~one) + uSq * + (((-768)*~one) + uSq * ((320*~one) + - (175*~one)*uSq))) + bigB = uSq/(1024*~one) * ((256*~one) + + uSq * (((-128)*~one) + + uSq * ((74*~one) - (47*~one)*uSq))) + deltaDelta = + bigB * sinDelta * (cos2DeltaM + + bigB/_4 * (cosDelta * (_2 * cos2DeltaM^pos2 - _1) + - bigB/_6 * cos2DeltaM * (_4 * sinDelta^pos2 - _3) + * (_4 * cos2DeltaM - _3))) + s = b * bigA * (delta - deltaDelta) + alpha1 = atan2(cosU2 * sin lambda) (cosU1 * sinU2 - sinU1 * cosU2 * cos lambda) + alpha2 = atan2(cosU1 * sin lambda) (cosU1 * sinU2 * cos lambda - sinU1 * cosU2) + return (s, alpha1, alpha2) + where + f = flattening $ ellipsoid p1 + a = majorRadius $ ellipsoid p1 + b = minorRadius $ ellipsoid p1 + l = abs $ longitude p1 - longitude p2 + u1 = atan ((_1-f) * tan (latitude p1)) + u2 = atan ((_1-f) * tan (latitude p2)) + sinU1 = sin u1 + cosU1 = cos u1 + sinU2 = sin u2 + cosU2 = cos u2 + + nextLambda lambda = (lambda1, (cos2Alpha, delta, sinDelta, cosDelta, cos2DeltaM)) + where + sinLambda = sin lambda + cosLambda = cos lambda + sinDelta = sqrt((cosU2 * sinLambda) ^ pos2 + + (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) ^ pos2) + cosDelta = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda + delta = atan2 sinDelta cosDelta + sinAlpha = cosU1 * cosU2 * sinLambda / sinDelta + cos2Alpha = _1 - sinAlpha ^ pos2 + cos2DeltaM = if cos2Alpha == _0 + then _0 + else cosDelta - _2 * sinU1 * sinU2 / cos2Alpha + c = f/(16 *~ one) * cos2Alpha * (_4 + f * (_4 - _3 * cos2Alpha)) + lambda1 = l + (_1-c) * f * sinAlpha + * (delta + c * sinDelta + * (cos2DeltaM + c * cosDelta *(_2 * cos2DeltaM ^ pos2 - _1))) + lambdas = iterate (nextLambda . fst) (l, undefined) + converging ((l1,_),(l2,_)) = abs (l1 - l2) > (1e-14 *~ one) + + +-- | Add or subtract multiples of 2*pi so that for all @t@, @-pi < properAngle t < pi@. +properAngle :: Angle Double -> Angle Double +properAngle t + | r1 <= negate pi = r1 + pi2 + | r1 > pi = r1 - pi2 + | otherwise = r1 + where + pf :: Double -> (Int, Double) + pf = properFraction -- Shut up GHC warning about defaulting to Integer. + (_,r) = pf (t/pi2 /~ one) + r1 = (r *~ one) * pi2 + pi2 = pi * _2 + diff --git a/src/Geodetics/Grid.hs b/src/Geodetics/Grid.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Grid.hs @@ -0,0 +1,184 @@ +{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-} + +module Geodetics.Grid ( + -- ** Grid types + GridClass (..), + GridPoint (..), + GridOffset (..), + -- ** Grid operations + polarOffset, + offsetScale, + offsetNegate, + applyOffset, + offsetDistance, + offsetDistanceSq, + offsetBearing, + gridOffset, + -- ** Unsafe conversion + unsafeGridCoerce, + -- ** Utility functions for grid references + fromGridDigits, + toGridDigits +) where + +import Data.Char +import Data.Function +import Data.Monoid +import Geodetics.Altitude +import Geodetics.Geodetic +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + +-- | A Grid is a two-dimensional projection of the ellipsoid onto a plane. Any given type of grid can +-- usually be instantiated with parameters such as a tangential point or line, and these parameters +-- will include the terrestrial reference frame ("Ellipsoid" in this library) used as a foundation. +-- Hence conversion from a geodetic to a grid point requires the \"basis\" for the grid in question, +-- and grid points carry that basis with them because without it there is no defined relationship +-- between the grid points and terrestrial positions. +class GridClass r e | r->e where + fromGrid :: GridPoint r -> Geodetic e + toGrid :: r -> Geodetic e -> GridPoint r + gridEllipsoid :: r -> e + + +-- | A point on the specified grid. +data GridPoint r = GridPoint { + eastings, northings, altGP :: Length Double, + gridBasis :: r +} deriving (Show) + + +instance Eq (GridPoint r) where + p1 == p2 = + eastings p1 == eastings p2 && + northings p1 == northings p2 && + altGP p1 == altGP p2 + +instance HasAltitude (GridPoint g) where + altitude = altGP + setAltitude h gp = gp{altGP = h} + + + +-- | A vector relative to a point on a grid. +-- Operations that use offsets will only give +-- meaningful results if all the points come from the same grid. +-- +-- The monoid instance is the sum of offsets. +data GridOffset = GridOffset { + deltaEast, deltaNorth, deltaAltitude :: Length Double +} deriving (Eq, Show) + +instance Monoid GridOffset where + mempty = GridOffset _0 _0 _0 + mappend g1 g2 = GridOffset (deltaEast g1 + deltaEast g2) + (deltaNorth g1 + deltaNorth g2) + (deltaAltitude g1 + deltaAltitude g2) + +-- | An offset defined by a distance and a bearing to the right of North. +-- +-- There is no elevation parameter because we are using a plane to approximate an ellipsoid, +-- so elevation would not provide a useful result. If you want to work with elevations +-- then "rayPath" will give meaningful results. +polarOffset :: Length Double -> Angle Double -> GridOffset +polarOffset r d = GridOffset (r * sin d) (r * cos d) _0 + + +-- | Scale an offset by a scalar. +offsetScale :: Dimensionless Double -> GridOffset -> GridOffset +offsetScale s off = GridOffset (deltaEast off * s) + (deltaNorth off * s) + (deltaAltitude off * s) + +-- | Invert an offset. +offsetNegate :: GridOffset -> GridOffset +offsetNegate off = GridOffset (negate $ deltaEast off) + (negate $ deltaNorth off) + (negate $ deltaAltitude off) + + +-- Add an offset on to a point to get another point. +applyOffset :: GridOffset -> GridPoint g -> GridPoint g +applyOffset off p = GridPoint (eastings p + deltaEast off) + (northings p + deltaNorth off) + (altitude p + deltaAltitude off) + (gridBasis p) + + +-- | The distance represented by an offset. +offsetDistance :: GridOffset -> Length Double +offsetDistance = sqrt . offsetDistanceSq + + +-- | The square of the distance represented by an offset. +offsetDistanceSq :: GridOffset -> Area Double +offsetDistanceSq off = + deltaEast off ^ pos2 + deltaNorth off ^ pos2 + deltaAltitude off ^ pos2 + + +-- | The direction represented by an offset, as bearing to the right of North. +offsetBearing :: GridOffset -> Angle Double +offsetBearing off = atan2 (deltaEast off) (deltaNorth off) + + +-- | The offset required to move from p1 to p2. +gridOffset :: GridPoint g -> GridPoint g -> GridOffset +gridOffset p1 p2 = GridOffset (eastings p2 - eastings p1) + (northings p2 - northings p1) + (altitude p2 - altitude p1) + + +-- | Coerce a grid point of one type into a grid point of a different type, +-- but with the same easting, northing and altitude. This is unsafe because it +-- will produce a different position unless the two grids are actually equal. +-- +-- It should be used only to convert between distinguished grids (e.g. "UkNationalGrid") and +-- their equivalent numerical definitions. +unsafeGridCoerce :: b -> GridPoint a -> GridPoint b +unsafeGridCoerce base p = GridPoint (eastings p) (northings p) (altitude p) base + + + +-- | Convert a list of digits to a distance. The first argument is the size of the +-- grid square within which these digits specify a position. The first digit is +-- in units of one tenth of the grid square, the second one hundredth, and so on. +-- The first result is the lower limit of the result, and the second is the size +-- of the specified offset. +-- So for instance @fromGridDigits (100 *~ kilo meter) "237"@ will return +-- +-- > Just (23700 meters, 100 meters) +-- +-- If there are any non-digits in the string then the function returns @Nothing@. +fromGridDigits :: Length Double -> String -> Maybe (Length Double, Length Double) +fromGridDigits sq ds = if all isDigit ds then Just (d, p) else Nothing + where + n = length ds + d = sum $ zipWith (*) + (map ((*~ one) . fromIntegral . digitToInt) ds) + (tail $ iterate (/ (10 *~ one)) sq) + p = sq / ((10 *~ one) ** (fromIntegral n *~ one)) + +-- | Convert a distance into a digit string suitable for printing as part +-- of a grid reference. The result is the nearest position to the specified +-- number of digits, expressed as an integer count of squares and a string of digits. +-- If any arguments are invalid then @Nothing@ is returned. +toGridDigits :: + Length Double -- ^ Size of enclosing grid square. Must be at least 1 km. + -> Int -- ^ Number of digits to return. Must be positive. + -> Length Double -- ^ Offset to convert into grid. + -> Maybe (Integer, String) +toGridDigits sq n d = + if sq < (1 *~ kilo meter) || n < 0 || d < _0 + then Nothing + else + Just (sqs, pad) + where + p :: Integer + p = 10 P.^ n + unit :: Length Double + unit = sq / (fromIntegral p *~ one) + u = round ((d / unit) /~ one) + (sqs, d1) = u `divMod` p + s = show d1 + pad = if n == 0 then "" else replicate (n P.- length s) '0' ++ s + diff --git a/src/Geodetics/LatLongParser.hs b/src/Geodetics/LatLongParser.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/LatLongParser.hs @@ -0,0 +1,194 @@ +-- | The default reader for Geodetic ground positions is flexible but slow. If you are +-- going to read positions in a known format and performance matters then use one of +-- the more specialised parsers here. + +module Geodetics.LatLongParser ( + + degreesMinutesSeconds, + degreesMinutesSecondsUnits, + degreesDecimalMinutes, + degreesDecimalMinutesUnits, + dms7, + angle, + latitudeNS, + longitudeEW, + signedLatLong, + latLong +) where + +import Control.Applicative +import Control.Monad +import Data.Char +import Text.ParserCombinators.ReadP as P + + + +-- | Parse an unsigned Integer value. +natural :: ReadP Integer -- Beware arithmetic overflow of Int +natural = read <$> munch1 isDigit + +-- | Parse an unsigned decimal value with optional decimal places but no exponent. +decimal :: ReadP Double +decimal = do + str1 <- munch1 isDigit + option (read str1) $ do + str2 <- char '.' *> munch1 isDigit + return $ read $ str1 ++ '.' : str2 + + +-- | Read a character indicating the sign of a value. Returns either +1 or -1. +signChar :: (Num a) => + Char -- ^ Positive sign + -> Char -- ^ Negative sign + -> ReadP a +signChar pos neg = do + c <- char pos +++ char neg + return $ if c == pos then 1 else (-1) + + +-- | Parse a signed decimal value. +signedDecimal :: ReadP Double +signedDecimal = (*) <$> option 1 (signChar '+' '-') <*> decimal + +-- | Parse an unsigned angle written using degrees, minutes and seconds separated by spaces. +-- All except the last must be integers. +degreesMinutesSeconds :: ReadP Double +degreesMinutesSeconds = do + d <- fromIntegral <$> natural + guard $ d <= 360 + skipSpaces + ms <- option 0 $ do + m <- fromIntegral <$> natural + guard $ m < 60 + skipSpaces + s <- option 0 decimal + guard $ s < 60 + return $ m / 60 + s / 3600 + return $ d + ms + + +-- | Parse an unsigned angle written using degrees, minutes and seconds with units (° ' \"). At least one component must be specified. +degreesMinutesSecondsUnits :: ReadP Double +degreesMinutesSecondsUnits = do + (s, a) <- gather $ do + d <- fromIntegral <$> option 0 (natural <* char '°') + guard $ d <= 360 + skipSpaces + m <- fromIntegral <$> option 0 (natural <* char '\'') + guard $ m < 60 + skipSpaces + s <- option 0 (decimal <* char '"') + guard $ s < 60 + return $ d + m / 60 + s / 3600 + guard $ not $ null s -- Must specify at least one component. + return a + + +-- | Parse an unsigned angle written using degrees and decimal minutes. +degreesDecimalMinutes :: ReadP Double +degreesDecimalMinutes = do + d <- fromIntegral <$> natural + skipSpaces + guard $ d <= 360 -- Difference from degreesMinutesSeconds just to shut style checker up. + m <- option 0 decimal + guard $ m < 60 + return $ d + m/60 + + +-- | Parse an unsigned angle written using degrees and decimal minutes with units (° ') +degreesDecimalMinutesUnits :: ReadP Double +degreesDecimalMinutesUnits = do + (s, a) <- gather $ do + d <- fromIntegral <$> option 0 (natural <* char '°') + guard $ d <= 360 + m <- option 0 (decimal <* char '\'') + guard $ m < 60 + return $ d + m / 60 + guard $ not $ null s -- Must specify at least one component. + return a + + +-- | Parse an unsigned angle written in DDDMMSS.ss format. +-- Leading zeros on the degrees and decimal places on the seconds are optional +dms7 :: ReadP Double +dms7 = do + str <- munch1 isDigit + decs <- option "0" (char '.' *> munch1 isDigit) + let c = length str + (ds, rs) = splitAt (c-4) str + (ms,ss) = splitAt 2 rs + d = read ds + m = read ms + s = read $ ss ++ '.' : decs + guard $ c >= 5 && c <= 7 + guard $ m < 60 + guard $ s < 60 + return $ d + m / 60 + s / 3600 + + +-- | Parse an unsigned angle, either in decimal degrees or in degrees, minutes and seconds. +-- In the latter case the unit indicators are optional. +angle :: ReadP Double +angle = choice [ + decimal, + degreesMinutesSeconds, + degreesMinutesSecondsUnits, + degreesDecimalMinutes, + degreesDecimalMinutesUnits, + dms7 + ] + + +-- | Parse latitude as an unsigned angle followed by 'N' or 'S' +latitudeNS :: ReadP Double +latitudeNS = do + ul <- angle + guard $ ul <= 90 + skipSpaces + sgn <- signChar 'N' 'S' + return $ sgn * ul + + +-- | Parse longitude as an unsigned angle followed by 'E' or 'W'. +longitudeEW :: ReadP Double +longitudeEW = do + ul <- angle + guard $ ul <= 180 + skipSpaces + sgn <- signChar 'E' 'W' + return $ sgn * ul + + +-- | Parse latitude and longitude as two signed decimal numbers in that order, optionally separated by a comma. +-- Longitudes in the western hemisphere may be represented either by negative angles down to -180 +-- or by positive angles less than 360. +signedLatLong :: ReadP (Double, Double) +signedLatLong = do + lat <- signedDecimal + guard $ lat >= (-90) + guard $ lat <= 90 + skipSpaces + P.optional $ char ',' >> skipSpaces + long <- signedDecimal + guard $ long >= (-180) + guard $ long < 360 + return (lat, if long > 180 then long-180 else long) + + +-- | Parse latitude and longitude in any format. +latLong :: ReadP (Double, Double) +latLong = latLong1 +++ longLat +++ signedLatLong + where + latLong1 = do + lat <- latitudeNS + skipSpaces + P.optional $ char ',' >> skipSpaces + long <- longitudeEW + return (lat, long) + longLat = do + long <- longitudeEW + skipSpaces + P.optional $ char ',' >> skipSpaces + lat <- latitudeNS + return (lat, long) + diff --git a/src/Geodetics/Path.hs b/src/Geodetics/Path.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Path.hs @@ -0,0 +1,271 @@ +-- | The implementation assumes IEEE 754 arithmetic. + +module Geodetics.Path where + +import Control.Monad +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + + +-- | Lower and upper exclusive bounds within which a path is valid. +type PathValidity = (Length Double, Length Double) + +-- | A path is a parametric function of distance along the path. The result is the +-- position, and the direction of the path at that point as heading and elevation angles. +-- +-- A well-behaved path must be continuous and monotonic (that is, if the distance increases +-- then the result is further along the path) for all distances within its validity range. +-- Ideally the physical distance along the path from the origin to the result should be equal +-- to the argument. If this is not practical then a reasonably close approximation may be used, +-- such as a spheroid instead of an ellipsoid, provided that this is documented. +-- Outside its validity the path function may +-- return anything or bottom. +data Path e = Path { + pathFunc :: Length Double -> (Geodetic e, Angle Double, Angle Double), + pathValidity :: PathValidity + } + +-- | Convenience value for paths that are valid for all distances. +alwaysValid :: PathValidity +alwaysValid = (negate inf, inf) where + inf = (1.0 *~ meter) / (0 *~ one) -- Assumes IEE arithmetic. + + +-- | True if the path is valid at that distance. +pathValidAt :: Path e -> Length Double -> Bool +pathValidAt path d = d > x1 && d < x2 + where (x1,x2) = pathValidity path + + +-- | Find where a path meets a given condition using bisection. This is not the +-- most efficient algorithm, but it will always produce a result if there is one +-- within the initial bounds. If there is more than one result then an arbitrary +-- one will be returned. +-- +-- The initial bounds must return one GT or EQ value and one LT or EQ value. If they +-- do not then @Nothing@ is returned. +bisect :: + Path e + -> (Geodetic e -> Ordering) -- ^ Evaluation function. + -> Length Double -- ^ Required accuracy in terms of distance along the path. + -> Length Double -> Length Double -- ^ Initial bounds. + -> Maybe (Length Double) +bisect path f t b1 b2 = do + guard $ pathValidAt path b1 + guard $ pathValidAt path b2 + let r = pairResults b1 b2 + guard $ hasRoot r + bisect1 $ sortPair r + where + f' d = let (p, _, _) = pathFunc path d in f p + pairResults d1 d2 = ((d1, f' d1), (d2, f' d2)) + hasRoot (v1, v2) = snd v1 <= EQ && EQ <= snd v2 + sortPair (v1, v2) = if snd v1 <= snd v2 then (v1, v2) else (v2, v1) + bisect1 ((d1, r1), (d2, r2)) = + let d3 = (d1 + d2) / _2 + r3 = f' d3 + c1 = ((d1, r1), (d3, r3)) + c2 = ((d3, r3), (d2, r2)) + in if abs (d1 - d2) <= t + then return d3 + else bisect1 $ if hasRoot c1 then c1 else c2 + + +-- | Try to find the intersection point of two ground paths (i.e. ignoring +-- altitude). Returns the distance of the intersection point along each path +-- using a modified Newton-Raphson method. If the two paths are +-- fairly straight and not close to parallel then this will converge rapidly. +-- +-- The algorithm projects great-circle paths forwards using the bearing at the +-- estimate to find the estimated intersection, and then uses the distances to this +-- intersection as the next estimates. +-- +-- If either estimate departs from its path validity then @Nothing@ is returned. +intersect :: (Ellipsoid e) => + Length Double -> Length Double -- ^ Starting estimates. + -> Length Double -- ^ Required accuracy. + -> Int -- ^ Iteration limit. Returns @Nothing@ if this is reached. + -> Path e -> Path e -- ^ Paths to intersect. + -> Maybe (Length Double, Length Double) +intersect d1 d2 accuracy n path1 path2 + | not $ pathValidAt path1 d1 = Nothing + | not $ pathValidAt path2 d2 = Nothing + | n <= 0 = Nothing + | mag < (1e-15 *~ one) = Nothing + | mag3 (nv1 `cross3` nv2) * r <= accuracy = Just (d1, d2) + -- Assumes that sin (accuracy/r) == accuracy/r + | otherwise = + if abs d1a + abs d2a < abs d1b + abs d2b + then intersect (d1 + d1a) (d2 + d2a) accuracy (pred n) path1 path2 + else intersect (d1 + d1b) (d2 + d2b) accuracy (pred n) path1 path2 + where + (pt1, h1, _) = pathFunc path1 d1 + (pt2, h2, _) = pathFunc path2 d2 + vectors :: Angle Double -> Angle Double -> Angle Double + -> (Vec3 (Dimensionless Double), Vec3 (Dimensionless Double)) + vectors lat lon b = ( + -- Unit vector of normal to surface at (lat,lon) + (cosLat*cosLon, cosLat*sinLon, sinLat), + -- Normal of great circle defined by bearing b at (lat,lon) + (sinLon * cosB - sinLat * cosLon * sinB, + negate cosLon * cosB - sinLat * sinLon * sinB, + cosLat * sinB)) + where + sinLon = sin lon + sinLat = sin lat + cosLon = cos lon + cosLat = cos lat + sinB = sin b + cosB = cos b + mag3 (x,y,z) = sqrt $ x*x + y*y + z*z + (nv1, gc1) = vectors (latitude pt1) (longitude pt1) h1 + (nv2, gc2) = vectors (latitude pt2) (longitude pt2) h2 + nv3 = gc1 `cross3` gc2 -- Intersection of the great circles + mag = mag3 nv3 + nv3a = scale3 nv3 (_1 / mag) -- Scale to unit. See outer function for case when mag3 == 0 + nv3b = negate3 nv3a -- Antipodal result. Take the closest. + -- Find "nearest" intersection, defined as smaller of sum of distances to current points. + d1a = gcDist gc1 nv1 nv3a * r + d2a = gcDist gc2 nv2 nv3a * r + d1b = gcDist gc1 nv1 nv3b * r + d2b = gcDist gc2 nv2 nv3b * r + -- Signed angle between v1 and v2, + gcDist norm v1 v2 = + let c = v1 `cross3` v2 + in (if c `dot3` norm < _0 then negate else id) $ atan2 (mag3 c) (v1 `dot3` v2) + r = majorRadius $ ellipsoid pt1 + +{- Note on derivation + +The algorithm is a variant of the Newton-Raphson method, and shares its advantage +of rapid convergence in many useful cases. Each path has a current approximation to +the intersection, and the next approximation is computed by projecting both paths +along great circles from the current approximation and finding the point where those +great circles intersect. A spherical Earth is assumed for simplicity. + +The Great Circle calculations use a vector method rather than spherical trigonometry. +This avoids a lot of transcendental functions and also the singularities inherent in +polar coordinate systems. This implementation is based on formulae from +http://www.movable-type.co.uk/scripts/latlong-vectors.html, which in turn is based on +"A Non-singular Horizontal Position Representation" by Kenneth Gade, THE JOURNAL OF +NAVIGATION (2010), 63, 395–417. +http://www.navlab.net/Publications/A_Nonsingular_Horizontal_Position_Representation.pdf + + +"pt1" is the current approximation for the result on "path1". The vector "nv1" is the +unit vector pointing "pt1". "gc1" is the normal to the great circle tangent to +"path1" at "pt1". "nv2" and "gc2" are similarly derived from "path2". + +The intersection of the planes of "gc1" and "gc2" is given by their cross product, "nv3". +This is scaled to a unit vector "nv3a", and the other intersection of the great circles is +at "nv3b", opposite. The distances from nv1 and nv2 to nv3a and nv3b are computed using +the corresponding great-circle normals to determine whether the distances are positive or +negative along the paths. The nearest solution (defined in terms of great-circle distance) +is taken as the basis for the next approximation. +-} + +-- | A ray from a point heading in a straight line in 3 dimensions. +rayPath :: (Ellipsoid e) => + Geodetic e -- ^ Start point. + -> Angle Double -- ^ Bearing. + -> Angle Double -- ^ Elevation. + -> Path e +rayPath pt1 bearing elevation = Path ray alwaysValid + where + ray distance = (Geodetic lat long alt (ellipsoid pt1), bearing2, elevation2) + where + pt2' = pt1' `add3` (delta `scale3` distance) -- ECEF of result point. + (lat,long,alt) = earthToGeo (ellipsoid pt1) pt2' -- Geodetic of result point. + (dE,dN,dU) = transform3 (trans3 $ ecefMatrix lat long) delta -- Direction of ray at result point. + elevation2 = asin dU + bearing2 = if dE == _0 && dN == _0 then bearing else atan2 dE dN -- Allow for vertical elevation. + + ecefMatrix lat long = -- Transform matrix for vectors from (East, North, Up) to (X,Y,Z). + ((negate sinLong, negate cosLong*sinLat, cosLong*cosLat), + -- East X North X Up X + ( cosLong, negate sinLong*sinLat, sinLong*cosLat), + -- East Y North Y Up Y + ( _0 , cosLat , sinLat)) + -- East Z North Z Up Z + where + sinLong = sin long + cosLong = cos long + sinLat = sin lat + cosLat = cos lat + + direction = (sinB*cosE, cosB*cosE, sinE) -- Direction of ray in ENU + delta = transform3 (ecefMatrix (latitude pt1) (longitude pt1)) direction -- Convert to ECEF + pt1' = geoToEarth pt1 -- ECEF of origin point. + sinB = sin bearing + cosB = cos bearing + sinE = sin elevation + cosE = cos elevation + +-- | Rhumb line: path following a constant course. Also known as a loxodrome. +-- +-- The valid range stops a few arc-minutes short of the poles to ensure that the +-- polar singularities are not included. +-- +-- Based on *Practical Sailing Formulas for Rhumb-Line Tracks on an Oblate Earth* +-- by G.H. Kaplan, U.S. Naval Observatory. Except for points close to the poles +-- the approximation is accurate to within a few meters over 1000km. +rhumbPath :: (Ellipsoid e) => + Geodetic e -- ^ Start point. + -> Angle Double -- ^ Course. + -> Path e +rhumbPath pt course = Path rhumb validity + where + rhumb distance = (Geodetic lat (properAngle lon) _0 (ellipsoid pt), course, _0) + where + lat' = lat0 + distance * cosC / m0 -- Kaplan Eq 13. + lat = lat0 + (m0 / (a*(_1-e2))) * ((_1-_3*e2/_4)*(lat'-lat0) + + (_3*e2/_8)*(sin (_2*lat') - sin (_2*lat0))) + lon | abs cosC > 1e-7 *~ one + = lon0 + tanC * (q lat - q0) -- Kaplan Eq 16. + | otherwise + = lon0 + distance * sinC / latitudeRadius (ellipsoid pt) ((lat0 + lat')/_2) + validity + | cosC > _0 = ((negate pi/_2 - latitude pt) * b / cosC, (pi/_2 - latitude pt) * b / cosC) + | otherwise = ((pi/_2 - latitude pt) * b / cosC, (negate pi/_2 - latitude pt) * b / cosC) + q0 = q lat0 + q phi = log (tan (pi/_4+phi/_2)) + e * log ((_1-eSinPhi)/(_1+eSinPhi)) / _2 + where -- Factor out expression from Eq 16 of Kaplan + eSinPhi = e * sin phi + sinC = sin course + cosC = cos course + tanC = tan course + lat0 = latitude pt + lon0 = longitude pt + e2 = eccentricity2 $ ellipsoid pt + e = sqrt e2 + m0 = meridianRadius (ellipsoid pt) lat0 + a = majorRadius $ ellipsoid pt + b = minorRadius $ ellipsoid pt + + +-- | A path following the line of latitude around the Earth eastwards. +-- +-- This is equivalent to @rhumbPath pt (pi/2)@ +latitudePath :: (Ellipsoid e) => + Geodetic e -- ^ Start point. + -> Path e +latitudePath pt = Path line alwaysValid + where + line distance = (pt2, pi/_2, _0) + where + pt2 = Geodetic + (latitude pt) (longitude pt + distance / r) + _0 (ellipsoid pt) + r = latitudeRadius (ellipsoid pt) (latitude pt) + + +-- | A path from the specified point to the North Pole. Use negative distances +-- for the southward path. +-- +-- This is equivalent to @rhumbPath pt _0@ +longitudePath :: (Ellipsoid e) => + Geodetic e -- ^ Start point. + -> Path e +longitudePath pt = rhumbPath pt _0 diff --git a/src/Geodetics/Stereographic.hs b/src/Geodetics/Stereographic.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/Stereographic.hs @@ -0,0 +1,121 @@ +{- | +The following is based on equations in Section 1.4.7.1 in +OGP Surveying and Positioning Guidance Note number 7, part 2 – August 2006 + http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf +-} + +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} +module Geodetics.Stereographic ( + GridStereo (gridTangent, gridOrigin, gridScale), + mkGridStereo +) where + + +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Geodetics.Grid +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + + +-- | A stereographic projection with its origin at an arbitrary point on Earth, other than the poles. +data GridStereo e = GridStereo { + gridTangent :: Geodetic e, -- ^ Point where the plane of projection touches the ellipsoid. Often known as the Natural Origin. + gridOrigin :: GridOffset, -- ^ Grid position of the tangent point. Often known as the False Origin. + gridScale :: Dimensionless Double, -- ^ Scaling factor that balances the distortion between the center and the edges. + -- Should be slightly less than unity. + + -- Memoised parameters derived from the tangent point. + gridR :: Length Double, + gridN, gridC, gridSin, gridCos :: Dimensionless Double, + gridLatC :: Angle Double, + gridG, gridH :: Length Double + } deriving (Show) + +-- | Create a stereographic projection. The tangency point must not be one of the poles. +mkGridStereo :: (Ellipsoid e) => Geodetic e -> GridOffset -> Dimensionless Double -> GridStereo e +mkGridStereo tangent origin scale = GridStereo { + gridTangent = tangent, + gridOrigin = origin, + gridScale = scale, + gridR = r, + gridN = n, + gridC = c, + gridSin = sinLatC1, + gridCos = sqrt $ _1 - sinLatC1 * sinLatC1, + gridLatC = asin sinLatC1, + gridG = g, + gridH = h + } + where + -- The reference seems to use χO to refer to two slightly different values. + -- Here these will be called LatC0 and LatC1. + ellipse = ellipsoid tangent + op :: Num a => Quantity d a -> Quantity d a -- Values of longitude, tangent longitude, E and N + op = if latitude tangent < _0 then negate else id -- must be negated in the southern hemisphere. + lat0 = op $ latitude tangent + sinLat0 = sin lat0 + e2 = eccentricity2 ellipse + e = sqrt e2 + r = sqrt $ meridianRadius ellipse lat0 * primeVerticalRadius ellipse lat0 + n = sqrt $ _1 + ((e2 * cos lat0 ^ pos4)/(_1 - e2)) + s1 = (_1 + sinLat0) / (_1 - sinLat0) + s2 = (_1 - e * sinLat0) / (_1 + e * sinLat0) + w1 = (s1 * s2 ** e) ** n + sinLatC0 = (w1 - _1)/(w1 + _1) + c = ((n + sin lat0) * (_1 - sinLatC0)) / ((n - sin lat0) * (_1 + sinLatC0)) + w2 = c * w1 + sinLatC1 = (w2 - _1)/(w2 + _1) + g = _2 * r * scale * tan (pi/_4 - latC1/_2) + h = _4 * r * scale * tan latC1 + g + latC1 = asin sinLatC1 + + +instance (Ellipsoid e) => GridClass (GridStereo e) e where + toGrid grid geo = applyOffset (gridOrigin grid) $ GridPoint east north (geoAlt geo) grid + where + op :: Num a => Quantity d a -> Quantity d a -- Values of longitude, tangent longitude, E and N + op = if latitude (gridTangent grid) < _0 then negate else id -- must be negated in the southern hemisphere. + sinLatC = (w - _1)/(w + _1) + cosLatC = sqrt $ _1 - sinLatC * sinLatC + longC = gridN grid * (op (longitude geo) - long0) + long0 + w = gridC grid * (sA * sB ** e) ** gridN grid + sA = (_1+sinLat) / (_1 - sinLat) + sB = (_1 - e*sinLat) / (_1 + e*sinLat) + sinLat = sin $ op $ latitude geo + e = sqrt $ eccentricity2 $ ellipsoid geo + long0 = op $ longitude $ gridTangent grid + b = _1 + sinLatC * gridSin grid + cosLatC * gridCos grid * cos (longC - long0) + east = _2 * gridR grid * gridScale grid * cosLatC * sin (longC - long0) / b + north = _2 * gridR grid * gridScale grid * (sinLatC * gridCos grid - cosLatC * gridSin grid * cos (longC - long0)) / b + + fromGrid gp = + {- trace ( -- Remove comment brackets for debugging. + "fromGrid values:\n i = " ++ show i ++ "\n j = " ++ show j ++ + "\n longC = " ++ show longC ++ "\n long = " ++ show long ++ + "\n latC = " ++ show latC ++ + "\n lat1 = " ++ show lat1 ++ "\n latN = " ++ show latN ) $ -} + Geodetic (op latN) (op long) height $ gridEllipsoid grid + where + op :: Num a => Quantity d a -> Quantity d a -- Values of longitude, tangent longitude, E and N + op = if latitude (gridTangent grid) < _0 then negate else id -- must be negated in the southern hemisphere. + GridPoint east north height _ = applyOffset (offsetNegate $ gridOrigin grid) gp + east' = east + north' = north + grid = gridBasis gp + long0 = op $ longitude $ gridTangent grid + i = atan2 east' (gridH grid + north') + j = atan2 east' (gridG grid - north') - i + latC = gridLatC grid + _2 * atan2 (north' - east' * tan (j/_2)) (_2 * gridR grid * gridScale grid) + longC = j + _2 * i + long0 + sinLatC = sin latC + long = (longC - long0) / gridN grid + long0 + isoLat = log ((_1 + sinLatC) / (gridC grid * (_1 - sinLatC))) / (_2 * gridN grid) + lat1 = _2 * atan (exp isoLat) - pi/_2 + next lat = lat - (isoN - isoLat) * cos lat * (_1 - e2 * sin lat ^ pos2) / (_1 - e2) + where isoN = isometricLatitude (gridEllipsoid grid) lat + e2 = eccentricity2 $ gridEllipsoid grid + lats = iterate next lat1 + latN = snd $ head $ dropWhile (\(v1, v2) -> abs (v1-v2) > 0.01 *~ arcsecond) $ zip lats $ tail lats + + gridEllipsoid = ellipsoid . gridTangent diff --git a/src/Geodetics/TransverseMercator.hs b/src/Geodetics/TransverseMercator.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/TransverseMercator.hs @@ -0,0 +1,158 @@ +{-# LANGUAGE MultiParamTypeClasses, FlexibleInstances #-} + +module Geodetics.TransverseMercator( + GridTM (trueOrigin, falseOrigin, gridScale), + mkGridTM +) where + +import Data.Function +import Data.Monoid +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Geodetics.Grid +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + +-- | A Transverse Mercator projection gives an approximate mapping of the ellipsoid on to a 2-D grid. It models +-- a sheet curved around the ellipsoid so that it touches it at one north-south line (hence making it part of +-- a slightly elliptical cylinder). +data GridTM e = GridTM { + trueOrigin :: Geodetic e, + -- ^ A point on the line where the projection touches the ellipsoid (altitude is ignored). + falseOrigin :: GridOffset, + -- ^ The grid position of the true origin. Used to avoid negative coordinates over + -- the area of interest. The altitude gives a vertical offset from the ellipsoid. + gridScale :: Dimensionless Double, + -- ^ A scaling factor that balances the distortion between the east & west edges and the middle + -- of the projection. + + -- Remaining elements are memoised parameters computed from the ellipsoid underlying the true origin. + gridN1, gridN2, gridN3, gridN4 :: Dimensionless Double +} deriving (Show) + + +-- | Create a Transverse Mercator grid. +mkGridTM :: (Ellipsoid e) => + Geodetic e -- ^ True origin. + -> GridOffset -- ^ Vector from true origin to false origin. + -> Dimensionless Double -- ^ Scale factor. + -> GridTM e +mkGridTM origin offset sf = + GridTM {trueOrigin = origin, + falseOrigin = offset, + gridScale = sf, + gridN1 = _1 + n + (_5/_4) * n^pos2 + (_5/_4) * n^pos3, + gridN2 = _3 * n + _3 * n^pos2 + ((21*~one)/_8) * n^pos3, + gridN3 = ((15*~one)/_8) * (n^pos2 + n^pos3), + gridN4 = ((35*~one)/(24*~one)) * n^pos3 + } + where + f = flattening $ ellipsoid origin + n = f / (_2-f) -- Equivalent to (a-b)/(a+b) where b = (1-f)*a + + + + +-- | Equation C3 from reference [1]. +m :: (Ellipsoid e) => GridTM e -> Dimensionless Double -> Length Double +m grid lat = bF0 * (gridN1 grid * dLat + - gridN2 grid * sin dLat * cos sLat + + gridN3 grid * sin (_2 * dLat) * cos (_2 * sLat) + - gridN4 grid * sin (_3 * dLat) * cos (_3 * sLat)) + where + dLat = lat - latitude (trueOrigin grid) + sLat = lat + latitude (trueOrigin grid) + bF0 = minorRadius (gridEllipsoid grid) * gridScale grid + + +instance (Ellipsoid e) => GridClass (GridTM e) e where + fromGrid p = Geodetic + (lat' - east' ^ pos2 * tanLat / (_2 * rho * v) -- Term VII + + east' ^ pos4 * (tanLat / ((24 *~ one) * rho * v ^ pos3)) + * (_5 + _3 * tanLat ^ pos2 + eta2 - _9 * tanLat ^ pos2 * eta2) -- Term VIII + - east' * east' ^ pos5 * (tanLat / ((720 *~ one) * rho * v ^ pos5)) + * (61 *~ one + (90 *~ one) * tanLat ^ pos2 + (45 *~ one) * tanLat ^ pos4)) -- Term IX + (longitude (trueOrigin grid) + + east' / (cosLat * v) -- Term X + - (east' ^ pos3 / (_6 * cosLat * v ^ pos3)) * (v / rho + _2 * tanLat ^ pos2) -- Term XI + + (east' ^ pos5 / ((120 *~ one) * cosLat * v ^ pos5)) + * (_5 + (28 *~ one) * tanLat ^ pos2 + (24 *~ one) * tanLat ^ pos4) -- Term XII + - (east' ^ pos5 * east' ^ pos2 / ((5040 *~ one) * cosLat * v * v * v ^ pos5)) + * ((61 *~ one) + (662 *~ one) * tanLat ^ pos2 + (1320 *~ one) * tanLat ^ pos4 + (720 *~ one) * tanLat * tanLat ^ pos5)) -- Term XIIa + (0 *~ meter) (gridEllipsoid grid) + + + where + GridPoint east' north' _ _ = (falseOrigin grid) `applyOffset` p + lat' = fst $ head $ dropWhile ((> 0.01 *~ milli meter) . snd) + $ tail $ iterate next (latitude $ trueOrigin grid, 1 *~ meter) + where + next (phi, _) = let delta = north' - m grid phi in (phi + delta / aF0, delta) + -- head and tail are safe because iterate returns an infinite list. + + sinLat = sin lat' + cosLat = cos lat' + tanLat = tan lat' + sinLat2 = sinLat ^ pos2 + v = aF0 / sqrt (_1 - e2 * sinLat2) + rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one) + eta2 = v / rho - _1 + + + aF0 = majorRadius (gridEllipsoid grid) * gridScale grid + e2 = eccentricity2 $ gridEllipsoid grid + grid = gridBasis p + + toGrid grid geo = applyOffset (off `mappend` (offsetNegate $ falseOrigin grid)) $ + GridPoint (0 *~ metre) (0 *~ metre) (0 *~ metre) grid + where + v = aF0 / sqrt (_1 - e2 * sinLat2) + rho = aF0 * (_1 - e2) * (_1 - e2 * sinLat2) ** ((-1.5) *~ one) + eta2 = v / rho - _1 + off = GridOffset + (dLong * term_IV + + dLong ^ pos3 * term_V + + dLong ^ pos5 * term_VI) + (m grid lat + dLong ^ pos2 * term_II + + dLong ^ pos4 * term_III + + dLong * dLong ^ pos5 * term_IIIa) + (0 *~ meter) + -- Terms defined in [1]. + term_II = (v/_2) * sinLat * cosLat + term_III = (v/(24*~one)) * sinLat * cosLat ^ pos3 + * (_5 - tanLat ^ pos2 + _9 * eta2) + term_IIIa = (v/(720*~one)) * sinLat * cosLat ^ pos5 + * ((61 *~ one) - (58 *~ one) * tanLat ^ pos2 + tanLat ^ pos4) + term_IV = v * cosLat + term_V = (v/_6) * cosLat ^ pos3 * (v/rho - tanLat ^ pos2) + term_VI = (v/(120*~one)) * cosLat ^ pos5 + * (_5 - (18*~one) * tanLat ^ pos2 + + tanLat ^ pos4 + (14*~one) * eta2 + - (58*~one) * tanLat ^ pos2 * eta2) + {- + -- Trace message for debugging. Uncomment this code for easy access to intermediate values. + traceMsg = concat [ + "v = ", show v, "\n", + "rho = ", show rho, "\n", + "eta2 = ", show eta2, "\n", + "M = ", show $ m grid lat, "\n", + "I = ", show $ m grid lat + deltaNorth (falseOrigin grid), "\n", + "II = ", show term_II, "\n", + "III = ", show term_III, "\n", + "IIIa = ", show term_IIIa, "\n", + "IV = ", show term_IV, "\n", + "V = ", show term_V, "\n", + "VI = ", show term_VI, "\n"] + -} + -- Common subexpressions + lat = latitude geo + long = longitude geo + dLong = long - longitude (trueOrigin grid) + sinLat = sin lat + cosLat = cos lat + tanLat = tan lat + sinLat2 = sinLat ^ pos2 + aF0 = (majorRadius $ gridEllipsoid grid) * gridScale grid + e2 = eccentricity2 $ gridEllipsoid grid + + gridEllipsoid = ellipsoid . trueOrigin diff --git a/src/Geodetics/UK.hs b/src/Geodetics/UK.hs new file mode 100644 --- /dev/null +++ b/src/Geodetics/UK.hs @@ -0,0 +1,145 @@ +{-# LANGUAGE MultiParamTypeClasses #-} + +-- | Distinguished coordinate systems for the United Kingdom. +module Geodetics.UK ( + OSGB36 (..), + UkNationalGrid (..), + ukGrid, + fromUkGridReference, + toUkGridReference +) where + +import Control.Applicative +import Control.Monad +import Data.Array +import Data.Char +import Data.Monoid +import Geodetics.Geodetic +import Geodetics.Grid +import Geodetics.Ellipsoids +import Geodetics.TransverseMercator +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P + + + +-- | Ellipsoid definition for Great Britain. Airy 1830 offset from the centre of the Earth +-- and rotated slightly. +-- +-- The Helmert parameters are from the Ordnance Survey document +-- \"A Guide to Coordinate Systems in Great Britain\", which notes that it +-- can be in error by as much as 5 meters and should not be used in applications +-- requiring greater accuracy. A more precise conversion requires a large table +-- of corrections for historical inaccuracies in the triangulation of the UK. +data OSGB36 = OSGB36 deriving (Eq, Show) + +instance Ellipsoid OSGB36 where + majorRadius _ = 6377563.396 *~ meter + flatR _ = 299.3249646 *~ one + helmert _ = Helmert { + cX = 446.448 *~ meter, cY = (-125.157) *~ meter, cZ = 542.06 *~ meter, + helmertScale = (-20.4894) *~ one, + rX = 0.1502 *~ arcsecond, rY = 0.247 *~ arcsecond, rZ = 0.8421 *~ arcsecond } + + +-- | The UK National Grid is a Transverse Mercator projection with a true origin at +-- 49 degrees North, 2 degrees West on OSGB36, and a false origin 400km West and 100 km North of +-- the true origin. The scale factor is defined as @10**(0.9998268 - 1)@. +data UkNationalGrid = UkNationalGrid deriving (Eq, Show) + +instance GridClass UkNationalGrid OSGB36 where + toGrid _ = unsafeGridCoerce UkNationalGrid . toGrid ukGrid + fromGrid = fromGrid . unsafeGridCoerce ukGrid + gridEllipsoid _ = OSGB36 + + + +ukTrueOrigin :: Geodetic OSGB36 +ukTrueOrigin = Geodetic { + latitude = 49 *~ degree, + longitude = (-2) *~ degree, + geoAlt = 0 *~ meter, + ellipsoid = OSGB36 +} + +ukFalseOrigin :: GridOffset +ukFalseOrigin = GridOffset ((-400) *~ kilo meter) (100 *~ kilo meter) (0 *~ meter) + + +-- | Numerical definition of the UK national grid. +ukGrid :: GridTM OSGB36 +ukGrid = mkGridTM ukTrueOrigin ukFalseOrigin + ((10 *~ one) ** (0.9998268 *~ one - _1)) + + +-- | Size of a UK letter-pair grid square. +ukGridSquare :: Length Double +ukGridSquare = 100 *~ kilo meter + + +-- | Convert a grid reference to a position, if the reference is valid. +-- This returns the position of the south-west corner of the nominated +-- grid square and an offset to its centre. Altitude is set to zero. +fromUkGridReference :: String -> Maybe (GridPoint UkNationalGrid, GridOffset) +fromUkGridReference str = if length str < 2 then Nothing else do + let + c1:c2:ds = str + n = length ds + guard $ even n + let (dsE, dsN) = splitAt (n `div` 2) ds + (east, sq) <- fromGridDigits ukGridSquare dsE + (north, _) <- fromGridDigits ukGridSquare dsN + base <- fromUkGridLetters c1 c2 + let half = sq / (2 *~ one) + return (applyOffset (GridOffset east north (0 *~ meter)) base, + GridOffset half half (0 *~ meter)) + + + + +-- | The south west corner of the nominated grid square, if it is a legal square. +-- This function works for all pairs of letters except 'I' (as that is not used). +-- In practice only those pairs covering the UK are actually considered meaningful. +fromUkGridLetters :: Char -> Char -> Maybe (GridPoint UkNationalGrid) +fromUkGridLetters c1 c2 = applyOffset <$> (mappend <$> g1 <*> g2) <*> letterOrigin + where + letterOrigin = Just $ GridPoint ((-1000) *~ kilo meter) ((-500) *~ kilo meter) m0 UkNationalGrid + gridIndex c = + if inRange ('A', 'H') c then Just $ ord c P.- ord 'A' -- 'I' is not used. + else if inRange ('J', 'Z') c then Just $ ord c P.- ord 'B' + else Nothing + gridSquare c = do -- Maybe monad + g <- gridIndex c + let (y,x) = g `divMod` 5 + return (fromIntegral x *~ one, _4 - fromIntegral y *~ one) + g1 = do + (x,y) <- gridSquare c1 + return $ GridOffset (x * (500 *~ kilo meter)) (y * (500 *~ kilo meter)) m0 + g2 = do + (x,y) <- gridSquare c2 + return $ GridOffset (x * (100 *~ kilo meter)) (y * (100 *~ kilo meter)) m0 + m0 = 0 *~ meter + + +-- | Find the nearest UK grid reference point to a specified position. The Int argument is the number of +-- digits precision, so 2 for a 4-figure reference and 3 for a 6-figure reference, although any value +-- between 0 and 5 can be used (giving a 1 meter precision). +-- Altitude is ignored. If the result is outside the area defined by the two letter grid codes then +-- @Nothing@ is returned. +toUkGridReference :: Int -> GridPoint UkNationalGrid -> Maybe String +toUkGridReference n p + | n < 0 = error "toUkGridReference: precision argument must not be negative." + | otherwise = do + (gx, strEast) <- toGridDigits ukGridSquare n $ eastings p + 1000 *~ kilo meter + (gy, strNorth) <- toGridDigits ukGridSquare n $ northings p + 500 *~ kilo meter + let (gx1, gx2) = (fromIntegral gx) `divMod` 5 + (gy1, gy2) = (fromIntegral gy) `divMod` 5 + guard (gx1 < 5 && gy1 < 5) + let c1 = gridSquare gx1 gy1 + c2 = gridSquare gx2 gy2 + return $ c1 : c2 : strEast ++ strNorth + where + gridSquare x y = letters ! (4 P.- y, x) + letters :: Array (Int, Int) Char + letters = listArray ((0,0),(4,4)) $ ['A'..'H'] ++ ['J'..'Z'] + diff --git a/test/ArbitraryInstances.hs b/test/ArbitraryInstances.hs new file mode 100644 --- /dev/null +++ b/test/ArbitraryInstances.hs @@ -0,0 +1,254 @@ +{-# LANGUAGE FlexibleInstances #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +-- | Orphan "Arbitrary" and related instances for testing purposes. + +module ArbitraryInstances where + +import Control.Applicative +import Control.Monad +import Geodetics.Altitude +import Geodetics.Geodetic +import Geodetics.Grid +import Geodetics.Ellipsoids +import Geodetics.Path +import Geodetics.Stereographic as SG +import Geodetics.TransverseMercator as TM +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P +import Test.QuickCheck + + + +-- | Shrink using a dimension, so that shrunk values are round numbers in that dimension. +shrinkDimension :: (Fractional a, Arbitrary a) => Unit d a -> Quantity d a -> [Quantity d a] +shrinkDimension u v = (*~ u) <$> shrink (v /~ u) + +-- | Wrapper for arbitrary angles. +newtype Bearing = Bearing (Dimensionless Double) + +instance Show Bearing where + show (Bearing b) = "Bearing " ++ showAngle b + +instance Arbitrary Bearing where + arbitrary = Bearing <$> (*~ degree) <$> choose (-180,180) + shrink (Bearing b) = Bearing <$> shrinkDimension degree b + + +newtype Azimuth = Azimuth (Dimensionless Double) + +instance Show Azimuth where + show (Azimuth a) = "Azimuth " ++ showAngle a + +instance Arbitrary Azimuth where + arbitrary = Azimuth <$> (*~ degree) <$> choose (0,90) + shrink (Azimuth a) = Azimuth <$> shrinkDimension degree a + + +-- | Wrapper for arbitrary distances up to 10,000 km +newtype Distance = Distance (Length Double) deriving (Show) + +instance Arbitrary Distance where + arbitrary = Distance <$> (*~ kilo meter) <$> choose (0,10000) + shrink (Distance d) = Distance <$> shrinkDimension (kilo meter) d + + +-- | Wrapper for arbitrary distances up to 1,000 km +newtype Distance2 = Distance2 (Length Double) deriving (Show) + +instance Arbitrary Distance2 where + arbitrary = Distance2 <$> (*~ kilo meter) <$> choose (0,1000) + shrink (Distance2 d) = Distance2 <$> shrinkDimension (kilo meter) d + +-- | Wrapper for arbitrary altitudes up to 10 km +newtype Altitude = Altitude (Length Double) deriving (Show) + +instance Arbitrary Altitude where + arbitrary = Altitude <$> (*~ kilo meter) <$> choose (0,10) + shrink (Altitude h) = Altitude <$> shrinkDimension (kilo meter) h + + +-- | Wrapper for arbitrary dimensionless numbers (-10 .. 10) +newtype Scalar = Scalar (Dimensionless Double) deriving (Show) + +instance Arbitrary Scalar where + arbitrary = Scalar <$> (*~ one) <$> choose (-10,10) + shrink (Scalar s) = Scalar <$> shrinkDimension one s + + +-- | Wrapper for arbitrary grid references. +newtype GridRef = GridRef String deriving Show + +instance Arbitrary GridRef where + arbitrary = do + n <- choose (0,4) + c1 <- elements "HJNOST" -- General vicinity of UK + c2 <- elements $ ['A'..'H'] ++ ['J'..'Z'] + dx <- vectorOf n $ choose ('0','9') + dy <- vectorOf n $ choose ('0','9') + return $ GridRef $ c1 : c2 : (dx ++ dy) + shrink = shrinkNothing + + +-- | Generate in range +/- m. +genOffset :: Double -> Gen (Length Double) +genOffset d = (*~ meter) <$> choose (-d, d) + +genAlt :: Gen (Length Double) +genAlt = (*~ meter) <$> choose (0,10000) + + +genLatitude :: Gen (Dimensionless Double) +genLatitude = (*~ degree) <$> choose (-90,90) + +genLongitude :: Gen (Dimensionless Double) +genLongitude = (*~ degree) <$> choose (-180,180) + +genSeconds :: Gen (Dimensionless Double) +genSeconds = (*~ arcsecond) <$> choose (-10,10) + + +-- | Shrinking with the original value preserved. Used for shrinking records. See +-- http://stackoverflow.com/questions/14006005/idiomatic-way-to-shrink-a-record-in-quickcheck for details. +shrink' :: (Arbitrary a) => a -> [a] +shrink' x = x : shrink x + +-- | Shrink a quantity in the given units. +shrinkQuantity :: (Arbitrary a, Fractional a) => Unit d a -> Quantity d a -> [Quantity d a] +shrinkQuantity u q = map (*~ u) $ shrink' $ q /~ u + +shrinkLength :: (Arbitrary a, Fractional a) => Length a -> [Length a] +shrinkLength = shrinkQuantity meter + +shrinkUnit :: (Arbitrary a, Fractional a) => Dimensionless a -> [Dimensionless a] +shrinkUnit = shrinkQuantity one + +shrinkAngle :: (Arbitrary a, Floating a) => Dimensionless a -> [Dimensionless a] +shrinkAngle = shrinkQuantity degree + + +instance Arbitrary Helmert where + arbitrary = + Helmert <$> genOffset 300 <*> genOffset 300 <*> genOffset 300 <*> + ((*~ one) <$> choose (-5,10)) <*> + genSeconds <*> genSeconds <*> genSeconds + shrink h = + tail $ Helmert <$> shrinkLength (cX h) <*> shrinkLength (cY h) <*> shrinkLength (cZ h) <*> + shrinkUnit (helmertScale h) <*> + shrinkUnit (rX h) <*> shrinkUnit (rY h) <*> shrinkUnit (rZ h) + + +instance Arbitrary WGS84 where + arbitrary = return WGS84 + shrink = shrinkNothing + + +instance Arbitrary LocalEllipsoid where + arbitrary = + LocalEllipsoid <$> (("Local_" ++) <$> replicateM 3 (choose ('A','Z'))) <*> -- name + ((*~ meter) <$> choose (6378100, 6378400)) <*> -- majorRadius + ((*~ one) <$> choose (297,300)) <*> -- flatR + arbitrary -- helmert + shrink e = tail $ LocalEllipsoid (nameLocal e) (majorRadius e) (flatR e) <$> shrink' (helmert e) + + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (Geodetic e) where + arbitrary = + Geodetic <$> genLatitude <*> genLongitude <*> genOffset 1 <*> arbitrary + shrink g = + tail $ Geodetic <$> shrinkAngle (latitude g) <*> shrinkAngle (longitude g) <*> + shrinkLength (altitude g) <*> shrink' (ellipsoid g) + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridTM e)) where + arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary + shrink p = tail $ GridPoint <$> + shrinkLength (eastings p) <*> + shrinkLength (northings p) <*> + shrinkLength (altitude p) <*> + shrink' (gridBasis p) + + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridPoint (GridStereo e)) where + arbitrary = GridPoint <$> genOffset 100000 <*> genOffset 100000 <*> genOffset 1 <*> arbitrary + shrink p = tail $ GridPoint <$> + shrinkLength (eastings p) <*> + shrinkLength (northings p) <*> + shrinkLength (altitude p) <*> + shrink' (gridBasis p) + + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridTM e) where + arbitrary = mkGridTM <$> arbitrary <*> arbitrary <*> ((*~ one) <$> choose (0.95,1.0)) + shrink tm = tail $ mkGridTM <$> shrink' (trueOrigin tm) <*> shrink' (falseOrigin tm) <*> [TM.gridScale tm] + + +instance Arbitrary GridOffset where + arbitrary = GridOffset <$> genOffset 100000 <*> genOffset 100000 <*> genAlt + shrink d = tail $ GridOffset <$> + shrinkLength (deltaEast d) <*> shrinkLength (deltaNorth d) <*> shrinkLength (deltaAltitude d) + + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (GridStereo e) where + arbitrary = mkGridStereo <$> arbitrary <*> arbitrary <*> ((*~ one) <$> choose (0.95,1.0)) + shrink sg = tail $ mkGridStereo <$> shrink' (gridTangent sg) <*> shrink' (gridOrigin sg) <*> [SG.gridScale sg] + + +-- | Wrapper for arbitrary rays, along with creation parameters for printing and shrinking. +data Ray e = Ray (Geodetic e) (Angle Double) (Angle Double) + +instance (Ellipsoid e) => Show (Ray e) where + show (Ray p0 b e ) = "(Ray " ++ show p0 ++ ", " ++ showAngle b ++ ", " ++ showAngle e ++ ")" + +getRay :: (Ellipsoid e) => Ray e -> Path e +getRay (Ray p0 b e) = rayPath p0 b e + +instance (Ellipsoid e, Arbitrary e) => Arbitrary (Ray e) where + arbitrary = do + p0 <- arbitrary + b <- (*~ degree) <$> choose (-180,180) + e <- (*~ degree) <$> choose (0,90) + return $ Ray p0 b e + shrink (Ray p0 b e) = tail $ do + p0' <- shrink' p0 + b' <- shrinkAngle b + e' <- shrinkAngle e + return $ Ray p0' b' e' + + +-- | Two rhumb paths starting not more than 1000 km apart. +data RhumbPaths2 = RP2 { + rp2Point0 :: Geodetic WGS84, + rp2Bearing0 :: Bearing, + rp2Distance :: Distance2, + rp2Bearing1 :: Bearing, + rp2Bearing2 :: Bearing + } + +instance Show RhumbPaths2 where + show rp2 = show (pt1, Bearing b1) ++ show (pt2, Bearing b2) + where + (p1, p2) = mk2RhumbPaths rp2 + (pt1, b1, _) = pathFunc p1 (0 *~ meter) + (pt2, b2, _) = pathFunc p2 (0 *~ meter) + +instance Arbitrary RhumbPaths2 where + arbitrary = RP2 + <$> arbitrary `suchThat` ((< 70 *~ degree) . abs . latitude) + <*> arbitrary <*> arbitrary <*> arbitrary <*> arbitrary + shrink rp = + tail $ RP2 <$> + shrink' (rp2Point0 rp) <*> + shrink' (rp2Bearing0 rp) <*> + shrink' (rp2Distance rp) <*> + shrink' (rp2Bearing1 rp) <*> + shrink' (rp2Bearing2 rp) + +mk2RhumbPaths :: RhumbPaths2 -> (Path WGS84, Path WGS84) +mk2RhumbPaths (RP2 pt0 (Bearing b0) (Distance2 d) (Bearing b1) (Bearing b2)) = + (path1, path2) + where + path0 = rhumbPath pt0 b0 + path1 = rhumbPath pt0 b1 + (pt2, _, _) = pathFunc path0 d + path2 = rhumbPath pt2 b2 + diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -0,0 +1,390 @@ + +module Main where + +import Data.Maybe +import Data.Monoid +import Numeric.Units.Dimensional.Prelude +import qualified Prelude as P +import Test.Framework (Test, defaultMainWithOpts, testGroup) +import Test.Framework.Options (TestOptions, TestOptions'(..)) +import Test.Framework.Runners.Options (RunnerOptions, RunnerOptions'(..)) +import Test.Framework.Providers.HUnit +import Test.Framework.Providers.QuickCheck2 (testProperty) +import qualified Test.HUnit as HU +import Test.QuickCheck + +import ArbitraryInstances +import Geodetics.Altitude +import Geodetics.Ellipsoids +import Geodetics.Geodetic +import Geodetics.Grid +import Geodetics.Path +import Geodetics.Stereographic +import Geodetics.TransverseMercator +import Geodetics.UK + + +main :: IO () +main = do + let empty_test_opts = mempty :: TestOptions + let my_test_opts = empty_test_opts { + topt_maximum_generated_tests = Just 1000 + } + + let empty_runner_opts = mempty :: RunnerOptions + let my_runner_opts = empty_runner_opts { + ropt_test_options = Just my_test_opts + } + + defaultMainWithOpts tests my_runner_opts + +tests :: [Test] +tests = [ + testGroup "Geodetic" [ + testProperty "WGS84 and back" prop_WGS84_and_back, + testGroup "UK Points" $ map pointTest ukPoints], + testGroup "World lines" $ map worldLineTests worldLines, + testGroup "Grid" [ + testProperty "Grid Offset 1" prop_offset1, + testProperty "Grid Offset 2" prop_offset2, + testProperty "Grid Offset 3" prop_offset3, + testProperty "Grid 1" prop_grid1 ], + testGroup "UK" [ + testProperty "UK Grid 1" prop_ukGrid1, + testGroup "UK Grid 2" $ map ukGridTest2 ukSampleGrid, + testGroup "UK Grid 3" $ map ukGridTest3 ukSampleGrid, + testGroup "UK Grid 4" $ map ukGridTest4 ukSampleGrid, + testGroup "UK Grid 5" $ map ukGridTest5 ukSampleGrid + ], + testGroup "Stereographic" [ + testCase "toGrid north" $ HU.assertBool "" stereographicToGridN, + testCase "fromGrid north" $ HU.assertBool "" stereographicFromGridN, + testCase "toGrid south" $ HU.assertBool "" stereographicToGridS, + testCase "fromGrid south" $ HU.assertBool "" stereographicFromGridS, + testProperty "Stereographic round trip" prop_stereographic + ], + testGroup "Paths" [ + testProperty "Ray Path 1" prop_rayPath1, + testProperty "Ray Continuity" prop_rayContinuity, + testProperty "Ray Bisection" prop_rayBisect, + testProperty "Rhumb Continuity" prop_rhumbContinuity, + testProperty "Rhumb Intersection" prop_rhumbIntersect + ] + ] + + +-- | The positions are within 30 cm. +samePlace :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool +samePlace p1 p2 = geometricalDistance p1 p2 < 0.3 *~ meter + + +-- | The positions are within 10 m. +closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool +closeEnough p1 p2 = geometricalDistance p1 p2 < 10 *~ meter + + +-- | The angles are within 0.01 arcsec +sameAngle :: Angle Double -> Angle Double -> Bool +sameAngle v1 v2 = abs (properAngle (v1 - v2)) < 0.01 *~ arcsecond + +-- | The grid positions are within 1mm +sameGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool +sameGrid p1 p2 = check eastings && check northings && check altitude + where check f = f p1 - f p2 < 1 *~ milli meter + + +-- | Grid offsets are within 1mm. +sameOffset :: GridOffset -> GridOffset -> Bool +sameOffset go1 go2 = check deltaNorth && check deltaEast && check deltaAltitude + where check f = f go1 - f go2 < 1 *~ milli meter + + +-- | The grid X and Y are both within 1 meter +closeGrid :: (GridClass r e) => GridPoint r -> GridPoint r -> Bool +closeGrid p1 p2 = check eastings && check northings && check altitude + where check f = f p1 - f p2 < 1 *~ meter + +-- | Degrees, minutes and seconds into radians. +dms :: Int -> Int -> Double -> Dimensionless Double +dms d m s = fromIntegral d *~ degree + fromIntegral m *~ arcminute + s *~ arcsecond + +-- | Round-trip from local to WGS84 and back is identity (approximately) +prop_WGS84_and_back :: Geodetic LocalEllipsoid -> Bool +prop_WGS84_and_back p = samePlace p $ toLocal (ellipsoid p) $ toWGS84 p + + +-- | Sample pairs of points with bearings and distances. +-- The Oracle for these values is the @FORWARD@ program from +-- http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html +worldLines :: [(String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double)] +worldLines = [ + ("Ordinary", Geodetic (40*~degree) (30*~degree) _0 WGS84, Geodetic (30*~degree) (50*~degree) _0 WGS84, + 2128852.999*~meter, 115.19596706*~degree, 126.79044315*~degree), + ("Over Pole", Geodetic (60*~degree) (0*~degree) _0 WGS84, Geodetic (60*~degree) (180*~degree) _0 WGS84, + 6695785.820*~meter, 0*~degree, 180*~degree), + ("Equator to Pole", Geodetic (0*~degree) (0*~degree) _0 WGS84, Geodetic (90*~degree) (180*~degree) _0 WGS84, + 10001965.729*~meter, 0*~degree, 180*~degree)] + + +worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Length Double, Dimensionless Double, Dimensionless Double) -> Test +worldLineTests (str, g1, g2, d, a, b) = testCase str $ HU.assertBool "" $ ok $ groundDistance g1 g2 + where + ok Nothing = False + ok (Just (d1, a1, b1)) = + abs (d - d1) < 0.01 *~ meter + && abs (a - a1) < 0.01 *~ arcsecond + && abs (b - b1) < 0.01 *~ arcsecond + +-- | Sample points for UK tests. The oracle for these values is the script at +-- http://www.movable-type.co.uk/scripts/latlong-convert-coords.html, which uses +-- the same Helmert transform as this library. Hence the results should match to within 30 cm. +ukPoints :: [(String, Geodetic WGS84, Geodetic OSGB36)] +ukPoints = [ + ("Greenwich", Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) _0 WGS84, + Geodetic (dms 51 28 39.00) (dms 0 0 0) _0 OSGB36), + ("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) _0 WGS84, + Geodetic (dms 55 56 56.51) (dms (-3) (-11) (-57.61)) _0 OSGB36), + ("Lands End", Geodetic (dms 50 03 56.68) (dms (-5) (-42) (-51.20)) _0 WGS84, + Geodetic (dms 50 03 54.51) (dms (-5) (-42) (-47.87)) _0 OSGB36), + ("Gt. Yarmouth Pier",Geodetic (dms 52 36 29.33) (dms 1 44 27.79) _0 WGS84, + Geodetic (dms 52 36 27.84) (dms 1 44 34.52) _0 OSGB36), + ("Stanhope", Geodetic (dms 54 44 49.08) (dms (-2) 0 (-19.89)) _0 WGS84, + Geodetic (dms 54 44 48.71) (dms (-2) 0 (-14.41)) _0 OSGB36) ] + + + +-- Convert a named point into a test +pointTest :: (Ellipsoid e2) => (String, Geodetic WGS84, Geodetic e2) -> Test +pointTest (name, wgs84, local) = testCase name $ HU.assertBool "" $ samePlace wgs84 (toWGS84 local) + + +-- The negation of the sum of a list of offsets is equal to the sum of the negated items. +prop_offset1 :: [GridOffset] -> Bool +prop_offset1 offsets = sameOffset (offsetNegate $ mconcat offsets) (mconcat $ map offsetNegate offsets) + +-- A polar offset multiplied by a scalar is equal to an offset in the same direction with the length multiplied. +prop_offset2 :: Distance -> Bearing -> Scalar -> Bool +prop_offset2 (Distance d) (Bearing h) (Scalar s) = sameOffset go1 go2 + where + go1 = offsetScale s $ polarOffset d h + go2 = polarOffset (d * s) h + +-- | A polar offset has the offset distance and bearing of its arguments. +prop_offset3 :: GridOffset -> Bool +prop_offset3 delta = sameOffset delta0 + (polarOffset (offsetDistance delta0) (offsetBearing delta)) + where delta0 = delta {deltaAltitude = 0 *~ meter} + +-- | Given a grid point and an offset, applying the offset to the point gives a new point which +-- is offset from the first point by the argument offset. +prop_grid1 :: GridPoint (GridTM LocalEllipsoid) -> GridOffset -> Bool +prop_grid1 p d = sameOffset d $ p `gridOffset` applyOffset d p + + +-- | Converting a UK grid reference to a GridPoint and back is a null operation. +prop_ukGrid1 :: GridRef -> Bool +prop_ukGrid1 (GridRef str) = + str == + (fromJust $ toUkGridReference ((length str P.- 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str) + +-- | UK Grid Reference points. The oracle for these points was the +-- UK Grid Reference Finder (gridreferencefinder.com), retrieved on 26 Jan 2013. +ukSampleGrid :: [(String, GridPoint UkNationalGrid, Geodetic WGS84, String)] +ukSampleGrid = map convert [ + -- Grid Reference, X, Y, Latitude, Longitude, Description + ("SW3425625070", 134256, 025070, 50.066230, -5.7148278, "Lands End"), + ("TR3302139945", 633021, 139945, 51.111396, 1.3277159, "Dover Harbour"), + ("TQ3001980417", 530019, 180417, 51.507736, -0.12793230, "Nelsons Column"), + ("TA2542370644", 525423, 470644, 54.116376, -0.082668990, "Flamborough Lighthouse"), + ("NK1354745166", 413547, 845166, 57.496512, -1.7756310, "Peterhead harbour"), + ("ND3804872787", 338048, 972787, 58.638518, -3.0688688, "John O Groats"), + ("SC3915875189", 239158, 475189, 54.147275, -4.4641148, "Douglas Harbour"), + ("ST1922474591", 319224, 174591, 51.464505, -3.1641741, "Torchwood HQ"), + ("SK3520736502", 435207, 336502, 52.924784, -1.4777486, "Derby Cathedral")] + where + convert (grid, x, y, lat, long, desc) = + (grid, GridPoint (x *~ meter) (y *~ meter) (0 *~ meter) UkNationalGrid, + Geodetic (lat *~ degree) (long *~ degree) (0 *~ meter) WGS84, desc) + +type GridPointTest = (String, GridPoint UkNationalGrid, Geodetic WGS84, String) -> Test + +-- | Check that grid reference to grid point works for sample points. +ukGridTest2 :: GridPointTest +ukGridTest2 (gridRef, gp, _, name) = testCase name $ HU.assertBool "" + $ (fst $ fromJust $ fromUkGridReference gridRef) == gp + +-- | Check that grid point to grid reference works for sample points. +ukGridTest3 :: GridPointTest +ukGridTest3 (gridRef, gp, _, name) = testCase name $ HU.assertBool "" + $ toUkGridReference 5 gp == Just gridRef + +-- | Check that grid point to WGS84 works close enough for sample points. +ukGridTest4 :: GridPointTest +ukGridTest4 (_, gp, geo, name) = testCase name $ HU.assertBool "" + $ closeEnough geo $ toWGS84 $ fromGrid gp + +-- | Check that WGS84 to grid point works close enough for sample points. +ukGridTest5 :: GridPointTest +ukGridTest5 (_, gp, geo, name) = testCase name $ HU.assertBool "" + $ offsetDistance (gridOffset gp $ toGrid UkNationalGrid $ toLocal OSGB36 geo) < 1 *~ meter + + +-- | Worked example for UK Geodetic to GridPoint, taken from "A Guide to Coordinate Systems in Great Britain" [1] +ukTest :: Geodetic OSGB36 +ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) (0 *~ meter) OSGB36 + +{- + v = 6.3885023333E+06 + rho = 6.3727564399E+06 + eta2 = 2.4708136169E-03 + m = 4.0668829596E+05 + I = 3.0668829596E+05 + II = 1.5404079092E+06 + III = 1.5606875424E+05 + IIIa = -2.0671123011E+04 + IV = 3.8751205749E+06 + V = -1.7000078208E+05 + VI = -1.0134470432E+05 + E = 651409.903 m + N = 313177.270 m +-} + + +-- | Standard stereographic grid for point tests in the Northern Hemisphere. +stereoGridN :: GridStereo LocalEllipsoid +stereoGridN = mkGridStereo tangent origin (0.9999079 *~ one) + where + ellipse = LocalEllipsoid "Bessel 1841" (6377397.155 *~ metre) (299.15281 *~ one) mempty + tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) (0 *~ meter) ellipse + origin = GridOffset (155000 *~ metre) (463000 *~ metre) (0 *~ meter) + + +-- | Standard steregraphic grid for point tests in the Southern Hemisphere. +-- +-- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated. +stereoGridS :: GridStereo LocalEllipsoid +stereoGridS = mkGridStereo tangent origin (0.9999079 *~ one) + where + ellipse = LocalEllipsoid "Bessel 1841" (6377397.155 *~ metre) (299.15281 *~ one) mempty + tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) (0 *~ meter) ellipse + origin = GridOffset ((-155000) *~ metre) (463000 *~ metre) (0 *~ meter) + + +-- | Data for the stereographic tests taken from http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf +stereographicToGridN :: Bool +stereographicToGridN = sameGrid g1 g1' + where + p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN + g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN + g1' = toGrid stereoGridN p1 + +stereographicFromGridN :: Bool +stereographicFromGridN = samePlace p1 p1' + where + p1 = Geodetic (dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridN + g1 = GridPoint (196105.283 *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridN + p1' = fromGrid g1 + + +stereographicToGridS :: Bool +stereographicToGridS = sameGrid g1 g1' + where + p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridS + g1 = GridPoint ((-196105.283) *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridS + g1' = toGrid stereoGridS p1 + + +stereographicFromGridS :: Bool +stereographicFromGridS = samePlace p1 p1' + where + p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) (0 *~ meter) $ gridEllipsoid stereoGridS + g1 = GridPoint ((-196105.283) *~ meter) (557057.739 *~ meter) (0 *~ meter) stereoGridS + p1' = fromGrid g1 + + +-- | Check the round trip for a stereographic projection. +prop_stereographic :: GridPoint (GridStereo LocalEllipsoid) -> Property +prop_stereographic p = + let g = fromGrid p + r = toGrid (gridBasis p) g + in printTestCase ("p = " ++ show p ++ "\ng = " ++ show g ++ "\nr = " ++ show r) $ + closeGrid p r + + + +-- | A ray at distance zero returns its original arguments. +prop_rayPath1 :: Ray WGS84 -> Bool +prop_rayPath1 r@(Ray pt b e) = + samePlace pt pt1 && sameAngle b b1 && sameAngle e e1 + where (pt1,b1,e1) = pathFunc (getRay r) _0 + + +type ContinuityTest e = Geodetic e -> Bearing -> Azimuth -> Distance -> Distance -> Property + +type ContinuityTest1 e = Geodetic e -> Bearing -> Distance2 -> Distance2 -> Property + +-- | Many paths can be specified by a start point, bearing and azimuth, +-- and have the property that any (point,bearing,azimuth) triple on +-- the path will specify the same path with a distance offset. +prop_pathContinuity :: (Ellipsoid e) => + (Geodetic e -> Angle Double -> Angle Double -> Path e) -> ContinuityTest e +prop_pathContinuity pf pt0 (Bearing b0) (Azimuth a0) (Distance d1) (Distance d2) = + printTestCase (show ((pt2, Bearing b2, Azimuth a2), (pt3, Bearing b3, Azimuth a3))) $ + pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==> + closeEnough pt2 pt3 && sameAngle b2 b3 && sameAngle a2 a3 + where + path0 = pf pt0 b0 a0 + (pt1, b1, a1) = pathFunc path0 d1 + path1 = pf pt1 b1 a1 + (pt2, b2, a2) = pathFunc path1 d2 + (pt3, b3, a3) = pathFunc path0 (d1 + d2) -- Points 2 and 3 should be the same. + + +-- | For continuity testing of ground-based paths (azimuth & altitude always zero) +-- where lower accuracy is required. +prop_pathContinuity1 :: (Ellipsoid e) => (Geodetic e -> Angle Double -> Path e) -> ContinuityTest1 e +prop_pathContinuity1 pf pt0 (Bearing b0) (Distance2 d1) (Distance2 d2) = + printTestCase (show ((pt2, Bearing b2), (pt3, Bearing b3))) $ + pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==> + closeEnough pt2 pt3 && sameAngle b2 b3 + where + path0 = pf pt0 b0 + (pt1, b1, _) = pathFunc path0 d1 + path1 = pf pt1 b1 + (pt2, b2, _) = pathFunc path1 d2 + (pt3, b3, _) = pathFunc path0 (d1 + d2) -- Points 2 and 3 should be the same. + + +-- | A point on a ray will continue along the same ray, and hence give the same points. +prop_rayContinuity :: ContinuityTest WGS84 +prop_rayContinuity = prop_pathContinuity rayPath + + +-- | A ray bisected to an altitude will give that altitude. +-- This is a test of bisection rather than rays. +prop_rayBisect :: Ray WGS84 -> Altitude -> Bool +prop_rayBisect r (Altitude height) = + case bisect ray0 f (1 *~ centi meter) (0 *~ meter) (1000 *~ kilo meter) of + Nothing -> False + Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1 *~ centi meter + where + f g = compare (altitude g) height + ray0 = getRay r + + +-- | A point on a rhumb line will continue along the same rhumb. +prop_rhumbContinuity :: ContinuityTest1 WGS84 +prop_rhumbContinuity = prop_pathContinuity1 rhumbPath + + +-- | Two rhumb paths intersect at the same place. +prop_rhumbIntersect :: RhumbPaths2 -> Property +prop_rhumbIntersect rp = + case intersect _0 _0 (10.0 *~ centi meter) 100 path1 path2 of + Just (d1, d2) -> + let (pt1, _, _) = pathFunc path1 d1 + (pt2, _, _) = pathFunc path2 d2 + in printTestCase (show (pt1, pt2)) $ label "Intersection" $ samePlace pt1 pt2 + Nothing -> label "No intersection" True + where + (path1, path2) = mk2RhumbPaths rp