1 /***************************************************************************
2      testqgsdistancearea.cpp
3      --------------------------------------
4     Date                 : Tue 14 Aug 2012
5     Copyright            : (C) 2012 by Magnus Homann
6     Email                : magnus at homann dot se
7  ***************************************************************************
8  *                                                                         *
9  *   This program is free software; you can redistribute it and/or modify  *
10  *   it under the terms of the GNU General Public License as published by  *
11  *   the Free Software Foundation; either version 2 of the License, or     *
12  *   (at your option) any later version.                                   *
13  *                                                                         *
14  ***************************************************************************/
15 #include "qgstest.h"
16 #include <QFile>
17 #include <QTextStream>
18 #include <QObject>
19 #include <QString>
20 #include <QStringList>
21 #include <qgsapplication.h>
22 //header for class being tested
23 #include <qgsdistancearea.h>
24 #include <qgspoint.h>
25 #include "qgslogger.h"
26 #include "qgsgeometryfactory.h"
27 #include "qgsgeometry.h"
28 #include "qgis.h"
29 #include "qgsproject.h"
30 #include <memory>
31 
32 class TestQgsDistanceArea: public QObject
33 {
34 
35     Q_OBJECT
36   private slots:
37     void initTestCase();
38     void cleanupTestCase();
39     void basic();
40     void noEllipsoid();
41     void cache();
42     void test_distances();
43     void regression13601();
44     void collections();
45     void measureUnits();
46     void measureAreaAndUnits();
47     void emptyPolygon();
48     void regression14675();
49     void regression16820();
50 
51 };
52 
initTestCase()53 void TestQgsDistanceArea::initTestCase()
54 {
55   //
56   // Runs once before any tests are run
57   //
58   // init QGIS's paths - true means that all path will be inited from prefix
59   QgsApplication::init();
60   QgsApplication::initQgis();
61   QgsApplication::showSettings();
62 }
63 
cleanupTestCase()64 void TestQgsDistanceArea::cleanupTestCase()
65 {
66   QgsApplication::exitQgis();
67 }
68 
basic()69 void TestQgsDistanceArea::basic()
70 {
71   QgsPointXY p1( 1.0, 3.0 ), p2( -2.0, -1.0 );
72   QgsDistanceArea daA;
73   double resultA, resultB, resultC;
74 
75   daA.setEllipsoid( geoNone() );
76   resultA = daA.measureLine( p1, p2 );
77   QCOMPARE( resultA, 5.0 );
78 
79   // Now, on an ellipsoid. Always less?
80   daA.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:32442" ) ), QgsProject::instance()->transformContext() );
81   daA.setEllipsoid( QStringLiteral( "WGS84" ) );
82   resultA = daA.measureLine( p1, p2 );
83   QVERIFY( resultA < 5.0 );
84 
85   // Test copy constructor
86   QgsDistanceArea daB( daA );
87   resultB = daB.measureLine( p1, p2 );
88   QCOMPARE( resultA, resultB );
89 
90   // Different Ellipsoid
91   daB.setEllipsoid( QStringLiteral( "WGS72" ) );
92   resultB = daB.measureLine( p1, p2 );
93   QVERIFY( ! qFuzzyCompare( resultA, resultB ) );
94 
95   // Test assignment
96   std::shared_ptr<QgsDistanceArea> daC( new QgsDistanceArea );
97   *daC = daB;
98   resultC = daC->measureLine( p1, p2 );
99   QCOMPARE( resultB, resultC );
100 
101   // Use parameter setting of ellipsoid radii (from WGS72 )
102   daA.setEllipsoid( 6378135.0, 6378135.0 - ( 6378135.0 / 298.26 ) );
103   resultA = daA.measureLine( p1, p2 );
104   QCOMPARE( resultA, resultB );
105 }
106 
noEllipsoid()107 void TestQgsDistanceArea::noEllipsoid()
108 {
109   QgsDistanceArea da;
110   da.setEllipsoid( geoNone() );
111   QVERIFY( !da.willUseEllipsoid() );
112   QCOMPARE( da.ellipsoid(), geoNone() );
113 }
114 
cache()115 void TestQgsDistanceArea::cache()
116 {
117   // test that ellipsoid can be retrieved correctly from cache
118   QgsDistanceArea da;
119 
120   // warm cache
121   QVERIFY( da.setEllipsoid( QStringLiteral( "Ganymede2000" ) ) );
122   QVERIFY( da.willUseEllipsoid() );
123   // different definition under proj 6, neither is incorrect!
124 #if PROJ_VERSION_MAJOR>=6
125   QCOMPARE( da.ellipsoidSemiMajor(), 2632345.0 );
126   QCOMPARE( da.ellipsoidSemiMinor(), 2632345.0 );
127 #else
128   QCOMPARE( da.ellipsoidSemiMajor(), 2632400.0 );
129   QCOMPARE( da.ellipsoidSemiMinor(), 2632350.0 );
130   QCOMPARE( da.ellipsoidInverseFlattening(), 52648.0 );
131 #endif
132   QCOMPARE( da.ellipsoid(), QStringLiteral( "Ganymede2000" ) );
133 
134   // a second time, so ellipsoid is fetched from cache
135   QgsDistanceArea da2;
136   QVERIFY( da2.setEllipsoid( QStringLiteral( "Ganymede2000" ) ) );
137   QVERIFY( da2.willUseEllipsoid() );
138 #if PROJ_VERSION_MAJOR>=6
139   QCOMPARE( da2.ellipsoidSemiMajor(), 2632345.0 );
140   QCOMPARE( da2.ellipsoidSemiMinor(), 2632345.0 );
141 #else
142   QCOMPARE( da2.ellipsoidSemiMajor(), 2632400.0 );
143   QCOMPARE( da2.ellipsoidSemiMinor(), 2632350.0 );
144   QCOMPARE( da2.ellipsoidInverseFlattening(), 52648.0 );
145 #endif
146   QCOMPARE( da2.ellipsoid(), QStringLiteral( "Ganymede2000" ) );
147 
148   // using parameters
149   QgsDistanceArea da3;
150   QVERIFY( da3.setEllipsoid( QStringLiteral( "PARAMETER:2631400:2341350" ) ) );
151   QVERIFY( da3.willUseEllipsoid() );
152   QCOMPARE( da3.ellipsoidSemiMajor(), 2631400.0 );
153   QCOMPARE( da3.ellipsoidSemiMinor(), 2341350.0 );
154   QGSCOMPARENEAR( da3.ellipsoidInverseFlattening(), 9.07223, 0.00001 );
155   QCOMPARE( da3.ellipsoid(), QStringLiteral( "PARAMETER:2631400:2341350" ) );
156 
157   // again, to check parameters with cache
158   QgsDistanceArea da4;
159   QVERIFY( da4.setEllipsoid( QStringLiteral( "PARAMETER:2631400:2341350" ) ) );
160   QVERIFY( da4.willUseEllipsoid() );
161   QCOMPARE( da4.ellipsoidSemiMajor(), 2631400.0 );
162   QCOMPARE( da4.ellipsoidSemiMinor(), 2341350.0 );
163   QGSCOMPARENEAR( da4.ellipsoidInverseFlattening(), 9.07223, 0.00001 );
164   QCOMPARE( da4.ellipsoid(), QStringLiteral( "PARAMETER:2631400:2341350" ) );
165 
166   // invalid
167   QgsDistanceArea da5;
168   QVERIFY( !da5.setEllipsoid( QStringLiteral( "MyFirstEllipsoid" ) ) );
169   QVERIFY( !da5.willUseEllipsoid() );
170   QCOMPARE( da5.ellipsoid(), geoNone() );
171 
172   // invalid again, should be cached
173   QgsDistanceArea da6;
174   QVERIFY( !da6.setEllipsoid( QStringLiteral( "MyFirstEllipsoid" ) ) );
175   QVERIFY( !da6.willUseEllipsoid() );
176   QCOMPARE( da6.ellipsoid(), geoNone() );
177 }
178 
test_distances()179 void TestQgsDistanceArea::test_distances()
180 {
181   // Read the file of Geod data
182   // Column 0 (first) is latitude point 1
183   // Column 1 is longitude point 1
184   // Column 3 is latitude point 2
185   // Column 4 is longitude point 3
186   // Column 6 is the resulting distance in meters on the WGS84 ellipsoid
187   // Note: lat is north/south, so the QgsPointXY should be ( long, lat )
188   // See http://geographiclib.sourceforge.net/html/geodesic.html#testgeod
189 
190   // Set up DA
191   QgsDistanceArea myDa;
192   myDa.setSourceCrs( QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4030" ) ), QgsProject::instance()->transformContext() );
193   myDa.setEllipsoid( QStringLiteral( "WGS84" ) );
194 
195   QString myFileName = QStringLiteral( TEST_DATA_DIR ) + "/GeodTest-nano.dat";
196 
197   QFile myFile( myFileName );
198   if ( ! myFile.open( QIODevice::ReadOnly | QIODevice::Text ) )
199   {
200     QFAIL( "Couldn't open file" );
201     return;
202   }
203   QTextStream in( & myFile );
204   while ( !in.atEnd() )
205   {
206     QString line = in.readLine();
207     // Some test points (antipodal) does not converge with the chosen algorithm!
208     // See calcaulator at http://www.movable-type.co.uk/scripts/latlong-vincenty.html
209     // These are commented out.
210     if ( line[0] != '#' )
211     {
212       QStringList myLineList = line.split( ' ' ); // Split fields on space.
213       // Create points
214       QgsPointXY p1( myLineList[1].toDouble(), myLineList[0].toDouble() );
215       QgsPointXY p2( myLineList[4].toDouble(), myLineList[3].toDouble() );
216       double result = myDa.measureLine( p1, p2 );
217       // QgsDebugMsg( QStringLiteral( "Distance from %1 to %2 is %3" ).arg( p1.toString( 15 ) ).arg( p2.toString( 15 ) ).arg( result, 0, 'g', 15 ) );
218       // QgsDebugMsg( QStringLiteral( "Distance should be %1" ).arg( myLineList[6] ) );
219       // Check result is less than 0.5mm from expected.
220       QGSCOMPARENEAR( result, myLineList[6].toDouble(), 0.0005 );
221     }
222   }
223 
224 }
225 
regression13601()226 void TestQgsDistanceArea::regression13601()
227 {
228   //test regression #13601
229   QgsDistanceArea calc;
230   calc.setEllipsoid( QStringLiteral( "NONE" ) );
231   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3148" ) ), QgsProject::instance()->transformContext() );
232   QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((252000 1389000, 265000 1389000, 265000 1385000, 252000 1385000, 252000 1389000))" ) ).release() );
233   QGSCOMPARENEAR( calc.measureArea( geom ), 52000000, 0.0001 );
234 }
235 
collections()236 void TestQgsDistanceArea::collections()
237 {
238   //test measuring for collections
239   QgsDistanceArea myDa;
240   myDa.setSourceCrs( QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4030" ) ), QgsProject::instance()->transformContext() );
241   myDa.setEllipsoid( QStringLiteral( "WGS84" ) );
242 
243   //collection of lines, should be sum of line length
244   QgsGeometry lines( QgsGeometryFactory::geomFromWkt( QStringLiteral( "GeometryCollection( LineString(0 36.53, 5.76 -48.16), LineString(0 25.54, 24.20 36.70) )" ) ).release() );
245   double result = myDa.measureLength( lines );
246   QGSCOMPARENEAR( result, 12006159, 1 );
247   result = myDa.measureArea( lines );
248   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
249 
250   //collection of polygons
251   QgsGeometry polys( QgsGeometryFactory::geomFromWkt( QStringLiteral( "GeometryCollection( Polygon((0 36.53, 5.76 -48.16, 0 25.54, 0 36.53)), Polygon((10 20, 15 20, 15 10, 10 20)) )" ) ).release() );
252   result = myDa.measureArea( polys );
253   QGSCOMPARENEAR( result, 663136985074LL, 1 );
254   result = myDa.measureLength( polys );
255   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
256 
257   //mixed collection
258   QgsGeometry mixed( QgsGeometryFactory::geomFromWkt( QStringLiteral( "GeometryCollection( LineString(0 36.53, 5.76 -48.16), LineString(0 25.54, 24.20 36.70), Polygon((0 36.53, 5.76 -48.16, 0 25.54, 0 36.53)), Polygon((10 20, 15 20, 15 10, 10 20)) )" ) ).release() );
259   //measure area specifically
260   result = myDa.measureArea( mixed );
261   QGSCOMPARENEAR( result, 663136985075LL, 1 );
262   //measure length
263   result = myDa.measureLength( mixed );
264   QGSCOMPARENEAR( result, 12006159, 1 );
265 }
266 
measureUnits()267 void TestQgsDistanceArea::measureUnits()
268 {
269   //test regression #13610
270   QgsDistanceArea calc;
271   calc.setEllipsoid( QStringLiteral( "NONE" ) );
272   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:2272" ) ), QgsProject::instance()->transformContext() );
273   QgsUnitTypes::DistanceUnit units;
274   QgsPointXY p1( 1341683.9854275715, 408256.9562717728 );
275   QgsPointXY p2( 1349321.7807031618, 408256.9562717728 );
276 
277   double result = calc.measureLine( p1, p2 );
278   units = calc.lengthUnits();
279   //no OTF, result will be in CRS unit (feet)
280   QCOMPARE( units, QgsUnitTypes::DistanceFeet );
281   QGSCOMPARENEAR( result, 7637.7952755903825, 0.001 );
282 
283   calc.setEllipsoid( QStringLiteral( "WGS84" ) );
284   units = calc.lengthUnits();
285   result = calc.measureLine( p1, p2 );
286   //OTF, result will be in meters
287   QCOMPARE( units, QgsUnitTypes::DistanceMeters );
288   QGSCOMPARENEAR( result, 2328.0988253106957, 0.001 );
289 }
290 
measureAreaAndUnits()291 void TestQgsDistanceArea::measureAreaAndUnits()
292 {
293   QgsDistanceArea da;
294   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext() );
295   da.setEllipsoid( QStringLiteral( "NONE" ) );
296   QgsPolylineXY ring;
297   ring << QgsPointXY( 0, 0 )
298        << QgsPointXY( 1, 0 )
299        << QgsPointXY( 1, 1 )
300        << QgsPointXY( 2, 1 )
301        << QgsPointXY( 2, 2 )
302        << QgsPointXY( 0, 2 )
303        << QgsPointXY( 0, 0 );
304   QgsPolygonXY poly;
305   poly << ring;
306 
307   QgsGeometry polygon( QgsGeometry::fromPolygonXY( poly ) );
308 
309   // We check both the measured area AND the units, in case the logic regarding
310   // ellipsoids and units changes in future
311   double area = da.measureArea( polygon );
312   QgsUnitTypes::AreaUnit units = da.areaUnits();
313 
314   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
315 
316   QVERIFY( ( qgsDoubleNear( area, 3.0, 0.00000001 ) && units == QgsUnitTypes::AreaSquareDegrees )
317            || ( qgsDoubleNear( area, 37176087091.5, 0.1 ) && units == QgsUnitTypes::AreaSquareMeters ) );
318 
319   da.setEllipsoid( QStringLiteral( "WGS84" ) );
320   area = da.measureArea( polygon );
321   units = da.areaUnits();
322   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
323   // should always be in Meters Squared
324   QGSCOMPARENEAR( area, 36918093794.1, 0.1 );
325   QCOMPARE( units, QgsUnitTypes::AreaSquareMeters );
326 
327   // test converting the resultant area
328   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareMiles );
329   QGSCOMPARENEAR( area, 14254.155703, 0.001 );
330 
331   // now try with a source CRS which is in feet
332   ring.clear();
333   ring << QgsPointXY( 1850000, 4423000 )
334        << QgsPointXY( 1851000, 4423000 )
335        << QgsPointXY( 1851000, 4424000 )
336        << QgsPointXY( 1852000, 4424000 )
337        << QgsPointXY( 1852000, 4425000 )
338        << QgsPointXY( 1851000, 4425000 )
339        << QgsPointXY( 1850000, 4423000 );
340   poly.clear();
341   poly << ring;
342   polygon = QgsGeometry::fromPolygonXY( poly );
343 
344 #if PROJ_VERSION_MAJOR>=6
345   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "ESRI:102635" ) ), QgsProject::instance()->transformContext() );
346 #else
347   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:102635" ) ), QgsProject::instance()->transformContext() );
348 #endif
349   da.setEllipsoid( QStringLiteral( "NONE" ) );
350   // measurement should be in square feet
351   area = da.measureArea( polygon );
352   units = da.areaUnits();
353   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
354   QGSCOMPARENEAR( area, 2000000, 0.001 );
355   QCOMPARE( units, QgsUnitTypes::AreaSquareFeet );
356 
357   // test converting the resultant area
358   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareYards );
359   QGSCOMPARENEAR( area, 222222.2222, 0.001 );
360 
361   da.setEllipsoid( QStringLiteral( "WGS84" ) );
362   // now should be in Square Meters again
363   area = da.measureArea( polygon );
364   units = da.areaUnits();
365   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
366   QGSCOMPARENEAR( area, 185818.590966, 1.0 );
367   QCOMPARE( units, QgsUnitTypes::AreaSquareMeters );
368 
369   // test converting the resultant area
370   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareYards );
371   QgsDebugMsg( QStringLiteral( "measured %1 in sq yrds" ).arg( area ) );
372   QGSCOMPARENEAR( area, 222237.185213, 1.0 );
373 }
374 
emptyPolygon()375 void TestQgsDistanceArea::emptyPolygon()
376 {
377   QgsDistanceArea da;
378   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext() );
379   da.setEllipsoid( QStringLiteral( "WGS84" ) );
380 
381   //test that measuring an empty polygon doesn't crash
382   da.measurePolygon( QVector< QgsPointXY >() );
383 }
384 
regression14675()385 void TestQgsDistanceArea::regression14675()
386 {
387   //test regression #14675
388   QgsDistanceArea calc;
389   calc.setEllipsoid( QStringLiteral( "GRS80" ) );
390   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:2154" ) ), QgsProject::instance()->transformContext() );
391   QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((917593.5791854317067191 6833700.00807378999888897, 917596.43389983859378844 6833700.67099479306489229, 917599.53056440979707986 6833700.78673478215932846, 917593.5791854317067191 6833700.00807378999888897))" ) ).release() );
392   //lots of tolerance here - the formulas get quite unstable with small areas due to division by very small floats
393   QGSCOMPARENEAR( calc.measureArea( geom ), 0.833010, 0.03 );
394 }
395 
regression16820()396 void TestQgsDistanceArea::regression16820()
397 {
398   QgsDistanceArea calc;
399   calc.setEllipsoid( QStringLiteral( "WGS84" ) );
400   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:32634" ) ), QgsProject::instance()->transformContext() );
401   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))" ) ).release() );
402   //lots of tolerance here - the formulas get quite unstable with small areas due to division by very small floats
403   QGSCOMPARENEAR( calc.measureArea( geom ), 43.3280029296875, 0.2 );
404 }
405 
406 QGSTEST_MAIN( TestQgsDistanceArea )
407 #include "testqgsdistancearea.moc"
408 
409 
410 
411 
412