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