hcoord-1.0.0.0: src/UTMRef.hs
{-# 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'