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