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