diff --git a/geo-uk.cabal b/geo-uk.cabal --- a/geo-uk.cabal +++ b/geo-uk.cabal @@ -2,7 +2,7 @@ -- documentation, see http://haskell.org/cabal/users-guide/ name: geo-uk -version: 0.1.0.0 +version: 0.1.0.1 synopsis: High precision conversion between GPS and UK Grid description: Conversion between GPS and UKGrid Coordinates one approximate and one using table. homepage: https://github.com/tolysz/geo-uk @@ -24,20 +24,18 @@ exposed-modules: Data.Geo.OSGM02Tab Data.Geo.Type Data.Geo.Point - Data.Geo.Conv -- other-modules: -- other-extensions: if flag(embeded) cpp-options: -DEMBEDED_DATA build-depends: bzlib + , file-embed build-depends: base >=4.7 && < 5 - , file-embed , array , binary , bytestring , th-lift , template-haskell - , geodetics hs-source-dirs: src default-language: Haskell2010 ghc-options: -O2 diff --git a/src/Data/Geo/Conv.hs b/src/Data/Geo/Conv.hs deleted file mode 100644 --- a/src/Data/Geo/Conv.hs +++ /dev/null @@ -1,358 +0,0 @@ - -module Data.Geo.Conv where -{- - - ( UKOS(..) - , WGS84 (..) - , convWGS84toUKOS - , convUKOStoWGS84 - , gridToGps - , gpsToGrid - , toUKOST - , toUKOSTI - ) -where - -import Data.Char (toUpper, toLower) -import Prelude --- import System.Cmd --- import System.IO.Unsafe --- import System.IO.Capture -import Text.Printf - -import Data.Geo.OSGM02Tab -import Data.Geo.Type -{- -capt p x = do - xa <- (capture (system (p++" "++x)) "") - return $ fst xa --} - --- WGS84 -- EPSG:4326 --- UKOS --- OSGB36 --- ETRS89 - -data UKOS = UKOST String - -- ^ Coordinates as string - | UKOSTI String - -- ^ My altilery endoding for String - | UKOS Int Int - -- ^ as Integers - | UKOSD Double Double Double - -- ^ With height - | UKOSI Int Int Int - -data WGS84 = WGS84 -- GPS - Double - Double - Double - deriving Show - -data ETRS89 = ETRS89 - Double - Double - Double - deriving Show ---type ETRS89 = WGS84EN - -gpsToGrid la lo = ( round e , round n) - where - (UKOSD e n _ ) = convWGS84toUKOS ( WGS84 la lo 0 ) - -gridToGps e n = ( la, lo) - where - ( WGS84 la lo _ ) = convUKOStoWGS84 (UKOS e n ) - - -instance Show UKOS where - show (UKOS a b) = (show a)++" "++(show b) - show (UKOST a) = a - show (UKOSTI a) = a - show (UKOSD e n h) = (show e)++" "++(show n) ++" "++ (show h) - show (UKOSI e n h) = (td e)++" "++(td n) ++" "++ (td h) - where - td x = printf "%.1f"(((fromIntegral x) * 0.1) :: Double) - -toRad x = x * pi / 180 -toDeg x = x * 180 / pi - - -toPTPea (la,lo) = (ea,no) - where - (ETRS89 ea no _) = convWGS84toETRS89 $ WGS84 la lo 0 - -convWGS84toUKOS = convETRS89toUKOS . convWGS84toETRS89 - -r1 = round . (10 *) - - -diffETRS89toUKOS :: Double -> Double -> (Double, Double, Double) -diffETRS89toUKOS we wn - | we < 0 - || we > 700000 - || wn < 0 - || wn > 1250000 = (we, wn, 0) - | otherwise = (se, sn, sg) - where - ei = round (we / 1000) - ni = round (wn / 1000) - (PTPDD se0 sn0 sg0 ) = gb ! (ei , ni ) - (PTPDD se1 sn1 sg1 ) = gb ! (ei + 1, ni ) - (PTPDD se2 sn2 sg2 ) = gb ! (ei + 1, ni + 1) - (PTPDD se3 sn3 sg3 ) = gb ! (ei , ni + 1) - dx = we - 1000 * (fromIntegral ei) - dy = wn - 1000 * (fromIntegral ni) - t = dx / 1000 - u = dy / 1000 - se = (1-t)*(1-u)*se0 + t*(1-u)*se1 + t*u*se2 + (1-t)* u * se3 - sn = (1-t)*(1-u)*sn0 + t*(1-u)*sn1 + t*u*sn2 + (1-t)* u * sn3 - sg = (1-t)*(1-u)*sg0 + t*(1-u)*sg1 + t*u*sg2 + (1-t)* u * sg3 - -convUKOStoETRS89 :: UKOS -> ETRS89 -convUKOStoETRS89 a@(UKOST _) = convUKOStoETRS89 $ normalizeUKOS a -convUKOStoETRS89 a@(UKOSTI _) = convUKOStoETRS89 $ normalizeUKOS $ toUKOST a -convUKOStoETRS89 (UKOS a b) = convUKOStoETRS89 $ UKOSI a b 0 -convUKOStoETRS89 (UKOSI a b c) = convUKOStoETRS89 $ UKOSD (f a) (f b) (f c) - where - f = fromIntegral - -convUKOStoETRS89 (UKOSD e n a) = ETRS89 (e - de) (n - dn) (a + da) - where - (de, dn, da) = loop 0 0 - loop de1 dn1 = let - h@(den, dnn, dan) = diffETRS89toUKOS (e - de1) (n - dn1) - in - case (abs(de1 - den) + abs(dn1 - dnn) < 0.000000001 ) of - True -> h - False -> loop den dnn - - -convETRS89toUKOS :: ETRS89 -> UKOS -convETRS89toUKOS (ETRS89 we wn wa) = UKOSD e n h -- UKOSI (r1 e,r1 n,r1 h) - where - (se, sn, sg) = diffETRS89toUKOS we wn - e = we + se - n = wn + sn - h = wa - sg - - - -convETRS89toWGS84 :: ETRS89 -> WGS84 -convETRS89toWGS84 (ETRS89 ea no alt) = (WGS84 phi' lam' alt) - where - a = 6378137.000 - b = 6356752.31414 --- a = 6377563.396 --- b = 6356256.910 - f0 = 0.9996012717 - phi0 = toRad 49 - lam0 = toRad (-2) - e0 = 400000 - n0 = -100000 - - de = ea - e0 - de2 = de * de - - a2 = a**2 - b2 = b**2 - e2 = (a2 - b2) / a2 -- B1 - - phi = fixP ((no - n0) / a / f0 + phi0) - - fixP p = let d = (no - n0 -(calcM a b p)) - r = d / a / f0 + p - in - if abs(d)<0.00001 then - r - else - fixP r - - n = (a - b)/(a+b) -- B2 - n2 = n*n - n3 = n2*n - s = sin phi - s2 = s*s - - omes = 1 / sqrt (1 - e2 *s2) - - ν = a * f0 * omes --B3 - ρ = a * f0 * (1 - e2)*omes**3 --B4 - ν2 = 1 / (ν * ν) - η = ν / ρ - 1 -- B5 nu² --- phi' = phi --- lam' = 0 - t = tan phi - t2 = t * t - vii = t / (2 * ρ * ν) - - viii' = vii * ν2 / 12 - viii = viii' * (5 + 3 * t2 * (1 - 3 * η) + η ) - - ix' = viii' * ν2 / 30 - ix = ix' *(61 + 45 * t2 *(2 + t2)) --- phi' = ix - phi' = toDeg $ phi - de2 * ( vii - de2 * (viii - de2 * ix)) - - x = 1 / ( ν * (cos phi)) - xi' = x * ν2 / 6 - xi = xi' * (1 + η + 2 * t2) - xii' = xi' * ν2 / 20 - xii = xii' * (5 + 4*t2*(7 + 6 * t2)) - xiia' = xii' * ν2 / 42 - xiia = xiia' * ( 61 + t2 *( 662 + t2*(1320 + 720 * t2))) --- phi' = phi - lam' = toDeg $ lam0 + de * (x - de2 * ( xi - de2 * (xii - de2 * xiia))) - -convUKOStoWGS84 :: UKOS -> WGS84 -convUKOStoWGS84 = convETRS89toWGS84 . convUKOStoETRS89 - - -convWGS84toETRS89 :: WGS84 -> ETRS89 -convWGS84toETRS89 (WGS84 phi' lam' alt) = ETRS89 ea no alt - where - phi = toRad phi' - lam = toRad lam' - a = 6378137.000 - b = 6356752.3141 - - s = sin phi - s2 = s*s - - c = cos phi - c2 = c*c - c3 = c2 * c - c5 = c3 * c2 - - ta2 = s2/c2 - ta4 = ta2 * ta2 - - a2 = a**2 - b2 = b**2 - e2 = (a2 - b2) / a2 -- B1 - - f0 = 0.9996012717 - phi0 = toRad 49 - lam0 = toRad (-2) - e0 = 400000 - n0 = -100000 - - omes = 1 / sqrt (1 - e2 *s2) - - ν = a * f0 * omes --B3 - ρ = a * f0 * (1 - e2)*omes**3 --B4 - - η = ν / ρ - 1 -- B5 nu² - - m = calcM a b phi - - i = m + n0 - ii = ν/2 * s * c - iii = ν/24 * s * c3 * (5 - ta2 + 9*η) - iiia = ν/720 * s * c5 * (61 -58 * ta2 + ta4) - iv = ν * c - v = ν/6 * c3 * (1 + η - ta2) - vi = ν/120 * c5 *(5 - 18 * ta2 + ta4 + η * (14 - 58 * ta2)) - - mlam = lam - lam0 - mlam2 = mlam * mlam - - ea = e0 + mlam * (iv + mlam2 * (v + mlam2 * vi)) - no = i + mlam2 * (ii + mlam2 * (iii + iiia * mlam2)) - - -calcM a b phi = m - where --- a = 6378137.000 --- b = 6356752.3141 --- a = 6377563.396 --- b = 6356256.910 - f0 = 0.9996012717 - phi0 = toRad 49 ---- lam0 = toRad (-2) - - n = (a - b)/(a+b) -- B2 - n2 = n*n - n3 = n2*n - - - pphi = phi + phi0 - mphi = phi - phi0 - m1 = (1 + n + 5/4 * (n2+n3)) * mphi - m2 = (3*(n+n2) + 2.625*n3) * (sin mphi) * (cos pphi) - m3 = 1.875*(n2+n3) * (sin $ 2 * mphi) * (cos $ 2 * pphi) - m4 = 35/24*n3 * (sin $ 3 * mphi) * (cos $ 3 * pphi) - m = f0 * b * (m1 - m2 + m3 - m4) - - - --- convUKOStoWGS84' :: UKOS -> WGS84 --- convUKOStoWGS84' u = WGS84 la lo a --- where --- (lo:la:a:as) = map (\x->(read x )::Double) $ words (unsafePerformIO $ capt "gridToGPS" (show u)) - - -fZ a = reverse $ (take 5) $ (reverse (show $ a `mod` 100000)) ++ "00000" - - -toUKOSTI :: UKOS -> UKOS -toUKOSTI (UKOST (a1:b1:a2:a3:a4:a5:a6:b2:b3:b4:b5:b6:[])) = UKOSTI (toLower a1:toLower b1:a2:b2:a3:b3:a4:b4:a5:b5:a6:b6:[]) -toUKOSTI a = toUKOSTI $ toUKOST a - - -toUKOST :: UKOS -> UKOS -toUKOST a@(UKOST _) = a -toUKOST (UKOSTI (a1:b1:a2:b2:a3:b3:a4:b4:a5:b5:a6:b6:[])) = (UKOST (toUpper a1:toUpper b1:a2:a3:a4:a5:a6:b2:b3:b4:b5:b6:[])) -toUKOST (UKOSTI a) = toUKOST $ UKOSTI $ fx $ ll (length a) - where - ll 0 = "NT55" - ll 1 = "T55" - ll a = if (even a) then "55" else "5" - fx b = take 12 (a ++ b ++ (repeat '0')) - -toUKOST (UKOS a b) = UKOST $ (toLetters a b) ++ (fZ a) ++ (fZ b) -toUKOST (UKOSI a b c) = toUKOST $ UKOS a b -toUKOST (UKOSD a b c) = toUKOST $ UKOS (round a) (round b) - ---- TG 5140 1317 can also be expressed as 65140,31317 -toLetters :: Int -> Int -> String -toLetters e' n' = [x,y] -- ++ (show $ e'-e * 100000) ++ (show $ n'-n * 100000) - where - e = e' `div` 100000 - n = n' `div` 100000 - l1 = 19 - n - (19 - n) `mod`5 + (e+10) `div` 5 - l2 = ((19 -n)*5) `mod` 25 + e `mod` 5 - l1f | l1 > 7 = 1 - | otherwise = 0 - l2f | l2 > 7 = 1 - | otherwise = 0 - x = toEnum (65 + l1 + l1f) - y = toEnum (65 + l2 + l2f) - - -normalizeUKOS :: UKOS -> UKOS - -normalizeUKOS (UKOST tt) = UKOS e' n' - where - (a:b:t) = filter (/=' ') tt - ia = (fromEnum $ toUpper a) - 65 - ia' = if (ia > 7) then (ia - 1) else ia - ib = (fromEnum $ toUpper b) - 65 - ib' = if (ib > 7) then (ib - 1) else ib - e = ((ia' - 2) `mod` 5) * 5 + (ib' `mod` 5) - n = (19 - (ia' `div` 5) * 5) - (ib' `div` 5) - l = (length t) `div` 2 - c = case l of - 1->5000 - 2->500 - 3->50 - 4->5 - _->0 - (e1,n1) = case l of - 0 -> (0,0) - _ -> ((read $ take l t)::Int,(read $ drop l t):: Int) - e' = (10 ^ l * e + e1) * 10 ^ (5 -l) + c - n' = (10 ^ l * n + n1) * 10 ^ (5 -l) + c -normalizeUKOS a@(UKOSTI _ ) = normalizeUKOS $ toUKOST a -normalizeUKOS a = a --} diff --git a/src/Data/Geo/OSGM02Tab.hs b/src/Data/Geo/OSGM02Tab.hs --- a/src/Data/Geo/OSGM02Tab.hs +++ b/src/Data/Geo/OSGM02Tab.hs @@ -13,19 +13,9 @@ #if EMBEDED_DATA import qualified Codec.Compression.BZip as BZip import qualified Data.ByteString.Lazy.Char8 as L8 --- import qualified Data.ByteString.Lazy.Internal as L8 --- import qualified System.Environment as Env import Data.Binary (decode) import Data.FileEmbed --- import Data.Array --- import System.IO.Unsafe --- import Language.Haskell.TH.InlineIO --- import Language.Haskell.TH.Lift - import qualified Data.ByteString - --- myFile :: Data.ByteString.ByteString --- myFile = $(embedFile "dirName/fileName" gb :: Maybe ArrDD gb = Just (decodePTPDDFile . L8.fromStrict $ $(embedFile "data/GB.dat")) diff --git a/src/Data/Geo/Type.hs b/src/Data/Geo/Type.hs --- a/src/Data/Geo/Type.hs +++ b/src/Data/Geo/Type.hs @@ -18,9 +18,6 @@ instance Lift PTPDD where lift (PTPDD a b c) = [| PTPDD $(lift a) $(lift b) $(lift c) |] --- instance Lift Double where --- lift x = [| $(litE $ rationalL $ toRational x) :: Double |] - instance Binary PTPDD where put (PTPDD de dn da ) = put de >>