geodetics-2.0.0: test/Main.hs
{-# OPTIONS_GHC -fno-warn-orphans #-}
{-# OPTIONS_GHC -Wno-unrecognised-pragmas #-}
{-# HLINT ignore "Redundant bracket" #-}
{-# LANGUAGE ViewPatterns #-}
module Main where
import Control.Monad
import Data.Char
import Data.Either
import Data.Maybe
import Test.Hspec
import Test.Hspec.QuickCheck
import Test.HUnit (assertFailure)
import Test.QuickCheck
import Test.QuickCheck.Checkers (EqProp, eq, (=-=), unbatch)
import Test.QuickCheck.Classes (monoid)
import ArbitraryInstances
import Geodetics.Altitude
import Geodetics.Ellipsoids
import Geodetics.Geodetic
import Geodetics.Grid
import Geodetics.MGRS
import Geodetics.Path
import Geodetics.Stereographic
import Geodetics.TransverseMercator
import Geodetics.UK
import Geodetics.UTM
import LatLongParser (parserTests)
import Geodetics.PolarStereographic as PS
main :: IO ()
main = hspec $ do
describe "Geodetic" $ do
prop "WGS84 and back" prop_WGS84_and_back
prop "Zero ground distance" prop_zero_ground
describe "UK Points" $ mapM_ pointTest ukPoints
describe "World line" $ mapM_ worldLineTests worldLines
parserTests
describe "Grid" $ do
prop "Grid Offset 1" prop_offset1
prop "Grid Offset 2" prop_offset2
prop "Grid Offset 3" prop_offset3
prop "Grid 1" prop_grid1
describe "TransverseMercator" $ prop "fromGrid . toGrid == id" prop_tmGridInverse
describe "UK" $ do
prop "UK Grid 1" prop_ukGrid1
describe "UK Grid 2" $ mapM_ ukGridTest2 ukSampleGrid
describe "UK Grid 3" $ mapM_ ukGridTest3 ukSampleGrid
describe "UK Grid 4" $ mapM_ ukGridTest4 ukSampleGrid
describe "UK Grid 5" $ mapM_ ukGridTest5 ukSampleGrid
describe "UTM" $ do
prop "UTM Grid 1" prop_utmGridTest1
describe "UTM Grid 2" $ mapM_ utmGridTest2 utmSampleGrid
describe "UTM Grid 3" $ mapM_ utmGridTest3 utmSampleGrid
describe "UTM Grid 4" $ mapM_ utmGridTest4 utmSampleGrid
describe "UTM Grid 5" $ mapM_ utmGridTest5 utmSampleGrid
describe "MGRS" $ do
prop "MGRS Grid 1" prop_mgrs_gridTest1
prop "MGRS Grid 2" prop_mgrs_gridTest2
describe "MGRS Grid 3" $ mapM_ mgrsGridTest3 utmSampleGrid
describe "MGRS Grid 4" $ mapM_ mgrsGridTest4 utmSampleGrid
describe "UPS" $ do
describe "UPS Grid 4" $ mapM_ upsGridTest4 upsSampleGrid
describe "UPS Grid 5" $ mapM_ upsGridTest5 upsSampleGrid
describe "UPS Grid 6" $ mapM_ upsGridTest6 upsSampleGrid
describe "UPS Grid 7" $ mapM_ upsGridTest7 upsSampleGrid
describe "UPS Grid 8" $ mapM_ upsGridTest8 upsSampleGrid
describe "Stereographic" $ do
it "toGrid north" stereographicToGridN
it "fromGrid north" stereographicFromGridN
it "toGrid south" stereographicToGridS
it "fromGrid south" stereographicFromGridS
prop "Stereographic round trip" prop_stereographic
describe "Paths" $ do
prop "Ray Path 1" prop_rayPath1
prop "Ray Continuity" prop_rayContinuity
prop "Ray Bisection" prop_rayBisect
prop "Rhumb Continuity" prop_rhumbContinuity
prop "Rhumb Intersection" prop_rhumbIntersect
describe "GridOffset monoid" $ mapM_ (uncurry prop) $ unbatch $ monoid (mempty :: GridOffset)
describe "Helmert monoid" $ mapM_ (uncurry prop) $ unbatch $ monoid (mempty :: Helmert)
instance EqProp GridOffset where
(GridOffset a b c) =-= (GridOffset a' b' c') =
eq True $ a ≈ a' && b ≈ b' && c ≈ c'
where x ≈ y = abs (x - y) < 0.00001
instance EqProp Helmert where
(Helmert cX' cY' cZ' s rX' rY' rZ') =-= (Helmert cX'' cY'' cZ'' s' rX'' rY'' rZ'') =
eq True $ and [cX' ≈ cX'', cY' ≈ cY'', cZ' ≈ cZ'',
s ≈- s',
rX' ≈- rX'', rY' ≈- rY'', rZ' ≈- rZ'']
where x ≈ y = abs (x - y) < 0.00001
x ≈- y = abs (x - y) < 1 / (5 * 2) ^ _5
-- | The positions are within 30 cm.
samePlace :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Expectation
samePlace p1 p2 = expectTrue msg $ geometricalDistance p1 p2 < 0.3
where
msg = "location " <> show p2 <> " is > 30cm from expected " <> show p1
samePlace' :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool
samePlace' p1 p2 = geometricalDistance p1 p2 < 0.3
-- | The positions are within 10 m.
closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Expectation
closeEnough p1 p2 = expectTrue msg $ geometricalDistance p1 p2 < 10
where
msg = "location " <> show p2 <> " is > 10m from expected " <> show p1
closeEnough' :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool
closeEnough' p1 p2 = geometricalDistance p1 p2 < 10
-- | The angles are within 0.01 arcsec
sameAngle :: Double -> Double -> Expectation
sameAngle v1 v2 = expectTrue msg $ abs (properAngle (v1 - v2)) < 0.01 * arcsecond
where
msg = "expected angle " <> show (v1 / degree) <> ", got " <> show (v2 / degree)
sameAngle' :: Double -> Double -> Bool
sameAngle' v1 v2 = abs (properAngle (v1 - v2)) < 0.01 * arcsecond
-- | The grid positions are within 1mm
sameGrid :: (Show r) => GridPoint r -> GridPoint r -> Expectation
sameGrid p1 p2 = expectTrue msg $ check eastings && check northings && check altitude
where
msg = "expected " <> show p1 <> ", got " <> show p2
check f = f p1 - f p2 < 1e-3
-- | Grid offsets are within 1mm.
sameOffset :: GridOffset -> GridOffset -> Expectation
sameOffset go1 go2 = expectTrue msg $ check deltaNorth && check deltaEast && check deltaAltitude
where
msg = "expected " <> show go1 <> ", got " <> show go2
check f = f go1 - f go2 < 1e-3
-- | The grid X and Y are both within 1 meter
closeGrid :: (Show r) => GridPoint r -> GridPoint r -> Expectation
closeGrid p1 p2 = expectTrue msg $ check eastings && check northings && check altitude
where
msg = "expected " <> show p1 <> ", got " <> show p2
check f = f p1 - f p2 < 1
-- | Degrees, minutes and seconds into radians.
dms :: Int -> Int -> Double -> Double
dms d m s = fromIntegral d * degree + fromIntegral m * arcminute + s * arcsecond
-- | Round-trip from local to WGS84 and back is identity (approximately)
prop_WGS84_and_back :: Geodetic LocalEllipsoid -> Expectation
prop_WGS84_and_back p = samePlace p $ toLocal (ellipsoid p) $ toWGS84 p
-- | Test that for all points p, the ground distance from p to p is zero.
prop_zero_ground :: Geodetic WGS84 -> Bool
prop_zero_ground p =
case groundDistance p p of
Nothing -> False
Just (d, _, _) -> abs d < 1e-3
-- | Sample pairs of points with bearings and distances.
-- The Oracle for these values is the @FORWARD@ program from
-- <http://www.ngs.noaa.gov/TOOLS/Inv_Fwd/Inv_Fwd.html>
worldLines :: [(String, Geodetic WGS84, Geodetic WGS84, {-Length-} Double, {-Angle-} Double, {-Angle-} Double)]
worldLines = [
("Ordinary", Geodetic (40 * degree) (30 * degree) 0 WGS84, Geodetic (30 * degree) (50 * degree) 0 WGS84,
2128852.999, 115.19596706 * degree, 126.79044315 * degree),
("Over Pole", Geodetic (60 * degree) (0 * degree) 0 WGS84, Geodetic (60 * degree) (180 * degree) 0 WGS84,
6695785.820, 0 * degree, 180 * degree),
("Equator to Pole", Geodetic (0 * degree) (0 * degree) 0 WGS84, Geodetic (90 * degree) (180 * degree) 0 WGS84,
10001965.729, 0 * degree, 180 * degree)]
worldLineTests :: (String, Geodetic WGS84, Geodetic WGS84, Double, Double, Double) -> SpecWith (Arg Expectation)
worldLineTests (str, g1, g2, d, a, b) = it str $ ok $ groundDistance g1 g2
where
ok Nothing = False
ok (Just (d1, a1, b1)) =
abs (d - d1) < 0.01
&& abs (a - a1) < 0.01 * arcsecond
&& abs (b - b1) < 0.01 * arcsecond
-- | Sample points for UK tests. The oracle for these values is the script at
-- <http://www.movable-type.co.uk/scripts/latlong-convert-coords.html>, which uses
-- the same Helmert transform as this library. Hence the results should match to within 30 cm.
ukPoints :: [(String, Geodetic WGS84, Geodetic OSGB36)]
ukPoints = [
("Greenwich", Geodetic (dms 51 28 40.86) (dms 0 0 (-5.83)) 0 WGS84,
Geodetic (dms 51 28 39.00) (dms 0 0 0) 0 OSGB36),
("Edinburgh Castle", Geodetic (dms 55 56 56.30) (dms (-3) (-12) (-2.73)) 0 WGS84,
Geodetic (dms 55 56 56.51) (dms (-3) (-11) (-57.61)) 0 OSGB36),
("Lands End", Geodetic (dms 50 03 56.68) (dms (-5) (-42) (-51.20)) 0 WGS84,
Geodetic (dms 50 03 54.51) (dms (-5) (-42) (-47.87)) 0 OSGB36),
("Gt. Yarmouth Pier",Geodetic (dms 52 36 29.33) (dms 1 44 27.79) 0 WGS84,
Geodetic (dms 52 36 27.84) (dms 1 44 34.52) 0 OSGB36),
("Stanhope", Geodetic (dms 54 44 49.08) (dms (-2) 0 (-19.89)) 0 WGS84,
Geodetic (dms 54 44 48.71) (dms (-2) 0 (-14.41)) 0 OSGB36) ]
-- Convert a named point into a test
pointTest :: (Ellipsoid e2) => (String, Geodetic WGS84, Geodetic e2) -> SpecWith (Arg Expectation)
pointTest (testName, wgs84, local) = it testName $ wgs84 `samePlace` toWGS84 local
-- The negation of the sum of a list of offsets is equal to the sum of the negated items.
prop_offset1 :: [GridOffset] -> Expectation
prop_offset1 offsets = sameOffset (offsetNegate $ mconcat offsets) (mconcat $ map offsetNegate offsets)
-- A polar offset multiplied by a scalar is equal to an offset in the same direction with the length multiplied.
prop_offset2 :: Distance -> Bearing -> Scalar -> Expectation
prop_offset2 (Distance d) (Bearing h) (Scalar s) = sameOffset go1 go2
where
go1 = offsetScale s $ polarOffset d h
go2 = polarOffset (d * s) h
-- | A polar offset has the offset distance and bearing of its arguments.
prop_offset3 :: GridOffset -> Expectation
prop_offset3 delta = sameOffset delta0
(polarOffset (offsetDistance delta0) (offsetBearing delta))
where delta0 = delta {deltaAltitude = 0}
-- | Given a grid point and an offset, applying the offset to the point gives a new point which
-- is offset from the first point by the argument offset.
prop_grid1 :: GridPoint (GridTM LocalEllipsoid) -> GridOffset -> Expectation
prop_grid1 p d = sameOffset d $ p `gridOffset` applyOffset d p
-- | Check that using toGrid/fromGrid for TransverseMercator projection are inverses
-- | for negative latitudes near the coordinates 0,0
prop_tmGridInverse :: Expectation
prop_tmGridInverse =
let origin = Geodetic
{ latitude = 0 * degree
, longitude = 0 * degree
, geoAlt = 0
, ellipsoid = WGS84
}
g = mkGridTM origin mempty 1
testPoint = origin { latitude = (-1) * arcminute }
tp1 = toGrid g testPoint
tp2 = fromGrid tp1
in tp2 `closeEnough` testPoint
-- | Converting a UK grid reference to a GridPoint and back is a null operation.
prop_ukGrid1 :: UkGridRef -> Expectation
prop_ukGrid1 (UkGridRef str) =
str `shouldBe`
fromJust (toUkGridReference ((length str - 2) `div` 2) $ fst $ fromJust $ fromUkGridReference str)
-- | UK Grid Reference points. The oracle for these points was the
-- UK Grid Reference Finder (gridreferencefinder.com), retrieved on 26 Jan 2013.
ukSampleGrid :: [(String, GridPoint UkNationalGrid, Geodetic WGS84, String)]
ukSampleGrid = map convert [
-- Grid Reference, X, Y, Latitude, Longitude, Description
("SW3425625070", 134256, 025070, 50.066230, -5.7148278, "Lands End"),
("TR3302139945", 633021, 139945, 51.111396, 1.3277159, "Dover Harbour"),
("TQ3001980417", 530019, 180417, 51.507736, -0.12793230, "Nelsons Column"),
("TA2542370644", 525423, 470644, 54.116376, -0.082668990, "Flamborough Lighthouse"),
("NK1354745166", 413547, 845166, 57.496512, -1.7756310, "Peterhead harbour"),
("ND3804872787", 338048, 972787, 58.638518, -3.0688688, "John O Groats"),
("SC3915875189", 239158, 475189, 54.147275, -4.4641148, "Douglas Harbour"),
("ST1922474591", 319224, 174591, 51.464505, -3.1641741, "Torchwood HQ"),
("SK3520736502", 435207, 336502, 52.924784, -1.4777486, "Derby Cathedral"),
("TG5141013177", 651410, 313177, 52.657979 , 1.7160519, "Caister Water Tower"),
("TG2623802646", 626238, 302646, 52.574548 , 1.3373749, "Framingham")]
-- Caister and Framingham are taken from Ordnance Survey worked examples.
where
convert (grid, x, y, lat, long, desc) =
(grid, GridPoint x y 0 UkNationalGrid,
Geodetic (lat * degree) (long * degree) 0 WGS84, desc)
type UkGridPointTest = (String, GridPoint UkNationalGrid, Geodetic WGS84, String) -> SpecWith (Arg Expectation)
-- | Check that grid reference to grid point works for sample points.
ukGridTest2 :: UkGridPointTest
ukGridTest2 (gridRef, gp, _, testName) =
it testName $ fst (fromJust $ fromUkGridReference gridRef) `shouldBe` gp
-- | Check that grid point to grid reference works for sample points.
ukGridTest3 :: UkGridPointTest
ukGridTest3 (gridRef, gp, _, testName) =
it testName $ toUkGridReference 5 gp `shouldBe` Just gridRef
-- | Check that grid point to WGS84 works close enough for sample points.
ukGridTest4 :: UkGridPointTest
ukGridTest4 (_, gp, geo, testName) =
it testName $ geo `closeEnough` toWGS84 (fromGrid gp)
-- | Check that WGS84 to grid point works close enough for sample points.
ukGridTest5 :: UkGridPointTest
ukGridTest5 (_, gp, geo, testName) =
it testName $ gp `closeGrid` toGrid UkNationalGrid (toLocal OSGB36 geo)
-- | Worked example for UK Geodetic to GridPoint, taken from "A Guide to Coordinate Systems in Great Britain" [1]
ukTest :: Geodetic OSGB36
ukTest = Geodetic (dms 52 39 27.2531) (dms 1 43 4.5177) 0 OSGB36
{-
v = 6.3885023333E+06
rho = 6.3727564399E+06
eta2 = 2.4708136169E-03
m = 4.0668829596E+05
I = 3.0668829596E+05
II = 1.5404079092E+06
III = 1.5606875424E+05
IIIa = -2.0671123011E+04
IV = 3.8751205749E+06
V = -1.7000078208E+05
VI = -1.0134470432E+05
E = 651409.903 m
N = 313177.270 m
-}
-- | Check that a UTM grid reference round-trips to a gridpoint and back.
prop_utmGridTest1 :: UtmGridRef -> Expectation
prop_utmGridTest1 (UtmGridRef str) =
str `shouldBe`
toUtmGridReference Nothing True 0 (fromRight (error str) $ fromUtmGridReference str)
-- | Sample points for UTM, in both UTM and MGRS formats. The oracle for these points was the Montana
-- State University converter. http://rcn.montana.edu/resources/Converter.aspx.
-- Retrieved on 1st Feb 2025.
utmSampleGrid :: [(String, String, GridPoint UtmZone, Geodetic WGS84, String)]
utmSampleGrid = map convert [
-- UTM Reference, MGRS Reference, X, Y, Latitude, Longitude, Description
("30N 699304 5710208", "30U XC 99304 10208", 699304, 5710208, 51.5078064, -0.1279388, "Nelson's Column"),
("35S 379101 8017747", "35K LA 79101 17747", 379101, 8017747, -17.9249501, 25.8585071, "Victoria Falls"),
("27N 454366 7111715", "27W VM 54366 11715", 454366, 7111715, 64.1289075, -21.9373288, "Reykjavik Airport"),
("32N 297697 6700532", "32V KN 97697 00532", 297697, 6700532, 60.3904298, 5.3284215, "Bergen"),
("23S 683473 7460697", "23K PQ 83473 60697", 683473, 7460697, -22.9518122, -43.2105383, "Christ the Redeemer"),
("21S 439699 4272868", "21F VC 39699 72868", 439699, 4272868, -51.6919073, -57.8724006, "Port Stanley"),
("59S 461474 4919822", "59G MK 61474 19822", 461474, 4919822, -45.8740880, 170.5035807, "Dunedin"),
("31N 166022 0" , "31N AA 66022 00000", 166022, 0, 0.0 , 0.0 , "The Origin"),
("31S 166022 9999999", "31M AV 66022 99999", 166022, 9999999, -0.00000903, 0.0 , "1 meter south")]
where
convert (utm, mgrs, x, y, lat, long, desc) =
(utm, mgrs, GridPoint x y 0 zone, geo, desc)
where
geo = Geodetic (lat * degree) (long * degree) 0 WGS84
znum = fromJust $ utmZoneNumber geo
hemi = if lat < 0 then UtmSouth else UtmNorth
zone = fromJust $ mkUtmZone hemi znum
type UtmGridPointTest = (String, String, GridPoint UtmZone, Geodetic WGS84, String) -> SpecWith (Arg Expectation)
-- | Check that UTM reference to grid point works for sample points.
utmGridTest2 :: UtmGridPointTest
utmGridTest2 (gridRef, _, gp, _, testName) =
it testName $ fromUtmGridReference gridRef `shouldBe` Right gp
-- | Check that Grid point to UTM reference works for sample points.
utmGridTest3 :: UtmGridPointTest
utmGridTest3 (gridRef, _, gp, _, testName) =
it testName $ toUtmGridReference Nothing False 0 gp `shouldBe` gridRef
-- | Check that grid point to WGS84 works close enough for sample points.
utmGridTest4 :: UtmGridPointTest
utmGridTest4 (_, _, gp, geo, testName) =
it testName $ geo `closeEnough` fromGrid gp
-- | Check that WGS84 to grid point works close enough for sample points.
utmGridTest5 :: UtmGridPointTest
utmGridTest5 (_, _, gp, geo, testName) =
it testName $ gp `closeGrid` toGrid (fromJust $ utmZone geo) geo
-- | Check that a UTM grid point round-trips to MGRS and back, with spaces.
prop_mgrs_gridTest1 :: UtmGridRef -> Expectation
prop_mgrs_gridTest1 (UtmGridRef str) =
case fromUtmGridReference str of
Left msg -> assertFailure $ "gridTest1 bogus UTM ref: " <> str <> ". Messages = " <> show msg
Right (utmToMgrsPoint -> gp) ->
fromMgrsGridReference <$> toMgrsGridReference True 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0))
-- | Check that a UTM grid point round-trips to MGRS and back, without spaces.
prop_mgrs_gridTest2 :: UtmGridRef -> Expectation
prop_mgrs_gridTest2 (UtmGridRef str) =
case fromUtmGridReference str of
Left msg -> assertFailure $ "gridTest2 bogus UTM ref: " <> str <> ". Messages = " <> show msg
Right (utmToMgrsPoint -> gp) ->
fromMgrsGridReference <$> toMgrsGridReference False 5 gp `shouldBe` Just (Right (gp, GridOffset 0.5 0.5 0))
-- | Check that MGRS reference to grid point works for sample points.
mgrsGridTest3 :: UtmGridPointTest
mgrsGridTest3 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) =
it testName $ fromMgrsGridReference mgrs `shouldBe` Right (gp, GridOffset 0.5 0.5 0)
-- | Check that grid point to MGRS reference works for sample points.
mgrsGridTest4 :: UtmGridPointTest
mgrsGridTest4 (_, mgrs, (utmToMgrsPoint -> gp), _, testName) = do
it (testName <> " with spaces") $ toMgrsGridReference True 5 gp `shouldBe` Just mgrs
it (testName <> " without spaces") $ toMgrsGridReference False 5 gp `shouldBe` Just (filter (not . isSpace) mgrs)
-- | DMA TM 8358.2 can be found at https://apps.dtic.mil/sti/tr/pdf/ADA266497.pdf
-- Oracle for other points is https://geographiclib.sourceforge.io/C++/doc/GeoConvert.1.html
upsSampleGrid :: [(String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String)]
upsSampleGrid = map convert [
-- MGRS UPS, X, Y, Latitude, Longitude
("ZAH0000000000", "2000000 2000000", 2000000.00, 2000000.00, 90, 0, "North pole"),
("BAN0000000000", "2000000 2000000", 2000000.00, 2000000.00, -90, 0, "South pole"),
("ZAE0000077930", "2000000 1777930", 2000000.00, 1777930.73, 88, 0, "88N 0W"),
("YXH7793000000", "1777930 2000000", 1777930.73, 2000000.00, 88, -90, "88N, 90W"),
("BAQ0000022069", "2000000 2222069", 2000000.00, 2222069.27, -88, 0, "88S 0W"),
("AXN7793000000", "1777930 2000000", 1777930.73, 2000000.00, -88, -90, "88S, 90W"),
("YTM3012526773", "1530125 2426773", 1530125.78, 2426773.60, 84.28723389, -132.2479892, "DMA TM 8358.2 sample 1"),
("BCK2297997474", "2222979 1797474", 2222979.47, 1797474.90, -87.28733333, 132.2478619, "DMA TM 8358.2 sample 3") ]
-- DMA TM 8358.2 sample 2 is outside the UPS grid regions and is therefore not a valid test.
where
convert (mgrs, ups, x, y, lat, long, desc) = (mgrs, ups, GridPoint { eastings = x, northings = y, altGP = 0, gridBasis = gb }, geo, desc)
where
geo = Geodetic { latitude = (lat * degree), longitude = (long * degree), geoAlt = 0, ellipsoid = WGS84 }
gb = mkGridPolarStereographic
(if lat < 0 then SouthPole else NorthPole)
WGS84
(GridOffset { deltaEast = -(2000 * kilometer), deltaNorth = -(2000 * kilometer), deltaAltitude = 0 })
0.994
type UpsGridPointTest = (String, String, GridPoint (PolarStereographic WGS84), Geodetic WGS84, String) -> SpecWith (Arg Expectation)
-- | Check that the UPS grid point to WGS84 works close enough for sample points.
upsGridTest4 :: UpsGridPointTest
upsGridTest4 (_, _, gp, geo, testName) =
it testName $ geo `closeEnough` fromGrid gp
-- | Check that WGS84 to UPS grid point works close enough for sample points.
upsGridTest5 :: UpsGridPointTest
upsGridTest5 (_, _, gp, geo, testName) =
it testName $ gp `closeGrid` toGrid (gridBasis gp) geo
-- | Check that UPS and MGRS parsers agree.
upsGridTest6 :: UpsGridPointTest
upsGridTest6 (mgrs, ups, gp, _, testName) =
it testName $
fst <$> fromMgrsGridReference mgrs `shouldBe`
upsToMgrsPoint <$> fromUpsGridReference (PS.trueOrigin $ gridBasis gp) ups
-- | Check MGRS grid generator
upsGridTest7 :: UpsGridPointTest
upsGridTest7 (mgrs, _, gp, _, testName) =
it testName $ toMgrsGridReference False 5 (upsToMgrsPoint gp) `shouldBe` Just mgrs
-- | Check UPS grid generator
upsGridTest8 :: UpsGridPointTest
upsGridTest8 (_, ups, gp, _, testName) =
it testName $ toUpsGridReference Nothing False 0 gp `shouldBe` ups
-- | Standard stereographic grid for point tests in the Northern Hemisphere.
stereoGridN :: GridStereo LocalEllipsoid
stereoGridN = mkGridStereo tangent origin 0.9999079
where
ellipse = LocalEllipsoid "Bessel 1841" 6377397.155 299.15281 mempty
tangent = Geodetic (dms 52 9 22.178) (dms 5 23 15.500) 0 ellipse
origin = GridOffset 155000 463000 0
-- | Standard steregraphic grid for point tests in the Southern Hemisphere.
--
-- This is the same as stereoGridN but with the tangent latitude and the false origin northings negated.
stereoGridS :: GridStereo LocalEllipsoid
stereoGridS = mkGridStereo tangent origin 0.9999079
where
ellipse = LocalEllipsoid "Bessel 1841" 6377397.155 299.15281 mempty
tangent = Geodetic (negate $ dms 52 9 22.178) (dms 5 23 15.500) 0 ellipse
origin = GridOffset (-155000) 463000 0
-- | Data for the stereographic tests taken from
-- <http://ftp.stu.edu.tw/BSD/NetBSD/pkgsrc/distfiles/epsg-6.11/G7-2.pdf>
stereographicToGridN :: Expectation
stereographicToGridN = sameGrid g1 g1'
where
p1 = Geodetic (dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridN
g1 = GridPoint 196105.283 557057.739 0 stereoGridN
g1' = toGrid stereoGridN p1
stereographicFromGridN :: Expectation
stereographicFromGridN = samePlace p1 p1'
where
p1 = Geodetic (dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridN
g1 = GridPoint 196105.283 557057.739 0 stereoGridN
p1' = fromGrid g1
stereographicToGridS :: Expectation
stereographicToGridS = sameGrid g1 g1'
where
p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridS
g1 = GridPoint (-196105.283) 557057.739 0 stereoGridS
g1' = toGrid stereoGridS p1
stereographicFromGridS :: Expectation
stereographicFromGridS = samePlace p1 p1'
where
p1 = Geodetic (negate $ dms 53 0 0) (dms 6 0 0) 0 $ gridEllipsoid stereoGridS
g1 = GridPoint (-196105.283) 557057.739 0 stereoGridS
p1' = fromGrid g1
-- | Check the round trip for a stereographic projection.
prop_stereographic :: GridPoint (GridStereo LocalEllipsoid) -> Property
prop_stereographic p =
let g = fromGrid p
r = toGrid (gridBasis p) g
in counterexample ("p = " ++ show p ++ "\ng = " ++ show g ++ "\nr = " ++ show r) $
closeGrid p r
-- | A ray at distance zero returns its original arguments.
prop_rayPath1 :: Ray WGS84 -> Bool
prop_rayPath1 r@(Ray pt b e) =
samePlace' pt pt1 && sameAngle' b b1 && sameAngle' e e1
where (pt1,b1,e1) = pathFunc (getRay r) 0
type ContinuityTest e = Geodetic e -> Bearing -> Azimuth -> Distance -> Distance -> Property
type ContinuityTest1 e = Geodetic e -> Bearing -> Distance2 -> Distance2 -> Property
-- | Many paths can be specified by a start point, bearing and azimuth,
-- and have the property that any (point,bearing,azimuth) triple on
-- the path will specify the same path with a distance offset.
prop_pathContinuity :: (Ellipsoid e) =>
(Geodetic e -> Double -> Double -> Path e) -> ContinuityTest e
prop_pathContinuity pf pt0 (Bearing b0) (Azimuth a0) (Distance d1) (Distance d2) =
counterexample (show ((pt2, Bearing b2, Azimuth a2), (pt3, Bearing b3, Azimuth a3))) $
pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>
closeEnough' pt2 pt3 && sameAngle' b2 b3 && sameAngle' a2 a3
where
path0 = pf pt0 b0 a0
(pt1, b1, a1) = pathFunc path0 d1
path1 = pf pt1 b1 a1
(pt2, b2, a2) = pathFunc path1 d2
(pt3, b3, a3) = pathFunc path0 (d1 + d2) -- Points 2 and 3 should be the same.
-- | For continuity testing of ground-based paths (azimuth & altitude always zero)
-- where lower accuracy is required.
prop_pathContinuity1 :: (Ellipsoid e) => (Geodetic e -> Double -> Path e) -> ContinuityTest1 e
prop_pathContinuity1 pf pt0 (Bearing b0) (Distance2 d1) (Distance2 d2) =
counterexample (show ((pt2, Bearing b2), (pt3, Bearing b3))) $
pathValidAt path0 d1 && pathValidAt path0 d2 && pathValidAt path0 (d1+d2) ==>
closeEnough' pt2 pt3 && sameAngle' b2 b3
where
path0 = pf pt0 b0
(pt1, b1, _) = pathFunc path0 d1
path1 = pf pt1 b1
(pt2, b2, _) = pathFunc path1 d2
(pt3, b3, _) = pathFunc path0 (d1 + d2) -- Points 2 and 3 should be the same.
-- | A point on a ray will continue along the same ray, and hence give the same points.
prop_rayContinuity :: ContinuityTest WGS84
prop_rayContinuity = prop_pathContinuity rayPath
-- | A ray bisected to an altitude will give that altitude.
-- This is a test of bisection rather than rays.
prop_rayBisect :: Ray WGS84 -> Altitude -> Bool
prop_rayBisect r (Altitude height) =
case bisect ray0 f 1e-2 0 (1000 * kilometer) of
Nothing -> False
Just d -> let (g, _, _) = pathFunc ray0 d in abs (altitude g - height) < 1e-2
where
f g = compare (altitude g) height
ray0 = getRay r
-- | A point on a rhumb line will continue along the same rhumb.
prop_rhumbContinuity :: ContinuityTest1 WGS84
prop_rhumbContinuity = prop_pathContinuity1 rhumbPath
-- | Two rhumb paths intersect at the same place.
prop_rhumbIntersect :: RhumbPaths2 -> Property
prop_rhumbIntersect rp =
case intersect 0 0 0.1 100 path1 path2 of
Just (d1, d2) ->
let (pt1, _, _) = pathFunc path1 d1
(pt2, _, _) = pathFunc path2 d2
in counterexample (show (pt1, pt2)) $ label "Intersection" $ samePlace pt1 pt2
Nothing -> label "No intersection" True
where
(path1, path2) = mk2RhumbPaths rp
-- | Copied from `Test.Hspec.Expectations` source. It ought to be exported from there.
expectTrue :: HasCallStack => String -> Bool -> Expectation
expectTrue msg b = unless b (expectationFailure msg)