Quantum GIS API Documentation  1.7.5-Wroclaw
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 /* $Id$ */
16 
17 #include <cmath>
18 #include <sqlite3.h>
19 #include <QDir>
20 #include <QString>
21 #include <QLocale>
22 #include <QObject>
23 
24 #include "qgis.h"
25 #include "qgspoint.h"
26 #include "qgscoordinatetransform.h"
28 #include "qgsgeometry.h"
29 #include "qgsdistancearea.h"
30 #include "qgsapplication.h"
31 #include "qgslogger.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  mProjectionsEnabled = false;
46  setSourceCrs( GEOCRS_ID ); // WGS 84
47  setEllipsoid( "WGS84" );
48 }
49 
50 
52 {
53  delete mCoordTransform;
54 }
55 
56 
58 {
59  mProjectionsEnabled = flag;
60 }
61 
63 {
65  srcCRS.createFromSrsId( srsid );
66  mCoordTransform->setSourceCrs( srcCRS );
67 }
68 
70 {
72  srcCRS.createFromOgcWmsCrs( QString( "EPSG:%1" ).arg( epsgId ) );
73  mCoordTransform->setSourceCrs( srcCRS );
74 }
75 
76 void QgsDistanceArea::setSourceAuthId( QString authId )
77 {
79  srcCRS.createFromOgcWmsCrs( authId );
80  mCoordTransform->setSourceCrs( srcCRS );
81 }
82 
83 bool QgsDistanceArea::setEllipsoid( const QString& ellipsoid )
84 {
85  QString radius, parameter2;
86  //
87  // SQLITE3 stuff - get parameters for selected ellipsoid
88  //
89  sqlite3 *myDatabase;
90  const char *myTail;
91  sqlite3_stmt *myPreparedStatement;
92  int myResult;
93 
94  // Shortcut if ellipsoid is none.
95  if ( ellipsoid == "NONE" )
96  {
97  mEllipsoid = "NONE";
98  return true;
99  }
100 
101  //check the db is available
102  myResult = sqlite3_open( QgsApplication::srsDbFilePath().toUtf8().data(), &myDatabase );
103  if ( myResult )
104  {
105  QgsDebugMsg( QString( "Can't open database: %1" ).arg( sqlite3_errmsg( myDatabase ) ) );
106  // XXX This will likely never happen since on open, sqlite creates the
107  // database if it does not exist.
108  return false;
109  }
110  // Set up the query to retrieve the projection information needed to populate the ELLIPSOID list
111  QString mySql = "select radius, parameter2 from tbl_ellipsoid where acronym='" + ellipsoid + "'";
112  myResult = sqlite3_prepare( myDatabase, mySql.toUtf8(), mySql.toUtf8().length(), &myPreparedStatement, &myTail );
113  // XXX Need to free memory from the error msg if one is set
114  if ( myResult == SQLITE_OK )
115  {
116  if ( sqlite3_step( myPreparedStatement ) == SQLITE_ROW )
117  {
118  radius = QString(( char * )sqlite3_column_text( myPreparedStatement, 0 ) );
119  parameter2 = QString(( char * )sqlite3_column_text( myPreparedStatement, 1 ) );
120  }
121  }
122  // close the sqlite3 statement
123  sqlite3_finalize( myPreparedStatement );
124  sqlite3_close( myDatabase );
125 
126  // row for this ellipsoid wasn't found?
127  if ( radius.isEmpty() || parameter2.isEmpty() )
128  {
129  QgsDebugMsg( QString( "setEllipsoid: no row in tbl_ellipsoid for acronym '%1'" ).arg( ellipsoid ) );
130  return false;
131  }
132 
133  // get major semiaxis
134  if ( radius.left( 2 ) == "a=" )
135  mSemiMajor = radius.mid( 2 ).toDouble();
136  else
137  {
138  QgsDebugMsg( QString( "setEllipsoid: wrong format of radius field: '%1'" ).arg( radius ) );
139  return false;
140  }
141 
142  // get second parameter
143  // one of values 'b' or 'f' is in field parameter2
144  // second one must be computed using formula: invf = a/(a-b)
145  if ( parameter2.left( 2 ) == "b=" )
146  {
147  mSemiMinor = parameter2.mid( 2 ).toDouble();
149  }
150  else if ( parameter2.left( 3 ) == "rf=" )
151  {
152  mInvFlattening = parameter2.mid( 3 ).toDouble();
154  }
155  else
156  {
157  QgsDebugMsg( QString( "setEllipsoid: wrong format of parameter2 field: '%1'" ).arg( parameter2 ) );
158  return false;
159  }
160 
161  QgsDebugMsg( QString( "setEllipsoid: a=%1, b=%2, 1/f=%3" ).arg( mSemiMajor ).arg( mSemiMinor ).arg( mInvFlattening ) );
162 
163 
164  // get spatial ref system for ellipsoid
165  QString proj4 = "+proj=longlat +ellps=" + ellipsoid + " +no_defs";
167  destCRS.createFromProj4( proj4 );
168 
169  // set transformation from project CRS to ellipsoid coordinates
170  mCoordTransform->setDestCRS( destCRS );
171 
172  // precalculate some values for area calculations
173  computeAreaInit();
174 
176  return true;
177 }
178 
179 
181 {
182  if ( !geometry )
183  return 0.0;
184 
185  unsigned char* wkb = geometry->asWkb();
186  if ( !wkb )
187  return 0.0;
188 
189  unsigned char* ptr;
190  unsigned int wkbType;
191  double res, resTotal = 0;
192  int count, i;
193 
194  memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
195 
196  // measure distance or area based on what is the type of geometry
197  bool hasZptr = false;
198 
199  switch ( wkbType )
200  {
202  hasZptr = true;
203  case QGis::WKBLineString:
204  measureLine( wkb, &res, hasZptr );
205  QgsDebugMsg( "returning " + QString::number( res ) );
206  return res;
207 
209  hasZptr = true;
211  count = *(( int* )( wkb + 5 ) );
212  ptr = wkb + 9;
213  for ( i = 0; i < count; i++ )
214  {
215  ptr = measureLine( ptr, &res, hasZptr );
216  resTotal += res;
217  }
218  QgsDebugMsg( "returning " + QString::number( resTotal ) );
219  return resTotal;
220 
221  case QGis::WKBPolygon25D:
222  hasZptr = true;
223  case QGis::WKBPolygon:
224  measurePolygon( wkb, &res, 0, hasZptr );
225  QgsDebugMsg( "returning " + QString::number( res ) );
226  return res;
227 
229  hasZptr = true;
231  count = *(( int* )( wkb + 5 ) );
232  ptr = wkb + 9;
233  for ( i = 0; i < count; i++ )
234  {
235  ptr = measurePolygon( ptr, &res, 0, hasZptr );
236  resTotal += res;
237  }
238  QgsDebugMsg( "returning " + QString::number( resTotal ) );
239  return resTotal;
240 
241  default:
242  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
243  return 0;
244  }
245 }
246 
248 {
249  if ( !geometry )
250  return 0.0;
251 
252  unsigned char* wkb = geometry->asWkb();
253  if ( !wkb )
254  return 0.0;
255 
256  unsigned char* ptr;
257  unsigned int wkbType;
258  double res, resTotal = 0;
259  int count, i;
260 
261  memcpy( &wkbType, ( wkb + 1 ), sizeof( wkbType ) );
262 
263  // measure distance or area based on what is the type of geometry
264  bool hasZptr = false;
265 
266  switch ( wkbType )
267  {
269  case QGis::WKBLineString:
272  return 0.0;
273 
274  case QGis::WKBPolygon25D:
275  hasZptr = true;
276  case QGis::WKBPolygon:
277  measurePolygon( wkb, 0, &res, hasZptr );
278  QgsDebugMsg( "returning " + QString::number( res ) );
279  return res;
280 
282  hasZptr = true;
284  count = *(( int* )( wkb + 5 ) );
285  ptr = wkb + 9;
286  for ( i = 0; i < count; i++ )
287  {
288  ptr = measurePolygon( ptr, 0, &res, hasZptr );
289  resTotal += res;
290  }
291  QgsDebugMsg( "returning " + QString::number( resTotal ) );
292  return resTotal;
293 
294  default:
295  QgsDebugMsg( QString( "measure: unexpected geometry type: %1" ).arg( wkbType ) );
296  return 0;
297  }
298 }
299 
300 
301 unsigned char* QgsDistanceArea::measureLine( unsigned char* feature, double* area, bool hasZptr )
302 {
303  unsigned char *ptr = feature + 5;
304  unsigned int nPoints = *(( int* )ptr );
305  ptr = feature + 9;
306 
307  QList<QgsPoint> points;
308  double x, y;
309 
310  QgsDebugMsg( "This feature WKB has " + QString::number( nPoints ) + " points" );
311  // Extract the points from the WKB format into the vector
312  for ( unsigned int i = 0; i < nPoints; ++i )
313  {
314  x = *(( double * ) ptr );
315  ptr += sizeof( double );
316  y = *(( double * ) ptr );
317  ptr += sizeof( double );
318  if ( hasZptr )
319  {
320  // totally ignore Z value
321  ptr += sizeof( double );
322  }
323 
324  points.append( QgsPoint( x, y ) );
325  }
326 
327  *area = measureLine( points );
328  return ptr;
329 }
330 
331 double QgsDistanceArea::measureLine( const QList<QgsPoint>& points )
332 {
333  if ( points.size() < 2 )
334  return 0;
335 
336  double total = 0;
337  QgsPoint p1, p2;
338 
339  try
340  {
341  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
342  p1 = mCoordTransform->transform( points[0] );
343  else
344  p1 = points[0];
345 
346  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
347  {
348  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
349  {
350  p2 = mCoordTransform->transform( *i );
351  total += computeDistanceBearing( p1, p2 );
352  }
353  else
354  {
355  p2 = *i;
356  total += measureLine( p1, p2 );
357  }
358 
359  p1 = p2;
360  }
361 
362  return total;
363  }
364  catch ( QgsCsException &cse )
365  {
366  Q_UNUSED( cse );
367  QgsLogger::warning( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
368  return 0.0;
369  }
370 
371 }
372 
373 double QgsDistanceArea::measureLine( const QgsPoint& p1, const QgsPoint& p2 )
374 {
375  try
376  {
377  QgsPoint pp1 = p1, pp2 = p2;
378  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
379  {
380  pp1 = mCoordTransform->transform( p1 );
381  pp2 = mCoordTransform->transform( p2 );
382  return computeDistanceBearing( pp1, pp2 );
383  }
384  else
385  {
386  return sqrt(( p2.x() - p1.x() )*( p2.x() - p1.x() ) + ( p2.y() - p1.y() )*( p2.y() - p1.y() ) );
387  }
388  }
389  catch ( QgsCsException &cse )
390  {
391  Q_UNUSED( cse );
392  QgsLogger::warning( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
393  return 0.0;
394  }
395 }
396 
397 
398 unsigned char* QgsDistanceArea::measurePolygon( unsigned char* feature, double* area, double* perimeter, bool hasZptr )
399 {
400  // get number of rings in the polygon
401  unsigned int numRings = *(( int* )( feature + 1 + sizeof( int ) ) );
402 
403  if ( numRings == 0 )
404  return 0;
405 
406  // Set pointer to the first ring
407  unsigned char* ptr = feature + 1 + 2 * sizeof( int );
408 
409  QList<QgsPoint> points;
410  QgsPoint pnt;
411  double x, y;
412  if ( area )
413  *area = 0;
414  if ( perimeter )
415  *perimeter = 0;
416 
417  try
418  {
419  for ( unsigned int idx = 0; idx < numRings; idx++ )
420  {
421  int nPoints = *(( int* )ptr );
422  ptr += 4;
423 
424  // Extract the points from the WKB and store in a pair of
425  // vectors.
426  for ( int jdx = 0; jdx < nPoints; jdx++ )
427  {
428  x = *(( double * ) ptr );
429  ptr += sizeof( double );
430  y = *(( double * ) ptr );
431  ptr += sizeof( double );
432  if ( hasZptr )
433  {
434  // totally ignore Z value
435  ptr += sizeof( double );
436  }
437 
438  pnt = QgsPoint( x, y );
439 
440  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
441  {
442  pnt = mCoordTransform->transform( pnt );
443  }
444  points.append( pnt );
445  }
446 
447  if ( points.size() > 2 )
448  {
449  if ( area )
450  {
451  double areaTmp = computePolygonArea( points );
452  if ( idx == 0 )
453  {
454  // exterior ring
455  *area += areaTmp;
456  }
457  else
458  {
459  *area -= areaTmp; // interior rings
460  }
461  }
462 
463  if ( perimeter )
464  {
465  *perimeter += measureLine( points );
466  }
467  }
468 
469  points.clear();
470 
471  if ( !area )
472  {
473  break;
474  }
475  }
476  }
477  catch ( QgsCsException &cse )
478  {
479  Q_UNUSED( cse );
480  QgsLogger::warning( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area or perimeter." ) );
481  }
482 
483  return ptr;
484 }
485 
486 
487 double QgsDistanceArea::measurePolygon( const QList<QgsPoint>& points )
488 {
489 
490  try
491  {
492  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
493  {
494  QList<QgsPoint> pts;
495  for ( QList<QgsPoint>::const_iterator i = points.begin(); i != points.end(); ++i )
496  {
497  pts.append( mCoordTransform->transform( *i ) );
498  }
499  return computePolygonArea( pts );
500  }
501  else
502  {
503  return computePolygonArea( points );
504  }
505  }
506  catch ( QgsCsException &cse )
507  {
508  Q_UNUSED( cse );
509  QgsLogger::warning( QObject::tr( "Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
510  return 0.0;
511  }
512 }
513 
514 
515 double QgsDistanceArea::bearing( const QgsPoint& p1, const QgsPoint& p2 )
516 {
517  QgsPoint pp1 = p1, pp2 = p2;
518  double bearing;
519 
520  if ( mProjectionsEnabled && ( mEllipsoid != "NONE" ) )
521  {
522  pp1 = mCoordTransform->transform( p1 );
523  pp2 = mCoordTransform->transform( p2 );
524  computeDistanceBearing( pp1, pp2, &bearing );
525  }
526  else //compute simple planar azimuth
527  {
528  double dx = p2.x() - p1.x();
529  double dy = p2.y() - p1.y();
530  bearing = atan2( dx, dy );
531  }
532 
533  return bearing;
534 }
535 
536 
538 // distance calculation
539 
541  const QgsPoint& p1, const QgsPoint& p2,
542  double* course1, double* course2 )
543 {
544  if ( p1.x() == p2.x() && p1.y() == p2.y() )
545  return 0;
546 
547  // ellipsoid
548  double a = mSemiMajor;
549  double b = mSemiMinor;
550  double f = 1 / mInvFlattening;
551 
552  double p1_lat = DEG2RAD( p1.y() ), p1_lon = DEG2RAD( p1.x() );
553  double p2_lat = DEG2RAD( p2.y() ), p2_lon = DEG2RAD( p2.x() );
554 
555  double L = p2_lon - p1_lon;
556  double U1 = atan(( 1 - f ) * tan( p1_lat ) );
557  double U2 = atan(( 1 - f ) * tan( p2_lat ) );
558  double sinU1 = sin( U1 ), cosU1 = cos( U1 );
559  double sinU2 = sin( U2 ), cosU2 = cos( U2 );
560  double lambda = L;
561  double lambdaP = 2 * M_PI;
562 
563  double sinLambda = 0;
564  double cosLambda = 0;
565  double sinSigma = 0;
566  double cosSigma = 0;
567  double sigma = 0;
568  double alpha = 0;
569  double cosSqAlpha = 0;
570  double cos2SigmaM = 0;
571  double C = 0;
572  double tu1 = 0;
573  double tu2 = 0;
574 
575  int iterLimit = 20;
576  while ( qAbs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
577  {
578  sinLambda = sin( lambda );
579  cosLambda = cos( lambda );
580  tu1 = ( cosU2 * sinLambda );
581  tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
582  sinSigma = sqrt( tu1 * tu1 + tu2 * tu2 );
583  cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
584  sigma = atan2( sinSigma, cosSigma );
585  alpha = asin( cosU1 * cosU2 * sinLambda / sinSigma );
586  cosSqAlpha = cos( alpha ) * cos( alpha );
587  cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
588  C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
589  lambdaP = lambda;
590  lambda = L + ( 1 - C ) * f * sin( alpha ) *
591  ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
592  }
593 
594  if ( iterLimit == 0 )
595  return -1; // formula failed to converge
596 
597  double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
598  double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
599  double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
600  double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
601  B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
602  double s = b * A * ( sigma - deltaSigma );
603 
604  if ( course1 )
605  {
606  *course1 = atan2( tu1, tu2 );
607  }
608  if ( course2 )
609  {
610  // PI is added to return azimuth from P2 to P1
611  *course2 = atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
612  }
613 
614  return s;
615 }
616 
617 
618 
620 // stuff for measuring areas - copied from GRASS
621 // don't know how does it work, but it's working .)
622 // see G_begin_ellipsoid_polygon_area() in area_poly1.c
623 
624 double QgsDistanceArea::getQ( double x )
625 {
626  double sinx, sinx2;
627 
628  sinx = sin( x );
629  sinx2 = sinx * sinx;
630 
631  return sinx *( 1 + sinx2 *( m_QA + sinx2 *( m_QB + sinx2 * m_QC ) ) );
632 }
633 
634 
635 double QgsDistanceArea::getQbar( double x )
636 {
637  double cosx, cosx2;
638 
639  cosx = cos( x );
640  cosx2 = cosx * cosx;
641 
642  return cosx *( m_QbarA + cosx2 *( m_QbarB + cosx2 *( m_QbarC + cosx2 * m_QbarD ) ) );
643 }
644 
645 
647 {
648  double a2 = ( mSemiMajor * mSemiMajor );
649  double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
650  double e4, e6;
651 
652  m_TwoPI = M_PI + M_PI;
653 
654  e4 = e2 * e2;
655  e6 = e4 * e2;
656 
657  m_AE = a2 * ( 1 - e2 );
658 
659  m_QA = ( 2.0 / 3.0 ) * e2;
660  m_QB = ( 3.0 / 5.0 ) * e4;
661  m_QC = ( 4.0 / 7.0 ) * e6;
662 
663  m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
664  m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
665  m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
666  m_QbarD = ( 4.0 / 49.0 ) * e6;
667 
668  m_Qp = getQ( M_PI / 2 );
669  m_E = 4 * M_PI * m_Qp * m_AE;
670  if ( m_E < 0.0 ) m_E = -m_E;
671 }
672 
673 
674 double QgsDistanceArea::computePolygonArea( const QList<QgsPoint>& points )
675 {
676  double x1, y1, x2, y2, dx, dy;
677  double Qbar1, Qbar2;
678  double area;
679 
680  QgsDebugMsgLevel( "Ellipsoid: " + mEllipsoid, 3 );
681  if (( ! mProjectionsEnabled ) || ( mEllipsoid == "NONE" ) )
682  {
683  return computePolygonFlatArea( points );
684  }
685  int n = points.size();
686  x2 = DEG2RAD( points[n-1].x() );
687  y2 = DEG2RAD( points[n-1].y() );
688  Qbar2 = getQbar( y2 );
689 
690  area = 0.0;
691 
692  for ( int i = 0; i < n; i++ )
693  {
694  x1 = x2;
695  y1 = y2;
696  Qbar1 = Qbar2;
697 
698  x2 = DEG2RAD( points[i].x() );
699  y2 = DEG2RAD( points[i].y() );
700  Qbar2 = getQbar( y2 );
701 
702  if ( x1 > x2 )
703  while ( x1 - x2 > M_PI )
704  x2 += m_TwoPI;
705  else if ( x2 > x1 )
706  while ( x2 - x1 > M_PI )
707  x1 += m_TwoPI;
708 
709  dx = x2 - x1;
710  area += dx * ( m_Qp - getQ( y2 ) );
711 
712  if (( dy = y2 - y1 ) != 0.0 )
713  area += dx * getQ( y2 ) - ( dx / dy ) * ( Qbar2 - Qbar1 );
714  }
715  if (( area *= m_AE ) < 0.0 )
716  area = -area;
717 
718  /* kludge - if polygon circles the south pole the area will be
719  * computed as if it cirlced the north pole. The correction is
720  * the difference between total surface area of the earth and
721  * the "north pole" area.
722  */
723  if ( area > m_E ) area = m_E;
724  if ( area > m_E / 2 ) area = m_E - area;
725 
726  return area;
727 }
728 
729 double QgsDistanceArea::computePolygonFlatArea( const QList<QgsPoint>& points )
730 {
731  // Normal plane area calculations.
732  double area = 0.0;
733  int i, size;
734 
735  size = points.size();
736 
737  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
738  for ( i = 0; i < size; i++ )
739  {
740  // QgsDebugMsg("Area from point: " + (points[i]).toString(2));
741  // Using '% size', so that we always end with the starting point
742  // and thus close the polygon.
743  area = area + points[i].x() * points[( i+1 ) % size].y() - points[( i+1 ) % size].x() * points[i].y();
744  }
745  // QgsDebugMsg("Area from point: " + (points[i % size]).toString(2));
746  area = area / 2.0;
747  return qAbs( area ); // All areas are positive!
748 }
749 
750 QString QgsDistanceArea::textUnit( double value, int decimals, QGis::UnitType u, bool isArea, bool keepBaseUnit )
751 {
752  QString unitLabel;
753 
754  switch ( u )
755  {
756  case QGis::Meters:
757  if ( isArea )
758  {
759  if ( keepBaseUnit )
760  {
761  unitLabel = QObject::tr( " m2" );
762  }
763  else if ( qAbs( value ) > 1000000.0 )
764  {
765  unitLabel = QObject::tr( " km2" );
766  value = value / 1000000.0;
767  }
768  else if ( qAbs( value ) > 10000.0 )
769  {
770  unitLabel = QObject::tr( " ha" );
771  value = value / 10000.0;
772  }
773  else
774  {
775  unitLabel = QObject::tr( " m2" );
776  }
777  }
778  else
779  {
780  if ( keepBaseUnit || qAbs( value ) == 0.0 )
781  {
782  unitLabel = QObject::tr( " m" );
783  }
784  else if ( qAbs( value ) > 1000.0 )
785  {
786  unitLabel = QObject::tr( " km" );
787  value = value / 1000;
788  }
789  else if ( qAbs( value ) < 0.01 )
790  {
791  unitLabel = QObject::tr( " mm" );
792  value = value * 1000;
793  }
794  else if ( qAbs( value ) < 0.1 )
795  {
796  unitLabel = QObject::tr( " cm" );
797  value = value * 100;
798  }
799  else
800  {
801  unitLabel = QObject::tr( " m" );
802  }
803  }
804  break;
805  case QGis::Feet:
806  if ( isArea )
807  {
808  if ( keepBaseUnit || qAbs( value ) <= ( 528.0*528.0 ) )
809  {
810  unitLabel = QObject::tr( " sq ft" );
811  }
812  else
813  {
814  unitLabel = QObject::tr( " sq mile" );
815  value = value / ( 5280.0 * 5280.0 );
816  }
817  }
818  else
819  {
820  if ( qAbs( value ) <= 528.0 || keepBaseUnit )
821  {
822  if ( qAbs( value ) == 1.0 )
823  {
824  unitLabel = QObject::tr( " foot" );
825  }
826  else
827  {
828  unitLabel = QObject::tr( " feet" );
829  }
830  }
831  else
832  {
833  unitLabel = QObject::tr( " mile" );
834  value = value / 5280.0;
835  }
836  }
837  break;
838  case QGis::Degrees:
839  if ( isArea )
840  {
841  unitLabel = QObject::tr( " sq.deg." );
842  }
843  else
844  {
845  if ( qAbs( value ) == 1.0 )
846  unitLabel = QObject::tr( " degree" );
847  else
848  unitLabel = QObject::tr( " degrees" );
849  }
850  break;
851  case QGis::UnknownUnit:
852  unitLabel = QObject::tr( " unknown" );
853  default:
854  QgsDebugMsg( QString( "Error: not picked up map units - actual value = %1" ).arg( u ) );
855  };
856 
857 
858  return QLocale::system().toString( value, 'f', decimals ) + unitLabel;
859 }