Backport for measuring from trunk to 0.8. Not everything is ported, but the main functionality should at least be correct.

git-svn-id: http://svn.osgeo.org/qgis/branches/Release-0_8_0@6818 c8812cc2-4d05-0410-92ff-de0c093fc19c

Author: homann

src/core/qgsdistancearea.cpp

@@ -17,27 +17,34 @@
 #include <cmath>
 #include <sqlite3.h>
 #include <QDir>
-#include <QSettings>
+#include <QString>
+#include <QLocale>
+#include <QObject>
 
 #include "qgis.h"
 #include "qgspoint.h"
-#include "qgsproject.h"
 #include "qgscoordinatetransform.h"
 #include "qgsspatialrefsys.h"
 #include "qgsgeometry.h"
 #include "qgsdistancearea.h"
 #include "qgsapplication.h"
 #include "qgslogger.h"
 
+// MSVC compiler doesn't have defined M_PI in math.h
+#ifndef M_PI
+#define M_PI          3.14159265358979323846
+#endif
+
 #define DEG2RAD(x)    ((x)*M_PI/180)
 
 
 QgsDistanceArea::QgsDistanceArea()
 {
   // init with default settings
+  mProjectionsEnabled = FALSE;
   mCoordTransform = new QgsCoordinateTransform;
-  setDefaultEllipsoid();
-  setProjectAsSourceSRS();
+  setSourceSRS(GEOSRS_ID); // WGS 84
+  setEllipsoid("WGS84");
 }
 
 
@@ -47,27 +54,17 @@ QgsDistanceArea::~QgsDistanceArea()
 }
 
 
-void QgsDistanceArea::setSourceSRS(long srsid)
+void QgsDistanceArea::setProjectionsEnabled(bool flag)
 {
-  QgsSpatialRefSys srcSRS;
-  srcSRS.createFromSrsId(srsid);
-  mCoordTransform->setSourceSRS(srcSRS);
+  mProjectionsEnabled = flag;
 }
 
 
-void QgsDistanceArea::setProjectAsSourceSRS()
+void QgsDistanceArea::setSourceSRS(long srsid)
 {
-  // This function used to only get the /ProjectSRSID if on-the-fly
-  // projection was enabled (and used a default value in all other
-  // cases). However, even if it was not, a valid value for
-  // /ProjectSRSID is most likely available, and we now use it in all
-  // cases (as it gives correct distances and areas even when
-  // on-the-fly projection is turned off). The default of GEOSRS_ID
-  // is still applied if all else fails.
-
-  int srsid = QgsProject::instance()->readNumEntry("SpatialRefSys","/ProjectSRSID",GEOSRS_ID);
-
-  setSourceSRS(srsid);
+  QgsSpatialRefSys srcSRS;
+  srcSRS.createFromSrsId(srsid);
+  mCoordTransform->setSourceSRS(srcSRS);
 }
 
 
@@ -81,6 +78,14 @@ bool QgsDistanceArea::setEllipsoid(const QString& ellipsoid)
   const char   *myTail;
   sqlite3_stmt *myPreparedStatement;
   int           myResult;
+
+  // Shortcut if ellipsoid is none.
+  if (ellipsoid == "NONE")
+  {
+    mEllipsoid = "NONE";
+    return true;
+  }
+
   //check the db is available
   myResult = sqlite3_open(QString(QgsApplication::qgisUserDbFilePath()).latin1(), &myDatabase);
   if(myResult) 
@@ -170,26 +175,6 @@ bool QgsDistanceArea::setEllipsoid(const QString& ellipsoid)
 }
 
 
-bool QgsDistanceArea::setDefaultEllipsoid()
-{
-  QString defEll("WGS84");
-  QString ellKey("/qgis/measure/ellipsoid");
-  QSettings settings;
-  QString ellipsoid = settings.readEntry(ellKey, defEll);
-
-  // Somehow/sometimes the settings file can have a blank ellipsoid
-  // value. This is undesirable, so force a valid default value in
-  // that case, and fix the problem by writing a valid value.
-  if (ellipsoid.isEmpty())
-  {
-    ellipsoid = defEll;
-    settings.writeEntry(ellKey, ellipsoid);
-  }
-
-  return setEllipsoid(ellipsoid);
-}
-
-
 double QgsDistanceArea::measure(QgsGeometry* geometry)
 {
   unsigned char* wkb = geometry->wkbBuffer();
@@ -246,6 +231,7 @@ unsigned char* QgsDistanceArea::measureLine(unsigned char* feature, double* area
 
   std::vector<QgsPoint> points(nPoints);
 
+  QgsDebugMsg("This feature WKB has " + QString::number(nPoints) + " points");
   // Extract the points from the WKB format into the vector
   for (unsigned int i = 0; i < nPoints; ++i)
   {
@@ -265,18 +251,32 @@ double QgsDistanceArea::measureLine(const std::vector<QgsPoint>& points)
   if (points.size() < 2)
     return 0;
 
+  double total = 0;
+  QgsPoint p1, p2;
+  
   try
   {
-    double total = 0;
-    QgsPoint p1, p2;
-    p1 = mCoordTransform->transform(points[0]);
-
-    for (int i = 1; i < points.size(); i++)
+    if (mProjectionsEnabled && (mEllipsoid != "NONE"))
+      p1 = mCoordTransform->transform(points[0]);
+    else
+      p1 = points[0];
+    
+    for (std::vector<QgsPoint>::size_type i = 1; i < points.size(); i++)
     {
-      p2 = mCoordTransform->transform(points[i]);
-      total = computeDistanceBearing(p1,p2);
+      if (mProjectionsEnabled && (mEllipsoid != "NONE"))
+      {
+        p2 = mCoordTransform->transform(points[i]);
+        total += computeDistanceBearing(p1,p2);
+      }
+      else
+      {
+        p2 = points[i];
+        total += measureLine(p1,p2);
+      }
+  
       p1 = p2;
     }
+    
     return total;
   }
   catch (QgsCsException &cse)
@@ -291,9 +291,17 @@ double QgsDistanceArea::measureLine(const QgsPoint& p1, const QgsPoint& p2)
 {
   try
   {
-    QgsPoint pp1 = mCoordTransform->transform(p1);
-    QgsPoint pp2 = mCoordTransform->transform(p2);
-    return computeDistanceBearing(pp1, pp2);
+    QgsPoint pp1 = p1, pp2 = p2;
+    if (mProjectionsEnabled && (mEllipsoid != "NONE"))
+    {
+      pp1 = mCoordTransform->transform(p1);
+      pp2 = mCoordTransform->transform(p2);
+      return computeDistanceBearing(pp1, pp2);
+    }
+    else
+    {
+      return sqrt((p2.x()-p1.x())*(p2.x()-p1.x()) + (p2.y()-p1.y())*(p2.y()-p1.y()));
+    }
   }
   catch (QgsCsException &cse)
   {
@@ -328,14 +336,18 @@ unsigned char* QgsDistanceArea::measurePolygon(unsigned char* feature, double* a
 
       // Extract the points from the WKB and store in a pair of
       // vectors.
-      for (unsigned int jdx = 0; jdx < nPoints; jdx++)
+      for (int jdx = 0; jdx < nPoints; jdx++)
       {
         x = *((double *) ptr);
         ptr += sizeof(double);
         y = *((double *) ptr);
         ptr += sizeof(double);
 
-        points[jdx] = mCoordTransform->transform(QgsPoint(x,y));
+        points[jdx] = QgsPoint(x,y);
+        if (mProjectionsEnabled && (mEllipsoid != "NONE"))
+        {
+          points[jdx] = mCoordTransform->transform(points[jdx]);
+        }
       }
 
       if (points.size() > 2)
@@ -359,14 +371,22 @@ unsigned char* QgsDistanceArea::measurePolygon(unsigned char* feature, double* a
 
 double QgsDistanceArea::measurePolygon(const std::vector<QgsPoint>& points)
 {
+  
   try
   {
-    std::vector<QgsPoint> pts(points.size());
-    for (std::vector<QgsPoint>::size_type i = 0; i < points.size(); i++)
+    if (mProjectionsEnabled && (mEllipsoid != "NONE"))
     {
-      pts[i] = mCoordTransform->transform(points[i]);
+      std::vector<QgsPoint> pts(points.size());
+      for (std::vector<QgsPoint>::size_type i = 0; i < points.size(); i++)
+      {
+        pts[i] = mCoordTransform->transform(points[i]);
+      }
+      return computePolygonArea(pts);
+    }
+    else
+    {
+      return computePolygonArea(points);
     }
-    return computePolygonArea(pts);
   }
   catch (QgsCsException &cse)
   {
@@ -376,6 +396,20 @@ double QgsDistanceArea::measurePolygon(const std::vector<QgsPoint>& points)
 }
 
 
+double QgsDistanceArea::getBearing(const QgsPoint& p1, const QgsPoint& p2)
+{
+  QgsPoint pp1 = p1, pp2 = p2;
+  if (mProjectionsEnabled && (mEllipsoid != "NONE"))
+  {
+    pp1 = mCoordTransform->transform(p1);
+    pp2 = mCoordTransform->transform(p2);
+  }
+  
+  double bearing;
+  computeDistanceBearing(pp1, pp2, &bearing);
+  return bearing;
+}
+
 
 ///////////////////////////////////////////////////////////
 // distance calculation
@@ -512,6 +546,11 @@ double QgsDistanceArea::computePolygonArea(const std::vector<QgsPoint>& points)
   double Qbar1, Qbar2;
   double area;
 
+  QgsDebugMsg("Ellipsoid: " + mEllipsoid);
+  if ((! mProjectionsEnabled) || (mEllipsoid == "NONE"))
+  {
+    return computePolygonFlatArea(points);
+  }
   int n = points.size();
   x2 = DEG2RAD(points[n-1].x());
   y2 = DEG2RAD(points[n-1].y());
@@ -556,3 +595,131 @@ double QgsDistanceArea::computePolygonArea(const std::vector<QgsPoint>& points)
   return area;
 }
 
+double QgsDistanceArea::computePolygonFlatArea(const std::vector<QgsPoint>& points)
+{
+  // Normal plane area calculations.
+  double area = 0.0;
+  int i, size;
+  
+  size = points.size();
+
+  // QgsDebugMsg("New area calc, nr of points: " + QString::number(size));
+  for(i = 0; i < size; i++)
+  {
+    // QgsDebugMsg("Area from point: " + (points[i]).stringRep(2));
+    // Using '% size', so that we always end with the starting point
+    // and thus close the polygon.
+    area = area + points[i].x()*points[(i+1) % size].y() - points[(i+1) % size].x()*points[i].y();
+  }
+  // QgsDebugMsg("Area from point: " + (points[i % size]).stringRep(2));
+  area = area / 2.0;
+  return fabs(area); // All areas are positive!
+}
+
+QString QgsDistanceArea::textUnit(double value, int decimals, QGis::units u, bool isArea)
+{
+  QString unitLabel;
+
+
+  switch (u)
+  {
+  case QGis::METERS: 
+    if (isArea)
+    {
+      if (fabs(value) > 1000000.0)
+      {
+	unitLabel = QObject::tr(" km2");
+	value = value / 1000000.0;
+      }
+      else if (fabs(value) > 1000.0)
+      {
+	unitLabel = QObject::tr(" ha");
+	value = value / 10000.0;
+      }
+      else
+      {
+	unitLabel = QObject::tr(" m2");
+      }
+    }
+    else
+    {
+      if (fabs(value) == 0.0)
+      {
+        // Special case for pretty printing.
+	unitLabel=QObject::tr(" m"); 
+        
+      }
+      else if (fabs(value) > 1000.0)
+      {
+	unitLabel=QObject::tr(" km");
+	value = value/1000;
+      }
+      else if (fabs(value) < 0.01)
+      {
+	unitLabel=QObject::tr(" mm");
+	value = value*1000;
+      }
+      else if (fabs(value) < 0.1)
+      {
+	unitLabel=QObject::tr(" cm");
+	value = value*100;
+      }
+      else
+      {
+	unitLabel=QObject::tr(" m"); 
+      }
+    }
+    break;
+  case QGis::FEET:
+    if (isArea)
+    {
+      if (fabs(value) > (528.0*528.0))
+      {
+        unitLabel = QObject::tr(" sq mile");
+        value = value / (5280.0*5280.0);
+      }
+      else
+      {
+        unitLabel = QObject::tr(" sq ft");
+      }
+    }
+    else
+    {
+      if (fabs(value) > 528.0)
+      {
+        unitLabel = QObject::tr(" mile");
+        value = value / 5280.0;
+      }
+      else
+      {
+        if (fabs(value) == 1.0)
+          unitLabel=QObject::tr(" foot"); 
+        else
+          unitLabel=QObject::tr(" feet"); 
+      }
+    }
+    break;
+  case QGis::DEGREES:
+    if (isArea)
+    {
+      unitLabel = QObject::tr(" sq.deg.");
+    }
+    else
+    {
+      if (fabs(value) == 1.0)
+	unitLabel=QObject::tr(" degree"); 
+      else
+	unitLabel=QObject::tr(" degrees"); 
+    }
+    break;
+  case QGis::UNKNOWN:
+    unitLabel=QObject::tr(" unknown");
+  default: 
+    std::cout << "Error: not picked up map units - actual value = " 
+	      << u << std::endl;
+  };
+
+
+  return QLocale::system().toString(value, 'f', decimals) + unitLabel;
+
+}