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   const 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   QCOMPARE( da.ellipsoidSemiMajor(), 2632345.0 );
124   QCOMPARE( da.ellipsoidSemiMinor(), 2632345.0 );
125   QCOMPARE( da.ellipsoid(), QStringLiteral( "Ganymede2000" ) );
126 
127   // a second time, so ellipsoid is fetched from cache
128   QgsDistanceArea da2;
129   QVERIFY( da2.setEllipsoid( QStringLiteral( "Ganymede2000" ) ) );
130   QVERIFY( da2.willUseEllipsoid() );
131   QCOMPARE( da2.ellipsoidSemiMajor(), 2632345.0 );
132   QCOMPARE( da2.ellipsoidSemiMinor(), 2632345.0 );
133   QCOMPARE( da2.ellipsoid(), QStringLiteral( "Ganymede2000" ) );
134 
135   // using parameters
136   QgsDistanceArea da3;
137   QVERIFY( da3.setEllipsoid( QStringLiteral( "PARAMETER:2631400:2341350" ) ) );
138   QVERIFY( da3.willUseEllipsoid() );
139   QCOMPARE( da3.ellipsoidSemiMajor(), 2631400.0 );
140   QCOMPARE( da3.ellipsoidSemiMinor(), 2341350.0 );
141   QGSCOMPARENEAR( da3.ellipsoidInverseFlattening(), 9.07223, 0.00001 );
142   QCOMPARE( da3.ellipsoid(), QStringLiteral( "PARAMETER:2631400:2341350" ) );
143 
144   // again, to check parameters with cache
145   QgsDistanceArea da4;
146   QVERIFY( da4.setEllipsoid( QStringLiteral( "PARAMETER:2631400:2341350" ) ) );
147   QVERIFY( da4.willUseEllipsoid() );
148   QCOMPARE( da4.ellipsoidSemiMajor(), 2631400.0 );
149   QCOMPARE( da4.ellipsoidSemiMinor(), 2341350.0 );
150   QGSCOMPARENEAR( da4.ellipsoidInverseFlattening(), 9.07223, 0.00001 );
151   QCOMPARE( da4.ellipsoid(), QStringLiteral( "PARAMETER:2631400:2341350" ) );
152 
153   // invalid
154   QgsDistanceArea da5;
155   QVERIFY( !da5.setEllipsoid( QStringLiteral( "MyFirstEllipsoid" ) ) );
156   QVERIFY( !da5.willUseEllipsoid() );
157   QCOMPARE( da5.ellipsoid(), geoNone() );
158 
159   // invalid again, should be cached
160   QgsDistanceArea da6;
161   QVERIFY( !da6.setEllipsoid( QStringLiteral( "MyFirstEllipsoid" ) ) );
162   QVERIFY( !da6.willUseEllipsoid() );
163   QCOMPARE( da6.ellipsoid(), geoNone() );
164 }
165 
test_distances()166 void TestQgsDistanceArea::test_distances()
167 {
168   // Read the file of Geod data
169   // Column 0 (first) is latitude point 1
170   // Column 1 is longitude point 1
171   // Column 3 is latitude point 2
172   // Column 4 is longitude point 3
173   // Column 6 is the resulting distance in meters on the WGS84 ellipsoid
174   // Note: lat is north/south, so the QgsPointXY should be ( long, lat )
175   // See http://geographiclib.sourceforge.net/html/geodesic.html#testgeod
176 
177   // Set up DA
178   QgsDistanceArea myDa;
179   myDa.setSourceCrs( QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4030" ) ), QgsProject::instance()->transformContext() );
180   myDa.setEllipsoid( QStringLiteral( "WGS84" ) );
181 
182   const QString myFileName = QStringLiteral( TEST_DATA_DIR ) + "/GeodTest-nano.dat";
183 
184   QFile myFile( myFileName );
185   if ( ! myFile.open( QIODevice::ReadOnly | QIODevice::Text ) )
186   {
187     QFAIL( "Couldn't open file" );
188     return;
189   }
190   QTextStream in( & myFile );
191   while ( !in.atEnd() )
192   {
193     QString line = in.readLine();
194     // Some test points (antipodal) does not converge with the chosen algorithm!
195     // See calcaulator at http://www.movable-type.co.uk/scripts/latlong-vincenty.html
196     // These are commented out.
197     if ( line[0] != '#' )
198     {
199       QStringList myLineList = line.split( ' ' ); // Split fields on space.
200       // Create points
201       const QgsPointXY p1( myLineList[1].toDouble(), myLineList[0].toDouble() );
202       const QgsPointXY p2( myLineList[4].toDouble(), myLineList[3].toDouble() );
203       const double result = myDa.measureLine( p1, p2 );
204       // QgsDebugMsg( QStringLiteral( "Distance from %1 to %2 is %3" ).arg( p1.toString( 15 ) ).arg( p2.toString( 15 ) ).arg( result, 0, 'g', 15 ) );
205       // QgsDebugMsg( QStringLiteral( "Distance should be %1" ).arg( myLineList[6] ) );
206       // Check result is less than 0.5mm from expected.
207       QGSCOMPARENEAR( result, myLineList[6].toDouble(), 0.0005 );
208     }
209   }
210 
211 }
212 
regression13601()213 void TestQgsDistanceArea::regression13601()
214 {
215   //test regression #13601
216   QgsDistanceArea calc;
217   calc.setEllipsoid( QStringLiteral( "NONE" ) );
218   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3148" ) ), QgsProject::instance()->transformContext() );
219   const QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((252000 1389000, 265000 1389000, 265000 1385000, 252000 1385000, 252000 1389000))" ) ).release() );
220   QGSCOMPARENEAR( calc.measureArea( geom ), 52000000, 0.0001 );
221 }
222 
collections()223 void TestQgsDistanceArea::collections()
224 {
225   //test measuring for collections
226   QgsDistanceArea myDa;
227   myDa.setSourceCrs( QgsCoordinateReferenceSystem::fromOgcWmsCrs( QStringLiteral( "EPSG:4030" ) ), QgsProject::instance()->transformContext() );
228   myDa.setEllipsoid( QStringLiteral( "WGS84" ) );
229 
230   //collection of lines, should be sum of line length
231   const QgsGeometry lines( QgsGeometryFactory::geomFromWkt( QStringLiteral( "GeometryCollection( LineString(0 36.53, 5.76 -48.16), LineString(0 25.54, 24.20 36.70) )" ) ).release() );
232   double result = myDa.measureLength( lines );
233   QGSCOMPARENEAR( result, 12006159, 1 );
234   result = myDa.measureArea( lines );
235   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
236 
237   //collection of polygons
238 
239   const QgsGeometry poly1 = QgsGeometry::fromWkt( QStringLiteral( "Polygon((0 36.53, 5.76 -48.16, 0 25.54, 0 36.53))" ) );
240   result = myDa.measureArea( poly1 );
241   QGSCOMPARENEAR( result, 439881520607.079712, 1 );
242   result = myDa.measureLength( poly1 );
243   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
244   const QgsGeometry poly2 = QgsGeometry::fromWkt( QStringLiteral( "Polygon((10 20, 15 20, 15 10, 10 20))" ) );
245   result = myDa.measureArea( poly2 );
246   QGSCOMPARENEAR( result, 290350317025.906982, 1 );
247   result = myDa.measureLength( poly2 );
248   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
249 
250   const 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() );
251   result = myDa.measureArea( polys );
252   QGSCOMPARENEAR( result, 730231837632.98669, 1 );
253   result = myDa.measureLength( polys );
254   QGSCOMPARENEAR( result, 0, 4 * std::numeric_limits<double>::epsilon() );
255 
256   //mixed collection
257   const 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() );
258   //measure area specifically
259   result = myDa.measureArea( mixed );
260   QGSCOMPARENEAR( result, 730231837632.98669, 1 );
261   //measure length
262   result = myDa.measureLength( mixed );
263   QGSCOMPARENEAR( result, 12006159, 1 );
264 }
265 
measureUnits()266 void TestQgsDistanceArea::measureUnits()
267 {
268   //test regression #13610
269   QgsDistanceArea calc;
270   calc.setEllipsoid( QStringLiteral( "NONE" ) );
271   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:2272" ) ), QgsProject::instance()->transformContext() );
272   QgsUnitTypes::DistanceUnit units;
273   const QgsPointXY p1( 1341683.9854275715, 408256.9562717728 );
274   const QgsPointXY p2( 1349321.7807031618, 408256.9562717728 );
275 
276   double result = calc.measureLine( p1, p2 );
277   units = calc.lengthUnits();
278   //no OTF, result will be in CRS unit (feet)
279   QCOMPARE( units, QgsUnitTypes::DistanceFeet );
280   QGSCOMPARENEAR( result, 7637.7952755903825, 0.001 );
281 
282   calc.setEllipsoid( QStringLiteral( "WGS84" ) );
283   units = calc.lengthUnits();
284   result = calc.measureLine( p1, p2 );
285   //OTF, result will be in meters
286   QCOMPARE( units, QgsUnitTypes::DistanceMeters );
287   QGSCOMPARENEAR( result, 2328.0988253106957, 0.001 );
288 }
289 
measureAreaAndUnits()290 void TestQgsDistanceArea::measureAreaAndUnits()
291 {
292   QgsDistanceArea da;
293   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext() );
294   da.setEllipsoid( QStringLiteral( "NONE" ) );
295   QgsPolylineXY ring;
296   ring << QgsPointXY( 0, 0 )
297        << QgsPointXY( 1, 0 )
298        << QgsPointXY( 1, 1 )
299        << QgsPointXY( 2, 1 )
300        << QgsPointXY( 2, 2 )
301        << QgsPointXY( 0, 2 )
302        << QgsPointXY( 0, 0 );
303   QgsPolygonXY poly;
304   poly << ring;
305 
306   QgsGeometry polygon( QgsGeometry::fromPolygonXY( poly ) );
307 
308   // We check both the measured area AND the units, in case the logic regarding
309   // ellipsoids and units changes in future
310   double area = da.measureArea( polygon );
311   QgsUnitTypes::AreaUnit units = da.areaUnits();
312 
313   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
314 
315   QVERIFY( ( qgsDoubleNear( area, 3.0, 0.00000001 ) && units == QgsUnitTypes::AreaSquareDegrees )
316            || ( qgsDoubleNear( area, 37176087091.5, 0.1 ) && units == QgsUnitTypes::AreaSquareMeters ) );
317 
318   da.setEllipsoid( QStringLiteral( "WGS84" ) );
319   area = da.measureArea( polygon );
320   units = da.areaUnits();
321   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
322   // should always be in Meters Squared
323   QGSCOMPARENEAR( area, 36922805935.961571, 0.1 );
324   QCOMPARE( units, QgsUnitTypes::AreaSquareMeters );
325 
326   // test converting the resultant area
327   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareMiles );
328   QGSCOMPARENEAR( area, 14255.975071, 0.001 );
329 
330   // now try with a source CRS which is in feet
331   ring.clear();
332   ring << QgsPointXY( 1850000, 4423000 )
333        << QgsPointXY( 1851000, 4423000 )
334        << QgsPointXY( 1851000, 4424000 )
335        << QgsPointXY( 1852000, 4424000 )
336        << QgsPointXY( 1852000, 4425000 )
337        << QgsPointXY( 1851000, 4425000 )
338        << QgsPointXY( 1850000, 4423000 );
339   poly.clear();
340   poly << ring;
341   polygon = QgsGeometry::fromPolygonXY( poly );
342 
343   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "ESRI:102635" ) ), QgsProject::instance()->transformContext() );
344   da.setEllipsoid( QStringLiteral( "NONE" ) );
345   // measurement should be in square feet
346   area = da.measureArea( polygon );
347   units = da.areaUnits();
348   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
349   QGSCOMPARENEAR( area, 2000000, 0.001 );
350   QCOMPARE( units, QgsUnitTypes::AreaSquareFeet );
351 
352   // test converting the resultant area
353   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareYards );
354   QGSCOMPARENEAR( area, 222222.2222, 0.001 );
355 
356   da.setEllipsoid( QStringLiteral( "WGS84" ) );
357   // now should be in Square Meters again
358   area = da.measureArea( polygon );
359   units = da.areaUnits();
360   QgsDebugMsg( QStringLiteral( "measured %1 in %2" ).arg( area ).arg( QgsUnitTypes::toString( units ) ) );
361   QGSCOMPARENEAR( area, 185825.206903, 1.0 );
362   QCOMPARE( units, QgsUnitTypes::AreaSquareMeters );
363 
364   // test converting the resultant area
365   area = da.convertAreaMeasurement( area, QgsUnitTypes::AreaSquareYards );
366   QgsDebugMsg( QStringLiteral( "measured %1 in sq yrds" ).arg( area ) );
367   QGSCOMPARENEAR( area, 222245.097808, 1.0 );
368 }
369 
emptyPolygon()370 void TestQgsDistanceArea::emptyPolygon()
371 {
372   QgsDistanceArea da;
373   da.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext() );
374   da.setEllipsoid( QStringLiteral( "WGS84" ) );
375 
376   //test that measuring an empty polygon doesn't crash
377   da.measurePolygon( QVector< QgsPointXY >() );
378 }
379 
regression14675()380 void TestQgsDistanceArea::regression14675()
381 {
382   //test regression #14675
383   QgsDistanceArea calc;
384   calc.setEllipsoid( QStringLiteral( "GRS80" ) );
385   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:2154" ) ), QgsProject::instance()->transformContext() );
386   const QgsGeometry geom( QgsGeometryFactory::geomFromWkt( QStringLiteral( "Polygon ((917593.5791854317067191 6833700.00807378999888897, 917596.43389983859378844 6833700.67099479306489229, 917599.53056440979707986 6833700.78673478215932846, 917593.5791854317067191 6833700.00807378999888897))" ) ).release() );
387   QGSCOMPARENEAR( calc.measureArea( geom ), 0.861747, 0.001 );
388 }
389 
regression16820()390 void TestQgsDistanceArea::regression16820()
391 {
392   QgsDistanceArea calc;
393   calc.setEllipsoid( QStringLiteral( "WGS84" ) );
394   calc.setSourceCrs( QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:32634" ) ), QgsProject::instance()->transformContext() );
395   const 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() );
396   QGSCOMPARENEAR( calc.measureArea( geom ), 43.201092, 0.001 );
397 }
398 
399 QGSTEST_MAIN( TestQgsDistanceArea )
400 #include "testqgsdistancearea.moc"
401 
402 
403 
404 
405