Inaccuracies in transforming to/from TransverseMercator projection

[martyall]

I've looked in your test suite and I see evidence (and have tested myself) that moving back and forth Geodetic WGS84 -> GridStereo WGS84 -> Geodetic WGS84 gives the identity (up to epsilon) for points that are near the tangent point of the projection.

But there is no evidence in the tests that Geodetic WGS84 -> GridTM WGS84 -> Geodetic WGS84 gives the identity transformation for points near the origin, and I have had problems writing tests of my own for this. I am seeing that while longitude is unchanged, latitude is distorted in a way that I wouldn't expect for points being so close to the origin. Here is a minimal example of a test that fails -- I only introduce the Coordinates type as a convenience, but the first test shows that the transformation to/from Geodetic is faithful:

module LocalizationSpec where

import           Control.Monad                     (forM_, replicateM, unless)
import qualified Data.Random                       as DR
import           Geodetics.Ellipsoids
import           Geodetics.Geodetic
import           Geodetics.Grid
import           Geodetics.TransverseMercator
import           Numeric.Units.Dimensional.Prelude (degree, meter, (*~), (/~),
                                                    _1)
import           System.Random.MWC                 (createSystemRandom)
import           Test.Hspec


--------------------------------------------------------------------------------
-- Converting to/from the Geodetics
--------------------------------------------------------------------------------

data Coordinates = Coordinates
  { lat :: Double
  , lon :: Double
  } deriving (Eq)

instance Show Coordinates where
    show Coordinates{..} = show (lat,lon)

randomCoordinates :: Int -> IO [Coordinates]
randomCoordinates n = do
  mwc <- createSystemRandom
  flip DR.runRVar mwc $ replicateM n $ do
    x <- DR.uniform @Double (-1) 1
    y <- DR.uniform @Double (-1) 1
    return Coordinates
      { lat = x
      , lon = y
      }

-- convert from latitude/longitude to a geodetic
-- at 0 altitude
latLonToGeodetic
  :: Coordinates
  -> Geodetic WGS84
latLonToGeodetic Coordinates{lat,lon} =
  Geodetic
    { latitude = lat *~ degree
    , longitude = lon *~ degree
    , geoAlt = 0 *~ meter
    , ellipsoid = WGS84
    }

-- convert from a geodetic to lat/lon
geodeticToLatLon
  :: Geodetic WGS84
  -> Coordinates
geodeticToLatLon Geodetic{..} =
  Coordinates
    { lat = latitude /~ degree
    , lon = longitude /~ degree
    }

--------------------------------------------------------------------------------
-- test
--------------------------------------------------------------------------------

spec :: Spec
spec = describe "Localization Spec" $ do
  let n = 100
      --- this doesn't hold on edge cases near the 180/0 line,
      --- but we're near the equator in the examples
      isBasicallyEqual :: Coordinates -> Coordinates -> Bool
      isBasicallyEqual c1 c2 =
        abs (lat c1 - lat c2) < 1E-5 &&
        abs (lon c1 - lon c2) < 1E-5

  it "Can convert to and from Geodetics" $ do
      zone <- randomCoordinates n
      let f c = geodeticToLatLon (latLonToGeodetic c) `isBasicallyEqual` c
      forM_ zone $ \c -> do
        let transformed = geodeticToLatLon . latLonToGeodetic $ c
        unless (transformed `isBasicallyEqual` c) $
          fail $ "Transformation Failed (original, transformed):\n" <> show (c, transformed)

  let -- | The positions are within 50 m.
      closeEnough :: (Ellipsoid e) => Geodetic e -> Geodetic e -> Bool
      closeEnough p1 p2 = geometricalDistance p1 p2 < 50 *~ meter

  it "Can convert geodetic to and from Grid" $ do
      cs <- randomCoordinates n
      let origin = latLonToGeodetic $ Coordinates 0 0
          gridBasis = mkGridTM origin mempty _1
      forM_ cs $ \c -> do
        let g = latLonToGeodetic c
            transformed = fromGrid . toGrid gridBasis $ g
        unless (transformed `closeEnough` g) $ do
          let failMsg = mconcat
                [ "\n"
                , "Test input: " <> show c <> "\n"
                , "Original geodetic: " <> show g <> "\n"
                , "Transformed geodetic: " <> show transformed
                , "\n"
                ]
          fail failMsg

It fails with the following ouput:

  test/LocalizationSpec.hs:86:3: 
  1) Localization, Localization Spec, Can convert geodetic to and from Grid
       uncaught exception: IOException of type UserError
       user error (
       Test input: (-0.9842999545277058,-0.3790717325408717)
       Original geodetic: 0° 59' 3.48" S, 0° 22' 44.66" W, 0.0 m WGS84
       Transformed geodetic: 0° 58' 39.76" S, 0° 22' 44.66" W, 0.0 m WGS84
       )

However, if I change to sample only positive random numbers for the latitude, then there is no problem and the tests pass. Am I missing something about how this is supposed to work?

[martyall]

I would note that if I change to the GridStereo WGS84 then everything passes (as is indicated by your test suite). But I would have thought that surely projecting to and from should leave me closer than 50m from where I started given that everything is happening in a 2° x 2° box around (0,0)

I have also tried random numbers in [-0.1, 0.1] for latitude and longitude and have the same errors but this in not really my use case.

I have also tried taking only positive numbers in [0,10] and this produces no errors. Is there a problem with using negative numbers here?

[PaulJohnson]

This definitely looks like a bug. I'll see what I can do.

[PaulJohnson]

Finally nailed it. Sorry it's taken so long. An iterative approximation was terminating early because the delta was negative. See #19.