QGIS API Documentation  2.4.0-Chugiak
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
qgsdistancearea.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsdistancearea.cpp - Distance and area calculations on the ellipsoid
3  ---------------------------------------------------------------------------
4  Date : September 2005
5  Copyright : (C) 2005 by Martin Dobias
6  email : won.der at centrum.sk
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include <cmath>
17 #include <sqlite3.h>
18 #include <QDir>
19 #include <QString>
20 #include <QLocale>
21 #include <QObject>
22 
23 #include "qgis.h"
24 #include "qgspoint.h"
25 #include "qgscoordinatetransform.h"
27 #include "qgsgeometry.h"
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;
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 
76 {
78  mEllipsoid = origDA.mEllipsoid;
79  mSemiMajor = origDA.mSemiMajor;
80  mSemiMinor = origDA.mSemiMinor;
82  // Some calculations and trig. Should not be TOO time consuming.
83  // Alternatively we could copy the temp vars?
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  {
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();
200  }
201  else if ( parameter2.left( 3 ) == "rf=" )
202  {
203  mInvFlattening = parameter2.mid( 3 ).toDouble();
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  // precalculate some values for area calculations
236  computeAreaInit();
237 
239  return true;
240 }
241 
243 // Inverse flattening is calculated with invf = a/(a-b)
244 // Also, b = a-(a/invf)
245 bool QgsDistanceArea::setEllipsoid( double semiMajor, double semiMinor )
246 {
247  mEllipsoid = QString( "PARAMETER:%1:%2" ).arg( semiMajor ).arg( semiMinor );
248  mSemiMajor = semiMajor;
249  mSemiMinor = semiMinor;
251 
252  computeAreaInit();
253 
254  return true;
255 }
256 
258 {
259  if ( !geometry )
260  return 0.0;
261 
262  const unsigned char* wkb = geometry->asWkb();
263  if ( !wkb )
264  return 0.0;
265 
266  QgsConstWkbPtr wkbPtr( wkb + 1 );
267 
268  QGis::WkbType wkbType;
269  wkbPtr >> wkbType;
270 
271  double res, resTotal = 0;
272  int count, i;
273 
274  // measure distance or area based on what is the type of geometry
275  bool hasZptr = false;
276 
277  switch ( wkbType )
278  {
280  hasZptr = true;
281  case QGis::WKBLineString:
282  measureLine( wkb, &res, hasZptr );
283  QgsDebugMsg( "returning " + QString::number( res ) );
284  return res;
285 
287  hasZptr = true;
289  wkbPtr >> count;
290  for ( i = 0; i < count; i++ )
291  {
292  wkbPtr = measureLine( wkbPtr, &res, hasZptr );
293  resTotal += res;
294  }
295  QgsDebugMsg( "returning " + QString::number( resTotal ) );
296  return resTotal;
297 
298  case QGis::WKBPolygon25D:
299  hasZptr = true;
300  case QGis::WKBPolygon:
301  measurePolygon( wkb, &res, 0, hasZptr );
302  QgsDebugMsg( "returning " + QString::number( res ) );
303  return res;
304 
306  hasZptr = true;
308  wkbPtr >> count;
309  for ( i = 0; i < count; i++ )
310  {
311  wkbPtr = measurePolygon( wkbPtr, &res, 0, hasZptr );
312  if ( !wkbPtr )
313  {
314  QgsDebugMsg( "measurePolygon returned 0" );
315  break;
316  }
317  resTotal += res;
318  }
319  QgsDebugMsg( "returning " + QString::number( resTotal ) );
320  return resTotal;
321 
322  default:
323  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
324  return 0;
325  }
326 }
327 
329 {
330  if ( !geometry )
331  return 0.0;
332 
333  const unsigned char* wkb = geometry->asWkb();
334  if ( !wkb )
335  return 0.0;
336 
337  QgsConstWkbPtr wkbPtr( wkb + 1 );
338  QGis::WkbType wkbType;
339  wkbPtr >> wkbType;
340 
341  double res = 0.0, resTotal = 0.0;
342  int count, i;
343 
344  // measure distance or area based on what is the type of geometry
345  bool hasZptr = false;
346 
347  switch ( wkbType )
348  {
350  case QGis::WKBLineString:
353  return 0.0;
354 
355  case QGis::WKBPolygon25D:
356  hasZptr = true;
357  case QGis::WKBPolygon:
358  measurePolygon( wkb, 0, &res, hasZptr );
359  QgsDebugMsg( "returning " + QString::number( res ) );
360  return res;
361 
363  hasZptr = true;
365  wkbPtr >> count;
366  for ( i = 0; i < count; i++ )
367  {
368  wkbPtr = measurePolygon( wkbPtr, 0, &res, hasZptr );
369  if ( !wkbPtr )
370  {
371  QgsDebugMsg( "measurePolygon returned 0" );
372  break;
373  }
374  resTotal += res;
375  }
376  QgsDebugMsg( "returning " + QString::number( resTotal ) );
377  return resTotal;
378 
379  default:
380  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
381  return 0;
382  }
383 }
384 
385 
386 const unsigned char* QgsDistanceArea::measureLine( const unsigned char* feature, double* area, bool hasZptr )
387 {
388  QgsConstWkbPtr wkbPtr( feature + 1 + sizeof( int ) );
389  int nPoints;
390  wkbPtr >> nPoints;
391 
392  QList<QgsPoint> points;
393  double x, y;
394 
395  QgsDebugMsg( "This feature WKB has " + QString::number( nPoints ) + " points" );
396  // Extract the points from the WKB format into the vector
397  for ( int i = 0; i < nPoints; ++i )
398  {
399  wkbPtr >> x >> y;
400  if ( hasZptr )
401  {
402  // totally ignore Z value
403  wkbPtr += sizeof( double );
404  }
405 
406  points.append( QgsPoint( x, y ) );
407  }
408 
409  *area = measureLine( points );
410  return wkbPtr;
411 }
412 
413 double QgsDistanceArea::measureLine( const QList<QgsPoint> &points )
414 {
415  if ( points.size() < 2 )
416  return 0;
417 
418  double total = 0;
419  QgsPoint p1, p2;
420 
421  try
422  {
423  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
424  p1 = mCoordTransform->transform( points[0] );
425  else
426  p1 = points[0];
427 
428  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
429  {
430  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
431  {
432  p2 = mCoordTransform->transform( *i );
433  total += computeDistanceBearing( p1, p2 );
434  }
435  else
436  {
437  p2 = *i;
438  total += measureLine( p1, p2 );
439  }
440 
441  p1 = p2;
442  }
443 
444  return total;
445  }
446  catch ( QgsCsException &cse )
447  {
448  Q_UNUSED( cse );
449  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
450  return 0.0;
451  }
452 
453 }
454 
455 double QgsDistanceArea::measureLine( const QgsPoint &p1, const QgsPoint &p2 )
456 {
457  double result;
458 
459  try
460  {
461  QgsPoint pp1 = p1, pp2 = p2;
462 
463  QgsDebugMsgLevel( QString( "Measuring from %1 to %2" ).arg( p1.toString( 4 ) ).arg( p2.toString( 4 ) ), 3 );
464  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
465  {
466  QgsDebugMsgLevel( QString( "Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
467  QgsDebugMsgLevel( QString( "From proj4 : %1" ).arg( mCoordTransform->sourceCrs().toProj4() ), 4 );
468  QgsDebugMsgLevel( QString( "To proj4 : %1" ).arg( mCoordTransform->destCRS().toProj4() ), 4 );
469  pp1 = mCoordTransform->transform( p1 );
470  pp2 = mCoordTransform->transform( p2 );
471  QgsDebugMsgLevel( QString( "New points are %1 and %2, calculating..." ).arg( pp1.toString( 4 ) ).arg( pp2.toString( 4 ) ), 4 );
472  result = computeDistanceBearing( pp1, pp2 );
473  }
474  else
475  {
476  QgsDebugMsgLevel( "Cartesian calculation on canvas coordinates", 4 );
477  result = computeDistanceFlat( p1, p2 );
478  }
479  }
480  catch ( QgsCsException &cse )
481  {
482  Q_UNUSED( cse );
483  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
484  result = 0.0;
485  }
486  QgsDebugMsgLevel( QString( "The result was %1" ).arg( result ), 3 );
487  return result;
488 }
489 
490 
491 const unsigned char *QgsDistanceArea::measurePolygon( const unsigned char* feature, double* area, double* perimeter, bool hasZptr )
492 {
493  if ( !feature )
494  {
495  QgsDebugMsg( "no feature to measure" );
496  return 0;
497  }
498 
499  QgsConstWkbPtr wkbPtr( feature + 1 + sizeof( int ) );
500 
501  // get number of rings in the polygon
502  int numRings;
503  wkbPtr >> numRings;
504 
505  if ( numRings == 0 )
506  {
507  QgsDebugMsg( "no rings to measure" );
508  return 0;
509  }
510 
511  // Set pointer to the first ring
512  QList<QgsPoint> points;
513  QgsPoint pnt;
514  double x, y;
515  if ( area )
516  *area = 0;
517  if ( perimeter )
518  *perimeter = 0;
519 
520  try
521  {
522  for ( int idx = 0; idx < numRings; idx++ )
523  {
524  int nPoints;
525  wkbPtr >> nPoints;
526 
527  // Extract the points from the WKB and store in a pair of
528  // vectors.
529  for ( int jdx = 0; jdx < nPoints; jdx++ )
530  {
531  wkbPtr >> x >> y;
532  if ( hasZptr )
533  {
534  // totally ignore Z value
535  wkbPtr += sizeof( double );
536  }
537 
538  pnt = QgsPoint( x, y );
539 
540  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
541  {
542  pnt = mCoordTransform->transform( pnt );
543  }
544  points.append( pnt );
545  }
546 
547  if ( points.size() > 2 )
548  {
549  if ( area )
550  {
551  double areaTmp = computePolygonArea( points );
552  if ( idx == 0 )
553  {
554  // exterior ring
555  *area += areaTmp;
556  }
557  else
558  {
559  *area -= areaTmp; // interior rings
560  }
561  }
562 
563  if ( perimeter )
564  {
565  if ( idx == 0 )
566  {
567  // exterior ring
568  *perimeter += computeDistance( points );
569  }
570  }
571  }
572 
573  points.clear();
574  }
575  }
576  catch ( QgsCsException &cse )
577  {
578  Q_UNUSED( cse );
579  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
580  }
581 
582  return wkbPtr;
583 }
584 
585 
586 double QgsDistanceArea::measurePolygon( const QList<QgsPoint>& points )
587 {
588  try
589  {
590  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
591  {
592  QList<QgsPoint> pts;
593  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
594  {
595  pts.append( mCoordTransform->transform( *i ) );
596  }
597  return computePolygonArea( pts );
598  }
599  else
600  {
601  return computePolygonArea( points );
602  }
603  }
604  catch ( QgsCsException &cse )
605  {
606  Q_UNUSED( cse );
607  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
608  return 0.0;
609  }
610 }
611 
612 
613 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 )
614 {
615  QgsPoint pp1 = p1, pp2 = p2;
616  double bearing;
617 
618  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
619  {
620  pp1 = mCoordTransform->transform( p1 );
621  pp2 = mCoordTransform->transform( p2 );
622  computeDistanceBearing( pp1, pp2, &bearing );
623  }
624  else //compute simple planar azimuth
625  {
626  double dx = p2.x() - p1.x();
627  double dy = p2.y() - p1.y();
628  bearing = atan2( dx, dy );
629  }
630 
631  return bearing;
632 }
633 
634 
636 // distance calculation
637 
639  const QgsPoint& p1, const QgsPoint& p2,
640  double* course1, double* course2 )
641 {
642  if ( p1.x() == p2.x() && p1.y() == p2.y() )
643  return 0;
644 
645  // ellipsoid
646  double a = mSemiMajor;
647  double b = mSemiMinor;
648  double f = 1 / mInvFlattening;
649 
650  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
651  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
652 
653  double L = p2_lon - p1_lon;
654  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
655  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
656  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
657  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
658  double lambda = L;
659  double lambdaP = 2 * M_PI;
660 
661  double sinLambda = 0;
662  double cosLambda = 0;
663  double sinSigma = 0;
664  double cosSigma = 0;
665  double sigma = 0;
666  double alpha = 0;
667  double cosSqAlpha = 0;
668  double cos2SigmaM = 0;
669  double C = 0;
670  double tu1 = 0;
671  double tu2 = 0;
672 
673  int iterLimit = 20;
674  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
675  {
676  sinLambda = sin( lambda );
677  cosLambda = cos( lambda );
678  tu1 = ( cosU2 * sinLambda );
679  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
680  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
681  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
682  sigma = atan2( sinSigma, cosSigma );
683  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
684  cosSqAlpha = cos( alpha ) * cos( alpha );
685  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
686  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
687  lambdaP = lambda;
688  lambda = L + ( 1 - C ) * f * sin( alpha ) *
689  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
690  }
691 
692  if ( iterLimit == 0 )
693  return -1; // formula failed to converge
694 
695  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
696  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
697  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
698  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
699  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
700  double s = b * A * ( sigma - deltaSigma );
701 
702  if ( course1 )
703  {
704  *course1 = atan2( tu1, tu2 );
705  }
706  if ( course2 )
707  {
708  // PI is added to return azimuth from P2 to P1
709  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
710  }
711 
712  return s;
713 }
714 
716 {
717  return sqrt(( p2.x() - p1.x() ) * ( p2.x() - p1.x() ) + ( p2.y() - p1.y() ) * ( p2.y() - p1.y() ) );
718 }
719 
720 double QgsDistanceArea::computeDistance( const QList<QgsPoint>& points )
721 {
722  if ( points.size() < 2 )
723  return 0;
724 
725  double total = 0;
726  QgsPoint p1, p2;
727 
728  try
729  {
730  p1 = points[0];
731 
732  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
733  {
734  p2 = *i;
735  if ( mEllipsoidalMode && ( mEllipsoid != GEO_NONE ) )
736  {
737  total += computeDistanceBearing( p1, p2 );
738  }
739  else
740  {
741  total += computeDistanceFlat( p1, p2 );
742  }
743 
744  p1 = p2;
745  }
746 
747  return total;
748  }
749  catch ( QgsCsException &cse )
750  {
751  Q_UNUSED( cse );
752  QgsMessageLog::logMessage( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
753  return 0.0;
754  }
755 }
756 
757 
758 
760 // stuff for measuring areas - copied from GRASS
761 // don't know how does it work, but it's working .)
762 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
763 
764 double QgsDistanceArea::getQ( double x )
765 {
766  double sinx, sinx2;
767 
768  sinx = sin( x );
769  sinx2 = sinx * sinx;
770 
771  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
772 }
773 
774 
775 double QgsDistanceArea::getQbar( double x )
776 {
777  double cosx, cosx2;
778 
779  cosx = cos( x );
780  cosx2 = cosx * cosx;
781 
782  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
783 }
784 
785 
787 {
788  double a2 = ( mSemiMajor * mSemiMajor );
789  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
790  double e4, e6;
791 
792  m_TwoPI = M_PI + M_PI;
793 
794  e4 = e2 * e2;
795  e6 = e4 * e2;
796 
797  m_AE = a2 * ( 1 - e2 );
798 
799  m_QA = ( 2.0 / 3.0 ) * e2;
800  m_QB = ( 3.0 / 5.0 ) * e4;
801  m_QC = ( 4.0 / 7.0 ) * e6;
802 
803  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
804  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
805  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
806  m_QbarD = ( 4.0 / 49.0 ) * e6;
807 
808  m_Qp = getQ( M_PI / 2 );
809  m_E = 4 * M_PI * m_Qp * m_AE;
810  if ( m_E < 0.0 )
811  m_E = -m_E;
812 }
813 
814 
815 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint>& points )
816 {
817  double x1, y1, x2, y2, dx, dy;
818  double Qbar1, Qbar2;
819  double area;
820 
821  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
822  if (( ! mEllipsoidalMode ) || ( mEllipsoid == GEO_NONE ) )
823  {
824  return computePolygonFlatArea( points );
825  }
826  int n = points.size();
827  x2 = DEG2RAD( points[n-1].x() );
828  y2 = DEG2RAD( points[n-1].y() );
829  Qbar2 = getQbar( y2 );
830 
831  area = 0.0;
832 
833  for ( int i = 0; i < n; i++ )
834  {
835  x1 = x2;
836  y1 = y2;
837  Qbar1 = Qbar2;
838 
839  x2 = DEG2RAD( points[i].x() );
840  y2 = DEG2RAD( points[i].y() );
841  Qbar2 = getQbar( y2 );
842 
843  if ( x1 > x2 )
844  while ( x1 - x2 > M_PI )
845  x2 += m_TwoPI;
846  else if ( x2 > x1 )
847  while ( x2 - x1 > M_PI )
848  x1 += m_TwoPI;
849 
850  dx = x2 - x1;
851  area += dx * ( m_Qp - getQ( y2 ) );
852 
853  if (( dy = y2 - y1 ) != 0.0 )
854  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
855  }
856  if (( area *= m_AE ) < 0.0 )
857  area = -area;
858 
859  /* kludge - if polygon circles the south pole the area will be
860  * computed as if it cirlced the north pole. The correction is
861  * the difference between total surface area of the earth and
862  * the "north pole" area.
863  */
864  if ( area > m_E )
865  area = m_E;
866  if ( area > m_E / 2 )
867  area = m_E - area;
868 
869  return area;
870 }
871 
872 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint>& points )
873 {
874  // Normal plane area calculations.
875  double area = 0.0;
876  int i, size;
877 
878  size = points.size();
879 
880  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
881  for ( i = 0; i < size; i++ )
882  {
883  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
884  // Using '% size', so that we always end with the starting point
885  // and thus close the polygon.
886  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
887  }
888  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
889  area = area / 2.0;
890  return qAbs( area ); // All areas are positive!
891 }
892 
893 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
894 {
895  QString unitLabel;
896 
897  switch ( u )
898  {
899  case QGis::Meters:
900  if ( isArea )
901  {
902  if ( keepBaseUnit )
903  {
904  unitLabel = QObject::trUtf8( " m²" );
905  }
906  else if ( qAbs( value ) > 1000000.0 )
907  {
908  unitLabel = QObject::trUtf8( " km²" );
909  value = value / 1000000.0;
910  }
911  else if ( qAbs( value ) > 10000.0 )
912  {
913  unitLabel = QObject::tr( " ha" );
914  value = value / 10000.0;
915  }
916  else
917  {
918  unitLabel = QObject::trUtf8( " m²" );
919  }
920  }
921  else
922  {
923  if ( keepBaseUnit || qAbs( value ) == 0.0 )
924  {
925  unitLabel = QObject::tr( " m" );
926  }
927  else if ( qAbs( value ) > 1000.0 )
928  {
929  unitLabel = QObject::tr( " km" );
930  value = value / 1000;
931  }
932  else if ( qAbs( value ) < 0.01 )
933  {
934  unitLabel = QObject::tr( " mm" );
935  value = value * 1000;
936  }
937  else if ( qAbs( value ) < 0.1 )
938  {
939  unitLabel = QObject::tr( " cm" );
940  value = value * 100;
941  }
942  else
943  {
944  unitLabel = QObject::tr( " m" );
945  }
946  }
947  break;
948  case QGis::Feet:
949  if ( isArea )
950  {
951  if ( keepBaseUnit || qAbs( value ) <= 0.5*43560.0 )
952  {
953  // < 0.5 acre show sq ft
954  unitLabel = QObject::tr( " sq ft" );
955  }
956  else if ( qAbs( value ) <= 0.5*5280.0*5280.0 )
957  {
958  // < 0.5 sq mile show acre
959  unitLabel = QObject::tr( " acres" );
960  value /= 43560.0;
961  }
962  else
963  {
964  // above 0.5 acre show sq mi
965  unitLabel = QObject::tr( " sq mile" );
966  value /= 5280.0 * 5280.0;
967  }
968  }
969  else
970  {
971  if ( qAbs( value ) <= 528.0 || keepBaseUnit )
972  {
973  if ( qAbs( value ) == 1.0 )
974  {
975  unitLabel = QObject::tr( " foot" );
976  }
977  else
978  {
979  unitLabel = QObject::tr( " feet" );
980  }
981  }
982  else
983  {
984  unitLabel = QObject::tr( " mile" );
985  value /= 5280.0;
986  }
987  }
988  break;
989  case QGis::NauticalMiles:
990  if ( isArea )
991  {
992  unitLabel = QObject::tr( " sq. NM" );
993  }
994  else
995  {
996  unitLabel = QObject::tr( " NM" );
997  }
998  break;
999  case QGis::Degrees:
1000  if ( isArea )
1001  {
1002  unitLabel = QObject::tr( " sq.deg." );
1003  }
1004  else
1005  {
1006  if ( qAbs( value ) == 1.0 )
1007  unitLabel = QObject::tr( " degree" );
1008  else
1009  unitLabel = QObject::tr( " degrees" );
1010  }
1011  break;
1012  case QGis::UnknownUnit:
1013  unitLabel = QObject::tr( " unknown" );
1014  default:
1015  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
1016  };
1017 
1018 
1019  return QLocale::system().toString( value, 'f', decimals ) + unitLabel;
1020 }
1021 
1022 void QgsDistanceArea::convertMeasurement( double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea )
1023 {
1024  // Helper for converting between meters and feet and degrees and NauticalMiles...
1025  // The parameters measure and measureUnits are in/out
1026 
1027  if (( measureUnits == QGis::Degrees || measureUnits == QGis::Feet || measureUnits == QGis::NauticalMiles ) &&
1028  mEllipsoid != GEO_NONE &&
1030  {
1031  // Measuring on an ellipsoid returned meters. Force!
1032  measureUnits = QGis::Meters;
1033  QgsDebugMsg( "We're measuring on an ellipsoid or using projections, the system is returning meters" );
1034  }
1035  else if ( mEllipsoidalMode && mEllipsoid == GEO_NONE )
1036  {
1037  // Measuring in plane within the source CRS. Force its map units
1038  measureUnits = mCoordTransform->sourceCrs().mapUnits();
1039  QgsDebugMsg( "We're measuing on planimetric distance/area on given CRS, measured value is in CRS units" );
1040  }
1041 
1042  // Gets the conversion factor between the specified units
1043  double factorUnits = QGis::fromUnitToUnitFactor( measureUnits, displayUnits );
1044  if ( isArea )
1045  factorUnits *= factorUnits;
1046 
1047  QgsDebugMsg( QString( "Converting %1 %2" ).arg( QString::number( measure ), QGis::toLiteral( measureUnits ) ) );
1048  measure *= factorUnits;
1049  QgsDebugMsg( QString( "to %1 %2" ).arg( QString::number( measure ), QGis::toLiteral( displayUnits ) ) );
1050  measureUnits = displayUnits;
1051 }
1052 
const QgsCoordinateReferenceSystem & sourceCrs() const
double computePolygonFlatArea(const QList< QgsPoint > &points)
double getQbar(double x)
double computeDistance(const QList< QgsPoint > &points)
calculate distance with given coordinates (does not do a transform anymore)
~QgsDistanceArea()
Destructor.
#define QgsDebugMsg(str)
Definition: qgslogger.h:36
bool saveAsUserCRS(QString name)
Copied from QgsCustomProjectionDialog /// Please refactor into SQL handler !!! ///.
QgsCoordinateTransform * mCoordTransform
used for transforming coordinates from source CRS to ellipsoid's coordinates
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)
WkbType
Used for symbology operations.
Definition: qgis.h:53
bool setEllipsoid(const QString &ellipsoid)
sets ellipsoid by its acronym
double x() const
Definition: qgspoint.h:110
static void logMessage(QString message, QString tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
void _copy(const QgsDistanceArea &origDA)
Copy helper.
bool createFromOgcWmsCrs(QString theCrs)
Set up this CRS from the given OGC CRS.
double measurePerimeter(QgsGeometry *geometry)
measures perimeter of polygon
double measurePolygon(const QList< QgsPoint > &points)
measures polygon area
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:37
double measure(QgsGeometry *geometry)
general measurement (line distance or polygon area)
double computePolygonArea(const QList< QgsPoint > &points)
calculates area of polygon on ellipsoid algorithm has been taken from GRASS: gis/area_poly1.c
double bearing(const QgsPoint &p1, const QgsPoint &p2)
compute bearing - in radians
#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:384
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:121
QgsPoint transform(const QgsPoint p, TransformDirection direction=ForwardTransform) const
A class to represent a point geometry.
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)
double getQ(double x)
double mSemiMajor
ellipsoid parameters
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.
double measureLine(const QList< QgsPoint > &points)
measures line
const CORE_EXPORT QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
Definition: qgis.cpp:73
static double fromUnitToUnitFactor(QGis::UnitType fromUnit, QGis::UnitType toUnit)
Returns the conversion factor between the specified units.
Definition: qgis.cpp:123
Class for doing transforms between two map coordinate systems.
UnitType
Map units that qgis supports.
Definition: qgis.h:229
double y() const
Definition: qgspoint.h:118
void convertMeasurement(double &measure, QGis::UnitType &measureUnits, QGis::UnitType displayUnits, bool isArea)
Helper for conversion between physical units.
Custom exception class for Coordinate Reference System related exceptions.
bool mEllipsoidalMode
indicates whether we will transform coordinates
static QString toLiteral(QGis::UnitType unit)
Provides the canonical name of the type value.
Definition: qgis.cpp:113
double computeDistanceFlat(const QgsPoint &p1, const QgsPoint &p2)
uses flat / planimetric / Euclidean distance
const QgsCoordinateReferenceSystem & destCRS() const
double size
Definition: qgssvgcache.cpp:77
QString mEllipsoid
ellipsoid acronym (from table tbl_ellipsoids)
const unsigned char * asWkb() const
Returns the buffer containing this geometry in WKB format.
QgsDistanceArea()
Constructor.
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.
double computeDistanceBearing(const QgsPoint &p1, const QgsPoint &p2, double *course1=NULL, double *course2=NULL)
calculates distance from two points on ellipsoid based on inverse Vincenty's formulae ...
void setSourceAuthId(QString authid)
sets source spatial reference system by authid
#define tr(sourceText)