40#define DEG2RAD(x) ((x)*M_PI/180)
41#define RAD2DEG(r) (180.0 * (r) / M_PI)
42#define POW2(x) ((x)*(x))
49 mInvFlattening = -1.0;
58 : mCoordTransform( other.mCoordTransform )
59 , mEllipsoid( other.mEllipsoid )
60 , mSemiMajor( other.mSemiMajor )
61 , mSemiMinor( other.mSemiMinor )
62 , mInvFlattening( other.mInvFlattening )
69 mCoordTransform = other.mCoordTransform;
70 mEllipsoid = other.mEllipsoid;
71 mSemiMajor = other.mSemiMajor;
72 mSemiMinor = other.mSemiMinor;
73 mInvFlattening = other.mInvFlattening;
108 setFromParams( params );
118 mSemiMajor = semiMajor;
119 mSemiMinor = semiMinor;
120 mInvFlattening = mSemiMajor / ( mSemiMajor - mSemiMinor );
127double QgsDistanceArea::measure(
const QgsAbstractGeometry *geomV2, MeasureType type )
const
134 const int geomDimension = geomV2->
dimension();
135 if ( geomDimension <= 0 )
140 MeasureType measureType = type;
141 if ( measureType == Default )
143 measureType = ( geomDimension == 1 ? Length : Area );
149 if ( measureType == Length )
155 return geomV2->
area();
167 sum += measure( collection->
geometryN( i ), measureType );
172 if ( measureType == Length )
174 const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( geomV2 );
187 const QgsSurface *surface = qgsgeometry_cast<const QgsSurface *>( geomV2 );
192 QgsCurvePolygon *curvePolygon = qgsgeometry_cast<QgsCurvePolygon *>( surface );
218 return measure( geomV2, Area );
227 return measure( geomV2, Length );
236 if ( !geomV2 || geomV2->
dimension() < 2 )
247 QVector< const QgsSurface * > surfaces;
248 const QgsSurface *surf = qgsgeometry_cast<const QgsSurface *>( geomV2 );
251 surfaces.append( surf );
253 const QgsMultiSurface *multiSurf = qgsgeometry_cast<const QgsMultiSurface *>( geomV2 );
256 surfaces.reserve( ( surf ? 1 : 0 ) + multiSurf->
numGeometries() );
264 QVector<const QgsSurface *>::const_iterator surfaceIt = surfaces.constBegin();
265 for ( ; surfaceIt != surfaces.constEnd(); ++surfaceIt )
272 QgsCurvePolygon *curvePolygon = qgsgeometry_cast<QgsCurvePolygon *>( *surfaceIt );
279 length += measure( outerRing );
282 for (
int i = 0; i < nInnerRings; ++i )
300 QVector<QgsPointXY> linePoints;
301 curve->
points( linePointsV2 );
308 if ( points.size() < 2 )
318 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::measureLine()",
"Error creating geod_geodesic object" );
326 p1 = mCoordTransform.
transform( points[0] );
330 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
339 geod_inverse( mGeod.get(), p1.
y(), p1.
x(), p2.
y(), p2.
x(), &distance, &azimuth1, &azimuth2 );
356 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
370 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::measureLine()",
"Error creating geod_geodesic object" );
382 QgsDebugMsgLevel( QStringLiteral(
"Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
387 QgsDebugMsgLevel( QStringLiteral(
"New points are %1 and %2, calculating..." ).arg( pp1.
toString( 4 ), pp2.toString( 4 ) ), 4 );
391 geod_inverse( mGeod.get(), pp1.
y(), pp1.
x(), pp2.y(), pp2.x(), &result, &azimuth1, &azimuth2 );
395 QgsDebugMsgLevel( QStringLiteral(
"Cartesian calculation on canvas coordinates" ), 4 );
402 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate line length." ) );
426 p2 = p1.
project( distance, azimuth );
428 QgsDebugMsgLevel( QStringLiteral(
"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]" )
429 .arg( QString::number( distance,
'f', 7 ),
431 QString::number( result,
'f', 7 ),
432 mCoordTransform.
sourceCrs().
isGeographic() ? QStringLiteral(
"Geographic" ) : QStringLiteral(
"Cartesian" ),
440 .arg( QStringLiteral(
"SemiMajor[%1] SemiMinor[%2] InvFlattening[%3] " ).arg( QString::number( mSemiMajor,
'f', 7 ), QString::number( mSemiMinor,
'f', 7 ), QString::number( mInvFlattening,
'f', 7 ) ) ), 4 );
441 if ( projectedPoint )
449 const QgsPointXY &p1,
double distance,
double azimuth )
const
459 geod_direct( mGeod.get(), p1.
y(), p1.
x(),
RAD2DEG( azimuth ), distance, &lat2, &lon2, &azimuth2 );
469 p1.
setX( p1.
x() + 360 );
471 p2.
setX( p2.
x() + 360 );
474 double p1x = p1.
x() < 180 ? p1.
x() : p2.
x();
475 double p1y = p1.
x() < 180 ? p1.
y() : p2.
y();
476 double p2x = p1.
x() < 180 ? p2.
x() : p1.
x();
477 double p2y = p1.
x() < 180 ? p2.
y() : p1.
y();
485 fractionAlongLine = ( 180 - p1x ) / ( p2x - p1x );
487 fractionAlongLine = 1 - fractionAlongLine;
488 return p1y + ( 180 - p1x ) / ( p2x - p1x ) * ( p2y - p1y );
493 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::latitudeGeodesicCrossesAntimeridian()",
"Error creating geod_geodesic object" );
497 geod_geodesicline line;
498 geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
500 const double totalDist = line.s13;
501 double intersectionDist = line.s13;
506 while ( std::fabs( lon - 180.0 ) > 0.00000001 && iterations < 100 )
508 if ( iterations > 0 && std::fabs( p2x - p1x ) > 5 )
521 QgsDebugMsgLevel( QStringLiteral(
"Narrowed window to %1, %2 - %3, %4" ).arg( p1x ).arg( p1y ).arg( p2x ).arg( p2y ), 4 );
523 geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
524 intersectionDist = line.s13 * 0.5;
531 intersectionDist *= ( 180.0 - p1x ) / ( lon - p1x );
536 geod_position( &line, intersectionDist, &lat, &lon, &t );
542 QgsDebugMsgLevel( QStringLiteral(
"After %1 iterations lon is %2, lat is %3, dist from p1: %4" ).arg( iterations ).arg( lon ).arg( lat ).arg( intersectionDist ), 4 );
545 fractionAlongLine = intersectionDist / totalDist;
547 fractionAlongLine = 1 - fractionAlongLine;
563 std::unique_ptr< QgsMultiLineString > res = std::make_unique< QgsMultiLineString >();
566 const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( *part );
574 const std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >();
581 QVector< QgsPoint > newPoints;
589 for (
int i = 0; i < line->
numPoints(); i++ )
595 x = std::fmod( x, 360.0 );
606 if ( i > 0 && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
615 z = prevZ + ( p.
z() - prevZ ) * fract;
619 m = prevM + ( p.
m() - prevM ) * fract;
623 if ( prevLon < -120 )
628 QgsPoint newPoint( antiMeridianPoint );
634 if ( std::isfinite( newPoint.
x() ) && std::isfinite( newPoint.
y() ) )
636 newPoints << newPoint;
641 newPoints.reserve( line->
numPoints() - i + 1 );
648 if ( std::isfinite( antiMeridianPoint.
x() ) && std::isfinite( antiMeridianPoint.
y() ) )
652 newPoint.
setX( antiMeridianPoint.
x() );
653 newPoint.
setY( antiMeridianPoint.
y() );
654 newPoints << newPoint;
670 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform linestring. Unable to calculate break point." ) );
671 res->addGeometry( line->
clone() );
684 return QVector< QVector< QgsPointXY > >() << ( QVector< QgsPointXY >() << p1 << p2 );
690 return QVector< QVector< QgsPointXY > >();
700 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate geodesic line." ) );
701 return QVector< QVector< QgsPointXY > >();
704 geod_geodesicline line;
705 geod_inverseline( &line, mGeod.get(), pp1.y(), pp1.x(), pp2.y(), pp2.x(), GEOD_ALL );
706 const double totalDist = line.s13;
708 QVector< QVector< QgsPointXY > > res;
709 QVector< QgsPointXY > currentPart;
712 double prevLon = pp1.
x();
713 double prevLat = pp1.y();
714 bool lastRun =
false;
728 geod_position( &line, d, &lat, &lon, &t );
731 if ( breakLine && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
742 if ( prevLon < -120 )
747 if ( std::isfinite( p.
x() ) && std::isfinite( p.
y() ) )
765 if ( std::isfinite( p.
x() ) && std::isfinite( p.
y() ) )
791 if ( d >= totalDist )
817 curve->
points( linePointsV2 );
818 QVector<QgsPointXY> linePoints;
830 QVector<QgsPointXY> pts;
831 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
833 pts.append( mCoordTransform.
transform( *i ) );
835 return computePolygonArea( pts );
839 return computePolygonArea( points );
845 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate polygon area." ) );
863 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::bearing()",
"Error creating geod_geodesic object" );
870 geod_inverse( mGeod.get(), pp1.
y(), pp1.
x(), pp2.y(), pp2.x(), &distance, &azimuth1, &azimuth2 );
876 const double dx = p2.
x() - p1.
x();
877 const double dy = p2.
y() - p1.
y();
882 bearing = std::atan2( dx, dy );
888void QgsDistanceArea::computeAreaInit()
const
897 mGeod.reset(
new geod_geodesic() );
898 geod_init( mGeod.get(), mSemiMajor, 1 / mInvFlattening );
917double QgsDistanceArea::computePolygonArea(
const QVector<QgsPointXY> &points )
const
919 if ( points.isEmpty() )
927 return computePolygonFlatArea( points );
932 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::computePolygonArea()",
"Error creating geod_geodesic object" );
936 struct geod_polygon p;
937 geod_polygon_init( &p, 0 );
939 const bool isClosed = points.constFirst() == points.constLast();
944 int i = points.size();
945 while ( ( isClosed && --i ) || ( !isClosed && --i >= 0 ) )
946 geod_polygon_addpoint( mGeod.get(), &p, points.at( i ).y(), points.at( i ).x() );
949 double perimeter = 0;
950 geod_polygon_compute( mGeod.get(), &p, 0, 1, &area, &perimeter );
952 return std::fabs( area );
955double QgsDistanceArea::computePolygonFlatArea(
const QVector<QgsPointXY> &points )
const
961 size = points.size();
964 for ( i = 0; i < size; i++ )
969 area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
973 return std::fabs( area );
992 const double result = length * factorUnits;
993 QgsDebugMsgLevel( QStringLiteral(
"Converted length of %1 %2 to %3 %4" ).arg( length )
1006 const double result = area * factorUnits;
1007 QgsDebugMsgLevel( QStringLiteral(
"Converted area of %1 %2 to %3 %4" ).arg( area )
DistanceUnit
Units of distance.
@ SquareMeters
Square meters.
@ Reverse
Reverse/inverse transform (from destination to source)
Abstract base class for all geometries.
bool isMeasure() const
Returns true if the geometry contains m values.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual double perimeter() const
Returns the planar, 2-dimensional perimeter of the geometry.
virtual double length() const
Returns the planar, 2-dimensional length of the geometry.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual double area() const
Returns the planar, 2-dimensional area of the geometry.
This class represents a coordinate reference system (CRS).
QString toProj() const
Returns a Proj string representation of this CRS.
Qgis::DistanceUnit mapUnits
Contains information about the context in which a coordinate transform is executed.
Custom exception class for Coordinate Reference System related exceptions.
Curve polygon geometry type.
int numInteriorRings() const
Returns the number of interior rings contained with the curve polygon.
const QgsCurve * exteriorRing() const
Returns the curve polygon's exterior ring.
const QgsCurve * interiorRing(int i) const
Retrieves an interior ring from the curve polygon.
virtual QgsPolygon * surfaceToPolygon() const
Gets a polygon representation of this surface.
Abstract base class for curved geometry type.
virtual void points(QgsPointSequence &pt) const =0
Returns a list of points within the curve.
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.
A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
QgsDistanceArea & operator=(const QgsDistanceArea &other)
double latitudeGeodesicCrossesAntimeridian(const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine) const
Calculates the latitude at which the geodesic line joining p1 and p2 crosses the antimeridian (longit...
static QString formatDistance(double distance, int decimals, Qgis::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
QgsCoordinateReferenceSystem sourceCrs() const
Returns the source spatial reference system.
double measureArea(const QgsGeometry &geometry) const
Measures the area of a geometry.
double convertLengthMeasurement(double length, Qgis::DistanceUnit toUnits) const
Takes a length measurement calculated by this QgsDistanceArea object and converts it to a different d...
QVector< QVector< QgsPointXY > > geodesicLine(const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine=false) const
Calculates the geodesic line between p1 and p2, which represents the shortest path on the ellipsoid b...
double measurePerimeter(const QgsGeometry &geometry) const
Measures the perimeter of a polygon geometry.
double measureLength(const QgsGeometry &geometry) const
Measures the length of a geometry.
double bearing(const QgsPointXY &p1, const QgsPointXY &p2) const
Computes the bearing (in radians) between two points.
QString ellipsoid() const
Returns ellipsoid's acronym.
double measureLine(const QVector< QgsPointXY > &points) const
Measures the length of a line with multiple segments.
void setSourceCrs(const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context)
Sets source spatial reference system crs.
QgsGeometry splitGeometryAtAntimeridian(const QgsGeometry &geometry) const
Splits a (Multi)LineString geometry at the antimeridian (longitude +/- 180 degrees).
Qgis::DistanceUnit lengthUnits() const
Returns the units of distance for length calculations made by this object.
double measurePolygon(const QVector< QgsPointXY > &points) const
Measures the area of the polygon described by a set of points.
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...
bool setEllipsoid(const QString &ellipsoid)
Sets the ellipsoid by its acronym.
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.
double convertAreaMeasurement(double area, Qgis::AreaUnit toUnits) const
Takes an area measurement calculated by this QgsDistanceArea object and converts it to a different ar...
static QString formatArea(double area, int decimals, Qgis::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
Qgis::AreaUnit areaUnits() const
Returns the units of area for areal calculations made by this object.
static EllipsoidParameters ellipsoidParameters(const QString &ellipsoid)
Returns the parameters for the specified ellipsoid.
int numGeometries() const
Returns the number of geometries within the collection.
const QgsAbstractGeometry * geometryN(int n) const
Returns a const reference to a geometry from within the collection.
A geometry is the spatial representation of a feature.
QgsAbstractGeometry::const_part_iterator const_parts_begin() const
Returns STL-style const iterator pointing to the first part of the geometry.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
void convertToStraightSegment(double tolerance=M_PI/180., QgsAbstractGeometry::SegmentationToleranceType toleranceType=QgsAbstractGeometry::MaximumAngle)
Converts the geometry to straight line segments, if it is a curved geometry type.
QgsAbstractGeometry::const_part_iterator const_parts_end() const
Returns STL-style iterator pointing to the imaginary part after the last part of the geometry.
static void convertPointList(const QVector< QgsPointXY > &input, QgsPointSequence &output)
Upgrades a point list from QgsPointXY to QgsPoint.
Qgis::WkbType wkbType() const
Returns type of the geometry as a WKB type (point / linestring / polygon etc.)
Line string geometry type, with support for z-dimension and m-values.
bool isEmpty() const override
Returns true if the geometry is empty.
int numPoints() const override
Returns the number of points in the curve.
QgsPoint pointN(int i) const
Returns the specified point from inside the line string.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
static void logMessage(const QString &message, const QString &tag=QString(), Qgis::MessageLevel level=Qgis::MessageLevel::Warning, bool notifyUser=true)
Adds a message to the log instance (and creates it if necessary).
Multi surface geometry collection.
A class to represent a 2D point.
QgsPointXY project(double distance, double bearing) const
Returns a new point which corresponds to this point projected by a specified distance in a specified ...
QString toString(int precision=-1) const
Returns a string representation of the point (x, y) with a preset precision.
double distance(double x, double y) const
Returns the distance between this point and a specified x, y coordinate.
QString asWkt() const
Returns the well known text representation for the point (e.g.
void setX(double x)
Sets the x value of the point.
Point geometry type, with support for z-dimension and m-values.
void setY(double y)
Sets the point's y-coordinate.
bool addMValue(double mValue=0) override
Adds a measure to the geometry, initialized to a preset value.
void clear() override
Clears the geometry, ie reset it to a null geometry.
bool addZValue(double zValue=0) override
Adds a z-dimension to the geometry, initialized to a preset value.
void setX(double x)
Sets the point's x-coordinate.
static Q_INVOKABLE QString toString(Qgis::DistanceUnit unit)
Returns a translated string representing a distance unit.
static Q_INVOKABLE QString formatArea(double area, int decimals, Qgis::AreaUnit unit, bool keepBaseUnit=false)
Returns an area formatted as a friendly string.
static Q_INVOKABLE double fromUnitToUnitFactor(Qgis::DistanceUnit fromUnit, Qgis::DistanceUnit toUnit)
Returns the conversion factor between the specified distance units.
static Q_INVOKABLE QString formatDistance(double distance, int decimals, Qgis::DistanceUnit unit, bool keepBaseUnit=false)
Returns an distance formatted as a friendly string.
static Q_INVOKABLE Qgis::AreaUnit distanceToAreaUnit(Qgis::DistanceUnit distanceUnit)
Converts a distance unit to its corresponding area unit, e.g., meters to square meters.
static Qgis::GeometryType geometryType(Qgis::WkbType type)
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
static bool isCurvedType(Qgis::WkbType type)
Returns true if the WKB type is a curved type or can contain curved geometries.
CONSTLATIN1STRING geoNone()
Constant that holds the string representation for "No ellips/No CRS".
QString qgsDoubleToString(double a, int precision=17)
Returns a string representation of a double.
QVector< QgsPoint > QgsPointSequence
#define QgsDebugMsgLevel(str, level)
Contains parameters for an ellipsoid.
double semiMajor
Semi-major axis.
bool valid
Whether ellipsoid parameters are valid.
double semiMinor
Semi-minor axis.
QgsCoordinateReferenceSystem crs
Associated coordinate reference system.
double inverseFlattening
Inverse flattening.
bool useCustomParameters
Whether custom parameters alone should be used (semiMajor/semiMinor only)