diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -3,11 +3,11 @@ 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 +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 +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. @@ -18,19 +18,3 @@ 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. - -Release Notes -------------- - -Version 0.0.2: Tided up cabal file and removed spurious dependency on Parsec. - -Version 0.0.3: Updated for Haskell Platform 2014.2.0.0 and GHC 7.8.3. Fixed - some minor documentation issues. - -Version 0.0.4: Updated for Dimensional 1.0. - -Version 0.0.5: Fixed bug in Monoid instance for Helmert. Created Semigroup - instance for Helmert. - -Version 0.0.6: Prevent attempted building on GHC 7.8 (it doesn't work) - and fix the build on 7.10 with a conditional semigroups dependency diff --git a/changelog.md b/changelog.md new file mode 100644 --- /dev/null +++ b/changelog.md @@ -0,0 +1,17 @@ +Release Notes +------------- + +Version 0.0.2: Tided up cabal file and removed spurious dependency on Parsec. + +Version 0.0.3: Updated for Haskell Platform 2014.2.0.0 and GHC 7.8.3. Fixed + some minor documentation issues. + +Version 0.0.4: Updated for Dimensional 1.0. + +Version 0.0.5: Fixed bug in Monoid instance for Helmert. Created Semigroup + instance for Helmert. + +Version 0.0.6: Prevent attempted building on GHC 7.8 (it doesn't work) + and fix the build on 7.10 with a conditional semigroups dependency + +Version 0.1.0: Updated for Dimensional 1.3 and GHC 8.6. diff --git a/geodetics.cabal b/geodetics.cabal --- a/geodetics.cabal +++ b/geodetics.cabal @@ -1,27 +1,28 @@ name: geodetics -version: 0.0.6 +version: 0.1.0 cabal-version: >= 1.10 build-type: Simple author: Paul Johnson -data-files: - AddingProjections.txt, - LICENSE, - README.md, +data-files: + AddingProjections.txt, + LICENSE, + README.md, + changelog.md, ToDo.txt license: BSD3 -copyright: Paul Johnson 2015. +copyright: Paul Johnson 2018. synopsis: Terrestrial coordinate systems and geodetic 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 + 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.10.2 +tested-with: GHC==8.6.3 source-repository head type: git @@ -29,14 +30,13 @@ library hs-source-dirs: src - build-depends: - base >= 4.7 && < 5, - dimensional >= 1.0, - array >= 0.4 - if !impl(ghc>=8.0) - build-depends: semigroups >= 0.9 && < 0.19 + build-depends: + base >= 4.6 && < 5, + dimensional >= 1.3, + array >= 0.4, + semigroups >= 0.9 ghc-options: -Wall - exposed-modules: + exposed-modules: Geodetics.Altitude, Geodetics.Ellipsoids, Geodetics.Geodetic, @@ -55,16 +55,16 @@ build-depends: geodetics, base >= 4.6 && < 5, HUnit >= 1.2, - dimensional >= 0.13, + dimensional >= 1.3, QuickCheck >= 2.4, test-framework >= 0.4.1, test-framework-quickcheck2, test-framework-hunit, array >= 0.4, checkers - hs-source-dirs: + hs-source-dirs: test ghc-options: -Wall -rtsopts - other-modules: + other-modules: ArbitraryInstances Default-Language: Haskell2010 diff --git a/src/Geodetics/Ellipsoids.hs b/src/Geodetics/Ellipsoids.hs --- a/src/Geodetics/Ellipsoids.hs +++ b/src/Geodetics/Ellipsoids.hs @@ -1,6 +1,21 @@ -{-# LANGUAGE FlexibleContexts, TypeOperators, TypeFamilies #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE RoleAnnotations #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} -{- | An Ellipsoid is a reasonable best fit for the surface of the +{- | 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. @@ -43,7 +58,8 @@ import Data.Semigroup (Semigroup, (<>)) import Numeric.Units.Dimensional import Numeric.Units.Dimensional.Prelude -import Prelude () -- Numeric instances. +import qualified Numeric.Units.Dimensional.Dimensions.TypeLevel as T +-- import Prelude () -- Numeric instances. -- | 3d vector as @(X,Y,Z)@. @@ -55,7 +71,7 @@ -- | Multiply a vector by a scalar. scale3 :: (Num a) => - Vec3 (Quantity d a) -> Quantity d' a -> Vec3 (Quantity (d * d') a) + Vec3 (Quantity d a) -> Quantity d' a -> Vec3 (Quantity (d T.* d') a) scale3 (x,y,z) s = (x*s, y*s, z*s) @@ -70,7 +86,7 @@ -- | Multiply a matrix by a vector in the Dimensional type system. transform3 :: (Num a) => - Matrix3 (Quantity d a) -> Vec3 (Quantity d' a) -> Vec3 (Quantity (d*d') a) + Matrix3 (Quantity d a) -> Vec3 (Quantity d' a) -> Vec3 (Quantity (d T.* 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 @@ -78,7 +94,7 @@ -- | Inverse of a 3x3 matrix. invert3 :: (Fractional a) => - Matrix3 (Quantity d a) -> Matrix3 (Quantity ((d*d)/(d*d*d)) a) + Matrix3 (Quantity d a) -> Matrix3 (Quantity ((d T.* d)/(d T.* d T.* d)) a) invert3 ((x1,y1,z1), (x2,y2,z2), (x3,y3,z3)) = @@ -96,12 +112,12 @@ -- | Dot product of two vectors dot3 :: (Num a) => - Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Quantity (d1 * d2) a + Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Quantity (d1 T.* d2) a dot3 (x1,y1,z1) (x2,y2,z2) = x1*x2 + y1*y2 + z1*z2 -- | Cross product of two vectors cross3 :: (Num a) => - Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Vec3 (Quantity (d1 * d2) a) + Vec3 (Quantity d1 a) -> Vec3 (Quantity d2 a) -> Vec3 (Quantity (d1 T.* d2) a) cross3 (x1,y1,z1) (x2,y2,z2) = (y1*z2 - z1*y2, z1*x2 - x1*z2, x1*y2 - y1*x2) @@ -122,8 +138,8 @@ -- | The inverse of a Helmert transformation. inverseHelmert :: Helmert -> Helmert -inverseHelmert h = Helmert (negate $ cX h) (negate $ cY h) (negate $ cZ h) - (negate $ helmertScale h) +inverseHelmert h = Helmert (negate $ cX h) (negate $ cY h) (negate $ cZ h) + (negate $ helmertScale h) (negate $ rX h) (negate $ rY h) (negate $ rZ h) @@ -141,16 +157,16 @@ s = _1 + helmertScale h * (1e-6 *~ one) --- | An Ellipsoid is defined by the major radius and the inverse flattening (which define its shape), +-- | 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 +-- 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 @@ -159,7 +175,7 @@ -- ^ Inverse of the flattening. helmert :: a -> Helmert helmertToWSG84 :: a -> ECEF -> ECEF - -- ^ The Helmert transform that will convert a position wrt + -- ^ 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 @@ -168,9 +184,9 @@ -- | 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 +-- 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 @@ -179,15 +195,15 @@ 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. -- @@ -200,7 +216,7 @@ helmertLocal :: Helmert } deriving (Eq) instance Show LocalEllipsoid where - show = nameLocal + show = nameLocal instance Ellipsoid LocalEllipsoid where majorRadius = majorRadiusLocal @@ -226,38 +242,38 @@ 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 +-- | 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 +-- | 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. +-- | 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) +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 +-- | 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)