QGIS API Documentation  2.12.0-Lyon
qgsdistancearea.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdistancearea.cpp - Distance and area calculations on the ellipsoid
3  ---------------------------------------------------------------------------
4  Date : September 2005
5  Copyright : (C) 2005 by Martin Dobias
6  email : won.der at centrum.sk
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 
16 #include <cmath>
17 #include <sqlite3.h>
18 #include <QDir>
19 #include <QString>
20 #include <QLocale>
21 #include <QObject>
22 
23 #include "qgis.h"
24 #include "qgspoint.h"
25 #include "qgscoordinatetransform.h"
27 #include "qgsgeometry.h"
29 #include "qgsdistancearea.h"
30 #include "qgsapplication.h"
31 #include "qgslogger.h"
32 #include "qgsmessagelog.h"
33 #include "qgsmultisurfacev2.h"
34 #include "qgswkbptr.h"
35 #include "qgslinestringv2.h"
36 #include "qgspolygonv2.h"
37 #include "qgssurfacev2.h"
38 
39 // MSVC compiler doesn't have defined M_PI in math.h
40 #ifndef M_PI
41 #define M_PI 3.14159265358979323846
42 #endif
43 
44 #define DEG2RAD(x) ((x)*M_PI/180)
45 
46 
48 {
49  // init with default settings
50  mEllipsoidalMode = false;
51  mCoordTransform = new QgsCoordinateTransform;
52  setSourceCrs( GEOCRS_ID ); // WGS 84
54 }
55 
56 
59  : mCoordTransform( 0 )
60 {
61  _copy( origDA );
62 }
63 
65 {
66  delete mCoordTransform;
67 }
68 
71 {
72  if ( this == & origDA )
73  {
74  // Do not copy unto self
75  return *this;
76  }
77  _copy( origDA );
78  return *this;
79 }
80 
82 void QgsDistanceArea::_copy( const QgsDistanceArea & origDA )
83 {
84  mEllipsoidalMode = origDA.mEllipsoidalMode;
85  mEllipsoid = origDA.mEllipsoid;
86  mSemiMajor = origDA.mSemiMajor;
87  mSemiMinor = origDA.mSemiMinor;
88  mInvFlattening = origDA.mInvFlattening;
89  // Some calculations and trig. Should not be TOO time consuming.
90  // Alternatively we could copy the temp vars?
92  delete mCoordTransform;
93  mCoordTransform = new QgsCoordinateTransform( origDA.mCoordTransform->sourceCrs(), origDA.mCoordTransform->destCRS() );
94 }
95 
97 {
98  mEllipsoidalMode = flag;
99 }
100 
102 {
104  srcCRS.createFromSrsId( srsid );
105  mCoordTransform->setSourceCrs( srcCRS );
106 }
107 
109 {
110  mCoordTransform->setSourceCrs( srcCRS );
111 }
112 
114 {
116  srcCRS.createFromOgcWmsCrs( authId );
117  mCoordTransform->setSourceCrs( srcCRS );
118 }
119 
120 bool QgsDistanceArea::setEllipsoid( const QString& ellipsoid )
121 {
122  QString radius, parameter2;
123  //
124  // SQLITE3 stuff - get parameters for selected ellipsoid
125  //
126  sqlite3 *myDatabase;
127  const char *myTail;
128  sqlite3_stmt *myPreparedStatement;
129  int myResult;
130 
131  // Shortcut if ellipsoid is none.
132  if ( ellipsoid == GEO_NONE )
133  {
134  mEllipsoid = GEO_NONE;
135  return true;
136  }
137 
138  // Check if we have a custom projection, and set from text string.
139  // Format is "PARAMETER:<semi-major axis>:<semi minor axis>
140  // Numbers must be with (optional) decimal point and no other separators (C locale)
141  // Distances in meters. Flattening is calculated.
142  if ( ellipsoid.startsWith( "PARAMETER" ) )
143  {
144  QStringList paramList = ellipsoid.split( ":" );
145  bool semiMajorOk, semiMinorOk;
146  double semiMajor = paramList[1].toDouble( & semiMajorOk );
147  double semiMinor = paramList[2].toDouble( & semiMinorOk );
148  if ( semiMajorOk && semiMinorOk )
149  {
150  return setEllipsoid( semiMajor, semiMinor );
151  }
152  else
153  {
154  return false;
155  }
156  }
157 
158  // Continue with PROJ.4 list of ellipsoids.
159 
160  //check the db is available
161  myResult = sqlite3_open_v2( QgsApplication::srsDbFilePath().toUtf8().data(), &myDatabase, SQLITE_OPEN_READONLY, NULL );
162  if ( myResult )
163  {
164  QgsMessageLog::logMessage( QObject::tr( "Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
165  // XXX This will likely never happen since on open, sqlite creates the
166  // database if it does not exist.
167  return false;
168  }
169  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
170  QString mySql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + "'";
171  myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
172  // XXX Need to free memory from the error msg if one is set
173  if ( myResult == SQLITE_OK )
174  {
175  if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
176  {
177  radius = QString(( char * )sqlite3_column_text( myPreparedStatement, 0 ) );
178  parameter2 = QString(( char * )sqlite3_column_text( myPreparedStatement, 1 ) );
179  }
180  }
181  // close the sqlite3 statement
182  sqlite3_finalize( myPreparedStatement );
183  sqlite3_close( myDatabase );
184 
185  // row for this ellipsoid wasn't found?
186  if ( radius.isEmpty() || parameter2.isEmpty() )
187  {
188  QgsDebugMsg( QString( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
189  return false;
190  }
191 
192  // get major semiaxis
193  if ( radius.left( 2 ) == "a=" )
194  mSemiMajor = radius.mid( 2 ).toDouble();
195  else
196  {
197  QgsDebugMsg( QString( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
198  return false;
199  }
200 
201  // get second parameter
202  // one of values 'b' or 'f' is in field parameter2
203  // second one must be computed using formula: invf = a/(a-b)
204  if ( parameter2.left( 2 ) == "b=" )
205  {
206  mSemiMinor = parameter2.mid( 2 ).toDouble();
207  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
208  }
209  else if ( parameter2.left( 3 ) == "rf=" )
210  {
211  mInvFlattening = parameter2.mid( 3 ).toDouble();
212  mSemiMinor = mSemiMajor - ( mSemiMajor / mInvFlattening );
213  }
214  else
215  {
216  QgsDebugMsg( QString( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
217  return false;
218  }
219 
220  QgsDebugMsg( QString( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
221 
222 
223  // get spatial ref system for ellipsoid
224  QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
226  destCRS.createFromProj4( proj4 );
227  //TODO: createFromProj4 used to save to the user database any new CRS
228  // this behavior was changed in order to separate creation and saving.
229  // Not sure if it necessary to save it here, should be checked by someone
230  // familiar with the code (should also give a more descriptive name to the generated CRS)
231  if ( destCRS.srsid() == 0 )
232  {
233  QString myName = QString( " * %1 (%2)" )
234  .arg( QObject::tr( "Generated CRS", "A CRS automatically generated from layer info get this prefix for description" ),
235  destCRS.toProj4() );
236  destCRS.saveAsUserCRS( myName );
237  }
238  //
239 
240  // set transformation from project CRS to ellipsoid coordinates
241  mCoordTransform->setDestCRS( destCRS );
242 
243  mEllipsoid = ellipsoid;
244 
245  // precalculate some values for area calculations
246  computeAreaInit();
247 
248  return true;
249 }
250 
252 // Inverse flattening is calculated with invf = a/(a-b)
253 // Also, b = a-(a/invf)
254 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
255 {
256  mEllipsoid = QString( "PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
257  mSemiMajor = semiMajor;
258  mSemiMinor = semiMinor;
259  mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
260 
261  computeAreaInit();
262 
263  return true;
264 }
265 
266 double QgsDistanceArea::measure( const QgsAbstractGeometryV2* geomV2, MeasureType type ) const
267 {
268  if ( !geomV2 )
269  {
270  return 0.0;
271  }
272 
273  int geomDimension = geomV2->dimension();
274  if ( geomDimension <= 0 )
275  {
276  return 0.0;
277  }
278 
279  MeasureType measureType = type;
280  if ( measureType == Default )
281  {
282  measureType = ( geomDimension == 1 ? Length : Area );
283  }
284 
285  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
286  {
287  //no transform required
288  if ( measureType == Length )
289  {
290  return geomV2->length();
291  }
292  else
293  {
294  return geomV2->area();
295  }
296  }
297  else
298  {
299  //multigeom is sum of measured parts
300  const QgsGeometryCollectionV2* collection = dynamic_cast<const QgsGeometryCollectionV2*>( geomV2 );
301  if ( collection )
302  {
303  double sum = 0;
304  for ( int i = 0; i < collection->numGeometries(); ++i )
305  {
306  sum += measure( collection->geometryN( i ), measureType );
307  }
308  return sum;
309  }
310 
311  if ( measureType == Length )
312  {
313  const QgsCurveV2* curve = dynamic_cast<const QgsCurveV2*>( geomV2 );
314  if ( !curve )
315  {
316  return 0.0;
317  }
318 
319  QgsLineStringV2* lineString = curve->curveToLine();
320  double length = measureLine( lineString );
321  delete lineString;
322  return length;
323  }
324  else
325  {
326  const QgsSurfaceV2* surface = dynamic_cast<const QgsSurfaceV2*>( geomV2 );
327  if ( !surface )
328  return 0.0;
329 
330  QgsPolygonV2* polygon = surface->surfaceToPolygon();
331 
332  double area = 0;
333  const QgsCurveV2* outerRing = polygon->exteriorRing();
334  area += measurePolygon( outerRing );
335 
336  for ( int i = 0; i < polygon->numInteriorRings(); ++i )
337  {
338  const QgsCurveV2* innerRing = polygon->interiorRing( i );
339  area -= measurePolygon( innerRing );
340  }
341  delete polygon;
342  return area;
343  }
344  }
345 }
346 
347 double QgsDistanceArea::measure( const QgsGeometry *geometry ) const
348 {
349  if ( !geometry )
350  return 0.0;
351 
352  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
353  return measure( geomV2 );
354 }
355 
356 double QgsDistanceArea::measureArea( const QgsGeometry* geometry ) const
357 {
358  if ( !geometry )
359  return 0.0;
360 
361  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
362  return measure( geomV2, Area );
363 }
364 
365 double QgsDistanceArea::measureLength( const QgsGeometry* geometry ) const
366 {
367  if ( !geometry )
368  return 0.0;
369 
370  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
371  return measure( geomV2, Length );
372 }
373 
374 double QgsDistanceArea::measurePerimeter( const QgsGeometry* geometry ) const
375 {
376  if ( !geometry )
377  return 0.0;
378 
379  const QgsAbstractGeometryV2* geomV2 = geometry->geometry();
380  if ( !geomV2 || geomV2->dimension() < 2 )
381  {
382  return 0.0;
383  }
384 
385  if ( !mEllipsoidalMode || mEllipsoid == GEO_NONE )
386  {
387  return geomV2->perimeter();
388  }
389 
390  //create list with (single) surfaces
392  const QgsSurfaceV2* surf = dynamic_cast<const QgsSurfaceV2*>( geomV2 );
393  if ( surf )
394  {
395  surfaces.append( surf );
396  }
397  const QgsMultiSurfaceV2* multiSurf = dynamic_cast<const QgsMultiSurfaceV2*>( geomV2 );
398  if ( multiSurf )
399  {
400  surfaces.reserve(( surf ? 1 : 0 ) + multiSurf->numGeometries() );
401  for ( int i = 0; i < multiSurf->numGeometries(); ++i )
402  {
403  surfaces.append( static_cast<const QgsSurfaceV2*>( multiSurf->geometryN( i ) ) );
404  }
405  }
406 
407  double length = 0;
409  for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
410  {
411  if ( !*surfaceIt )
412  {
413  continue;
414  }
415 
416  QgsPolygonV2* poly = ( *surfaceIt )->surfaceToPolygon();
417  const QgsCurveV2* outerRing = poly->exteriorRing();
418  if ( outerRing )
419  {
420  length += measure( outerRing );
421  }
422  int nInnerRings = poly->numInteriorRings();
423  for ( int i = 0; i < nInnerRings; ++i )
424  {
425  length += measure( poly->interiorRing( i ) );
426  }
427  delete poly;
428  }
429  return length;
430 }
431 
432 double QgsDistanceArea::measureLine( const QgsCurveV2* curve ) const
433 {
434  if ( !curve )
435  {
436  return 0.0;
437  }
438 
439  QList<QgsPointV2> linePointsV2;
440  QList<QgsPoint> linePoints;
441  curve->points( linePointsV2 );
442  QgsGeometry::convertPointList( linePointsV2, linePoints );
443  return measureLine( linePoints );
444 }
445 
446 double QgsDistanceArea::measureLine( const QList<QgsPoint> &points ) const
447 {
448  if ( points.size() < 2 )
449  return 0;
450 
451  double total = 0;
452  QgsPoint p1, p2;
453 
454  try
455  {
456  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
457  p1 = mCoordTransform->transform( points[0] );
458  else
459  p1 = points[0];
460 
461  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
462  {
463  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
464  {
465  p2 = mCoordTransform->transform( *i );
466  total += computeDistanceBearing( p1, p2 );
467  }
468  else
469  {
470  p2 = *i;
471  total += measureLine( p1, p2 );
472  }
473 
474  p1 = p2;
475  }
476 
477  return total;
478  }
479  catch ( QgsCsException &cse )
480  {
481  Q_UNUSED( cse );
482  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
483  return 0.0;
484  }
485 
486 }
487 
488 double QgsDistanceArea::measureLine( const QgsPoint &p1, const QgsPoint &p2 ) const
489 {
490  QGis::UnitType units;
491  return measureLine( p1, p2, units );
492 }
493 
494 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2, QGis::UnitType& units ) const
495 {
496  double result;
497  units = mCoordTransform->sourceCrs().mapUnits();
498 
499  try
500  {
501  QgsPoint pp1 = p1, pp2 = p2;
502 
503  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ), p2.toString( 4 ) ), 3 );
504  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
505  {
506  units = QGis::Meters;
507  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
508  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform->sourceCrs().toProj4() ), 4 );
509  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform->destCRS().toProj4() ), 4 );
510  pp1 = mCoordTransform->transform( p1 );
511  pp2 = mCoordTransform->transform( p2 );
512  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ), pp2.toString( 4 ) ), 4 );
513  result = computeDistanceBearing( pp1, pp2 );
514  }
515  else
516  {
517  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
518  result = computeDistanceFlat( p1, p2 );
519  }
520  }
521  catch ( QgsCsException &cse )
522  {
523  Q_UNUSED( cse );
524  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
525  result = 0.0;
526  }
527  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
528  return result;
529 }
530 
531 
532 const unsigned char *QgsDistanceArea::measurePolygon( const unsigned char* feature, double* area, double* perimeter, bool hasZptr ) const
533 {
534  if ( !feature )
535  {
536  QgsDebugMsg( "no feature to measure" );
537  return 0;
538  }
539 
540  QgsConstWkbPtr wkbPtr( feature + 1 + sizeof( int ) );
541 
542  // get number of rings in the polygon
543  int numRings;
544  wkbPtr >> numRings;
545 
546  if ( numRings == 0 )
547  {
548  QgsDebugMsg( "no rings to measure" );
549  return 0;
550  }
551 
552  // Set pointer to the first ring
553  QList<QgsPoint> points;
554  QgsPoint pnt;
555  double x, y;
556  if ( area )
557  *area = 0;
558  if ( perimeter )
559  *perimeter = 0;
560 
561  try
562  {
563  for ( int idx = 0; idx < numRings; idx++ )
564  {
565  int nPoints;
566  wkbPtr >> nPoints;
567 
568  // Extract the points from the WKB and store in a pair of
569  // vectors.
570  for ( int jdx = 0; jdx < nPoints; jdx++ )
571  {
572  wkbPtr >> x >> y;
573  if ( hasZptr )
574  {
575  // totally ignore Z value
576  wkbPtr += sizeof( double );
577  }
578 
579  pnt = QgsPoint( x, y );
580 
581  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
582  {
583  pnt = mCoordTransform->transform( pnt );
584  }
585  points.append( pnt );
586  }
587 
588  if ( points.size() > 2 )
589  {
590  if ( area )
591  {
592  double areaTmp = computePolygonArea( points );
593  if ( idx == 0 )
594  {
595  // exterior ring
596  *area += areaTmp;
597  }
598  else
599  {
600  *area -= areaTmp; // interior rings
601  }
602  }
603 
604  if ( perimeter )
605  {
606  if ( idx == 0 )
607  {
608  // exterior ring
609  *perimeter += computeDistance( points );
610  }
611  }
612  }
613 
614  points.clear();
615  }
616  }
617  catch ( QgsCsException &cse )
618  {
619  Q_UNUSED( cse );
620  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
621  }
622 
623  return wkbPtr;
624 }
625 
626 double QgsDistanceArea::measurePolygon( const QgsCurveV2* curve ) const
627 {
628  if ( !curve )
629  {
630  return 0.0;
631  }
632 
633  QList<QgsPointV2> linePointsV2;
634  curve->points( linePointsV2 );
635  QList<QgsPoint> linePoints;
636  QgsGeometry::convertPointList( linePointsV2, linePoints );
637  return measurePolygon( linePoints );
638 }
639 
640 
642 {
643  try
644  {
645  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
646  {
647  QList<QgsPoint> pts;
648  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
649  {
650  pts.append( mCoordTransform->transform( *i ) );
651  }
652  return computePolygonArea( pts );
653  }
654  else
655  {
656  return computePolygonArea( points );
657  }
658  }
659  catch ( QgsCsException &cse )
660  {
661  Q_UNUSED( cse );
662  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
663  return 0.0;
664  }
665 }
666 
667 
668 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 ) const
669 {
670  QgsPoint pp1 = p1, pp2 = p2;
671  double bearing;
672 
673  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
674  {
675  pp1 = mCoordTransform->transform( p1 );
676  pp2 = mCoordTransform->transform( p2 );
677  computeDistanceBearing( pp1, pp2, &bearing );
678  }
679  else //compute simple planar azimuth
680  {
681  double dx = p2.x() - p1.x();
682  double dy = p2.y() - p1.y();
683  bearing = atan2( dx, dy );
684  }
685 
686  return bearing;
687 }
688 
689 
691 // distance calculation
692 
694  const QgsPoint& p1, const QgsPoint& p2,
695  double* course1, double* course2 ) const
696 {
697  if ( p1.x() == p2.x() && p1.y() == p2.y() )
698  return 0;
699 
700  // ellipsoid
701  double a = mSemiMajor;
702  double b = mSemiMinor;
703  double f = 1 / mInvFlattening;
704 
705  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
706  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
707 
708  double L = p2_lon - p1_lon;
709  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
710  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
711  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
712  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
713  double lambda = L;
714  double lambdaP = 2 * M_PI;
715 
716  double sinLambda = 0;
717  double cosLambda = 0;
718  double sinSigma = 0;
719  double cosSigma = 0;
720  double sigma = 0;
721  double alpha = 0;
722  double cosSqAlpha = 0;
723  double cos2SigmaM = 0;
724  double C = 0;
725  double tu1 = 0;
726  double tu2 = 0;
727 
728  int iterLimit = 20;
729  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
730  {
731  sinLambda = sin( lambda );
732  cosLambda = cos( lambda );
733  tu1 = ( cosU2 * sinLambda );
734  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
735  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
736  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
737  sigma = atan2( sinSigma, cosSigma );
738  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
739  cosSqAlpha = cos( alpha ) * cos( alpha );
740  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
741  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
742  lambdaP = lambda;
743  lambda = L + ( 1 - C ) * f * sin( alpha ) *
744  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
745  }
746 
747  if ( iterLimit == 0 )
748  return -1; // formula failed to converge
749 
750  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
751  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
752  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
753  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
754  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
755  double s = b * A * ( sigma - deltaSigma );
756 
757  if ( course1 )
758  {
759  *course1 = atan2( tu1, tu2 );
760  }
761  if ( course2 )
762  {
763  // PI is added to return azimuth from P2 to P1
764  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
765  }
766 
767  return s;
768 }
769 
770 double QgsDistanceArea::computeDistanceFlat( const QgsPoint& p1, const QgsPoint& p2 ) const
771 {
772  return sqrt(( p2.x() - p1.x() ) * ( p2.x() - p1.x() ) + ( p2.y() - p1.y() ) * ( p2.y() - p1.y() ) );
773 }
774 
776 {
777  if ( points.size() < 2 )
778  return 0;
779 
780  double total = 0;
781  QgsPoint p1, p2;
782 
783  try
784  {
785  p1 = points[0];
786 
787  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
788  {
789  p2 = *i;
790  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
791  {
792  total += computeDistanceBearing( p1, p2 );
793  }
794  else
795  {
796  total += computeDistanceFlat( p1, p2 );
797  }
798 
799  p1 = p2;
800  }
801 
802  return total;
803  }
804  catch ( QgsCsException &cse )
805  {
806  Q_UNUSED( cse );
807  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
808  return 0.0;
809  }
810 }
811 
812 
813 
815 // stuff for measuring areas - copied from GRASS
816 // don't know how does it work, but it's working .)
817 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
818 
819 double QgsDistanceArea::getQ( double x ) const
820 {
821  double sinx, sinx2;
822 
823  sinx = sin( x );
824  sinx2 = sinx * sinx;
825 
826  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
827 }
828 
829 
830 double QgsDistanceArea::getQbar( double x ) const
831 {
832  double cosx, cosx2;
833 
834  cosx = cos( x );
835  cosx2 = cosx * cosx;
836 
837  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
838 }
839 
840 
842 {
843  //don't try to perform calculations if no ellipsoid
844  if ( mEllipsoid == GEO_NONE )
845  {
846  return;
847  }
848 
849  double a2 = ( mSemiMajor * mSemiMajor );
850  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
851  double e4, e6;
852 
853  m_TwoPI = M_PI + M_PI;
854 
855  e4 = e2 * e2;
856  e6 = e4 * e2;
857 
858  m_AE = a2 * ( 1 - e2 );
859 
860  m_QA = ( 2.0 / 3.0 ) * e2;
861  m_QB = ( 3.0 / 5.0 ) * e4;
862  m_QC = ( 4.0 / 7.0 ) * e6;
863 
864  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
865  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
866  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
867  m_QbarD = ( 4.0 / 49.0 ) * e6;
868 
869  m_Qp = getQ( M_PI / 2 );
870  m_E = 4 * M_PI * m_Qp * m_AE;
871  if ( m_E < 0.0 )
872  m_E = -m_E;
873 }
874 
875 
877 {
878  double x1, y1, x2, y2, dx, dy;
879  double Qbar1, Qbar2;
880  double area;
881 
882  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
883  if (( ! mEllipsoidalMode ) || ( mEllipsoid == GEO_NONE ) )
884  {
885  return computePolygonFlatArea( points );
886  }
887  int n = points.size();
888  x2 = DEG2RAD( points[n-1].x() );
889  y2 = DEG2RAD( points[n-1].y() );
890  Qbar2 = getQbar( y2 );
891 
892  area = 0.0;
893 
894  for ( int i = 0; i < n; i++ )
895  {
896  x1 = x2;
897  y1 = y2;
898  Qbar1 = Qbar2;
899 
900  x2 = DEG2RAD( points[i].x() );
901  y2 = DEG2RAD( points[i].y() );
902  Qbar2 = getQbar( y2 );
903 
904  if ( x1 > x2 )
905  while ( x1 - x2 > M_PI )
906  x2 += m_TwoPI;
907  else if ( x2 > x1 )
908  while ( x2 - x1 > M_PI )
909  x1 += m_TwoPI;
910 
911  dx = x2 - x1;
912  area += dx * ( m_Qp - getQ( y2 ) );
913 
914  if (( dy = y2 - y1 ) != 0.0 )
915  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
916  }
917  if (( area *= m_AE ) < 0.0 )
918  area = -area;
919 
920  /* kludge - if polygon circles the south pole the area will be
921  * computed as if it cirlced the north pole. The correction is
922  * the difference between total surface area of the earth and
923  * the "north pole" area.
924  */
925  if ( area > m_E )
926  area = m_E;
927  if ( area > m_E / 2 )
928  area = m_E - area;
929 
930  return area;
931 }
932 
934 {
935  // Normal plane area calculations.
936  double area = 0.0;
937  int i, size;
938 
939  size = points.size();
940 
941  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
942  for ( i = 0; i < size; i++ )
943  {
944  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
945  // Using '% size', so that we always end with the starting point
946  // and thus close the polygon.
947  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
948  }
949  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
950  area = area / 2.0;
951  return qAbs( area ); // All areas are positive!
952 }
953 
954 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
955 {
956  QString unitLabel;
957 
958  switch ( u )
959  {
960  case QGis::Meters:
961  if ( isArea )
962  {
963  if ( keepBaseUnit )
964  {
965  unitLabel = QObject::trUtf8( " m²" );
966  }
967  else if ( qAbs( value ) > 1000000.0 )
968  {
969  unitLabel = QObject::trUtf8( " km²" );
970  value = value / 1000000.0;
971  }
972  else if ( qAbs( value ) > 10000.0 )
973  {
974  unitLabel = QObject::tr( " ha" );
975  value = value / 10000.0;
976  }
977  else
978  {
979  unitLabel = QObject::trUtf8( " m²" );
980  }
981  }
982  else
983  {
984  if ( keepBaseUnit || qAbs( value ) == 0.0 )
985  {
986  unitLabel = QObject::tr( " m" );
987  }
988  else if ( qAbs( value ) > 1000.0 )
989  {
990  unitLabel = QObject::tr( " km" );
991  value = value / 1000;
992  }
993  else if ( qAbs( value ) < 0.01 )
994  {
995  unitLabel = QObject::tr( " mm" );
996  value = value * 1000;
997  }
998  else if ( qAbs( value ) < 0.1 )
999  {
1000  unitLabel = QObject::tr( " cm" );
1001  value = value * 100;
1002  }
1003  else
1004  {
1005  unitLabel = QObject::tr( " m" );
1006  }
1007  }
1008  break;
1009  case QGis::Feet:
1010  if ( isArea )
1011  {
1012  if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
1013  {
1014  // < 0.5 acre show sq ft
1015  unitLabel = QObject::tr( " sq ft" );
1016  }
1017  else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
1018  {
1019  // < 0.5 sq mile show acre
1020  unitLabel = QObject::tr( " acres" );
1021  value /= 43560.0;
1022  }
1023  else
1024  {
1025  // above 0.5 acre show sq mi
1026  unitLabel = QObject::tr( " sq mile" );
1027  value /= 5280.0 * 5280.0;
1028  }
1029  }
1030  else
1031  {
1032  if ( qAbs( value ) <= 528.0 || keepBaseUnit )
1033  {
1034  if ( qAbs( value ) == 1.0 )
1035  {
1036  unitLabel = QObject::tr( " foot" );
1037  }
1038  else
1039  {
1040  unitLabel = QObject::tr( " feet" );
1041  }
1042  }
1043  else
1044  {
1045  unitLabel = QObject::tr( " mile" );
1046  value /= 5280.0;
1047  }
1048  }
1049  break;
1050  case QGis::NauticalMiles:
1051  if ( isArea )
1052  {
1053  unitLabel = QObject::tr( " sq. NM" );
1054  }
1055  else
1056  {
1057  unitLabel = QObject::tr( " NM" );
1058  }
1059  break;
1060  case QGis::Degrees:
1061  if ( isArea )
1062  {
1063  unitLabel = QObject::tr( " sq.deg." );
1064  }
1065  else
1066  {
1067  if ( qAbs( value ) == 1.0 )
1068  unitLabel = QObject::tr( " degree" );
1069  else
1070  unitLabel = QObject::tr( " degrees" );
1071  }
1072  break;
1073  case QGis::UnknownUnit:
1074  unitLabel = QObject::tr( " unknown" );
1075  //intentional fall-through
1076  default:
1077  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
1078  }
1079 
1080  return QLocale::system().toString( value, 'f', decimals ) + unitLabel;
1081 }
1082 
1083 void QgsDistanceArea::convertMeasurement( double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea ) const
1084 {
1085  // Helper for converting between meters and feet and degrees and NauticalMiles...
1086  // The parameters measure and measureUnits are in/out
1087 
1088  if (( measureUnits == QGis::Degrees || measureUnits == QGis::Feet || measureUnits == QGis::NauticalMiles ) &&
1089  mEllipsoid != GEO_NONE &&
1090  mEllipsoidalMode )
1091  {
1092  // Measuring on an ellipsoid returned meters. Force!
1093  measureUnits = QGis::Meters;
1094  QgsDebugMsg( "We're measuring on an ellipsoid or using projections, the system is returning meters" );
1095  }
1096  else if ( mEllipsoidalMode && mEllipsoid == GEO_NONE )
1097  {
1098  // Measuring in plane within the source CRS. Force its map units
1099  measureUnits = mCoordTransform->sourceCrs().mapUnits();
1100  QgsDebugMsg( "We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1101  }
1102 
1103  // Gets the conversion factor between the specified units
1104  double factorUnits = QGis::fromUnitToUnitFactor( measureUnits, displayUnits );
1105  if ( isArea )
1106  factorUnits *= factorUnits;
1107 
1108  QgsDebugMsg( QString( "Converting %1 %2" ).arg( QString::number( measure ), QGis::toLiteral( measureUnits ) ) );
1109  measure *= factorUnits;
1110  QgsDebugMsg( QString( "to %1 %2" ).arg( QString::number( measure ), QGis::toLiteral( displayUnits ) ) );
1111  measureUnits = displayUnits;
1112 }
1113 
double bearing(const QgsPoint &p1, const QgsPoint &p2) const
compute bearing - in radians
const QgsCoordinateReferenceSystem & sourceCrs() const
void clear()
virtual QgsPolygonV2 * surfaceToPolygon() const =0
void convertMeasurement(double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea) const
Helper for conversion between physical units.
QString toString(qlonglong i) const
~QgsDistanceArea()
Destructor.
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
long srsid() const
Get the SrsId - if possible.
UnitType
Map units that qgis supports.
Definition: qgis.h:147
void setSourceCrs(const QgsCoordinateReferenceSystem &theCRS)
QgsAbstractGeometryV2 * geometry() const
Returns the underlying geometry store.
virtual double length() const
Returns the length of the geometry.
QStringList split(const QString &sep, SplitBehavior behavior, Qt::CaseSensitivity cs) const
virtual void points(QList< QgsPointV2 > &pt) const =0
Returns a list of points within the curve.
void reserve(int alloc)
void setSourceCrs(long srsid)
sets source spatial reference system (by QGIS CRS)
void computeAreaInit()
precalculates some values (must be called always when changing ellipsoid)
#define DEG2RAD(x)
QgsCurveV2 * exteriorRing() const
Abstract base class for all geometries.
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:76
QgsPoint transform(const QgsPoint &p, TransformDirection direction=ForwardTransform) const
Transform the point from Source Coordinate System to Destination Coordinate System If the direction i...
bool setEllipsoid(const QString &ellipsoid)
sets ellipsoid by its acronym
int length() const
double toDouble(bool *ok) const
QString tr(const char *sourceText, const char *disambiguation, int n)
double measurePolygon(const QList< QgsPoint > &points) const
measures polygon area
double x() const
Get the x value of the point.
Definition: qgspoint.h:126
QString trUtf8(const char *sourceText, const char *disambiguation, int n)
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:75
int size() const
QLocale system()
Polygon geometry type.
Definition: qgspolygonv2.h:29
bool createFromOgcWmsCrs(QString theCrs)
Set up this CRS from the given OGC CRS.
Q_DECL_DEPRECATED double measure(const QgsGeometry *geometry) const
General measurement (line distance or polygon area)
QString number(int n, int base)
bool createFromSrsId(const long theSrsId)
Set up this srs by fetching the appropriate information from the sqlite backend.
double computeDistance(const QList< QgsPoint > &points) const
calculate distance with given coordinates (does not do a transform anymore)
void append(const T &value)
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
double computePolygonArea(const QList< QgsPoint > &points) const
calculates area of polygon on ellipsoid algorithm has been taken from GRASS: gis/area_poly1.c
Line string geometry type.
double computeDistanceBearing(const QgsPoint &p1, const QgsPoint &p2, double *course1=NULL, double *course2=NULL) const
calculates distance from two points on ellipsoid based on inverse Vincenty's formulae ...
static void convertPointList(const QList< QgsPoint > &input, QList< QgsPointV2 > &output)
Upgrades a point list from QgsPoint to QgsPointV2.
bool isEmpty() const
double computePolygonFlatArea(const QList< QgsPoint > &points) const
bool startsWith(const QString &s, Qt::CaseSensitivity cs) const
#define M_PI
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
Definition: qgis.h:336
static void logMessage(const QString &message, const QString &tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
double measureLine(const QList< QgsPoint > &points) const
measures line
Multi surface geometry collection.
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:126
double measurePerimeter(const QgsGeometry *geometry) const
measures perimeter of polygon
bool saveAsUserCRS(const QString &name)
Save the proj4-string as a custom CRS.
A class to represent a point.
Definition: qgspoint.h:63
iterator end()
double measureArea(const QgsGeometry *geometry) const
Measures the area of a geometry.
struct sqlite3 sqlite3
QgsDistanceArea & operator=(const QgsDistanceArea &origDA)
Assignment operator.
static QString textUnit(double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit=false)
General purpose distance and area calculator.
void setDestCRS(const QgsCoordinateReferenceSystem &theCRS)
int numGeometries() const
Returns the number of geometries within the collection.
double measureLength(const QgsGeometry *geometry) const
Measures the length of a geometry.
virtual QgsLineStringV2 * curveToLine() const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
QString mid(int position, int n) const
const QString & ellipsoid() const
returns ellipsoid's acronym
QgsPolygonV2 * surfaceToPolygon() const override
void setSourceAuthId(const QString &authid)
sets source spatial reference system by authid
Class for storing a coordinate reference system (CRS)
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual double area() const
Returns the area of the geometry.
static double fromUnitToUnitFactor(QGis::UnitType fromUnit, QGis::UnitType toUnit)
Returns the conversion factor between the specified units.
Definition: qgis.cpp:187
Class for doing transforms between two map coordinate systems.
const QgsAbstractGeometryV2 * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
QString left(int n) const
double y() const
Get the y value of the point.
Definition: qgspoint.h:134
static QString srsDbFilePath()
Returns the path to the srs.db file.
Custom exception class for Coordinate Reference System related exceptions.
const_iterator constEnd() const
static QString toLiteral(QGis::UnitType unit)
Provides the canonical name of the type value.
Definition: qgis.cpp:165
int numInteriorRings() const
virtual double perimeter() const
Returns the perimeter of the geometry.
const_iterator constBegin() const
Abstract base class for curved geometry type.
Definition: qgscurvev2.h:32
const QgsCoordinateReferenceSystem & destCRS() const
QString arg(qlonglong a, int fieldWidth, int base, const QChar &fillChar) const
QgsDistanceArea()
Constructor.
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2) const
uses flat / planimetric / Euclidean distance
iterator begin()
QgsCurveV2 * interiorRing(int i) const
void setEllipsoidalMode(bool flag)
sets whether coordinates must be projected to ellipsoid before measuring
bool createFromProj4(const QString &theProjString)
Set up this srs by passing it a proj4 style formatted string.
QString toProj4() const
Get the Proj Proj4 string representation of this srs.
QByteArray toUtf8() const
QGis::UnitType mapUnits() const
Get the units that the projection is in.