sharpetronics/QGIS
1
0
Fork 0
mirror of https://github.com/qgis/QGIS.git synced 2025-04-17 00:04:02 -04:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4828 from nyalldawson/fix_16820

Browse Source 
Fix incorrect area calculation in corner cases (fix #16820)
This commit is contained in:
Nyall Dawson 2017-07-10 10:06:41 +10:00 committed by GitHub
commit 61bb2ac0b3
2 changed files with 40 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
src/core/qgsdistancearea.cpp
View File

@ -678,6 +678,14 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
double Qbar1, Qbar2; double Qbar1, Qbar2;
double area; 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 ); QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
if ( !willUseEllipsoid() ) if ( !willUseEllipsoid() )
{ {
@ -708,11 +716,28 @@ double QgsDistanceArea::computePolygonArea( const QList<QgsPointXY> &points ) co
x1 += m_TwoPI; x1 += m_TwoPI;
dx = x2 - x1; dx = x2 - x1;
dy = y2 - y1;
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 ) ); area += dx * ( m_Qp - getQ( y2 ) );
dy = y2 - y1; /* original:
if ( !qgsDoubleNear( dy, 0.0 ) ) * area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 ); */
}
} }
if ( ( area *= m_AE ) < 0.0 ) if ( ( area *= m_AE ) < 0.0 )
area = -area; area = -area;

11
tests/src/core/testqgsdistancearea.cpp
View File

@ -46,6 +46,7 @@ class TestQgsDistanceArea: public QObject
void measureAreaAndUnits(); void measureAreaAndUnits();
void emptyPolygon(); void emptyPolygon();
void regression14675(); void regression14675();
void regression16820();
}; };
@ -379,6 +380,16 @@ void TestQgsDistanceArea::regression14675()
QGSCOMPARENEAR( calc.measureArea( geom ), 0.83301, 0.02 ); 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 ) QGSTEST_MAIN( TestQgsDistanceArea )
#include "testqgsdistancearea.moc" #include "testqgsdistancearea.moc"