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