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" );
324 p1 = mCoordTransform.
transform( points[0] );
328 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
337 geod_inverse( mGeod.get(), p1.
y(), p1.
x(), p2.
y(), p2.
x(), &distance, &azimuth1, &azimuth2 );
360 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::measureLine()",
"Error creating geod_geodesic object" );
370 QgsDebugMsgLevel( QStringLiteral(
"Ellipsoidal calculations is enabled, using ellipsoid %1" ).arg( mEllipsoid ), 4 );
375 QgsDebugMsgLevel( QStringLiteral(
"New points are %1 and %2, calculating..." ).arg( pp1.
toString( 4 ), pp2.toString( 4 ) ), 4 );
379 geod_inverse( mGeod.get(), pp1.
y(), pp1.
x(), pp2.y(), pp2.x(), &result, &azimuth1, &azimuth2 );
383 QgsDebugMsgLevel( QStringLiteral(
"Cartesian calculation on canvas coordinates" ), 4 );
408 p2 = p1.
project( distance, azimuth );
410 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]" )
411 .arg( QString::number( distance,
'f', 7 ),
413 QString::number( result,
'f', 7 ),
414 mCoordTransform.
sourceCrs().
isGeographic() ? QStringLiteral(
"Geographic" ) : QStringLiteral(
"Cartesian" ),
422 .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 );
423 if ( projectedPoint )
431 const QgsPointXY &p1,
double distance,
double azimuth )
const
441 geod_direct( mGeod.get(), p1.
y(), p1.
x(),
RAD2DEG( azimuth ), distance, &lat2, &lon2, &azimuth2 );
451 p1.
setX( p1.
x() + 360 );
453 p2.
setX( p2.
x() + 360 );
456 double p1x = p1.
x() < 180 ? p1.
x() : p2.
x();
457 double p1y = p1.
x() < 180 ? p1.
y() : p2.
y();
458 double p2x = p1.
x() < 180 ? p2.
x() : p1.
x();
459 double p2y = p1.
x() < 180 ? p2.
y() : p1.
y();
467 fractionAlongLine = ( 180 - p1x ) / ( p2x - p1x );
469 fractionAlongLine = 1 - fractionAlongLine;
470 return p1y + ( 180 - p1x ) / ( p2x - p1x ) * ( p2y - p1y );
475 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::latitudeGeodesicCrossesAntimeridian()",
"Error creating geod_geodesic object" );
479 geod_geodesicline line;
480 geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
482 const double totalDist = line.s13;
483 double intersectionDist = line.s13;
488 while ( std::fabs( lon - 180.0 ) > 0.00000001 && iterations < 100 )
490 if ( iterations > 0 && std::fabs( p2x - p1x ) > 5 )
503 QgsDebugMsgLevel( QStringLiteral(
"Narrowed window to %1, %2 - %3, %4" ).arg( p1x ).arg( p1y ).arg( p2x ).arg( p2y ), 4 );
505 geod_inverseline( &line, mGeod.get(), p1y, p1x, p2y, p2x, GEOD_ALL );
506 intersectionDist = line.s13 * 0.5;
513 intersectionDist *= ( 180.0 - p1x ) / ( lon - p1x );
518 geod_position( &line, intersectionDist, &lat, &lon, &t );
524 QgsDebugMsgLevel( QStringLiteral(
"After %1 iterations lon is %2, lat is %3, dist from p1: %4" ).arg( iterations ).arg( lon ).arg( lat ).arg( intersectionDist ), 4 );
527 fractionAlongLine = intersectionDist / totalDist;
529 fractionAlongLine = 1 - fractionAlongLine;
545 std::unique_ptr< QgsMultiLineString > res = std::make_unique< QgsMultiLineString >();
548 const QgsLineString *line = qgsgeometry_cast< const QgsLineString * >( *part );
556 const std::unique_ptr< QgsLineString > l = std::make_unique< QgsLineString >();
563 QVector< QgsPoint > newPoints;
571 for (
int i = 0; i < line->
numPoints(); i++ )
577 x = std::fmod( x, 360.0 );
588 if ( i > 0 && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
597 z = prevZ + ( p.
z() - prevZ ) * fract;
601 m = prevM + ( p.
m() - prevM ) * fract;
605 if ( prevLon < -120 )
610 QgsPoint newPoint( antiMeridianPoint );
616 if ( std::isfinite( newPoint.
x() ) && std::isfinite( newPoint.
y() ) )
618 newPoints << newPoint;
623 newPoints.reserve( line->
numPoints() - i + 1 );
630 if ( std::isfinite( antiMeridianPoint.
x() ) && std::isfinite( antiMeridianPoint.
y() ) )
634 newPoint.
setX( antiMeridianPoint.
x() );
635 newPoint.
setY( antiMeridianPoint.
y() );
636 newPoints << newPoint;
652 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform linestring. Unable to calculate break point." ) );
653 res->addGeometry( line->
clone() );
666 return QVector< QVector< QgsPointXY > >() << ( QVector< QgsPointXY >() << p1 << p2 );
672 return QVector< QVector< QgsPointXY > >();
682 QgsMessageLog::logMessage( QObject::tr(
"Caught a coordinate system exception while trying to transform a point. Unable to calculate geodesic line." ) );
683 return QVector< QVector< QgsPointXY > >();
686 geod_geodesicline line;
687 geod_inverseline( &line, mGeod.get(), pp1.y(), pp1.x(), pp2.y(), pp2.x(), GEOD_ALL );
688 const double totalDist = line.s13;
690 QVector< QVector< QgsPointXY > > res;
691 QVector< QgsPointXY > currentPart;
694 double prevLon = pp1.
x();
695 double prevLat = pp1.y();
696 bool lastRun =
false;
710 geod_position( &line, d, &lat, &lon, &t );
713 if ( breakLine && ( ( prevLon < -120 && lon > 120 ) || ( prevLon > 120 && lon < -120 ) ) )
724 if ( prevLon < -120 )
729 if ( std::isfinite( p.
x() ) && std::isfinite( p.
y() ) )
747 if ( std::isfinite( p.
x() ) && std::isfinite( p.
y() ) )
773 if ( d >= totalDist )
799 curve->
points( linePointsV2 );
800 QVector<QgsPointXY> linePoints;
810 QVector<QgsPointXY> pts;
811 for ( QVector<QgsPointXY>::const_iterator i = points.constBegin(); i != points.constEnd(); ++i )
813 pts.append( mCoordTransform.
transform( *i ) );
815 return computePolygonArea( pts );
819 return computePolygonArea( points );
836 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::bearing()",
"Error creating geod_geodesic object" );
843 geod_inverse( mGeod.get(), pp1.
y(), pp1.
x(), pp2.y(), pp2.x(), &distance, &azimuth1, &azimuth2 );
849 const double dx = p2.
x() - p1.
x();
850 const double dy = p2.
y() - p1.
y();
855 bearing = std::atan2( dx, dy );
861void QgsDistanceArea::computeAreaInit()
const
870 mGeod.reset(
new geod_geodesic() );
871 geod_init( mGeod.get(), mSemiMajor, 1 / mInvFlattening );
890double QgsDistanceArea::computePolygonArea(
const QVector<QgsPointXY> &points )
const
892 if ( points.isEmpty() )
900 return computePolygonFlatArea( points );
905 Q_ASSERT_X(
static_cast<bool>( mGeod ),
"QgsDistanceArea::computePolygonArea()",
"Error creating geod_geodesic object" );
909 struct geod_polygon p;
910 geod_polygon_init( &p, 0 );
912 const bool isClosed = points.constFirst() == points.constLast();
917 int i = points.size();
918 while ( ( isClosed && --i ) || ( !isClosed && --i >= 0 ) )
919 geod_polygon_addpoint( mGeod.get(), &p, points.at( i ).y(), points.at( i ).x() );
922 double perimeter = 0;
923 geod_polygon_compute( mGeod.get(), &p, 0, 1, &area, &perimeter );
925 return std::fabs( area );
928double QgsDistanceArea::computePolygonFlatArea(
const QVector<QgsPointXY> &points )
const
934 size = points.size();
937 for ( i = 0; i < size; i++ )
942 area = area + points[i].x() * points[( i + 1 ) % size].y() - points[( i + 1 ) % size].x() * points[i].y();
946 return std::fabs( area );
965 const double result = length * factorUnits;
966 QgsDebugMsgLevel( QStringLiteral(
"Converted length of %1 %2 to %3 %4" ).arg( length )
979 const double result = area * factorUnits;
980 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)