Skip to content

Commit

Permalink
Fix incorrect area calculation in corner cases (fix #16820)
Browse files Browse the repository at this point in the history
In certain circumstances very proximal nodes could cause instability
in the ellipsoidal area calculation.

Port (slightly tweaked) fix from grass changeset 71167 for same issue,
and add a unit test
  • Loading branch information
nyalldawson committed Jul 9, 2017
1 parent dae0bf8 commit d850393
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 4 deletions.
33 changes: 29 additions & 4 deletions src/core/qgsdistancearea.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -678,6 +678,14 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
double Qbar1, Qbar2;
double area;

/* GRASS comment: threshold for dy, should be between 1e-4 and 1e-7
* QGIS note: while the grass comment states that thres should be between 1e-4->1e-7,
* a value of 1e-7 caused TestQgsDistanceArea::regression14675() to regress
* The maximum threshold possible which permits regression14675() to pass
* was found to be ~0.7e-7.
*/
const double thresh = 0.7e-7;

QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
if ( !willUseEllipsoid() )
{
Expand Down Expand Up @@ -708,11 +716,28 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
x1 += m_TwoPI;

dx = x2 - x1;
area += dx * ( m_Qp - getQ( y2 ) );

dy = y2 - y1;
if ( !qgsDoubleNear( dy, 0.0 ) )
area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
if ( qAbs( dy ) > thresh )
{
/* account for different latitudes y1, y2 */
area += dx * ( m_Qp - ( Qbar2 - Qbar1 ) / dy );
}
else
{
/* latitudes y1, y2 are (nearly) identical */

/* if y2 becomes similar to y1, i.e. y2 -> y1
* Qbar2 - Qbar1 -> 0 and dy -> 0
* (Qbar2 - Qbar1) / dy -> ?
* (Qbar2 - Qbar1) / dy should approach Q((y1 + y2) / 2)
* Metz 2017
*/
area += dx * ( m_Qp - getQ( y2 ) );

/* original:
* area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
*/
}
}
if ( ( area *= m_AE ) < 0.0 )
area = -area;
Expand Down
11 changes: 11 additions & 0 deletions tests/src/core/testqgsdistancearea.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,7 @@ class TestQgsDistanceArea: public QObject
void measureAreaAndUnits();
void emptyPolygon();
void regression14675();
void regression16820();

};

Expand Down Expand Up @@ -379,6 +380,16 @@ void TestQgsDistanceArea::regression14675()
QGSCOMPARENEAR( calc.measureArea( geom ), 0.83301, 0.02 );
}

void TestQgsDistanceArea::regression16820()
{
QgsDistanceArea calc;
calc.setEllipsoid( QStringLiteral( "WGS84" ) );
calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:32634" ) ) );
QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((110250.54038314701756462 5084495.57398066483438015, 110243.46975068224128336 5084507.17200060561299324, 110251.23908144699817058 5084506.68309532757848501, 110251.2394439501222223 5084506.68307251576334238, 110250.54048078990308568 5084495.57553235255181789, 110250.54038314701756462 5084495.57398066483438015))" ) ) );
//lots of tolerance here - the formulas get quite unstable with small areas due to division by very small floats
QGSCOMPARENEAR( calc.measureArea( geom ), 43.183369, 0.02 );
}

QGSTEST_MAIN( TestQgsDistanceArea )
#include "testqgsdistancearea.moc"

Expand Down

0 comments on commit d850393

Please sign in to comment.