hcoord (empty) → 1.0.0.0
raw patch · 16 files changed
+2111/−0 lines, 16 filesdep +HUnitdep +basedep +data-defaultsetup-changed
Dependencies added: HUnit, base, data-default, hcoord, ieee754, mtl, regex-pcre
Files
- LICENSE +30/−0
- README.md +103/−0
- Setup.hs +2/−0
- app/Main.hs +4/−0
- hcoord.cabal +57/−0
- src/DMS.hs +40/−0
- src/Datum.hs +370/−0
- src/ECEFRef.hs +62/−0
- src/Ellipsoid.hs +115/−0
- src/IrishRef.hs +209/−0
- src/LatLng.hs +212/−0
- src/MGRSRef.hs +295/−0
- src/MathExtensions.hs +21/−0
- src/OSRef.hs +395/−0
- src/UTMRef.hs +177/−0
- test/Spec.hs +19/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Daniele Francesconi (c) 2016++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * 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.++ * Neither the name of Daniele Francesconi nor the names of other+ 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+OWNER 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.
+ README.md view
@@ -0,0 +1,103 @@+# HCoord + +[](https://hackage.haskell.org/package/hcoord) +[](https://travis-ci.org/danfran/hcoord) + +Note that this repo is a _translation_ of the 2006 [Jcoord](http://www.jstott.me.uk/jcoord) project to Haskell. + +## Original Project Homepage ## +Conversion between latitude/longitude and grid references is not as simple as one would expect. Indeed, the maths involved is rather complex and not really for the faint-hearted! I have created a Java package which provides a simple API to allow conversion between OSGB (Ordnance Survey of Great Britain) grid references, UTM (Universal Transverse Mercator) references and latitude/longitude. OSGB Grid references can be returned as a standard six-figure grid reference using the 100km grid square letters (e.g. SK619847). + +_Note: it is now possible to obtain an eight-figure grid reference as from [release 1.2](https://github.com/xni06/JCoord/releases/tag/1.2)._ + +The Jcoord script also contains a method to allow the calculation of the surface distance between two points of latitude/longitude. + +Looking for the same thing but with JavaScript or PHP? Translations of this package with the same functionality are available for JavaScript (see JScoord) and PHP (see PHPcoord). + +Information on converting between UTM and latitude/longitude is available at www.posc.org. + +The Ordnance Survey have published a concise guide on converting between OSGB and latitude/longitude which is available at www.gps.gov.uk. + +Bear in mind that Jcoord is still in development, so check back often for updates. To contact me about Jcoord, send an e-mail to jcoord@jstott.me.uk. + +``` +-------------------------------------------------------------------------- + Jcoord + readme.txt + + (c) 2006 Jonathan Stott + + Created on 11-Feb-2006 +-------------------------------------------------------------------------- + +Jcoord is a collection of Java classes that provide functions for handling +various co-ordinate systems and converting between them. Currently, OSGB +(Ordnance Survey of Great Britain) grid references, UTM (Universal Transverse +Mercator) references and latitude/longitude are supported. A function is +also provided to find the surface distance between two points of latitude +and longitude. + +When using the OSGB conversions, the majority of applications use the +WGS84 datum rather than the OSGB36 datum. Conversions between the two +data were added in v1.1 - the conversions should be accurate to within +5m or so. If accuracy is not important (i.e. to within 200m or so), +then it isn't necessary to perform the conversions. + +Examples of how to use the classes in Jcoord can be found in the +uk.me.jstott.jcoord.Test class + +See http://www.jstott.me.uk/jcoord/ for latest releases and information. + + +DISCLAIMER + +Accuracy of the co-ordinate conversions contained within the Jcoord +package is not guaranteed. Use of the conversions is entirely at your +own risk and I cannot be held responsible for any consequences of +errors created by the conversions. I do not recommend using the package +for mission-critical applications. + + +LICENSING + +This software product is available under the GNU General Public License +(GPL). Terms of the GPL can be read at http://www.jstott.me.uk/gpl/. +Any commercial use requires the purchase of a license - contact me at +jcoord@jstott.me.uk for details. + + +VERSION HISTORY + + 1.1 - + - Added Junit unit test suite in the uk.me.jstott.jcoord.junit package. + - Added NotDefinedOnUTMGridException which is thrown whenever an attempt is + made to convert a LatLng object to a UTMRef object where the latitude falls + outside the defined UTM grid. + - Can format a latitude/longitude into a String in degrees-minutes-seconds + format. + - Added a class for Military Grid Reference System (MGRS) to allow for + conversrions between MGRS and UTM and latitude/longitude. + - Documentation (javadoc) improvements. + - LatLng.getLat() has been deprecated. Use LatLng.getLatitude() instead. + - LatLng.getLng() has been deprecated. Use LatLng.getLongitude() instead. + - LatLng constructor now validates parameters and throws an + IllegalArgumentException if either the latitude or longitude are invalid. + - UTMRef(double,double,char,int) constructor has been deprecated and + replaced with UTMRef(int,char,double,double). This is to make the order of + the parameters more sensible. + - Fixed slight error with values for converting between WGS84 and OSGB36. + - New method: LatLng.distanceMiles(LatLng) to calculate the distance between + two points of latitude/longitude in miles. + - New constructors in LatLng to allow for specifying a datum and/or the + latitude and longitude as degrees/minutes/seconds. + - Added IrishRef class to allow conversion between latitude and longitude and + Irish National Grid References. + - Added new uk.me.jstott.jcoord.datum package with a number of data that + all extend uk.me.jstott.jcoord.datum.Datum. + - Added new uk.me.jstott.jcoord.ellipsoid package with a number of definitions + of reference ellipsoids that all extend + uk.me.jstott.jcoord.ellipsoid.Ellipsoid. + - Fixed slight errors in OSRef to LatLng conversions. + 1.0 - 11 Feb 2006 + - Initial version created from PHPcoord v2.2. +```
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ app/Main.hs view
@@ -0,0 +1,4 @@+module Main where++main :: IO ()+main = undefined
+ hcoord.cabal view
@@ -0,0 +1,57 @@+name: hcoord+version: 1.0.0.0+synopsis: Easily convert between latitude/longitude, UTM and OSGB.+description: Please see README.md+homepage: https://github.com/danfran/hcoord#readme+license: BSD3+license-file: LICENSE+author: Daniele Francesconi+maintainer: dfrancesconi12@gmail.com+copyright: copyleft+category: geography+build-type: Simple+extra-source-files: README.md+cabal-version: >=1.10++library+ hs-source-dirs: src+ exposed-modules:+ Datum+ DMS+ ECEFRef+ Ellipsoid+ IrishRef+ LatLng+ MathExtensions+ MGRSRef+ OSRef+ UTMRef+ build-depends: base >= 4 && < 5+ , mtl+ , regex-pcre+ default-language: Haskell2010++executable hcoord-exe+ hs-source-dirs: app+ main-is: Main.hs+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ build-depends: base+ , hcoord+ default-language: Haskell2010++test-suite hcoord-test+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Spec.hs+ build-depends: base+ , HUnit >= 1.2 && < 1.6+ , hcoord+ , mtl+ , ieee754+ , data-default+ ghc-options: -threaded -rtsopts -with-rtsopts=-N+ default-language: Haskell2010++source-repository head+ type: git+ location: https://github.com/danfran/hcoord
+ src/DMS.hs view
@@ -0,0 +1,40 @@+{-# LANGUAGE CPP #-}+module DMS where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Datum+import LatLng++data LatitudeDMS = North DMS -- ^ Latitude is north of the equator.+ | South DMS -- ^ Latitude is south of the equator.++data LongitudeDMS = East DMS -- ^ Longitude is east of the prime meridian.+ | West DMS -- ^ Longitude is west of the prime meridian.++data DMS = DMS { degrees :: Double+ , minutes :: Double+ , seconds :: Double+ }+++toLatLng :: LatitudeDMS -> LongitudeDMS -> Datum -> Except String LatLng+toLatLng lat lng dtm = do+ lt <- withExcept (const "Invalid latitude") (evalLatitude lat)+ ln <- withExcept (const "Invalid longitude") (evalLongitude lng)+ pure LatLng { latitude = lt, longitude = ln, height = 0, datum = dtm }+ where evalLatitude :: LatitudeDMS -> Except String Double+ evalLatitude (North p) = dmsToLatLngPoint p 1+ evalLatitude (South p) = dmsToLatLngPoint p (-1)++ evalLongitude :: LongitudeDMS -> Except String Double+ evalLongitude (East p) = dmsToLatLngPoint p 1+ evalLongitude (West p) = dmsToLatLngPoint p (-1)++ dmsToLatLngPoint :: DMS -> Double -> Except String Double+ dmsToLatLngPoint DMS { degrees = d, minutes = m, seconds = s } cardinal+ | d < 0.0 || m < 0.0 || s < 0.0 || d > 90.0 || m >= 60.0 || s >= 60.0 = throwError "Invalid point"+ | otherwise = pure (cardinal * (d + m / 60.0 + s / 3600.0))
+ src/Datum.hs view
@@ -0,0 +1,370 @@+module Datum where++import Ellipsoid++{- |+ To represent a set of parameters for describing a particular datum,+ including a name, the reference ellipsoid used and the seven parameters+ required to translate co-ordinates in this datum to the WGS84 datum.+-}+data Datum = Datum { name :: String -- ^ The name of this Datum.+ , ellipsoid :: Ellipsoid -- ^ The reference ellipsoid associated with this Datum.+ {-|+ Translation along the x-axis for use in 7-parameter Helmert+ transformations. This value should be used to convert a co-ordinate in a+ given datum to the WGS84 datum.+ -}+ , dx :: Double+ {-|+ Translation along the y-axis for use in 7-parameter Helmert+ transformations. This value should be used to convert a co-ordinate in a+ given datum to the WGS84 datum.+ -}+ , dy :: Double+ {-|+ Translation along the z-axis for use in 7-parameter Helmert+ transformations. This value should be used to convert a co-ordinate in a+ given datum to the WGS84 datum.+ -}+ , dz :: Double+ {-|+ Scale factor for use in 7-parameter Helmert transformations. This value+ should be used to convert a co-ordinate in a given datum to the WGS84+ datum.+ -}+ , ds :: Double+ {-|+ Rotation about the x-axis for use in 7-parameter Helmert transformations.+ This value should be used to convert a co-ordinate in a given datum to the+ WGS84 datum.+ -}+ , rx :: Double+ {-|+ Rotation about the y-axis for use in 7-parameter Helmert transformations.+ This value should be used to convert a co-ordinate in a given datum to the+ WGS84 datum.+ -}+ , ry :: Double+ {-|+ Rotation about the z-axis for use in 7-parameter Helmert transformations.+ This value should be used to convert a co-ordinate in a given datum to the+ WGS84 datum.+ -}+ , rz :: Double+ } deriving (Eq, Show)++-- instance Show Datum where+-- show Datum{name = n, ellipsoid = ell, dx = dx, dy = dy, dz = dz,+-- ds = ds, rx = rx, ry = ry, rz = rz} =+-- n ++ " " ++ show ell ++ " dx=" ++ show dx ++ " dy=" ++ show dy +++-- " dz=" ++ show dz ++ " ds=" ++ show ds ++ " rx=" ++ show rx ++ " ry=" ++ show ry ++ " rz=" ++ show rz+++-- | Pre-determined data:++etrf89Datum :: Datum+etrf89Datum = Datum { name = "European Terrestrial Reference Frame (ETRF89)"+ , ellipsoid = wgs84Ellipsoid+ , dx = 0.0+ , dy = 0.0+ , dz = 0.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }++ireland1965Datum :: Datum+ireland1965Datum = Datum { name = "Ireland 1965"+ , ellipsoid = modifiedAiryEllipsoid+ , dx = 482.53+ , dy = -130.596+ , dz = 564.557+ , ds = 8.15+ , rx = -1.042+ , ry = -0.214+ , rz = -0.631+ }++osgb36Datum :: Datum+osgb36Datum = Datum { name = "Ordnance Survey of Great Britain 1936 (OSGB36)"+ , ellipsoid = airy1830Ellipsoid+ , dx = 446.448+ , dy = -125.157+ , dz = 542.06+ , ds = -20.4894+ , rx = 0.1502+ , ry = 0.2470+ , rz = 0.8421+ }++wgs84Datum :: Datum+wgs84Datum = Datum { name = "World Geodetic System 1984 (WGS84)"+ , ellipsoid = wgs84Ellipsoid+ , dx = 0.0+ , dy = 0.0+ , dz = 0.0+ , ds = 1.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }++nad27GreenlandDatum :: Datum+nad27GreenlandDatum = Datum { name = "North American Datum 1927 (NAD27) - Greenland"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 11.0+ , dy = 114.0+ , dz = 195.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27WesternUSDatum :: Datum+nad27WesternUSDatum = Datum { name = "North American Datum 1927 (NAD27) - Western US"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -8.0+ , dy = 159.0+ , dz = 175.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27AlaskaDatum :: Datum+nad27AlaskaDatum = Datum { name = "North American Datum 1927 (NAD27) - Alaska"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -5.0+ , dy = 135.0+ , dz = 172.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CentralAmericaDatum :: Datum+nad27CentralAmericaDatum = Datum { name = "North American Datum 1927 (NAD27) - Central America"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 0.0+ , dy = 125.0+ , dz = 194.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27SanSalvadorDatum :: Datum+nad27SanSalvadorDatum = Datum { name = "North American Datum 1927 (NAD27) - San Salvador"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 1.0+ , dy = 140.0+ , dz = 165.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27AlbertaBritishColumbiaDatum :: Datum+nad27AlbertaBritishColumbiaDatum = Datum { name = "North American Datum 1927 (NAD27) - Alberta and British Columbia"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -7.0+ , dy = 162.0+ , dz = 188.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanadaEastDatum :: Datum+nad27CanadaEastDatum = Datum { name = "North American Datum 1927 (NAD27) - Canada East"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -22.0+ , dy = 160.0+ , dz = 190.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27ContiguousUSDatum :: Datum+nad27ContiguousUSDatum = Datum { name = "North American Datum 1927 (NAD27) - Contiguous United States"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -8.0+ , dy = 160.0+ , dz = 176.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27MexicoDatum :: Datum+nad27MexicoDatum = Datum { name = "North American Datum 1927 (NAD27) - Mexico"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -12.0+ , dy = 130.0+ , dz = 190.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27BahamasDatum :: Datum+nad27BahamasDatum = Datum { name = "North American Datum 1927 (NAD27) - Bahamas"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -4.0+ , dy = 154.0+ , dz = 178.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanadaNWTerritoryDatum :: Datum+nad27CanadaNWTerritoryDatum = Datum { name = "North American Datum 1927 (NAD27) - Canada NW Territory"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 4.0+ , dy = 159.0+ , dz = 188.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanadaManitobaOntarioDatum :: Datum+nad27CanadaManitobaOntarioDatum = Datum { name = "North American Datum 1927 (NAD27) - Canada Manitoba/Ontario"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -9.0+ , dy = 157.0+ , dz = 184.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanadaDatum :: Datum+nad27CanadaDatum = Datum { name = "North American Datum 1927 (NAD27) - Canada"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -10.0+ , dy = 158.0+ , dz = 187.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanadaYukonDatum :: Datum+nad27CanadaYukonDatum = Datum { name = "North American Datum 1927 (NAD27) - Canada Yukon"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -7.0+ , dy = 139.0+ , dz = 181.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CubaDatum :: Datum+nad27CubaDatum = Datum { name = "North American Datum 1927 (NAD27) - Cuba"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -9.0+ , dy = 152.0+ , dz = 178.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27AleutianWestDatum :: Datum+nad27AleutianWestDatum = Datum { name = "North American Datum 1927 (NAD27) - Aleutian West"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 2.0+ , dy = 204.0+ , dz = 105.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27AleutianEastDatum :: Datum+nad27AleutianEastDatum = Datum { name = "North American Datum 1927 (NAD27) - Aleutian East"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -2.0+ , dy = 152.0+ , dz = 149.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CanalZoneDatum :: Datum+nad27CanalZoneDatum = Datum { name = "North American Datum 1927 (NAD27) - Canal Zone"+ , ellipsoid = clarke1866Ellipsoid+ , dx = 0.0+ , dy = 125.0+ , dz = 201.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27EasternUSDatum :: Datum+nad27EasternUSDatum = Datum { name = "North American Datum 1927 (NAD27) - Eastern US"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -9.0+ , dy = 161.0+ , dz = 179.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }+++nad27CaribbeanDatum :: Datum+nad27CaribbeanDatum = Datum { name = "North American Datum 1927 (NAD27) - Caribbean"+ , ellipsoid = clarke1866Ellipsoid+ , dx = -3.0+ , dy = 142.0+ , dz = 183.0+ , ds = 0.0+ , rx = 0.0+ , ry = 0.0+ , rz = 0.0+ }
+ src/ECEFRef.hs view
@@ -0,0 +1,62 @@+module ECEFRef where++import Datum+import Ellipsoid+import qualified LatLng as L+import MathExtensions++{- |+ ECEF (earth-centred, earth-fixed) Cartesian co-ordinates are used to define a+ point in three-dimensional space. ECEF co-ordinates are defined relative to+ an x-axis (the intersection of the equatorial plane and the plane defined by+ the prime meridian), a y-axis (at 90° to the x-axis and its intersection+ with the equator) and a z-axis (intersecting the North Pole). All the axes+ intersect at the point defined by the centre of mass of the Earth.+-}++data ECEFRef = ECEFRef { x :: Double+ , y :: Double+ , z :: Double+ , datum :: Datum+ } deriving (Show)+++-- | Create a new earth-centred, earth-fixed reference from the given latitude and longitude.+toECEFRef :: L.LatLng -> ECEFRef+toECEFRef (L.LatLng latitude longitude height datum) = do+ let+ el = ellipsoid datum+ a = semiMajorAxis el+ f = flattening el+ eSquared = f * (2 - f)+ phi = toRadians latitude+ nphi = a / sqrt (1 - eSquared * sinSquared phi)+ lambda = toRadians longitude++ ECEFRef { x = (nphi + height) * cos phi * cos lambda+ , y = (nphi + height) * cos phi * sin lambda+ , z = (nphi * (1 - eSquared) + height) * sin phi+ , datum = datum+ }++-- | Convert this ECEFRef object to a LatLng.+toLatLng :: ECEFRef -> L.LatLng+toLatLng (ECEFRef x y z datum) = do+ let+ el = ellipsoid datum+ a = semiMajorAxis el+ b = semiMinorAxis el+ e2Squared = (a / b) ** 2 - 1+ f = flattening el+ eSquared = f * (2 - f)+ p = sqrt (x ** 2 + y ** 2)+ theta = atan (z * a / (p * b))++ phi = atan $ (z + (e2Squared * b * sinCubed theta ))+ / (p - eSquared * a * cosCubed theta)+ lambda = atan2 y x++ nphi = a / sqrt(1 - eSquared * sinSquared phi )+ h = (p / cos phi ) - nphi++ L.LatLng (toDegrees phi) (toDegrees lambda) h wgs84Datum
+ src/Ellipsoid.hs view
@@ -0,0 +1,115 @@+module Ellipsoid where++-- | To represent a reference ellipsoid.+data Ellipsoid =+ Ellipsoid { semiMajorAxis :: Double -- ^ Semi major axis.+ , semiMinorAxis :: Double -- ^ Semi minor axis.+ , eccentricitySquared :: Double -- ^ Eccentricity squared.+ , flattening :: Double -- ^ Flattening.+ } deriving (Eq)++-- | Parameters to build Ellipsoid.+data EllipsoidParams = SemiMajMinAxises Double Double -- ^ Semi major and semi minor axises.+ | SemiMajAxisEccentricitySquared Double Double -- ^ Semi major axis and eccentric squared.++instance Show Ellipsoid where+ show (Ellipsoid sMajAx sMinAx _ _) =+ "[semi-major axis = " ++ show sMajAx ++ ", semi-minor axis = " ++ show sMinAx ++ "]"+++-- | Create an Ellipsoid with the given parameters.+mkEllipsoid :: EllipsoidParams -- ^ Ellipsoid parameters.+ -> Ellipsoid -- ^ The created Ellipsoid.++mkEllipsoid (SemiMajMinAxises sMajAx sMinAx) =+ Ellipsoid { semiMajorAxis = sMajAx+ , semiMinorAxis = sMinAx+ , eccentricitySquared = 1 - (sMinAx / sMajAx) ** 2 :: Double+ , flattening = 1 - sMinAx / sMajAx+ }++mkEllipsoid (SemiMajAxisEccentricitySquared sMajAx eccSq) =+ Ellipsoid { semiMajorAxis = sMajAx+ , semiMinorAxis = sMinAx+ , eccentricitySquared = eccSq+ , flattening = 1 - sMinAx / sMajAx+ }+ where sMinAx = sqrt (sMajAx * sMajAx * (1 - eccSq))+++-- | Pre-determined ellipsoids:++airy1830Ellipsoid :: Ellipsoid+airy1830Ellipsoid = mkEllipsoid (SemiMajMinAxises 6377563.396 6356256.909)++australianNational1966Ellipsoid :: Ellipsoid+australianNational1966Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378160.0 6356774.719)++bessel1841Ellipsoid :: Ellipsoid+bessel1841Ellipsoid = mkEllipsoid (SemiMajMinAxises 6377397.155 6356078.9629)++clarke1866Ellipsoid :: Ellipsoid+clarke1866Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378206.4 6356583.8)++clarke1880Ellipsoid :: Ellipsoid+clarke1880Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378249.145 6356514.8696)++everest1830Ellipsoid :: Ellipsoid+everest1830Ellipsoid = mkEllipsoid (SemiMajMinAxises 6377276.34518 6356075.41511)++fischer1960Ellipsoid :: Ellipsoid+fischer1960Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378166.0 6356784.284)++fischer1968Ellipsoid :: Ellipsoid+fischer1968Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378150.0 6356768.337)++grs67Ellipsoid :: Ellipsoid+grs67Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378160.0 6356774.51609)++grs75Ellipsoid :: Ellipsoid+grs75Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378140.0 6356755.288)++grs80Ellipsoid :: Ellipsoid+grs80Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378137 6356752.3141)++hayford1910Ellipsoid :: Ellipsoid+hayford1910Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378388.0 6356911.946)++helmert1906Ellipsoid :: Ellipsoid+helmert1906Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378200.0 6356818.17)++hough1956Ellipsoid :: Ellipsoid+hough1956Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378270.0 6356794.34)++iers1989Ellipsoid :: Ellipsoid+iers1989Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378136.0 6356751.302)++internationalEllipsoid :: Ellipsoid+internationalEllipsoid = mkEllipsoid (SemiMajMinAxises 6378388 6356911.9462)++krassovsky1940Ellipsoid :: Ellipsoid+krassovsky1940Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378245.0 6356863.019)++modifiedAiryEllipsoid :: Ellipsoid+modifiedAiryEllipsoid = mkEllipsoid (SemiMajAxisEccentricitySquared 6377340.189 0.00667054015)++modifiedEverestEllipsoid :: Ellipsoid+modifiedEverestEllipsoid = mkEllipsoid (SemiMajMinAxises 6377304.063 6356103.039)++newInternational1967Ellipsoid :: Ellipsoid+newInternational1967Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378157.5 6356772.2)++southAmerican1969Ellipsoid :: Ellipsoid+southAmerican1969Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378160.0 6356774.7192)++wgs60Ellipsoid :: Ellipsoid+wgs60Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378165.0 6356783.287)++wgs66Ellipsoid :: Ellipsoid+wgs66Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378145.0 6356759.770)++wgs72Ellipsoid :: Ellipsoid+wgs72Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378135 6356750.5)++wgs84Ellipsoid :: Ellipsoid+wgs84Ellipsoid = mkEllipsoid (SemiMajMinAxises 6378137 6356752.3142)
+ src/IrishRef.hs view
@@ -0,0 +1,209 @@+{-# LANGUAGE CPP #-}+{-|+ To represent an Irish National Grid reference.++ Projection: Transverse Mercator+ Reference ellipsoid: Modified Airy+ Units: metres+ Origin: 53°30'N, 8°W+ False co-ordinates of origin: 200000m east, 250000m north+-}+module IrishRef where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Data.Char+import Datum+import Ellipsoid+import qualified LatLng as L+import MathExtensions++data IrishRef = IrishRef { easting :: Double -- ^ The easting in metres relative to the origin of the British National Grid.+ , northing :: Double -- ^ The northing in metres relative to the origin of the British National Grid.+ , datum :: Datum+ } deriving (Show)++scaleFactor, falseOriginLatitude, falseOriginLongitude, falseOriginEasting, falseOriginNorthing :: Double+scaleFactor = 1.000035+falseOriginLatitude = 53.5+falseOriginLongitude = -8.0+falseOriginEasting = 200000.0+falseOriginNorthing = 250000.0++{-|+ Create a new Ordnance Survey grid reference using a given easting and+ northing. The easting and northing must be in metres and must be relative+ to the origin of the British National Grid.+-}+mkIrishRef :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 400000.0.+ -> Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than or equal to 500000.0.+ -> Except String IrishRef -- ^ Throws an exception if either the easting or the northing are invalid.+mkIrishRef e n = do+ est <- withExcept (const "Invalid easting") (evalEasting e)+ nrt <- withExcept (const "Invalid northing") (evalNorthing n)+ pure IrishRef { easting = est, northing = nrt, datum = ireland1965Datum }+++{-|+ Take a string formatted as a six-figure OS grid reference (e.g. "TG514131")+ and create a new OSRef object that represents that grid reference. The+ first character must be H, N, S, O or T. The second character can be any+ uppercase character from A through Z excluding I.+-}+mkIrishRef' :: String -- ^ A String representing a six-figure Ordnance Survey grid reference in the form XY123456.+ -> Except String IrishRef -- ^ Throws an exception if ref is not of the form XY123456.+mkIrishRef' ref = do+ -- TODO 2006-02-05 : check format+ let+ east = (read (take 3 $ drop 1 ref) :: Int) * 100+ north = (read (take 3 $ drop 4 ref) :: Int) * 100+ firstCh = ord (ref !! 0)+ ch = if firstCh > 73 then firstCh - 1 else firstCh -- Adjust for no I+ nx = ((ch - 65) `mod` 5) * 100000+ ny = (4 - floor ((fromIntegral ch - 65) / 5)) * 100000++ mkIrishRef (fromIntegral (east + nx)) (fromIntegral (north + ny))+++-- | Create an IrishRef object from the given latitude and longitude.+mkIrishRef'' :: L.LatLng -> Except String IrishRef+mkIrishRef'' (L.LatLng latitude longitude height _) = do+ let+ n0 = falseOriginNorthing+ e0 = falseOriginEasting+ phi0 = toRadians falseOriginLatitude+ lambda0 = toRadians falseOriginLongitude++ el = ellipsoid ireland1965Datum+ a = scaleFactor * semiMajorAxis el+ b = scaleFactor * semiMinorAxis el+ eSquared = eccentricitySquared el++ phi = toRadians latitude+ lambda = toRadians longitude+ n = (a - b) / (a + b)++ va = a * (1 - eSquared * sinSquared phi) ** (-0.5)+ rho = a * (1 - eSquared) * (1 - eSquared * sinSquared phi) ** (-1.5)+ etaSquared = va / rho - 1++ m = b * ((1 + n + 1.25 * n ** 2 + 1.25 * n ** 3) * (phi - phi0)+ - (3 * n + 3 * n ** 2 + 21.0 / 8.0 * n ** 3) * sin (phi - phi0) * cos (phi + phi0)+ + (15.0 / 8.0 * n ** 2 + 15.0 / 8.0 * n ** 3) * sin(2.0 * (phi - phi0)) * cos(2.0 * (phi + phi0))+ - (35.0 / 24.0 * n ** 3) * sin(3.0 * (phi - phi0)) * cos(3.0 * (phi + phi0)))++ i = m + n0+ ii = va / 2.0 * sin phi * cos phi+ iii = va / 24.0 * sin phi * cos phi ** 3 * (5.0 - tanSquared phi + 9.0 * etaSquared)+ iiia = va / 720.0 * sin phi * cos phi ** 5 * (61.0 - 58.0 * tanSquared phi + tan phi ** 4)+ iv = va * cos(phi)+ v = va / 6.0 * cos phi ** 3 * (va / rho - tanSquared phi)+ vi = va / 120.0 * cos phi ** 5.0 * (5.0 - 18.0 * tanSquared phi + tan phi ** 4.0+ + 14 * etaSquared - 58 * tanSquared phi * etaSquared)++ north = i + ii * (lambda - lambda0) ** 2+ + iii * (lambda - lambda0) ** 4+ + iiia * (lambda - lambda0) ** 6++ east = e0 + iv * (lambda - lambda0)+ + v * (lambda - lambda0) ** 3+ + vi * (lambda - lambda0) ** 5++ mkIrishRef east north+++{-|+ Return a String representation of this Irish grid reference using the+ six-figure notation in the form X123456+-}+toSixFigureString :: IrishRef -> String+toSixFigureString (IrishRef easting northing datum) = do+ let+ hundredkmE = floor (easting / 100000)+ hundredkmN = floor (northing / 100000)++ charOffset = 4 - hundredkmN+ i = 65 + 5 * charOffset + hundredkmE+ index = if (i >= 73) then i + 1 else i++ e = floor ((easting - 100000 * fromIntegral hundredkmE) / 100)+ n = floor ((northing - 100000 * fromIntegral hundredkmN) / 100)++ chr index : compose e ++ compose n++ where compose :: Int -> String+ compose x = (if (x < 100) then "0" else "")+ ++ (if (x < 10) then "0" else "")+ ++ show x+++{-|+ Convert this Irish grid reference to a latitude/longitude pair using the+ Ireland 1965 datum. Note that, the LatLng object may need to be converted to the+ WGS84 datum depending on the application.+-}+toLatLng :: IrishRef+ -> Except String L.LatLng -- ^ To represent Irish grid reference using the Ireland 1965 datum.+toLatLng (IrishRef easting northing datum) = do+ let+ n0 = falseOriginNorthing+ e0 = falseOriginEasting+ phi0 = toRadians falseOriginLatitude+ lambda0 = toRadians falseOriginLongitude++ el = ellipsoid datum+ a = semiMajorAxis el+ b = semiMinorAxis el+ eSquared = eccentricitySquared el++ n = (a - b) / (a + b)+ phiPrime = calcPhiPrime ((northing - n0) / (a * scaleFactor) + phi0) a b n phi0 n0 northing++ va = a * scaleFactor * (1 - eSquared * sinSquared phiPrime) ** (-0.5)+ rho = a * scaleFactor * (1 - eSquared) * (1 - eSquared * sinSquared phiPrime) ** (-1.5)+ etaSquared = va / rho - 1++ vii = tan phiPrime / (2 * rho * va)+ viii = tan phiPrime / (24 * rho * va ** 3)+ * (5 + 3 * tanSquared phiPrime + etaSquared - 9 * tanSquared phiPrime * etaSquared)+ ix = tan phiPrime / (720 * rho * va ** 5)+ * (61 + 90 * tanSquared phiPrime + 45 * tanSquared phiPrime ** 2)+ x = sec phiPrime / va+ xi = sec phiPrime / (6 * va ** 3) * (va / rho + 2 * tanSquared phiPrime)+ xii = sec phiPrime / (120 * va ** 5)+ * (5 + 28 * tanSquared phiPrime + 24 * tanSquared phiPrime ** 2)+ xiia = sec phiPrime / (5040 * va ** 7)+ * (61 + 662 * tanSquared phiPrime + 1320 * tanSquared phiPrime ** 2+ + 720 * tanSquared phiPrime ** 3)++ phi = phiPrime - vii * (easting - e0) ** 2 + viii * (easting - e0) ** 4.0 - ix * (easting - e0) ** 6.0+ lambda = lambda0 + x * (easting - e0) - xi * (easting - e0) ** 3.0 + xii * (easting - e0) ** 5.0 - xiia * (easting - e0) ** 7.0+ L.mkLatLng (toDegrees phi) (toDegrees lambda) 0 wgs84Datum+ where calcPhiPrime :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double+ calcPhiPrime phi a b n phi0 n0 northing = do+ let m = b * scaleFactor+ * ((1 + n + 1.25 * n ** 2 + 1.25 * n ** 3) * (phi - phi0)+ - (3 * n + 3 * n ** 2 + 21.0 / 8.0 * n ** 3) * sin (phi - phi0) * cos (phi + phi0)+ + (15.0 / 8.0 * n ** 2 + 15.0 / 8.0 * n ** 3) * sin(2.0 * (phi - phi0)) * cos(2.0 * (phi + phi0))+ - (35.0 / 24.0 * n ** 3) * sin(3.0 * (phi - phi0)) * cos(3.0 * (phi + phi0)))+ delta = northing - n0 - m+ phiPrime = phi + delta / (a * scaleFactor)+ if (delta >= 0.001) then calcPhiPrime phiPrime a b n phi0 n0 northing+ else phiPrime+++-- | Validate the easting.+evalEasting :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 400000.0.+ -> Except String Double -- ^ Throws an Exception if the easting is invalid.+evalEasting e | e < 0.0 || e >= 400000.0 = throwError ("Easting (" ++ show e ++ ") is invalid. Must be greater than or equal to 0.0 and less than 400000.0.")+ | otherwise = pure (e)+++-- | Validate the northing.+evalNorthing :: Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than or equal to 500000.0.+ -> Except String Double -- ^ Throws an Exception if the northing is invalid.+evalNorthing n | n < 0.0 || n > 500000.0 = throwError ("Northing (" ++ show n ++ ") is invalid. Must be greather than or equal to 0.0 and less than or equal to 500000.0.")+ | otherwise = pure (n)
+ src/LatLng.hs view
@@ -0,0 +1,212 @@+{-# LANGUAGE CPP #-}+-- | To represent a latitude/longitude pair based on a particular datum.+module LatLng where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Datum+import Ellipsoid+import MathExtensions+++data LatLng = LatLng { latitude :: Double+ , longitude :: Double+ , height :: Double+ , datum :: Datum+ } deriving (Show)+++mkLatLng :: Double -> Double -> Double -> Datum -> Except String LatLng+mkLatLng lat lng h dtm = do+ lt <- withExcept (const $ "Latitude (" ++ show lat ++ ") is invalid. Must be between -90.0 and 90.0 inclusive.") (validateLatitude lat)+ ln <- withExcept (const $ "Longitude (" ++ show lng ++ ") is invalid. Must be between -180.0 and 180.0 inclusive.") (validateLongitude lng)+ pure LatLng { latitude = lt, longitude = ln, height = h, datum = dtm }+ where+ validateLatitude :: Double -> Except String Double+ validateLatitude l+ | l < -90.0 || l > 90.0 = throwError "Invalid latitude"+ | otherwise = pure l+ validateLongitude :: Double -> Except String Double+ validateLongitude l+ | l < -180.0 || l > 180.0 = throwError "Invalid longitude"+ | otherwise = pure l+++calcPhiN' :: Double -> Double -> Double -> Double -> Double -> Double+calcPhiN' zB eSquared2 a2 p phiN = do+ let v = a2 / sqrt (1 - eSquared2 * sinSquared phiN)+ atan $ (zB + eSquared2 * v * sin phiN) / p+++-- | Convert this LatLng from the OSGB36 datum to the WGS84 datum using an approximate Helmert transformation.+toWGS84 :: LatLng -> LatLng+toWGS84 (LatLng latitude longitude height datum) = do+ let+ a = semiMajorAxis airy1830Ellipsoid+ eSquared = eccentricitySquared airy1830Ellipsoid++ phi = toRadians latitude+ lambda = toRadians longitude+ v = a / sqrt(1 - eSquared * sinSquared phi)+ x = v * cos phi * cos lambda+ y = v * cos phi * sin lambda+ z = (1 - eSquared) * v * sin phi++ tx = 446.448 :: Double+ -- ty : Incorrect value in v1.0 (-124.157). Corrected in v1.1.+ ty = -125.157 :: Double+ tz = 542.060 :: Double+ s = -0.0000204894 :: Double+ rx = toRadians (0.00004172222) :: Double+ ry = toRadians (0.00006861111) :: Double+ rz = toRadians (0.00023391666) :: Double++ xB = tx + x * (1 + s) + (-rx) * y + ry * z+ yB = ty + rz * x + y * (1 + s) + (-rx) * z+ zB = tz + (-ry) * x + rx * y + z * (1 + s)++ a2 = semiMajorAxis wgs84Ellipsoid+ eSquared2 = eccentricitySquared wgs84Ellipsoid++ p = sqrt (xB ** 2 + yB ** 2)+ phiN = atan(zB / (p * (1 - eSquared2)))++ calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN+ LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum+++toDatum :: LatLng -> Datum -> LatLng+toDatum ll@(LatLng latitude longitude height datum) newDatum+ | datum /= wgs84Datum && newDatum /= wgs84Datum = mkToDatum ll newDatum 1+ | datum == wgs84Datum && newDatum /= wgs84Datum = mkToDatum ll newDatum (-1)+ | datum /= wgs84Datum && newDatum == wgs84Datum = ll+ | otherwise = mkToDatum ll newDatum 1+ where mkToDatum :: LatLng -> Datum -> Double -> LatLng+ mkToDatum ll@(LatLng latitude longitude height datum) newDatum invert = do+ let+ a = semiMajorAxis $ ellipsoid $ datum+ eSquared = eccentricitySquared $ ellipsoid $ datum++ phi = toRadians latitude+ lambda = toRadians longitude+ v = a / sqrt(1 - eSquared * sinSquared phi)+ x = (v + height) * cos phi * cos lambda+ y = (v + height) * cos phi * sin lambda+ z = ((1 - eSquared) * v + height) * sin phi++ dx2 = invert * dx newDatum -- 446.448+ dy2 = invert * dy newDatum -- -125.157+ dz2 = invert * dz newDatum -- 542.060+ ds2 = invert * ds newDatum / 1000000.0 -- -0.0000204894+ rx2 = invert * toRadians (rx newDatum / 3600.0) -- toRadians(0.00004172222)+ ry2 = invert * toRadians (ry newDatum / 3600.0) -- toRadians(0.00006861111)+ rz2 = invert * toRadians (rz newDatum / 3600.0) -- toRadians(0.00023391666)++ sc = 1 + ds2+ xB = dx2 + sc * (x + (-rx2) * y + ry2 * z)+ yB = dy2 + sc * (rz2 * x + y + (-rx2) * z)+ zB = dz2 + sc * ((-ry2) * x + rx2 * y + z)++ a2 = semiMajorAxis $ ellipsoid $ newDatum+ eSquared2 = eccentricitySquared $ ellipsoid $ newDatum++ p = sqrt(xB ** 2 + yB ** 2)+ phiN = atan (zB / (p * (1 - eSquared2)))++ calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN+ LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum+++-- | Convert this LatLng from the WGS84 datum to the OSGB36 datum using an approximate Helmert transformation.+toOSGB36 :: LatLng -> LatLng+toOSGB36 (LatLng latitude longitude height datum) = do+ let+ a = semiMajorAxis wgs84Ellipsoid+ eSquared = eccentricitySquared wgs84Ellipsoid++ phi = toRadians latitude+ lambda = toRadians longitude+ v = a / sqrt(1 - eSquared * sinSquared phi)+ x = v * cos phi * cos lambda+ y = v * cos phi * sin lambda+ z = (1 - eSquared) * v * sin phi++ tx = -446.448 :: Double+ -- ty : Incorrect value in v1.0 (-124.157). Corrected in v1.1.+ ty = 125.157 :: Double+ tz = -542.060 :: Double+ s = 0.0000204894 :: Double+ rx = toRadians (-0.00004172222) :: Double+ ry = toRadians (-0.00006861111) :: Double+ rz = toRadians (-0.00023391666) :: Double++ xB = tx + x * (1 + s) + (-rx) * y + ry * z+ yB = ty + rz * x + y * (1 + s) + (-rx) * z+ zB = tz + (-ry) * x + rx * y + z * (1 + s)++ a2 = semiMajorAxis airy1830Ellipsoid+ eSquared2 = eccentricitySquared airy1830Ellipsoid++ p = sqrt (xB ** 2 + yB ** 2)+ phiN = atan(zB / (p * (1 - eSquared2)))++ calcPhiN = last $ take 10 $ iterate (calcPhiN' zB eSquared2 a2 p) phiN+ LatLng (toDegrees calcPhiN) (toDegrees $ atan (yB / xB)) height datum+++-- | Calculate the surface distance in kilometres from this LatLngPoint to the given LatLngPoint.+distance :: LatLng -> LatLng -> Double+distance (LatLng latitudeFrom longitudeFrom _ _) (LatLng latitudeTo longitudeTo _ _) = do+ let+ er = 6366.707 :: Double+ latFrom = toRadians latitudeFrom+ latTo = toRadians latitudeTo+ lngFrom = toRadians longitudeFrom+ lngTo = toRadians longitudeTo+ acos (sin latFrom * sin latTo + cos latFrom * cos latTo * cos (lngTo - lngFrom)) * er+++-- | Calculate the surface distance in miles from this LatLngPoint to the given LatLngPoint.+distanceMiles :: LatLng -> LatLng -> Double+distanceMiles from to = (distance from to) / 1.609344+++latitudeDegrees, longitudeDegrees :: LatLng -> Int+latitudeDegrees (LatLng l _ _ _) = calcDegrees l+longitudeDegrees (LatLng _ l _ _) = calcDegrees l++calcDegrees :: Double -> Int+calcDegrees l = do+ let deg = floor l+ minx = (l - fromIntegral deg) :: Double+ case l of _ | l < 0 && minx /= 0 -> (deg :: Int) + 1+ | otherwise -> deg :: Int+++latitudeMinutes, longitudeMinutes :: LatLng -> Int+latitudeMinutes (LatLng l _ _ _) = calcMinutes l+longitudeMinutes (LatLng _ l _ _) = calcMinutes l++calcMinutes :: Double -> Int+calcMinutes l = do+ let deg = floor l+ minx = (l - fromIntegral deg) :: Double+ min | l < 0 && minx /= 0 = 1 - minx+ | otherwise = minx+ floor (60 * min) :: Int+++latitudeSeconds, longitudeSeconds :: LatLng -> Double+latitudeSeconds (LatLng l _ _ _) = calcSeconds l+longitudeSeconds (LatLng _ l _ _) = calcSeconds l++calcSeconds :: Double -> Double+calcSeconds l = do+ let deg = floor l+ minx = (l - fromIntegral deg) :: Double+ min | l < 0 && minx /= 0 = 1 - minx+ | otherwise = minx+ 60 * (60 * min - fromIntegral (floor (60 * min)))
+ src/MGRSRef.hs view
@@ -0,0 +1,295 @@+{-# LANGUAGE CPP #-}+{-|+ To represent a Military Grid Reference System (MGRS) reference.++ Military Grid Reference System (MGRS)++ The Military Grid Reference System (MGRS) is an extension of the Universal+ Transverse Mercator (UTM) reference system. An MGRS reference is made from 5+ parts:++ UTM Longitude Zone+ This is a number indicating which UTM longitude zone the reference falls+ into. Zones are numbered from 1 (starting at 180°W) through 60. Each zone+ is 6° wide.++ UTM Latitude Zone+ Latitude is split into regions that are 8° high, starting at 80°S.+ Latitude zones are lettered using C through X, but omitting I and O as they+ can easily be confused with the numbers 1 and 0.++ 100,000m Square identification+ Each UTM zone is treated as a square 100,000m to a side. The 50,000m easting+ is centred on the centre-point of the UTM zone. 100,000m squares are+ identified using two characters - one to identify the row and one to identify+ the column.+ Row identifiers use the characters A through V (omitting I and O again). The+ sequence is repeated every 2,000,000m from the equator. If the UTM longitude+ zone is odd, then the lettering is advanced by five characters to start at F.+ Column identifiers use the characters A through Z (again omitting I and O).++ Easting and northing+ Each 100,000m grid square is further divided into smaller squares+ representing 1m, 10m, 100m, 1,000m and 10,000m precision. The easting and+ northing are given using the numeric row and column reference of the square,+ starting at the bottom-left corner of the square.++ MGRS Reference Example+ 18SUU8362601432 is an example of an MGRS reference. '18' is the UTM longitude+ zone, 'S' is the UTM latitude zone, 'UU' is the 100,000m square+ identification, 83626 is the easting reference to 1m precision and 01432 is+ the northing reference to 1m precision.++ MGRSRef+ Methods are provided to query an MGRSRef object for its+ parameters. As MGRS references are related to UTM references, a+ {@link MGRSRef#toUTMRef() toUTMRef()} method is provided to+ convert an MGRSRef object into a {@link UTMRef} + object. The reverse conversion can be made using the+ {@link #MGRSRef(UTMRef) MGRSRef(UTMRef)} constructor.+ MGRSRef objects can be converted to {@link LatLng LatLng}+ objects using the {@link MGRSRef#toLatLng() toLatLng()} method. The reverse+ conversion is made using the {@link LatLng#toMGRSRef() LatLng.toMGRSRef()} method.+ Some MGRS references use the Bessel 1841 ellipsoid rather than the Geodetic+ Reference System 1980 (GRS 1980), International or World Geodetic System 1984+ (WGS84) ellipsoids. Use the constructors with the optional boolean parameter+ to be able to specify whether your MGRS reference uses the Bessel 1841+ ellipsoid. Note that no automatic determination of the correct ellipsoid to+ use is made.++ Important note: There is currently no support for MGRS references in+ polar regions north of 84°N and south of 80°S. There is also no+ account made for UTM zones with slightly different sizes to normal around+ Svalbard and Norway.+-}+module MGRSRef where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Data.Char+import Data.Fixed+import Text.Regex.PCRE++import Datum+import qualified UTMRef as U++data MGRSRef = MGRSRef { easting :: Int+ , northing :: Int+ , eastingId :: Char+ , northingId :: Char+ , utmZoneNumber :: Int+ , utmZoneChar :: Char+ , precision :: Precision+ , isBessel :: Bool+ , datum :: Datum+ } deriving (Eq)+++data Precision = M1 -- ^ precision of 1m+ | M10 -- ^ precision of 10m+ | M100 -- ^ precision of 100m+ | M1000 -- ^ precision of 1000m (1km)+ | M10000 -- ^ precision of 10000m (10km)+ deriving (Eq)+++instance Enum Precision where+ fromEnum M1 = 1+ fromEnum M10 = 10+ fromEnum M100 = 100+ fromEnum M1000 = 1000+ fromEnum M10000 = 10000+ toEnum 1 = M1+ toEnum 10 = M10+ toEnum 100 = M100+ toEnum 1000 = M1000+ toEnum 10000 = M10000+++-- | Northing characters+northingIds :: [Char]+northingIds = [ 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'J', 'K'+ , 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'U', 'V' ]+++instance Show MGRSRef where+ show m@(MGRSRef _ _ _ _ _ _ precision _ _) = showWithPrecision m precision++-- | Return a String representation of this MGRS reference to 1m, 10m, 100m, 1000m or 10000m precision.+showWithPrecision :: MGRSRef+ -> Precision -- ^ One of MGRSRef.PRECISION_1M, MGRSRef.PRECISION_10M, MGRSRef.PRECISION_100M, MGRSRef.PRECISION_1000M, MGRSRef.PRECISION_10000M.+ -> String -- ^ ~ String representation of this MGRS reference to the required precision.+showWithPrecision (MGRSRef easting northing eastingId northingId utmZoneNumber utmZoneChar _ _ _) precision = do+ let+ p = fromEnum precision+ eastingR = easting `div` p+ northingR = northing `div` p+ padding = case precision of+ M1 -> 5+ M10 -> 4+ M100 -> 3+ M1000 -> 2+ M10000 -> 1+ eastingRs = padWithZeros (show eastingR) padding+ northingRs = padWithZeros (show northingR) padding+ utmZonePadding = if (utmZoneNumber < 10) then "0" else ""+ utmZonePadding ++ show utmZoneNumber ++ [utmZoneChar] ++ [eastingId] ++ [northingId] ++ eastingRs ++ northingRs+ where padWithZeros :: String -> Int -> String+ padWithZeros text padding =+ if (padding - length text > 0) then padWithZeros ('0' : text) padding+ else text+++{-|+ Create a new MGRS reference object from the given UTM reference. It is+ assumed that this MGRS reference represents a point using the GRS 1980,+ International or WGS84 ellipsoids. It is assumed that the UTMRef object is+ valid.+-}+toMGRSRef :: U.UTMRef -- ^ UTM reference.+ -> Bool -- ^ True if the parameters represent an MGRS reference using the Bessel 1841 ellipsoid; false is the parameters represent an MGRS reference using the GRS 1980, International or WGS84 ellipsoids.+ -> MGRSRef+toMGRSRef (U.UTMRef easting northing latZone lngZone _) isBessel = do+ let+ set = (lngZone - 1) `mod` 6 + 1++ eId = (\x -> x + (if x > 14 then 1 else 0))+ . (\x -> x + (if x > 8 then 1 else 0))+ $ (floor (easting / 100000.0) + 8 * ((set - 1) `mod` 3))+ eIdC = chr $ eId + 64++ nId = (\x -> x - (if x > 19 then 20 else 0))+ . (\x -> x + (if isBessel then 10 else 0))+ . (\x -> x + (if even set then 5 else 0))+ $ floor ((northing `mod'` 2000000) / 100000.0)+ nIdC = northingIds !! nId++ MGRSRef (round(easting) `mod` 100000) (round(northing) `mod` 100000) eIdC nIdC lngZone latZone M1 isBessel wgs84Datum+++mkMGRSRef :: Int -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 100000.+ -> Int -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than or equal to 500000.0.+ -> Char -- ^ The character representing the 100,000km easting square.+ -> Char -- ^ The character representing the 100,000km northing square.+ -> Int -- ^ The UTM zone number representing the longitude.+ -> Char -- ^ The UTM zone character representing the latitude.+ -> Precision -- ^ The precision of the given easting and northing.+ -> Bool -- ^ True if the parameters represent an MGRS reference using the Bessel 1841 ellipsoid; False is the parameters represent an MGRS reference using the GRS 1980, International or WGS84 ellipsoids.+ -> Except String MGRSRef -- ^ Throws an exception if any of the given parameters are invalid. Note that the parameters are only checked for the range of values that they can take on. Being able to create an MGRSRef object does not necessarily imply that the reference is valid.+mkMGRSRef e n eId nId un uc p b = do+ est <- withExcept (const "Invalid easting") (evalEasting e)+ nrt <- withExcept (const "Invalid northing") (evalNorthing n)+ estId <- withExcept (const "Invalid eastingId") (evalEastingId eId)+ nrtId <- withExcept (const "Invalid northingId") (evalNorthingId nId)+ uzn <- withExcept (const "Invalid utmZoneNumber") (evalUtmZoneNumber un)+ uzc <- withExcept (const "Invalid utmZoneChar") (evalUtmZoneChar uc)+ pure MGRSRef { easting = est+ , northing = nrt+ , eastingId = estId+ , northingId = nrtId+ , utmZoneNumber = uzn+ , utmZoneChar = uzc+ , precision = p+ , isBessel = b+ , datum = wgs84Datum }+++ where evalEasting :: Int -> Except String Int+ evalEasting e | e < 0 || e > 99999 = throwError ("Invalid easting (" ++ show e ++ ")")+ | otherwise = pure (e)++ evalNorthing :: Int -> Except String Int+ evalNorthing n | n < 0 || n > 99999 = throwError ("Invalid northing (" ++ show n ++ ")")+ | otherwise = pure (n)++ evalUtmZoneNumber :: Int -> Except String Int+ evalUtmZoneNumber u | u < 1 || u > 60 = throwError ("Invalid utmZoneNumber (" ++ show u ++ ")")+ | otherwise = pure (u)++ evalUtmZoneChar :: Char -> Except String Char+ evalUtmZoneChar u | u < 'A' || u > 'Z' = throwError ("Invalid utmZoneChar (" ++ show u ++ ")")+ | otherwise = pure (u)++ evalEastingId :: Char -> Except String Char+ evalEastingId e | e < 'A' || e > 'Z' || e == 'I' || e == 'O' = throwError ("Invalid eastingId (" ++ show e ++ ")")+ | otherwise = pure (e)++ evalNorthingId :: Char -> Except String Char+ evalNorthingId n | n < 'A' || n > 'Z' || n == 'I' || n == 'O' = throwError ("Invalid northingId (" ++ show n ++ ")")+ | otherwise = pure (n)+++{-|+ Create a new MGRS reference object from the given String. Must be correctly+ formatted otherwise an IllegalArgumentException will be thrown.+-}+mkMGRSRef' :: String -- ^ A String to create an MGRS reference from.+ -> Except String MGRSRef -- ^ Throws an exception if the given String is not correctly formatted.+mkMGRSRef' ref = do+ groups <- withExcept (const "Invalid easting") (evalRefMatch ref)+ let+ group = head groups+ en = group !! 5+ enlh = length en `div` 2+ p = 10 ^ (5 - enlh)+ pure MGRSRef { easting = p * (read (take enlh en) :: Int)+ , northing = p * (read (drop enlh en) :: Int)+ , eastingId = (group !! 3) !! 0+ , northingId = (group !! 4) !! 0+ , utmZoneNumber = read (group !! 1) :: Int+ , utmZoneChar = (group !! 2) !! 0+ , precision = toEnum p+ , isBessel = False+ , datum = wgs84Datum }+++ where evalRefMatch :: String -> Except String [[String]]+ evalRefMatch ref = do+ let groups = ref =~ "(\\d{1,2})((?![IO])[C-X])((?![IO])[A-Z])((?![IO])[A-Z])(\\d{2,10})" :: [[String]]+ case groups of+ [] -> throwError ("No matches")+ [_,_,_,_,_,en]+ | odd $ length en -> throwError ("Invalid ref")+ | otherwise -> pure groups+ _ -> pure groups+++-- | Convert this MGRS reference to an equivelent UTM reference.+toUTMRef :: MGRSRef -> Except String U.UTMRef+toUTMRef (MGRSRef easting northing eastingId northingId utmZoneNumber utmZoneChar _ _ _) = do+ let+ set = (utmZoneNumber - 1) `mod` 6 + 1+ isOffset = even set++ e = (\x -> if (x > 9) then x - 1 else x)+ . (\x -> if (x > 15) then x - 1 else x)+ $ ord eastingId - 65+ ex = (easting + ((e `mod` 8 + 1) * 100000)) `mod` 1000000++ n = (\x -> if (x < 0) then x + 16 else x)+ . (\x -> if isOffset then x - 5 else x)+ . (\x -> if (x > 9) then x - 1 else x)+ . (\x -> if (x > 15) then x - 1 else x)+ $ ord northingId - 64++ nn = case utmZoneChar of+ 'Q'+ | (not(isOffset) && northingId < 'T') || (isOffset && (northingId < 'C' || northingId > 'E')) -> 2000000+ | otherwise -> 0+ 'R' -> 2000000+ 'S'+ | (not(isOffset) && northingId < 'R') || (isOffset && northingId > 'E') -> 4000000+ | otherwise -> 2000000+ 'T' -> 4000000+ 'U'+ | (not(isOffset) && northingId < 'P') || (isOffset && northingId > 'U') -> 6000000+ | otherwise -> 4000000+ 'V' -> 6000000+ 'W' -> 6000000+ 'X' -> 8000000++ nx = nn + (100000 * (n - 1)) + northing+ U.mkUTMRef (fromIntegral ex) (fromIntegral nx) utmZoneChar utmZoneNumber
+ src/MathExtensions.hs view
@@ -0,0 +1,21 @@+module MathExtensions where++toRadians :: Double -> Double+toRadians d = d / 180 * pi++toDegrees :: Double -> Double+toDegrees a = a * 180.0 / pi++sinSquared, sinCubed :: Double -> Double+sinSquared phi = sin phi ** 2+sinCubed phi = sin phi ** 3++cosSquared, cosCubed :: Double -> Double+cosSquared phi = cos phi ** 2+cosCubed phi = cos phi ** 3++tanSquared :: Double -> Double+tanSquared phi = tan phi ** 2++sec :: Double -> Double+sec x = 1 / cos x
+ src/OSRef.hs view
@@ -0,0 +1,395 @@+{-# LANGUAGE CPP #-}+{-|+ To represent an Ordnance Survey of Great Britain (OSGB) grid reference.++ British National Grid+ Projection: Transverse Mercator+ Reference ellipsoid: Airy 1830+ Units: metres+ Origin: 49°N, 2°W+ False co-ordinates of origin: 400000m east, -100000m north+ + A full reference includes a two-character code identifying a particular+ 100,000m grid square. The table below shows how the two-character 100,000m+ grid squares are identified. The bottom left corner is at the false origin of+ the grid. Squares without values fall outside the boundaries of the British+ National Grid.++ km+ 0+ 100+ 200+ 300+ 400+ 500+ 600+ 700+ + + 1200+ HL+ HM+ HN+ HO+ HP+ JL+ JM+ + + + 1100+ HQ+ HR+ HS+ HT+ HU+ JQ+ JR+ + + + 1000+ HV+ HW+ HX+ HY+ HZ+ JV+ JW+ + + + 900+ NA+ NB+ NC+ ND+ NE+ OA+ OB+ + + + 800+ NF+ NG+ NH+ NJ+ NK+ OF+ OG+ OH+ + + 700+ NL+ NM+ NN+ NO+ NP+ OL+ OM+ ON+ + + 600+ NQ+ NR+ NS+ NT+ NU+ OQ+ OR+ OS+ + + 500+ + NW+ NX+ NY+ NZ+ OV+ OW+ OX+ + + 400+ + SB+ SC+ SD+ SE+ TA+ TB+ TC+ + + 300+ + SG+ SH+ SJ+ SK+ TF+ TG+ TH+ + + 200+ + SM+ SN+ SO+ SP+ TL+ TM+ TN+ + + 100+ SQ+ SR+ SS+ ST+ SU+ TQ+ TR+ TS+ + + 0+ SV+ SW+ SX+ SY+ SZ+ TV+ TW+ + + Within each 100,000m square, the grid is further subdivided into 1000m+ squares. These 1km squares are shown on Ordnance Survey 1:25000 and 1:50000+ mapping as the main grid. To reference a 1km square, give the easting and+ then the northing, e.g. TR2266. In this example, TR represents the 100,000m+ square, 22 represents the easting (in km) and 66 represents the northing (in+ km). This is commonly called a four-figure grid reference.+ + When providing local references, the 2 characters representing the 100,000m+ square are often omitted.+-}+module OSRef where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Data.Char+import Datum+import Ellipsoid+import qualified LatLng as L+import MathExtensions++data OSRef = OSRef { easting :: Double -- ^ The easting in metres relative to the origin of the British National Grid.+ , northing :: Double -- ^ The northing in metres relative to the origin of the British National Grid.+ , datum :: Datum+ } deriving (Eq, Show)++data Precision = SixDigits -- ^ A six-figure representation this OSGB grid reference i.e XY123456+ | EightDigits -- ^ A eight-figure representation this OSGB grid reference i.e XY12345678+++scaleFactor, falseOriginLatitude, falseOriginLongitude, falseOriginEasting, falseOriginNorthing :: Double+scaleFactor = 0.9996012717 -- OSGB_F0+falseOriginLatitude = 49.0+falseOriginLongitude = -2.0+falseOriginEasting = 400000.0+falseOriginNorthing = -100000.0+++{-|+ Create a new Ordnance Survey grid reference using a given easting and+ northing. The easting and northing must be in metres and must be relative+ to the origin of the British National Grid.+-}+mkOSRef :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 800000.0.+ -> Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than 1400000.0.+ -> Except String OSRef -- ^ Throws an exception if either the easting or the northing are invalid.+mkOSRef e n = do+ est <- withExcept (const "Invalid easting") (evalEasting e)+ nrt <- withExcept (const "Invalid northing") (evalNorthing n)+ pure OSRef { easting = est, northing = nrt, datum = osgb36Datum }+++{-|+ Take a string formatted as a six-figure OS grid reference (e.g. "TG514131")+ and create a new OSRef object that represents that grid reference. The+ first character must be H, N, S, O or T. The second character can be any+ uppercase character from A through Z excluding I.+-}+mkOSRef' :: String -- ^ a String representing a six-figure Ordnance Survey grid reference in the form XY123456+ -> Except String OSRef -- ^ Throws an exception if ref is not of the form XY123456.+mkOSRef' ref = do+ -- TODO 2006-02-05 : check format+ let coords = findCoords (fromLabel ref 2) (fromLabel ref 5) (ref !! 0)+ where+ fromLabel :: String -> Int -> Int+ fromLabel r p = (read (take 3 $ drop p r) :: Int) * 100++ findCoords :: Int -> Int -> Char -> (Int, Int)+ findCoords e n c+ | c == 'H' = (e, n + 1000000)+ | c == 'N' = (e, n + 500000)+ | c == 'O' = (e + 500000, n + 500000)+ | c == 'T' = (e + 500000, n)++ let char2Ord = ord $ ref !! 1+ let c2 = if (char2Ord > 73) then char2Ord - 66 else char2Ord - 65++ let nx = c2 `mod` 5 * 100000+ let ny = (4 - floor ((fromIntegral c2 :: Double) / 5)) * 100000++ mkOSRef (fromIntegral (fst coords + nx)) (fromIntegral (snd coords + ny))+++{-|+ Convert latitude and longitude into an OSGB (Ordnance Survey of Great Britain) grid reference.+-}+toOSRef :: L.LatLng -> Except String OSRef+toOSRef (L.LatLng latitude longitude _ _) = do+ let osgb_f0 = 0.9996012717 :: Double+ n0 = -100000.0 :: Double+ e0 = 400000.0 :: Double+ phi0 = toRadians 49.0+ lambda0 = toRadians (-2.0)++ a = semiMajorAxis airy1830Ellipsoid+ b = semiMinorAxis airy1830Ellipsoid+ eSquared = eccentricitySquared airy1830Ellipsoid++ phi = toRadians latitude+ lambda = toRadians longitude++ n = (a - b) / (a + b)+ vc = a * osgb_f0 * (1.0 - eSquared * sinSquared phi) ** (-0.5)+ rho = a * osgb_f0 * (1.0 - eSquared) * (1.0 - eSquared * sinSquared phi) ** (-1.5)+ etaSquared = vc / rho - 1.0++ cp = cos phi+ sp = sin phi+ ts = tan phi ** 2++ m = (b * osgb_f0) * (((1.0 + n + (1.25 * n * n) + (1.25 * n * n * n)) * (phi - phi0))+ - (((3 * n) + (3 * n * n) + (2.625 * n * n * n)) * sin(phi - phi0) * cos(phi + phi0))+ + (((1.875 * n * n) + (1.875 * n * n * n)) * sin(2.0 * (phi - phi0)) * cos(2.0 * (phi + phi0)))+ - (((35.0 / 24.0) * n * n * n) * sin(3.0 * (phi - phi0)) * cos(3.0 * (phi + phi0))))+ i = m + n0+ ii = (vc / 2.0) * sp * cp+ iii = (vc / 24.0) * sp * (cp ** 3) * (5.0 - ts + 9.0 * etaSquared)+ iiia = (vc / 720.0) * sp * (cp ** 5) * (61.0 - 58.0 * ts + ts ** 2)+ iv = vc * cp+ v = (vc / 6.0) * (cp ** 3) * (vc / rho - ts)+ vi = (vc / 120.0) * (cp ** 5.0) * (5.0 - 18.0 * ts + ts ** 2 + 14 * etaSquared - 58 * ts * etaSquared)++ mkOSRef+ (e0 + iv * (lambda - lambda0) + v * (lambda - lambda0) ** 3 + vi * (lambda - lambda0) ** 5)+ (i + ii * (lambda - lambda0) ** 2 + iii * (lambda - lambda0) ** 4 + iiia * (lambda - lambda0) ** 6)+++getOsRefWithPrecisionOf :: Precision -> OSRef -> String+getOsRefWithPrecisionOf SixDigits (OSRef e n _) = evalOsRef 100 e n+getOsRefWithPrecisionOf EightDigits (OSRef e n _) = evalOsRef 10 e n++evalOsRef :: Double -> Double -> Double -> String+evalOsRef precision easting northing = do+ let+ hundredkmE = floor (easting / 100000)+ hundredkmN = floor (northing / 100000)++ firstLetter+ | hundredkmN < 5 = if (hundredkmE < 5) then 'S' else 'T'+ | hundredkmN < 10 = if (hundredkmE < 5) then 'N' else 'O'+ | otherwise = 'H'++ i = 85 - 5 * (hundredkmN `mod` 5) + (hundredkmE `mod` 5)+ secondLetter = chr $ if (i >= 73) then i + 1 else i++ e = floor ((easting - 100000 * fromIntegral hundredkmE) / precision)+ n = floor ((northing - 100000 * fromIntegral hundredkmN) / precision)++ firstLetter : secondLetter : compose e ++ compose n++ where compose :: Int -> String+ compose x = (if (x < 100) then "0" else "")+ ++ (if (x < 10) then "0" else "")+ ++ show x+++{-|+ Convert this OSGB grid reference to a latitude/longitude pair using the+ OSGB36 datum. Note that, the LatLng object may need to be converted to the+ WGS84 datum depending on the application.+-}+toLatLng :: OSRef+ -> Except String L.LatLng -- ^ To represent OSGB grid reference using the OSGB36 datum.+toLatLng (OSRef easting northing datum) = do+ let+ n0 = falseOriginNorthing+ e0 = falseOriginEasting+ phi0 = toRadians falseOriginLatitude+ lambda0 = toRadians falseOriginLongitude++ el = ellipsoid datum+ a = semiMajorAxis el+ b = semiMinorAxis el+ eSquared = eccentricitySquared el++ n = (a - b) / (a + b)+ phiPrime = calcPhiPrime ((northing - n0) / (a * scaleFactor) + phi0) a b n phi0 n0 northing++ va = a * scaleFactor * (1 - eSquared * sinSquared phiPrime) ** (-0.5)+ rho = a * scaleFactor * (1 - eSquared) * (1 - eSquared * sinSquared phiPrime) ** (-1.5)+ etaSquared = va / rho - 1++ vii = tan phiPrime / (2 * rho * va)+ viii = tan phiPrime / (24 * rho * va ** 3)+ * (5 + 3 * tanSquared phiPrime + etaSquared - 9 * tanSquared phiPrime * etaSquared)+ ix = tan phiPrime / (720 * rho * va ** 5)+ * (61 + 90 * tanSquared phiPrime + 45 * tanSquared phiPrime ** 2)+ x = sec phiPrime / va+ xi = sec phiPrime / (6 * va ** 3) * (va / rho + 2 * tanSquared phiPrime)+ xii = sec phiPrime / (120 * va ** 5)+ * (5 + 28 * tanSquared phiPrime + 24 * tanSquared phiPrime ** 2)+ xiia = sec phiPrime / (5040 * va ** 7)+ * (61 + 662 * tanSquared phiPrime + 1320 * tanSquared phiPrime ** 2+ + 720 * tanSquared phiPrime ** 3)++ phi = phiPrime - vii * (easting - e0) ** 2 + viii * (easting - e0) ** 4.0 - ix * (easting - e0) ** 6.0+ lambda = lambda0 + x * (easting - e0) - xi * (easting - e0) ** 3.0 + xii * (easting - e0) ** 5.0 - xiia * (easting - e0) ** 7.0+ L.mkLatLng (toDegrees phi) (toDegrees lambda) 0 wgs84Datum+ where calcPhiPrime :: Double -> Double -> Double -> Double -> Double -> Double -> Double -> Double+ calcPhiPrime phi a b n phi0 n0 northing = do+ let m = b * scaleFactor+ * ((1 + n + 1.25 * n ** 2 + 1.25 * n ** 3) * (phi - phi0)+ - (3 * n + 3 * n ** 2 + 21.0 / 8.0 * n ** 3) * sin (phi - phi0) * cos (phi + phi0)+ + (15.0 / 8.0 * n ** 2 + 15.0 / 8.0 * n ** 3) * sin(2.0 * (phi - phi0)) * cos(2.0 * (phi + phi0))+ - (35.0 / 24.0 * n ** 3) * sin(3.0 * (phi - phi0)) * cos(3.0 * (phi + phi0)))+ delta = northing - n0 - m+ phiPrime = phi + delta / (a * scaleFactor)+ if (delta >= 0.001) then calcPhiPrime phiPrime a b n phi0 n0 northing+ else phiPrime+++-- | Validate the easting.+evalEasting :: Double -- ^ The easting in metres. Must be greater than or equal to 0.0 and less than 800000.0.+ -> Except String Double -- ^ Throws an Exception if the easting is invalid.+evalEasting e | e < 0.0 || e >= 800000.0 = throwError ("Easting (" ++ show e ++ ") is invalid. Must be greater than or equal to 0.0 and less than 800000.0.")+ | otherwise = pure (e)+++-- | Validate the northing.+evalNorthing :: Double -- ^ The northing in metres. Must be greater than or equal to 0.0 and less than 1400000.0.+ -> Except String Double -- ^ Throws an Exception if the northing is invalid.+evalNorthing n | n < 0.0 || n >= 1400000.0 = throwError ("Northing (" ++ show n ++ ") is invalid. Must be greather than or equal to 0.0 and less than or equal to 1400000.0.")+ | otherwise = pure (n)
+ src/UTMRef.hs view
@@ -0,0 +1,177 @@+{-# LANGUAGE CPP #-}+-- | To represent a Universal Transverse Mercator (UTM) reference.+module UTMRef where++#if __GLASGOW_HASKELL__ < 710+import Control.Applicative+#endif++import Control.Monad.Except+import Datum+import Ellipsoid+import LatLng+import MathExtensions++data UTMRef = UTMRef { easting :: Double -- ^ Easting+ , northing :: Double -- ^ Northing+ , latZone :: Char -- ^ Latitude zone character+ , lngZone :: Int -- ^ Longitude zone number+ , datum :: Datum+ } deriving (Show)++{-|+ Create a new UTM reference object. Checks are made to make sure that the+ given parameters are roughly valid, but the checks are not exhaustive with+ regards to the easting value. Catching an exception does not necessarily+ mean that the UTM reference is well-formed. This is because that valid+ values for the easting vary depending on the latitude.+-}+mkUTMRef :: Double -- ^ The easting in metres+ -> Double -- ^ The northing in metres+ -> Char -- ^ The latitude zone character+ -> Int -- ^ The longitude zone number+ -> Except String UTMRef -- ^ If any of the parameters are invalid. Be careful that a valid+ -- value for the easting does not necessarily mean that the UTM+ -- reference is well-formed. The current checks do not take into+ -- account the varying range of valid values for the easting for+ -- different latitudes.+mkUTMRef easting northing latZone lngZone+ | lngZone < 1 || lngZone > 60 = throwError $ "Longitude zone (" ++ show lngZone ++ ") is not defined on the UTM grid."+ | latZone < 'C' || latZone > 'X' = throwError $ "Latitude zone (" ++ show latZone ++ ") is not defined on the UTM grid."+ | easting < 0.0 || easting > 1000000.0 = throwError $ "Easting (" ++ show easting ++ ") is not defined on the UTM grid."+ | northing < 0.0 || northing > 10000000.0 = throwError $ "Northing (" ++ show northing ++ ") is not defined on the UTM grid."+ | otherwise = pure $ UTMRef easting northing latZone lngZone wgs84Datum+++-- | Convert this UTM reference to a latitude and longitude.+toLatLng :: UTMRef -> Except String LatLng+toLatLng (UTMRef east north ltz lnz datum) = do+ let+ utm_f0 = 0.9996 :: Double+ el = ellipsoid datum+ a = semiMajorAxis el+ eSquared = eccentricitySquared el++ ePrimeSquared = eSquared / (1 - eSquared)+ es1 = sqrt (1 - eSquared)+ e1 = (1 - es1) / (1 + es1)+ x = east - 500000.0+ y = if (ltz < 'N') then north - 10000000.0 else north+ longitudeOrigin = fromIntegral $ 6 * lnz - 183++ m = y / utm_f0+ mu = m / (a * (1 - eSquared / 4 - 3 / 64 * eSquared ** 2 - 5 / 256 * eSquared ** 3))++ phi1Rad = mu + (1.5 * e1 - 27 / 32 * e1 ** 3) * sin(2 * mu)+ + (21 / 16 * e1 ** 2 - 55 / 32 * e1 ** 4)+ * sin(4 * mu) + (151 / 96 * e1 ** 3) * sin(6 * mu)++ n = a / sqrt(1 - eSquared * sinSquared phi1Rad)+ t = tanSquared phi1Rad+ c = ePrimeSquared * cosSquared phi1Rad+ r = a * (1 - eSquared) / (1 - eSquared * sinSquared phi1Rad) ** 1.5+ d = x / (n * utm_f0)++ latitude = (phi1Rad - (n * tan phi1Rad / r)+ * (d * d / 2 - (5 + 3 * t + 10 * c - 4 * c ** 2 - 9 * ePrimeSquared)+ * d ** 4 / 24 + (61 + 90 * t + 298 * c + 45 * t ** 2+ - 252 * ePrimeSquared - 3 * c ** 2) * d ** 6 / 720)) * 180 / pi++ longitude = longitudeOrigin+ + ((d - (1 + 2 * t + c) / 6 * d ** 3+ + (5 - 2 * c + 28 * t - 3 * c ** 2 + 8 * ePrimeSquared + 24 * t ** 2)+ * d ** 5 / 120) / cos phi1Rad) * 180 / pi++ mkLatLng latitude longitude 0 wgs84Datum+++{-|+ Convert latitude and longitude to a UTM reference.+ If an attempt is made to convert a LatLng that falls outside the area covered by the UTM grid.+ The UTM grid is only defined for latitudes south of 84°N and north of 80°S.+-}+toUTMRef :: LatLng -> Except String UTMRef+toUTMRef (LatLng latitude longitude _ datum) =+ do+ lt <- withExcept (const $ "Latitude (" ++ show latitude ++ ") falls outside the UTM grid.") (validateLatitude latitude)+ let+ utm_f0 = 0.9996 :: Double+ a = semiMajorAxis wgs84Ellipsoid+ eSquared = eccentricitySquared wgs84Ellipsoid++ pifr = pi / 180.0 :: Double+ latitudeRad = lt * pifr+ ln = if longitude == 180.0 then (-180.0) else longitude+ longitudeRad = ln * pifr+ longitudeZone = evalLongitudeZone lt ln++ longitudeOrigin = fromIntegral $ 6 * longitudeZone - 183 :: Double+ longitudeOriginRad = longitudeOrigin * pifr++ utmZone = getUTMLatitudeZoneLetter lt++ ePrimeSquared = eSquared / (1 - eSquared)++ n = a / sqrt (1 - eSquared * sin latitudeRad * sin latitudeRad)+ t = tan latitudeRad ** 2+ c = ePrimeSquared * cos latitudeRad ** 2+ aa = cos latitudeRad * (longitudeRad - longitudeOriginRad)++ eFoured = eSquared * eSquared+ eSixed = eFoured * eSquared++ m = a * ((1 - eSquared / 4 - 3 * eFoured / 64 - 5 * eSixed / 256) * latitudeRad+ - (3 * eSquared / 8 + 3 * eFoured / 32 + 45 * eSixed / 1024) * sin (2 * latitudeRad)+ + (15 * eFoured / 256 + 45 * eSixed / 1024) * sin (4 * latitudeRad)+ - (35 * eSixed / 3072) * sin (6 * latitudeRad))++ utmEasting = (utm_f0 * n * (aa + (1 - t + c) * aa ** 3 / 6+ + (5 - 18 * t + t * t + 72 * c - 58 * ePrimeSquared) * aa ** 5.0 / 120) + 500000.0)++ -- Adjust for the southern hemisphere+ utmNorthingAdj = if (lt < 0) then 10000000.0 else 0.0+ utmNorthing = (utm_f0 * (m + n * tan latitudeRad * (aa * aa / 2 + (5 - t + 9 * c + 4 * c * c) * aa ** 4.0 / 24+ + (61 - 58 * t + t * t + 600 * c - 330 * ePrimeSquared) * aa ** 6.0 / 720))) + utmNorthingAdj++ mkUTMRef utmEasting utmNorthing utmZone longitudeZone+ where+ validateLatitude :: Double -> Except String Double+ validateLatitude lat+ | lat < -80.0 || lat > 84.0 = throwError "Invalid latitude"+ | lat == 180.0 = pure (-180.0)+ | otherwise = pure lat++ evalLongitudeZone :: Double -> Double -> Int+ evalLongitudeZone lat lng+ | lat >= 56.0 && lat < 64.0 && lng >= 3.0 && lng < 12.0 = 32 -- Special zone for Norway+ | lat >= 72.0 && lat < 84.0 && lng >= 0.0 && lng < 9.0 = 31 -- Special zones for Svalbard+ | lat >= 72.0 && lat < 84.0 && lng >= 9.0 && lng < 21.0 = 33+ | lat >= 72.0 && lat < 84.0 && lng >= 21.0 && lng < 33.0 = 35+ | lat >= 72.0 && lat < 84.0 && lng >= 33.0 && lng < 42.0 = 37+ | otherwise = fromIntegral $ floor (lng / 6.0 + 30.0) + 1 :: Int+++-- | Work out the UTM latitude zone from the latitude.+getUTMLatitudeZoneLetter :: Double -> Char+getUTMLatitudeZoneLetter latitude+ | 84 >= latitude && latitude >= 72 = 'X'+ | 72 > latitude && latitude >= 64 = 'W'+ | 64 > latitude && latitude >= 56 = 'V'+ | 56 > latitude && latitude >= 48 = 'U'+ | 48 > latitude && latitude >= 40 = 'T'+ | 40 > latitude && latitude >= 32 = 'S'+ | 32 > latitude && latitude >= 24 = 'R'+ | 24 > latitude && latitude >= 16 = 'Q'+ | 16 > latitude && latitude >= 8 = 'P'+ | 8 > latitude && latitude >= 0 = 'N'+ | 0 > latitude && latitude >= -8 = 'M'+ | -8 > latitude && latitude >= -16 = 'L'+ | -16 > latitude && latitude >= -24 = 'K'+ | -24 > latitude && latitude >= -32 = 'J'+ | -32 > latitude && latitude >= -40 = 'H'+ | -40 > latitude && latitude >= -48 = 'G'+ | -48 > latitude && latitude >= -56 = 'F'+ | -56 > latitude && latitude >= -64 = 'E'+ | -64 > latitude && latitude >= -72 = 'D'+ | -72 > latitude && latitude >= -80 = 'C'+ | otherwise = 'Z'
+ test/Spec.hs view
@@ -0,0 +1,19 @@+import Control.Monad+import Test.HUnit+import DMS.Tests (dmsTests)+import ECEFRef.Tests (ecefrefTests)+import IrishRef.Tests (irishrefTests)+import LatLng.Tests (latlngTests)+import MGRSRef.Tests (mgrsrefTests)+import OSRef.Tests (osrefTests)+import UTMRef.Tests (utmrefTests)++main :: IO ()+main = forM_ [ dmsTests+ , ecefrefTests+ , irishrefTests+ , latlngTests+ , mgrsrefTests+ , osrefTests+ , utmrefTests+ ] runTestTT