added measureLineProjected and computeSpheroidProject to QgsDistanceArea

python/core/qgsdistancearea.sip

@@ -146,6 +146,22 @@ class QgsDistanceArea
      */
     double measureLine( const QgsPoint& p1, const QgsPoint& p2 ) const;
 
+    /**
+     * calculates distance from one point with distance in meters and azimuth (direction)
+     * based on inverse Vincenty's formula when Geographic
+     * otherwise QgsPoint.project() will be called after QgsUnitTypes::fromUnitToUnitFactor() has been applied to the distance
+     * note: usage:
+     *  QgsDistanceArea myDa;
+      * myDa.setSourceCrs( layer->crs() );
+     * \since QGIS 3.0
+     * @param p1 start point [can be Cartesian or Geographic]
+     * @param distance expected to be in meters
+     * @param azimuth - azimuth in radians. [default M_PI/2 - East of p1]
+     * @param projectedPoint calculated projected point
+     * @return distance in mapUnits
+     */
+    double measureLineProjected( const QgsPoint &p1, double distance = 1, double azimuth = M_PI / 2, QgsPoint *projectedPoint /Out/ = nullptr ) const;
+
     /** Returns the units of distance for length calculations made by this object.
      * @note added in QGIS 2.14
      * @see areaUnits()
@@ -228,6 +244,25 @@ class QgsDistanceArea
     double computeDistanceBearing( const QgsPoint& p1, const QgsPoint& p2,
                                    double* course1 = 0, double* course2 = 0 ) const;
 
+    /**
+     * Given a location, an azimuth and a distance, computes the
+     * location of the projected point. Based on Vincenty's formula
+     * for the geodetic direct problem as described in "Geocentric
+     * Datum of Australia Technical Manual", Chapter 4. Tested against:
+     * http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
+     * and
+     * http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
+     * @note code (and documentation) taken from rttopo project
+     * https://git.osgeo.org/gogs/rttopo/librttopo
+     * - spheroid_project.spheroid_project(...)
+     * \since QGIS 3.0
+     * @param p1 - location of first Geographic (lat,long) point.
+     * @param distance - distance in meters. [default 1 meter]
+     * @param azimuth - azimuth in radians. [default M_PI/2 - East of p1]
+     * @return p2 - location of projected point.
+     */
+    QgsPoint computeSpheroidProject( const QgsPoint &p1, double distance = 1, double azimuth = M_PI / 2 ) const;
+
     //! uses flat / planimetric / Euclidean distance
     double computeDistanceFlat( const QgsPoint& p1, const QgsPoint& p2 ) const;
 

src/core/qgsdistancearea.cpp

@@ -43,6 +43,8 @@
 #endif
 
 #define DEG2RAD(x)    ((x)*M_PI/180)
+#define RAD2DEG(r) (180.0 * (r) / M_PI)
+#define POW2(x) ((x)*(x))
 
 
 QgsDistanceArea::QgsDistanceArea()
@@ -78,9 +80,12 @@ void QgsDistanceArea::_copy( const QgsDistanceArea &origDA )
   mSemiMajor = origDA.mSemiMajor;
   mSemiMinor = origDA.mSemiMinor;
   mInvFlattening = origDA.mInvFlattening;
-  // Some calculations and trig. Should not be TOO time consuming.
-  // Alternatively we could copy the temp vars?
-  computeAreaInit();
+  if ( ( mSemiMajor > 0 ) && ( mSemiMinor > 0 ) )
+  {
+    // Some calculations and trig. Should not be TOO time consuming.
+    // Alternatively we could copy the temp vars?
+    computeAreaInit();
+  }
   mCoordTransform = origDA.mCoordTransform;
 }
 
@@ -537,6 +542,92 @@ double QgsDistanceArea::measureLine( const QgsPoint &p1, const QgsPoint &p2, Qgs
   return result;
 }
 
+double QgsDistanceArea::measureLineProjected( const QgsPoint &p1, double distance, double azimuth, QgsPoint *projectedPoint ) const
+{
+  double result = 0.0;
+  QgsPoint p2;
+  if ( geographic() )
+  {
+    if ( mEllipsoid != GEO_NONE )
+    {
+      p2 = computeSpheroidProject( p1, distance, azimuth );
+      result = p1.distance( p2 );
+    }
+  }
+  else // cartesian coordinates
+  {
+    result = distance; // Avoid rounding errors when using meters [return as sent]
+    if ( sourceCrs().mapUnits() != QgsUnitTypes::DistanceMeters )
+    {
+      distance = ( distance * QgsUnitTypes::fromUnitToUnitFactor( QgsUnitTypes::DistanceMeters, sourceCrs().mapUnits() ) );
+      result = p1.distance( p2 );
+    }
+    p2 = p1.project( distance, azimuth );
+  }
+  if ( projectedPoint )
+  {
+    *projectedPoint = QgsPoint( p2 );
+  }
+  return result;
+}
+
+QgsPoint QgsDistanceArea::computeSpheroidProject(
+  const QgsPoint &p1, double distance, double azimuth ) const
+{
+  // ellipsoid
+  double a = mSemiMajor;
+  double b = mSemiMinor;
+  double f = 1 / mInvFlattening;
+  double radians_lat = DEG2RAD( p1.y() );
+  double radians_long = DEG2RAD( p1.x() );
+  double b2 = POW2( b ); // spheroid_mu2
+  double omf = 1 - f;
+  double tan_u1 = omf * tan( radians_lat );
+  double u1 = atan( tan_u1 );
+  double sigma, last_sigma, delta_sigma, two_sigma_m;
+  double sigma1, sin_alpha, alpha, cos_alphasq;
+  double u2, A, B;
+  double lat2, lambda, lambda2, C, omega;
+  int i = 0;
+  if ( azimuth < 0.0 )
+  {
+    azimuth = azimuth + M_PI * 2.0;
+  }
+  if ( azimuth > ( M_PI * 2.0 ) )
+  {
+    azimuth = azimuth - M_PI * 2.0;
+  }
+  sigma1 = atan2( tan_u1, cos( azimuth ) );
+  sin_alpha = cos( u1 ) * sin( azimuth );
+  alpha = asin( sin_alpha );
+  cos_alphasq = 1.0 - POW2( sin_alpha );
+  u2 = POW2( cos( alpha ) ) * ( POW2( a ) - b2 ) / b2; // spheroid_mu2
+  A = 1.0 + ( u2 / 16384.0 ) * ( 4096.0 + u2 * ( -768.0 + u2 * ( 320.0 - 175.0 * u2 ) ) );
+  B = ( u2 / 1024.0 ) * ( 256.0 + u2 * ( -128.0 + u2 * ( 74.0 - 47.0 * u2 ) ) );
+  sigma = ( distance / ( b * A ) );
+  do
+  {
+    two_sigma_m = 2.0 * sigma1 + sigma;
+    delta_sigma = B * sin( sigma ) * ( cos( two_sigma_m ) + ( B / 4.0 ) * ( cos( sigma ) * ( -1.0 + 2.0 * POW2( cos( two_sigma_m ) ) - ( B / 6.0 ) * cos( two_sigma_m ) * ( -3.0 + 4.0 * POW2( sin( sigma ) ) ) * ( -3.0 + 4.0 * POW2( cos( two_sigma_m ) ) ) ) ) );
+    last_sigma = sigma;
+    sigma = ( distance / ( b * A ) ) + delta_sigma;
+    i++;
+  }
+  while ( i < 999 && qAbs( ( last_sigma - sigma ) / sigma ) > 1.0e-9 );
+
+  lat2 = atan2( ( sin( u1 ) * cos( sigma ) + cos( u1 ) * sin( sigma ) *
+                  cos( azimuth ) ), ( omf * sqrt( POW2( sin_alpha ) +
+                                      POW2( sin( u1 ) * sin( sigma ) - cos( u1 ) * cos( sigma ) *
+                                          cos( azimuth ) ) ) ) );
+  lambda = atan2( ( sin( sigma ) * sin( azimuth ) ), ( cos( u1 ) * cos( sigma ) -
+                  sin( u1 ) * sin( sigma ) * cos( azimuth ) ) );
+  C = ( f / 16.0 ) * cos_alphasq * ( 4.0 + f * ( 4.0 - 3.0 * cos_alphasq ) );
+  omega = lambda - ( 1.0 - C ) * f * sin_alpha * ( sigma + C * sin( sigma ) *
+          ( cos( two_sigma_m ) + C * cos( sigma ) * ( -1.0 + 2.0 * POW2( cos( two_sigma_m ) ) ) ) );
+  lambda2 = radians_long + omega;
+  return QgsPoint( RAD2DEG( lambda2 ), RAD2DEG( lat2 ) );
+}
+
 QgsUnitTypes::DistanceUnit QgsDistanceArea::lengthUnits() const
 {
   return willUseEllipsoid() ? QgsUnitTypes::DistanceMeters : mCoordTransform.sourceCrs().mapUnits();

src/core/qgsdistancearea.h

@@ -222,6 +222,23 @@ class CORE_EXPORT QgsDistanceArea
      */
     double measureLine( const QgsPoint &p1, const QgsPoint &p2, QgsUnitTypes::DistanceUnit &units ) const;
 
+    /**
+     * Calculates distance from one point with distance in meters and azimuth (direction)
+     * When the sourceCrs is Geographic, computeSpheroidProject will be called
+     * otherwise QgsPoint.project() will be called after QgsUnitTypes::fromUnitToUnitFactor() has been applied to the distance
+     * \note:
+     *  The input Point must be in the CoordinateReferenceSystem being used
+     * \since QGIS 3.0
+     * \param p1 start point [can be Cartesian or Geographic]
+     * \param distance must be in meters
+     * \param azimuth - azimuth in radians. [default M_PI/2 - East of p1]
+     * \param projectedPoint calculated projected point
+     * \return distance in mapUnits
+     * \see sourceCrs()
+     * \see computeSpheroidProject()
+     */
+    double measureLineProjected( const QgsPoint &p1, double distance = 1, double azimuth = M_PI / 2, QgsPoint *projectedPoint = nullptr ) const;
+
     /** Returns the units of distance for length calculations made by this object.
      * \since QGIS 2.14
      * \see areaUnits()
@@ -304,6 +321,30 @@ class CORE_EXPORT QgsDistanceArea
     double computeDistanceBearing( const QgsPoint &p1, const QgsPoint &p2,
                                    double *course1 = nullptr, double *course2 = nullptr ) const;
 
+    /**
+     * Given a location, an azimuth and a distance, computes the
+     * location of the projected point. Based on Vincenty's formula
+     * for the geodetic direct problem as described in "Geocentric
+     * Datum of Australia Technical Manual", Chapter 4.
+     * \note code (and documentation) taken from rttopo project
+     * https://git.osgeo.org/gogs/rttopo/librttopo
+     * - spheroid_project.spheroid_project(...)
+     * \since QGIS 3.0
+     * \param p1 - location of first Geographic (lat,long) point as degrees.
+     * \param distance - distance in meters. [default 1 meter]
+     * \param azimuth - azimuth in radians. [default M_PI/2 - East of p1]
+     * \return p2 - location of projected point as degrees.
+     */
+
+    /*
+     *  From original rttopo documentation:
+     *  Tested against:
+     *   http://mascot.gdbc.gov.bc.ca/mascot/util1b.html
+     *  and
+     *   http://www.ga.gov.au/nmd/geodesy/datums/vincenty_direct.jsp
+     */
+    QgsPoint computeSpheroidProject( const QgsPoint &p1, double distance = 1, double azimuth = M_PI / 2 ) const;
+
     //! uses flat / planimetric / Euclidean distance
     double computeDistanceFlat( const QgsPoint &p1, const QgsPoint &p2 ) const;