37 #define DEG2RAD(x) ((x)*M_PI/180) 38 #define RAD2DEG(r) (180.0 * (r) / M_PI) 39 #define POW2(x) ((x)*(x)) 46 mInvFlattening = -1.0;
80 setFromParams( params );
90 mSemiMajor = semiMajor;
91 mSemiMinor = semiMinor;
92 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
99 double QgsDistanceArea::measure(
const QgsAbstractGeometry *geomV2, MeasureType type )
const 107 if ( geomDimension <= 0 )
112 MeasureType measureType = type;
113 if ( measureType == Default )
115 measureType = ( geomDimension == 1 ? Length : Area );
121 if ( measureType == Length )
127 return geomV2->
area();
139 sum += measure( collection->
geometryN( i ), measureType );
144 if ( measureType == Length )
186 return measure( geomV2, Area );
195 return measure( geomV2, Length );
204 if ( !geomV2 || geomV2->
dimension() < 2 )
215 QVector< const QgsSurface * > surfaces;
219 surfaces.append( surf );
224 surfaces.reserve( ( surf ? 1 : 0 ) + multiSurf->
numGeometries() );
227 surfaces.append( static_cast<const QgsSurface *>( multiSurf->
geometryN( i ) ) );
232 QVector<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
233 for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
244 length += measure( outerRing );
247 for (
int i = 0; i < nInnerRings; ++i )
264 QVector<QgsPointXY> linePoints;
265 curve->
points( linePointsV2 );
272 if ( points.size() < 2 )
281 p1 = mCoordTransform.
transform( points[0] );
285 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
290 total += computeDistanceBearing( p1, p2 );
306 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
323 QgsDebugMsgLevel( QString(
"Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
328 QgsDebugMsgLevel( QString(
"New points are %1 and %2, calculating..." ).arg( pp1.
toString( 4 ), pp2.toString( 4 ) ), 4 );
329 result = computeDistanceBearing( pp1, pp2 );
340 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
364 p2 = p1.
project( distance, azimuth );
366 QgsDebugMsgLevel( QString(
"Converted distance of %1 %2 to %3 distance %4 %5, using azimuth[%6] from point[%7] to point[%8] sourceCrs[%9] mEllipsoid[%10] isGeographic[%11] [%12]" )
367 .arg( QString::number( distance,
'f', 7 ) )
369 .arg( QString::number( result,
'f', 7 ) )
370 .arg( ( ( mCoordTransform.
sourceCrs().
isGeographic() ) == 1 ? QString(
"Geographic" ) : QString(
"Cartesian" ) ) )
378 .arg( QString(
"SemiMajor[%1] SemiMinor[%2] InvFlattening[%3] " ).arg( QString::number( mSemiMajor,
'f', 7 ) ).arg( QString::number( mSemiMinor,
'f', 7 ) ).arg( QString::number( mInvFlattening,
'f', 7 ) ) ), 4 );
379 if ( projectedPoint )
394 const QgsPointXY &p1,
double distance,
double azimuth )
const 397 double a = mSemiMajor;
398 double b = mSemiMinor;
399 double f = 1 / mInvFlattening;
400 if ( ( ( a < 0 ) && ( b < 0 ) ) ||
401 ( ( p1.
x() < -180.0 ) || ( p1.
x() > 180.0 ) || ( p1.
y() < -85.05115 ) || ( p1.
y() > 85.05115 ) ) )
407 double radians_lat =
DEG2RAD( p1.
y() );
408 double radians_long =
DEG2RAD( p1.
x() );
409 double b2 =
POW2( b );
411 double tan_u1 = omf * std::tan( radians_lat );
412 double u1 = std::atan( tan_u1 );
413 double sigma, last_sigma, delta_sigma, two_sigma_m;
414 double sigma1, sin_alpha, alpha, cos_alphasq;
416 double lat2, lambda, lambda2, C, omega;
420 azimuth = azimuth + M_PI * 2.0;
422 if ( azimuth > ( M_PI * 2.0 ) )
424 azimuth = azimuth - M_PI * 2.0;
426 sigma1 = std::atan2( tan_u1, std::cos( azimuth ) );
427 sin_alpha = std::cos( u1 ) * std::sin( azimuth );
428 alpha = std::asin( sin_alpha );
429 cos_alphasq = 1.0 -
POW2( sin_alpha );
430 u2 =
POW2( std::cos( alpha ) ) * (
POW2( a ) - b2 ) / b2;
431 A = 1.0 + ( u2 / 16384.0 ) * ( 4096.0 + u2 * ( -768.0 + u2 * ( 320.0 - 175.0 * u2 ) ) );
432 B = ( u2 / 1024.0 ) * ( 256.0 + u2 * ( -128.0 + u2 * ( 74.0 - 47.0 * u2 ) ) );
433 sigma = ( distance / ( b * A ) );
436 two_sigma_m = 2.0 * sigma1 + sigma;
437 delta_sigma = B * std::sin( sigma ) * ( std::cos( two_sigma_m ) + ( B / 4.0 ) * ( std::cos( sigma ) * ( -1.0 + 2.0 *
POW2( std::cos( two_sigma_m ) ) - ( B / 6.0 ) * std::cos( two_sigma_m ) * ( -3.0 + 4.0 *
POW2( std::sin( sigma ) ) ) * ( -3.0 + 4.0 *
POW2( std::cos( two_sigma_m ) ) ) ) ) );
439 sigma = ( distance / ( b * A ) ) + delta_sigma;
442 while ( i < 999 && std::fabs( ( last_sigma - sigma ) / sigma ) > 1.0e-9 );
444 lat2 = std::atan2( ( std::sin( u1 ) * std::cos( sigma ) + std::cos( u1 ) * std::sin( sigma ) *
445 std::cos( azimuth ) ), ( omf * std::sqrt(
POW2( sin_alpha ) +
446 POW2( std::sin( u1 ) * std::sin( sigma ) - std::cos( u1 ) * std::cos( sigma ) *
447 std::cos( azimuth ) ) ) ) );
448 lambda = std::atan2( ( std::sin( sigma ) * std::sin( azimuth ) ), ( std::cos( u1 ) * std::cos( sigma ) -
449 std::sin( u1 ) * std::sin( sigma ) * std::cos( azimuth ) ) );
450 C = ( f / 16.0 ) * cos_alphasq * ( 4.0 + f * ( 4.0 - 3.0 * cos_alphasq ) );
451 omega = lambda - ( 1.0 - C ) * f * sin_alpha * ( sigma + C * std::sin( sigma ) *
452 ( std::cos( two_sigma_m ) + C * std::cos( sigma ) * ( -1.0 + 2.0 *
POW2( std::cos( two_sigma_m ) ) ) ) );
453 lambda2 = radians_long + omega;
476 curve->
points( linePointsV2 );
477 QVector<QgsPointXY> linePoints;
489 QVector<QgsPointXY> pts;
490 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
492 pts.append( mCoordTransform.
transform( *i ) );
494 return computePolygonArea( pts );
498 return computePolygonArea( points );
504 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
519 computeDistanceBearing( pp1, pp2, &bearing );
523 double dx = p2.
x() - p1.
x();
524 double dy = p2.
y() - p1.
y();
525 bearing = std::atan2( dx, dy );
535 double QgsDistanceArea::computeDistanceBearing(
537 double *course1,
double *course2 )
const 543 double a = mSemiMajor;
544 double b = mSemiMinor;
545 double f = 1 / mInvFlattening;
550 double L = p2_lon - p1_lon;
551 double U1 = std::atan( ( 1 - f ) * std::tan( p1_lat ) );
552 double U2 = std::atan( ( 1 - f ) * std::tan( p2_lat ) );
553 double sinU1 = std::sin( U1 ), cosU1 = std::cos( U1 );
554 double sinU2 = std::sin( U2 ), cosU2 = std::cos( U2 );
556 double lambdaP = 2 * M_PI;
558 double sinLambda = 0;
559 double cosLambda = 0;
564 double cosSqAlpha = 0;
565 double cos2SigmaM = 0;
571 while ( std::fabs( lambda - lambdaP ) > 1e-12 && --iterLimit > 0 )
573 sinLambda = std::sin( lambda );
574 cosLambda = std::cos( lambda );
575 tu1 = ( cosU2 * sinLambda );
576 tu2 = ( cosU1 * sinU2 - sinU1 * cosU2 * cosLambda );
577 sinSigma = std::sqrt( tu1 * tu1 + tu2 * tu2 );
578 cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
579 sigma = std::atan2( sinSigma, cosSigma );
580 alpha = std::asin( cosU1 * cosU2 * sinLambda / sinSigma );
581 cosSqAlpha = std::cos( alpha ) * std::cos( alpha );
582 cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
583 C = f / 16 * cosSqAlpha * ( 4 + f * ( 4 - 3 * cosSqAlpha ) );
585 lambda = L + ( 1 - C ) * f * std::sin( alpha ) *
586 ( sigma + C * sinSigma * ( cos2SigmaM + C * cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) ) );
589 if ( iterLimit == 0 )
592 double uSq = cosSqAlpha * ( a * a - b * b ) / ( b * b );
593 double A = 1 + uSq / 16384 * ( 4096 + uSq * ( -768 + uSq * ( 320 - 175 * uSq ) ) );
594 double B = uSq / 1024 * ( 256 + uSq * ( -128 + uSq * ( 74 - 47 * uSq ) ) );
595 double deltaSigma = B * sinSigma * ( cos2SigmaM + B / 4 * ( cosSigma * ( -1 + 2 * cos2SigmaM * cos2SigmaM ) -
596 B / 6 * cos2SigmaM * ( -3 + 4 * sinSigma * sinSigma ) * ( -3 + 4 * cos2SigmaM * cos2SigmaM ) ) );
597 double s = b * A * ( sigma - deltaSigma );
601 *course1 = std::atan2( tu1, tu2 );
606 *course2 = std::atan2( cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda ) + M_PI;
617 double QgsDistanceArea::getQ(
double x )
const 621 sinx = std::sin( x );
624 return sinx * ( 1 + sinx2 * ( m_QA + sinx2 * ( m_QB + sinx2 * m_QC ) ) );
628 double QgsDistanceArea::getQbar(
double x )
const 632 cosx = std::cos( x );
635 return cosx * ( m_QbarA + cosx2 * ( m_QbarB + cosx2 * ( m_QbarC + cosx2 * m_QbarD ) ) );
639 void QgsDistanceArea::computeAreaInit()
647 double a2 = ( mSemiMajor * mSemiMajor );
648 double e2 = 1 - ( a2 / ( mSemiMinor * mSemiMinor ) );
651 m_TwoPI = M_PI + M_PI;
656 m_AE = a2 * ( 1 - e2 );
658 m_QA = ( 2.0 / 3.0 ) * e2;
659 m_QB = ( 3.0 / 5.0 ) * e4;
660 m_QC = ( 4.0 / 7.0 ) * e6;
662 m_QbarA = -1.0 - ( 2.0 / 3.0 ) * e2 - ( 3.0 / 5.0 ) * e4 - ( 4.0 / 7.0 ) * e6;
663 m_QbarB = ( 2.0 / 9.0 ) * e2 + ( 2.0 / 5.0 ) * e4 + ( 4.0 / 7.0 ) * e6;
664 m_QbarC = - ( 3.0 / 25.0 ) * e4 - ( 12.0 / 35.0 ) * e6;
665 m_QbarD = ( 4.0 / 49.0 ) * e6;
667 m_Qp = getQ( M_PI_2 );
668 m_E = 4 * M_PI * m_Qp * m_AE;
690 double QgsDistanceArea::computePolygonArea(
const QVector<QgsPointXY> &points )
const 692 if ( points.isEmpty() )
701 double x1, y1, x2, y2, dx, dy;
708 const double thresh = 1e-6;
713 return computePolygonFlatArea( points );
715 int n = points.size();
716 x2 =
DEG2RAD( points[n - 1].x() );
717 y2 =
DEG2RAD( points[n - 1].y() );
718 Qbar2 = getQbar( y2 );
722 for (
int i = 0; i < n; i++ )
730 Qbar2 = getQbar( y2 );
733 while ( x1 - x2 > M_PI )
736 while ( x2 - x1 > M_PI )
741 if ( std::fabs( dy ) > thresh )
744 area += dx * ( m_Qp - ( Qbar2 - Qbar1 ) / dy );
756 area += dx * ( m_Qp - getQ( ( y1 + y2 ) / 2.0 ) );
763 if ( ( area *= m_AE ) < 0.0 )
773 if ( area > m_E / 2 )
779 double QgsDistanceArea::computePolygonFlatArea(
const QVector<QgsPointXY> &points )
const 785 size = points.size();
788 for ( i = 0; i < size; i++ )
793 area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
797 return std::fabs( area );
816 double result = length * factorUnits;
817 QgsDebugMsgLevel( QString(
"Converted length of %1 %2 to %3 %4" ).arg( length )
830 double result = area * factorUnits;
const QgsCurve * interiorRing(int i) const
QgsUnitTypes::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)
static Q_INVOKABLE QString toString(QgsUnitTypes::DistanceUnit unit)
Returns a translated string representing a distance unit.
bool isNull() const
Returns true if the geometry is null (ie, contains no underlying geometry accessible via geometry() )...
QgsPolygon * surfaceToPolygon() const override
Get a polygon representation of this surface.
A class to represent a 2D point.
QString toProj4() const
Returns a Proj4 string representation of this CRS.
A geometry is the spatial representation of a feature.
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
double convertLengthMeasurement(double length, QgsUnitTypes::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
QgsUnitTypes::DistanceUnit mapUnits
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
double convertAreaMeasurement(double area, QgsUnitTypes::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
Contains parameters for an ellipsoid.
double measurePolygon(const QVector< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
Compare two doubles (but allow some difference)
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::Warning)
add a message to the instance (and create it if necessary)
Multi surface geometry collection.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
const QString GEO_NONE
Constant that holds the string representation for "No ellips/No CRS".
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
QgsPointXY computeSpheroidProject(const QgsPointXY &p1, double distance=1, double azimuth=M_PI_2) const
Given a location, an azimuth and a distance, computes the location of the projected point...
bool willUseEllipsoid() const
Returns whether calculations will use the ellipsoid.
virtual QgsPolygon * surfaceToPolygon() const =0
Get a polygon representation of this surface.
static Q_INVOKABLE QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
int numInteriorRings() const
virtual double length() const
Returns the length of the geometry.
bool valid
Whether ellipsoid parameters are valid.
#define QgsDebugMsgLevel(str, level)
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
const long GEOCRS_ID
Magic number for a geographic coord sys in QGIS srs.db tbl_srs.srs_id.
T qgsgeometry_cast(const QgsAbstractGeometry *geom)
QString description() const
Returns the descriptive name of the CRS, e.g., "WGS 84" or "GDA 94 / Vicgrid94".
QgsPointXY project(double distance, double bearing) const
Returns a new point which corresponds to this point projected by a specified distance in a specified ...
virtual double area() const
Returns the area of the geometry.
bool isGeographic() const
Returns whether the CRS is a geographic CRS (using lat/lon coordinates)
Abstract base class for curved geometry type.
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
Abstract base class for all geometries.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
Contains information about the context in which a coordinate transform is executed.
const QgsCurve * exteriorRing() const
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
Point geometry type, with support for z-dimension and m-values.
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
int numGeometries() const
Returns the number of geometries within the collection.
double semiMinor
Semi-minor axis.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
DistanceUnit
Units of distance.
QVector< QgsPoint > QgsPointSequence
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
static void convertPointList(const QVector< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPoint.
QgsUnitTypes::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
QString ellipsoid() const
Returns ellipsoid's acronym.
static Q_INVOKABLE QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
double measureLineProjected(const QgsPointXY &p1, double distance=1, double azimuth=M_PI_2, QgsPointXY *projectedPoint=nullptr) const
Calculates the distance from one point with distance in meters and azimuth (direction) When the sourc...
Line string geometry type, with support for z-dimension and m-values.
virtual double perimeter() const
Returns the perimeter of the geometry.
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve...
This class represents a coordinate reference system (CRS).
double inverseFlattening
Inverse flattening.
static QString formatArea(double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
Custom exception class for Coordinate Reference System related exceptions.
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
static Q_INVOKABLE double fromUnitToUnitFactor(QgsUnitTypes::DistanceUnit fromUnit, QgsUnitTypes::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
QgsDistanceArea()
Constructor.
static Q_INVOKABLE QgsUnitTypes::AreaUnit distanceToAreaUnit(QgsUnitTypes::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
double semiMajor
Semi-major axis.
static QgsCoordinateReferenceSystem fromSrsId(long srsId)
Creates a CRS from a specified QGIS SRS ID.
QString asWkt() const
Return the well known text representation for the point (e.g.
static QString formatDistance(double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.