Branch data Line data Source code
1 : : /***************************************************************************
2 : : qgsdistancearea.h - Distance and area calculations on the ellipsoid
3 : : ---------------------------------------------------------------------------
4 : : Date : September 2005
5 : : Copyright : (C) 2005 by Martin Dobias
6 : : email : won.der at centrum.sk
7 : : ***************************************************************************
8 : : * *
9 : : * This program is free software; you can redistribute it and/or modify *
10 : : * it under the terms of the GNU General Public License as published by *
11 : : * the Free Software Foundation; either version 2 of the License, or *
12 : : * (at your option) any later version. *
13 : : * *
14 : : ***************************************************************************/
15 : :
16 : : #ifndef QGSDISTANCEAREA_H
17 : : #define QGSDISTANCEAREA_H
18 : :
19 : : #include "qgis_core.h"
20 : : #include <QVector>
21 : : #include <QReadWriteLock>
22 : : #include "qgscoordinatetransform.h"
23 : : #include "qgsunittypes.h"
24 : : #include "qgsellipsoidutils.h"
25 : :
26 : : class QgsGeometry;
27 : : class QgsAbstractGeometry;
28 : : class QgsCurve;
29 : : class geod_geodesic;
30 : :
31 : : /**
32 : : * \ingroup core
33 : : * \brief A general purpose distance and area calculator, capable of performing ellipsoid based calculations.
34 : : *
35 : : * Measurements can either be performed on existing QgsGeometry objects, or using
36 : : * lists of points.
37 : : *
38 : : * If a valid ellipsoid() has been set for the QgsDistanceArea, all calculations will be
39 : : * performed using ellipsoidal algorithms (e.g. using Vincenty's formulas). If no
40 : : * ellipsoid has been set, all calculations will be performed using Cartesian
41 : : * formulas only. The behavior can be determined by calling willUseEllipsoid().
42 : : *
43 : : * In order to perform accurate calculations, the source coordinate reference system
44 : : * of all measured geometries must first be specified using setSourceCrs().
45 : : *
46 : : * Usually, the measurements returned by QgsDistanceArea are in meters. If no valid
47 : : * ellipsoid is set, then the units may not be meters. The units can be retrieved
48 : : * by calling lengthUnits() and areaUnits().
49 : : *
50 : : * Internally, the GeographicLib library is used to calculate all ellipsoid based measurements.
51 : : */
52 : : class CORE_EXPORT QgsDistanceArea
53 : : {
54 : : public:
55 : :
56 : : //! Constructor
57 : : QgsDistanceArea();
58 : : ~QgsDistanceArea();
59 : :
60 : : //! Copy constructor
61 : : QgsDistanceArea( const QgsDistanceArea &other );
62 : : QgsDistanceArea &operator=( const QgsDistanceArea &other );
63 : :
64 : : /**
65 : : * Returns whether calculations will use the ellipsoid. Calculations will only use the
66 : : * ellipsoid if a valid ellipsoid() has been set.
67 : : * \see ellipsoid()
68 : : * \since QGIS 2.14
69 : : */
70 : : bool willUseEllipsoid() const;
71 : :
72 : : /**
73 : : * Sets source spatial reference system \a crs.
74 : : * \see sourceCrs()
75 : : * \since QGIS 2.2
76 : : */
77 : : void setSourceCrs( const QgsCoordinateReferenceSystem &crs, const QgsCoordinateTransformContext &context );
78 : :
79 : : /**
80 : : * Returns the source spatial reference system.
81 : : * \see setSourceCrs()
82 : : * \see ellipsoidCrs()
83 : : */
84 : 0 : QgsCoordinateReferenceSystem sourceCrs() const { return mCoordTransform.sourceCrs(); }
85 : :
86 : : /**
87 : : * Returns the ellipsoid (destination) spatial reference system.
88 : : * \see sourceCrs()
89 : : * \see ellipsoid()
90 : : * \since QGIS 3.6
91 : : */
92 : : QgsCoordinateReferenceSystem ellipsoidCrs() const { return mCoordTransform.destinationCrs(); }
93 : :
94 : : /**
95 : : * Sets the \a ellipsoid by its acronym. Known ellipsoid acronyms can be
96 : : * retrieved using QgsEllipsoidUtils::acronyms().
97 : : * Calculations will only use the ellipsoid if a valid ellipsoid has been set.
98 : : * \returns TRUE if ellipsoid was successfully set
99 : : * \see ellipsoid()
100 : : * \see willUseEllipsoid()
101 : : */
102 : : bool setEllipsoid( const QString &ellipsoid );
103 : :
104 : : /**
105 : : * Sets ellipsoid by supplied radii. Calculations will only use the ellipsoid if
106 : : * a valid ellipsoid been set.
107 : : * \returns TRUE if ellipsoid was successfully set
108 : : * \see ellipsoid()
109 : : * \see willUseEllipsoid()
110 : : */
111 : : bool setEllipsoid( double semiMajor, double semiMinor );
112 : :
113 : : /**
114 : : * Returns ellipsoid's acronym. Calculations will only use the
115 : : * ellipsoid if a valid ellipsoid has been set.
116 : : * \see setEllipsoid()
117 : : * \see willUseEllipsoid()
118 : : * \see ellipsoidCrs()
119 : : */
120 : : QString ellipsoid() const { return mEllipsoid; }
121 : :
122 : : /**
123 : : * Returns the ellipsoid's semi major axis.
124 : : * \see ellipsoid()
125 : : * \see ellipsoidSemiMinor()
126 : : * \see ellipsoidInverseFlattening()
127 : : */
128 : : double ellipsoidSemiMajor() const { return mSemiMajor; }
129 : :
130 : : /**
131 : : * Returns ellipsoid's semi minor axis.
132 : : * \see ellipsoid()
133 : : * \see ellipsoidSemiMajor()
134 : : * \see ellipsoidInverseFlattening()
135 : : */
136 : : double ellipsoidSemiMinor() const { return mSemiMinor; }
137 : :
138 : : /**
139 : : * Returns ellipsoid's inverse flattening.
140 : : * The inverse flattening is calculated with invf = a/(a-b).
141 : : * \see ellipsoid()
142 : : * \see ellipsoidSemiMajor()
143 : : * \see ellipsoidSemiMinor()
144 : : */
145 : : double ellipsoidInverseFlattening() const { return mInvFlattening; }
146 : :
147 : : /**
148 : : * Measures the area of a geometry.
149 : : * \param geometry geometry to measure
150 : : * \returns area of geometry. For geometry collections, non surface geometries will be ignored. The units for the
151 : : * returned area can be retrieved by calling areaUnits().
152 : : * \see measureLength()
153 : : * \see measurePerimeter()
154 : : * \see areaUnits()
155 : : * \since QGIS 2.12
156 : : */
157 : : double measureArea( const QgsGeometry &geometry ) const;
158 : :
159 : : /**
160 : : * Measures the length of a geometry.
161 : : * \param geometry geometry to measure
162 : : * \returns length of geometry. For geometry collections, non curve geometries will be ignored. The units for the
163 : : * returned distance can be retrieved by calling lengthUnits().
164 : : * \see lengthUnits()
165 : : * \see measureArea()
166 : : * \see measurePerimeter()
167 : : * \since QGIS 2.12
168 : : */
169 : : double measureLength( const QgsGeometry &geometry ) const;
170 : :
171 : : /**
172 : : * Measures the perimeter of a polygon geometry.
173 : : * \param geometry geometry to measure
174 : : * \returns perimeter of geometry. For geometry collections, any non-polygon geometries will be ignored. The units for the
175 : : * returned perimeter can be retrieved by calling lengthUnits().
176 : : * \see lengthUnits()
177 : : * \see measureArea()
178 : : * \see measurePerimeter()
179 : : * \since QGIS 2.12
180 : : */
181 : : double measurePerimeter( const QgsGeometry &geometry ) const;
182 : :
183 : : /**
184 : : * Measures the length of a line with multiple segments.
185 : : * \param points list of points in line
186 : : * \returns length of line. The units for the returned length can be retrieved by calling lengthUnits().
187 : : * \see lengthUnits()
188 : : */
189 : : double measureLine( const QVector<QgsPointXY> &points ) const;
190 : :
191 : : /**
192 : : * Measures the distance between two points.
193 : : * \param p1 start of line
194 : : * \param p2 end of line
195 : : * \returns distance between points. The units for the returned distance can be retrieved by calling lengthUnits().
196 : : * \see lengthUnits()
197 : : */
198 : : double measureLine( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
199 : :
200 : : /**
201 : : * Calculates the distance from one point with distance in meters and azimuth (direction)
202 : : * When the sourceCrs() is geographic, computeSpheroidProject() will be called
203 : : * otherwise QgsPoint.project() will be called after QgsUnitTypes::fromUnitToUnitFactor() has been applied to the distance
204 : : * \param p1 start point [can be Cartesian or Geographic]
205 : : * \param distance must be in meters
206 : : * \param azimuth - azimuth in radians, clockwise from North
207 : : * \param projectedPoint calculated projected point
208 : : * \return distance in mapUnits
209 : : * \see sourceCrs()
210 : : * \see computeSpheroidProject()
211 : : * \note The input Point must be in the coordinate reference system being used
212 : : * \since QGIS 3.0
213 : : */
214 : : double measureLineProjected( const QgsPointXY &p1, double distance = 1, double azimuth = M_PI_2, QgsPointXY *projectedPoint SIP_OUT = nullptr ) const;
215 : :
216 : : /**
217 : : * Returns the units of distance for length calculations made by this object.
218 : : * \see areaUnits()
219 : : * \since QGIS 2.14
220 : : */
221 : : QgsUnitTypes::DistanceUnit lengthUnits() const;
222 : :
223 : : /**
224 : : * Returns the units of area for areal calculations made by this object.
225 : : * \see lengthUnits()
226 : : * \since QGIS 2.14
227 : : */
228 : : QgsUnitTypes::AreaUnit areaUnits() const;
229 : :
230 : : /**
231 : : * Measures the area of the polygon described by a set of points.
232 : : */
233 : : double measurePolygon( const QVector<QgsPointXY> &points ) const;
234 : :
235 : : /**
236 : : * Computes the bearing (in radians) between two points.
237 : : */
238 : : double bearing( const QgsPointXY &p1, const QgsPointXY &p2 ) const;
239 : :
240 : : /**
241 : : * Returns an distance formatted as a friendly string.
242 : : * \param distance distance to format
243 : : * \param decimals number of decimal places to show
244 : : * \param unit unit of distance
245 : : * \param keepBaseUnit set to FALSE to allow conversion of large distances to more suitable units, e.g., meters to
246 : : * kilometers
247 : : * \returns formatted distance string
248 : : * \see formatArea()
249 : : * \since QGIS 2.16
250 : : */
251 : : static QString formatDistance( double distance, int decimals, QgsUnitTypes::DistanceUnit unit, bool keepBaseUnit = false );
252 : :
253 : : /**
254 : : * Returns an area formatted as a friendly string.
255 : : * \param area area to format
256 : : * \param decimals number of decimal places to show
257 : : * \param unit unit of area
258 : : * \param keepBaseUnit set to FALSE to allow conversion of large areas to more suitable units, e.g., square meters to
259 : : * square kilometers
260 : : * \returns formatted area string
261 : : * \see formatDistance()
262 : : * \since QGIS 2.14
263 : : */
264 : : static QString formatArea( double area, int decimals, QgsUnitTypes::AreaUnit unit, bool keepBaseUnit = false );
265 : :
266 : : /**
267 : : * Takes a length measurement calculated by this QgsDistanceArea object and converts it to a
268 : : * different distance unit.
269 : : * \param length length value calculated by this class to convert. It is assumed that the length
270 : : * was calculated by this class, ie that its unit of length is equal to lengthUnits().
271 : : * \param toUnits distance unit to convert measurement to
272 : : * \returns converted distance
273 : : * \see convertAreaMeasurement()
274 : : * \since QGIS 2.14
275 : : */
276 : : double convertLengthMeasurement( double length, QgsUnitTypes::DistanceUnit toUnits ) const;
277 : :
278 : : /**
279 : : * Takes an area measurement calculated by this QgsDistanceArea object and converts it to a
280 : : * different areal unit.
281 : : * \param area area value calculated by this class to convert. It is assumed that the area
282 : : * was calculated by this class, ie that its unit of area is equal to areaUnits().
283 : : * \param toUnits area unit to convert measurement to
284 : : * \returns converted area
285 : : * \see convertLengthMeasurement()
286 : : * \since QGIS 2.14
287 : : */
288 : : double convertAreaMeasurement( double area, QgsUnitTypes::AreaUnit toUnits ) const;
289 : :
290 : : /**
291 : : * Given a location, an azimuth and a distance, computes the
292 : : * location of the projected point.
293 : : *
294 : : * \param p1 - location of first geographic (latitude/longitude) point as degrees.
295 : : * \param distance - distance in meters.
296 : : * \param azimuth - azimuth in radians, clockwise from North
297 : : * \return p2 - location of projected point as longitude/latitude.
298 : : *
299 : : * \since QGIS 3.0
300 : : */
301 : : QgsPointXY computeSpheroidProject( const QgsPointXY &p1, double distance = 1, double azimuth = M_PI_2 ) const;
302 : :
303 : : /**
304 : : * Calculates the geodesic line between \a p1 and \a p2, which represents the shortest path on the
305 : : * ellipsoid between these two points.
306 : : *
307 : : * The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
308 : : *
309 : : * \a p1 and \a p2 must be in the sourceCrs() of this QgsDistanceArea object. The returned line
310 : : * will also be in this same CRS.
311 : : *
312 : : * The \a interval parameter gives the maximum distance between points on the computed line.
313 : : * This argument is always specified in meters. A shorter distance results in a denser line,
314 : : * at the cost of extra computing time.
315 : : *
316 : : * If the geodesic line crosses the antimeridian (+/- 180 degrees longitude) and \a breakLine is TRUE, then
317 : : * the line will be split into two parts, broken at the antimeridian. In this case the function
318 : : * will return two lines, corresponding to the portions at either side of the antimeridian.
319 : : *
320 : : * \since QGIS 3.6
321 : : */
322 : : QVector<QVector<QgsPointXY> > geodesicLine( const QgsPointXY &p1, const QgsPointXY &p2, double interval, bool breakLine = false ) const;
323 : :
324 : : /**
325 : : * Calculates the latitude at which the geodesic line joining \a p1 and \a p2 crosses
326 : : * the antimeridian (longitude +/- 180 degrees).
327 : : *
328 : : * The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
329 : : *
330 : : * \a p1 and \a p2 must be in the ellipsoidCrs() of this QgsDistanceArea object. The returned latitude
331 : : * will also be in this same CRS.
332 : : *
333 : : * \param p1 Starting point, in ellipsoidCrs()
334 : : * \param p2 Ending point, in ellipsoidCrs()
335 : : * \param fractionAlongLine will be set to the fraction along the geodesic line joining \a p1 to \a p2 at which the antimeridian crossing occurs.
336 : : *
337 : : * \returns the latitude at which the geodesic crosses the antimeridian
338 : : * \see splitGeometryAtAntimeridian()
339 : : *
340 : : * \since QGIS 3.6
341 : : */
342 : : double latitudeGeodesicCrossesAntimeridian( const QgsPointXY &p1, const QgsPointXY &p2, double &fractionAlongLine SIP_OUT ) const;
343 : :
344 : : /**
345 : : * Splits a (Multi)LineString \a geometry at the antimeridian (longitude +/- 180 degrees).
346 : : * The returned geometry will always be a multi-part geometry.
347 : : *
348 : : * Whenever line segments in the input geometry cross the antimeridian, they will be
349 : : * split into two segments, with the latitude of the breakpoint being determined using a geodesic
350 : : * line connecting the points either side of this segment.
351 : : *
352 : : * The ellipsoid settings defined on this QgsDistanceArea object will be used during the calculations.
353 : : *
354 : : * \a geometry must be in the sourceCrs() of this QgsDistanceArea object. The returned geometry
355 : : * will also be in this same CRS.
356 : : *
357 : : * If \a geometry contains M or Z values, these will be linearly interpolated for the new vertices
358 : : * created at the antimeridian.
359 : : *
360 : : * \note Non-(Multi)LineString geometries will be returned unchanged.
361 : : *
362 : : * \see latitudeGeodesicCrossesAntimeridian()
363 : : * \since QGIS 3.6
364 : : */
365 : : QgsGeometry splitGeometryAtAntimeridian( const QgsGeometry &geometry ) const;
366 : :
367 : : private:
368 : :
369 : : /**
370 : : * Calculates area of polygon on ellipsoid
371 : : * algorithm has been taken from GRASS: gis/area_poly1.c
372 : : */
373 : : double computePolygonArea( const QVector<QgsPointXY> &points ) const;
374 : :
375 : : double computePolygonFlatArea( const QVector<QgsPointXY> &points ) const;
376 : :
377 : : /**
378 : : * Precalculates some values
379 : : * (must be called always when changing ellipsoid)
380 : : */
381 : : void computeAreaInit() const;
382 : :
383 : : void setFromParams( const QgsEllipsoidUtils::EllipsoidParameters ¶ms );
384 : :
385 : : enum MeasureType
386 : : {
387 : : Default,
388 : : Area,
389 : : Length
390 : : };
391 : :
392 : : //! used for transforming coordinates from source CRS to ellipsoid's coordinates
393 : : QgsCoordinateTransform mCoordTransform;
394 : :
395 : : //! ellipsoid acronym (from table tbl_ellipsoids)
396 : : QString mEllipsoid;
397 : :
398 : : //! ellipsoid parameters
399 : : double mSemiMajor, mSemiMinor, mInvFlattening;
400 : :
401 : : mutable std::unique_ptr< geod_geodesic > mGeod;
402 : :
403 : : // utility functions for polygon area measurement
404 : :
405 : : double measure( const QgsAbstractGeometry *geomV2, MeasureType type = Default ) const;
406 : : double measureLine( const QgsCurve *curve ) const;
407 : : double measurePolygon( const QgsCurve *curve ) const;
408 : :
409 : : };
410 : :
411 : : #endif
412 : :
|
