{-# 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)