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