Skip to content
Permalink
Browse files

Avoid loss of precision when reading ellipsoids from proj

  • Loading branch information
nyalldawson committed Mar 27, 2019
1 parent 583f3d3 commit 163f046945bde1073d707d40d244c3d3c545a043
Showing with 34 additions and 20 deletions.
  1. +7 −11 src/core/qgsellipsoidutils.cpp
  2. +27 −9 tests/src/python/test_qgsellipsoidutils.py
@@ -357,22 +357,18 @@ QList<QgsEllipsoidUtils::EllipsoidDefinition> QgsEllipsoidUtils::definitions()
def.description = QStringLiteral( "%1 (%2:%3)" ).arg( name, authority, code );

double semiMajor, semiMinor, invFlattening;
int semiMinorComputed;
int semiMinorComputed = 0;
if ( proj_ellipsoid_get_parameters( context, ellipsoid, &semiMajor, &semiMinor, &semiMinorComputed, &invFlattening ) )
{
def.parameters.semiMajor = semiMajor;
if ( semiMinor > 0 )
{
def.parameters.semiMinor = semiMinor;
def.parameters.inverseFlattening = def.parameters.semiMajor / ( def.parameters.semiMajor - def.parameters.semiMinor );
def.parameters.semiMinor = semiMinor;
def.parameters.inverseFlattening = invFlattening;
if ( !semiMinorComputed )
def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +b=%2 +no_defs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.semiMinor, 0, 'g', 17 ) );
}
else
{
def.parameters.inverseFlattening = invFlattening;
def.parameters.semiMinor = def.parameters.semiMajor * ( 1 - def.parameters.inverseFlattening );
else if ( !qgsDoubleNear( def.parameters.inverseFlattening, 0.0 ) )
def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +rf=%2 +no_defs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ).arg( def.parameters.inverseFlattening, 0, 'g', 17 ) );
}
else
def.parameters.crs = QgsCoordinateReferenceSystem::fromProj4( QStringLiteral( "+proj=longlat +a=%1 +no_defs" ).arg( def.parameters.semiMajor, 0, 'g', 17 ) );
}
else
{
@@ -13,7 +13,6 @@
__revision__ = '$Format:%H$'

import qgis # NOQA
import math
from qgis.core import (QgsEllipsoidUtils,
QgsProjUtils)
from qgis.testing import start_app, unittest
@@ -40,25 +39,44 @@ def testParams(self):
if QgsProjUtils.projVersionMajor() < 6:
self.assertEqual(params.crs.authid(), 'EPSG:4030')
else:
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=6378137 +b=6356752.3142451793 +no_defs')
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=6378137 +rf=298.25722356300003 +no_defs')

for i in range(2):
params = QgsEllipsoidUtils.ellipsoidParameters("Ganymede2000")
self.assertTrue(params.valid)
self.assertEqual(params.semiMajor, 2632400.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(params.semiMinor, 2632350.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(params.inverseFlattening, 52648.0 if QgsProjUtils.projVersionMajor() < 6 else math.inf)
self.assertEqual(params.inverseFlattening, 52648.0 if QgsProjUtils.projVersionMajor() < 6 else 0)
self.assertFalse(params.useCustomParameters)
self.assertEqual(params.crs.authid(), '')
if QgsProjUtils.projVersionMajor() < 6:
self.assertEqual(params.crs.authid(), '')
else:
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=2632345 +no_defs')

if QgsProjUtils.projVersionMajor() >= 6:
params = QgsEllipsoidUtils.ellipsoidParameters("ESRI:107916")
self.assertTrue(params.valid)
self.assertEqual(params.semiMajor, 2632400.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(params.semiMinor, 2632350.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(params.inverseFlattening, 52648.0 if QgsProjUtils.projVersionMajor() < 6 else math.inf)
self.assertEqual(params.semiMajor, 2632345.0)
self.assertEqual(params.semiMinor, 2632345.0)
self.assertEqual(params.inverseFlattening, 0)
self.assertFalse(params.useCustomParameters)
self.assertEqual(params.crs.authid(), '')
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=2632345 +no_defs')

params = QgsEllipsoidUtils.ellipsoidParameters("EPSG:7001")
self.assertTrue(params.valid)
self.assertEqual(params.semiMajor, 6377563.396)
self.assertEqual(params.semiMinor, 6356256.909237285)
self.assertEqual(params.inverseFlattening, 299.3249646)
self.assertFalse(params.useCustomParameters)
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=6377563.3959999997 +rf=299.32496459999999 +no_defs')

params = QgsEllipsoidUtils.ellipsoidParameters("EPSG:7008")
self.assertTrue(params.valid)
self.assertEqual(params.semiMajor, 6378206.4)
self.assertEqual(params.semiMinor, 6356583.8)
self.assertEqual(params.inverseFlattening, 294.9786982138982)
self.assertFalse(params.useCustomParameters)
self.assertEqual(params.crs.toProj4(), '+proj=longlat +a=6378206.4000000004 +b=6356583.7999999998 +no_defs')

# using parameters
for i in range(2):
@@ -89,7 +107,7 @@ def testDefinitions(self):
self.assertTrue(gany_defs.parameters.valid)
self.assertEqual(gany_defs.parameters.semiMajor, 2632400.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(gany_defs.parameters.semiMinor, 2632350.0 if QgsProjUtils.projVersionMajor() < 6 else 2632345.0)
self.assertEqual(gany_defs.parameters.inverseFlattening, 52648.0 if QgsProjUtils.projVersionMajor() < 6 else math.inf)
self.assertEqual(gany_defs.parameters.inverseFlattening, 52648.0 if QgsProjUtils.projVersionMajor() < 6 else 0.0)
self.assertFalse(gany_defs.parameters.useCustomParameters)
self.assertEqual(gany_defs.parameters.crs.authid(), '')

0 comments on commit 163f046

Please sign in to comment.
You can’t perform that action at this time.