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