QGIS API Documentation 3.99.0-Master (b927df884fe)
qgsgeos.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsgeos.cpp
3 -------------------------------------------------------------------
4Date : 22 Sept 2014
5Copyright : (C) 2014 by Marco Hugentobler
6email : marco.hugentobler at sourcepole dot com
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#include "qgsgeos.h"
17#include "qgsabstractgeometry.h"
19#include "qgsgeometryfactory.h"
20#include "qgslinestring.h"
21#include "qgsmulticurve.h"
22#include "qgsmultilinestring.h"
23#include "qgsmultipoint.h"
24#include "qgsmultipolygon.h"
25#include "qgslogger.h"
26#include "qgspolygon.h"
30#include <limits>
31#include <cstdio>
32
33#define DEFAULT_QUADRANT_SEGMENTS 8
34
35#define CATCH_GEOS(r) \
36 catch (QgsGeosException &) \
37 { \
38 return r; \
39 }
40
41#define CATCH_GEOS_WITH_ERRMSG(r) \
42 catch (QgsGeosException &e) \
43 { \
44 if ( errorMsg ) \
45 { \
46 *errorMsg = e.what(); \
47 } \
48 return r; \
49 }
50
52
53static void throwQgsGeosException( const char *fmt, ... )
54{
55 va_list ap;
56 char buffer[1024];
57
58 va_start( ap, fmt );
59 vsnprintf( buffer, sizeof buffer, fmt, ap );
60 va_end( ap );
61
62 QString message = QString::fromUtf8( buffer );
63
64#ifdef _MSC_VER
65 // stupid stupid MSVC, *SOMETIMES* raises it's own exception if we throw QgsGeosException, resulting in a crash!
66 // see https://github.com/qgis/QGIS/issues/22709
67 // if you want to test alternative fixes for this, run the testqgsexpression.cpp test suite - that will crash
68 // and burn on the "line_interpolate_point point" test if a QgsGeosException is thrown.
69 // TODO - find a real fix for the underlying issue
70 try
71 {
72 throw QgsGeosException( message );
73 }
74 catch ( ... )
75 {
76 // oops, msvc threw an exception when we tried to throw the exception!
77 // just throw nothing instead (except your mouse at your monitor)
78 }
79#else
80 throw QgsGeosException( message );
81#endif
82}
83
84
85static void printGEOSNotice( const char *fmt, ... )
86{
87#if defined(QGISDEBUG)
88 va_list ap;
89 char buffer[1024];
90
91 va_start( ap, fmt );
92 vsnprintf( buffer, sizeof buffer, fmt, ap );
93 va_end( ap );
94#else
95 Q_UNUSED( fmt )
96#endif
97}
98
99//
100// QgsGeosContext
101//
102
103#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
104thread_local QgsGeosContext QgsGeosContext::sGeosContext;
105#else
106QThreadStorage< QgsGeosContext * > QgsGeosContext::sGeosContext;
107#endif
108
109
111{
112#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=5 )
113 mContext = GEOS_init_r();
114 GEOSContext_setNoticeHandler_r( mContext, printGEOSNotice );
115 GEOSContext_setErrorHandler_r( mContext, throwQgsGeosException );
116#else
117 mContext = initGEOS_r( printGEOSNotice, throwQgsGeosException );
118#endif
119}
120
122{
123#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=5 )
124 GEOS_finish_r( mContext );
125#else
126 finishGEOS_r( mContext );
127#endif
128}
129
130GEOSContextHandle_t QgsGeosContext::get()
131{
132#if defined(USE_THREAD_LOCAL) && !defined(Q_OS_WIN)
133 return sGeosContext.mContext;
134#else
135 GEOSContextHandle_t gContext = nullptr;
136 if ( sGeosContext.hasLocalData() )
137 {
138 gContext = sGeosContext.localData()->mContext;
139 }
140 else
141 {
142 sGeosContext.setLocalData( new QgsGeosContext() );
143 gContext = sGeosContext.localData()->mContext;
144 }
145 return gContext;
146#endif
147}
148
149//
150// geos
151//
152
153void geos::GeosDeleter::operator()( GEOSGeometry *geom ) const
154{
155 GEOSGeom_destroy_r( QgsGeosContext::get(), geom );
156}
157
158void geos::GeosDeleter::operator()( const GEOSPreparedGeometry *geom ) const
159{
160 GEOSPreparedGeom_destroy_r( QgsGeosContext::get(), geom );
161}
162
163void geos::GeosDeleter::operator()( GEOSBufferParams *params ) const
164{
165 GEOSBufferParams_destroy_r( QgsGeosContext::get(), params );
166}
167
168void geos::GeosDeleter::operator()( GEOSCoordSequence *sequence ) const
169{
170 GEOSCoordSeq_destroy_r( QgsGeosContext::get(), sequence );
171}
172
173
175
176
178 : QgsGeometryEngine( geometry )
179 , mGeos( nullptr )
180 , mPrecision( precision )
181{
182 cacheGeos( flags );
183}
184
186{
188 GEOSGeom_destroy_r( QgsGeosContext::get(), geos );
189 return g;
190}
191
193{
194 QgsGeometry g( QgsGeos::fromGeos( geos.get() ) );
195 return g;
196}
197
198std::unique_ptr<QgsAbstractGeometry> QgsGeos::makeValid( Qgis::MakeValidMethod method, bool keepCollapsed, QString *errorMsg ) const
199{
200 if ( !mGeos )
201 {
202 return nullptr;
203 }
204
205 GEOSContextHandle_t context = QgsGeosContext::get();
206
207#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<10
208 if ( method != Qgis::MakeValidMethod::Linework )
209 throw QgsNotSupportedException( QObject::tr( "The structured method to make geometries valid requires a QGIS build based on GEOS 3.10 or later" ) );
210
211 if ( keepCollapsed )
212 throw QgsNotSupportedException( QObject::tr( "The keep collapsed option for making geometries valid requires a QGIS build based on GEOS 3.10 or later" ) );
213 geos::unique_ptr geos;
214 try
215 {
216 geos.reset( GEOSMakeValid_r( context, mGeos.get() ) );
217 }
218 CATCH_GEOS_WITH_ERRMSG( nullptr )
219#else
220
221 GEOSMakeValidParams *params = GEOSMakeValidParams_create_r( context );
222 switch ( method )
223 {
225 GEOSMakeValidParams_setMethod_r( context, params, GEOS_MAKE_VALID_LINEWORK );
226 break;
227
229 GEOSMakeValidParams_setMethod_r( context, params, GEOS_MAKE_VALID_STRUCTURE );
230 break;
231 }
232
233 GEOSMakeValidParams_setKeepCollapsed_r( context,
234 params,
235 keepCollapsed ? 1 : 0 );
236
237 geos::unique_ptr geos;
238 try
239 {
240 geos.reset( GEOSMakeValidWithParams_r( context, mGeos.get(), params ) );
241 GEOSMakeValidParams_destroy_r( context, params );
242 }
243 catch ( QgsGeosException &e )
244 {
245 if ( errorMsg )
246 {
247 *errorMsg = e.what();
248 }
249 GEOSMakeValidParams_destroy_r( context, params );
250 return nullptr;
251 }
252#endif
253
254 return fromGeos( geos.get() );
255}
256
257geos::unique_ptr QgsGeos::asGeos( const QgsGeometry &geometry, double precision, Qgis::GeosCreationFlags flags )
258{
259 if ( geometry.isNull() )
260 {
261 return nullptr;
262 }
263
264 return asGeos( geometry.constGet(), precision, flags );
265}
266
268{
269 if ( geometry.isNull() )
270 {
272 }
273 if ( !newPart )
274 {
276 }
277
278 std::unique_ptr< QgsAbstractGeometry > geom = fromGeos( newPart );
279 return QgsGeometryEditUtils::addPart( geometry.get(), std::move( geom ) );
280}
281
283{
284 mGeos.reset();
285 mGeosPrepared.reset();
287}
288
290{
291 if ( mGeosPrepared )
292 {
293 // Already prepared
294 return;
295 }
296 if ( mGeos )
297 {
298 mGeosPrepared.reset( GEOSPrepare_r( QgsGeosContext::get(), mGeos.get() ) );
299 }
300}
301
302void QgsGeos::cacheGeos( Qgis::GeosCreationFlags flags ) const
303{
304 if ( mGeos )
305 {
306 // Already cached
307 return;
308 }
309 if ( !mGeometry )
310 {
311 return;
312 }
313
314 mGeos = asGeos( mGeometry, mPrecision, flags );
315}
316
317QgsAbstractGeometry *QgsGeos::intersection( const QgsAbstractGeometry *geom, QString *errorMsg, const QgsGeometryParameters &parameters ) const
318{
319 return overlay( geom, OverlayIntersection, errorMsg, parameters ).release();
320}
321
322QgsAbstractGeometry *QgsGeos::difference( const QgsAbstractGeometry *geom, QString *errorMsg, const QgsGeometryParameters &parameters ) const
323{
324 return overlay( geom, OverlayDifference, errorMsg, parameters ).release();
325}
326
327std::unique_ptr<QgsAbstractGeometry> QgsGeos::clip( const QgsRectangle &rect, QString *errorMsg ) const
328{
329 if ( !mGeos || rect.isNull() || rect.isEmpty() )
330 {
331 return nullptr;
332 }
333
334 try
335 {
336 geos::unique_ptr opGeom( GEOSClipByRect_r( QgsGeosContext::get(), mGeos.get(), rect.xMinimum(), rect.yMinimum(), rect.xMaximum(), rect.yMaximum() ) );
337 return fromGeos( opGeom.get() );
338 }
339 catch ( QgsGeosException &e )
340 {
341 logError( QStringLiteral( "GEOS" ), e.what() );
342 if ( errorMsg )
343 {
344 *errorMsg = e.what();
345 }
346 return nullptr;
347 }
348}
349
350void QgsGeos::subdivideRecursive( const GEOSGeometry *currentPart, int maxNodes, int depth, QgsGeometryCollection *parts, const QgsRectangle &clipRect, double gridSize ) const
351{
352 GEOSContextHandle_t context = QgsGeosContext::get();
353 int partType = GEOSGeomTypeId_r( context, currentPart );
354 if ( qgsDoubleNear( clipRect.width(), 0.0 ) && qgsDoubleNear( clipRect.height(), 0.0 ) )
355 {
356 if ( partType == GEOS_POINT )
357 {
358 parts->addGeometry( fromGeos( currentPart ).release() );
359 return;
360 }
361 else
362 {
363 return;
364 }
365 }
366
367 if ( partType == GEOS_MULTILINESTRING || partType == GEOS_MULTIPOLYGON || partType == GEOS_GEOMETRYCOLLECTION )
368 {
369 int partCount = GEOSGetNumGeometries_r( context, currentPart );
370 for ( int i = 0; i < partCount; ++i )
371 {
372 subdivideRecursive( GEOSGetGeometryN_r( context, currentPart, i ), maxNodes, depth, parts, clipRect, gridSize );
373 }
374 return;
375 }
376
377 if ( depth > 50 )
378 {
379 parts->addGeometry( fromGeos( currentPart ).release() );
380 return;
381 }
382
383 int vertexCount = GEOSGetNumCoordinates_r( context, currentPart );
384 if ( vertexCount == 0 )
385 {
386 return;
387 }
388 else if ( vertexCount < maxNodes )
389 {
390 parts->addGeometry( fromGeos( currentPart ).release() );
391 return;
392 }
393
394 // chop clipping rect in half by longest side
395 double width = clipRect.width();
396 double height = clipRect.height();
397 QgsRectangle halfClipRect1 = clipRect;
398 QgsRectangle halfClipRect2 = clipRect;
399 if ( width > height )
400 {
401 halfClipRect1.setXMaximum( clipRect.xMinimum() + width / 2.0 );
402 halfClipRect2.setXMinimum( halfClipRect1.xMaximum() );
403 }
404 else
405 {
406 halfClipRect1.setYMaximum( clipRect.yMinimum() + height / 2.0 );
407 halfClipRect2.setYMinimum( halfClipRect1.yMaximum() );
408 }
409
410 if ( height <= 0 )
411 {
412 halfClipRect1.setYMinimum( halfClipRect1.yMinimum() - std::numeric_limits<double>::epsilon() );
413 halfClipRect2.setYMinimum( halfClipRect2.yMinimum() - std::numeric_limits<double>::epsilon() );
414 halfClipRect1.setYMaximum( halfClipRect1.yMaximum() + std::numeric_limits<double>::epsilon() );
415 halfClipRect2.setYMaximum( halfClipRect2.yMaximum() + std::numeric_limits<double>::epsilon() );
416 }
417 if ( width <= 0 )
418 {
419 halfClipRect1.setXMinimum( halfClipRect1.xMinimum() - std::numeric_limits<double>::epsilon() );
420 halfClipRect2.setXMinimum( halfClipRect2.xMinimum() - std::numeric_limits<double>::epsilon() );
421 halfClipRect1.setXMaximum( halfClipRect1.xMaximum() + std::numeric_limits<double>::epsilon() );
422 halfClipRect2.setXMaximum( halfClipRect2.xMaximum() + std::numeric_limits<double>::epsilon() );
423 }
424
425 geos::unique_ptr clipPart1( GEOSClipByRect_r( context, currentPart, halfClipRect1.xMinimum(), halfClipRect1.yMinimum(), halfClipRect1.xMaximum(), halfClipRect1.yMaximum() ) );
426 geos::unique_ptr clipPart2( GEOSClipByRect_r( context, currentPart, halfClipRect2.xMinimum(), halfClipRect2.yMinimum(), halfClipRect2.xMaximum(), halfClipRect2.yMaximum() ) );
427
428 ++depth;
429
430 if ( clipPart1 )
431 {
432 if ( gridSize > 0 )
433 {
434 clipPart1.reset( GEOSIntersectionPrec_r( context, mGeos.get(), clipPart1.get(), gridSize ) );
435 }
436 subdivideRecursive( clipPart1.get(), maxNodes, depth, parts, halfClipRect1, gridSize );
437 }
438 if ( clipPart2 )
439 {
440 if ( gridSize > 0 )
441 {
442 clipPart2.reset( GEOSIntersectionPrec_r( context, mGeos.get(), clipPart2.get(), gridSize ) );
443 }
444 subdivideRecursive( clipPart2.get(), maxNodes, depth, parts, halfClipRect2, gridSize );
445 }
446}
447
448std::unique_ptr<QgsAbstractGeometry> QgsGeos::subdivide( int maxNodes, QString *errorMsg, const QgsGeometryParameters &parameters ) const
449{
450 if ( !mGeos )
451 {
452 return nullptr;
453 }
454
455 // minimum allowed max is 8
456 maxNodes = std::max( maxNodes, 8 );
457
458 std::unique_ptr< QgsGeometryCollection > parts = QgsGeometryFactory::createCollectionOfType( mGeometry->wkbType() );
459 try
460 {
461 subdivideRecursive( mGeos.get(), maxNodes, 0, parts.get(), mGeometry->boundingBox(), parameters.gridSize() );
462 }
463 CATCH_GEOS_WITH_ERRMSG( nullptr )
464
465 return std::move( parts );
466}
467
468QgsAbstractGeometry *QgsGeos::combine( const QgsAbstractGeometry *geom, QString *errorMsg, const QgsGeometryParameters &parameters ) const
469{
470 return overlay( geom, OverlayUnion, errorMsg, parameters ).release();
471}
472
473QgsAbstractGeometry *QgsGeos::combine( const QVector<QgsAbstractGeometry *> &geomList, QString *errorMsg, const QgsGeometryParameters &parameters ) const
474{
475 std::vector<geos::unique_ptr> geosGeometries;
476 geosGeometries.reserve( geomList.size() );
477 for ( const QgsAbstractGeometry *g : geomList )
478 {
479 if ( !g )
480 continue;
481
482 geosGeometries.emplace_back( asGeos( g, mPrecision ) );
483 }
484
485 GEOSContextHandle_t context = QgsGeosContext::get();
486 geos::unique_ptr geomUnion;
487 try
488 {
489 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
490 if ( parameters.gridSize() > 0 )
491 {
492 geomUnion.reset( GEOSUnaryUnionPrec_r( context, geomCollection.get(), parameters.gridSize() ) );
493 }
494 else
495 {
496 geomUnion.reset( GEOSUnaryUnion_r( context, geomCollection.get() ) );
497 }
498 }
499 CATCH_GEOS_WITH_ERRMSG( nullptr )
500
501 std::unique_ptr< QgsAbstractGeometry > result = fromGeos( geomUnion.get() );
502 return result.release();
503}
504
505QgsAbstractGeometry *QgsGeos::combine( const QVector<QgsGeometry> &geomList, QString *errorMsg, const QgsGeometryParameters &parameters ) const
506{
507 std::vector<geos::unique_ptr> geosGeometries;
508 geosGeometries.reserve( geomList.size() );
509 for ( const QgsGeometry &g : geomList )
510 {
511 if ( g.isNull() )
512 continue;
513
514 geosGeometries.emplace_back( asGeos( g.constGet(), mPrecision ) );
515 }
516
517 GEOSContextHandle_t context = QgsGeosContext::get();
518 geos::unique_ptr geomUnion;
519 try
520 {
521 geos::unique_ptr geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geosGeometries );
522
523 if ( parameters.gridSize() > 0 )
524 {
525 geomUnion.reset( GEOSUnaryUnionPrec_r( context, geomCollection.get(), parameters.gridSize() ) );
526 }
527 else
528 {
529 geomUnion.reset( GEOSUnaryUnion_r( context, geomCollection.get() ) );
530 }
531
532 }
533 CATCH_GEOS_WITH_ERRMSG( nullptr )
534
535 std::unique_ptr< QgsAbstractGeometry > result = fromGeos( geomUnion.get() );
536 return result.release();
537}
538
539QgsAbstractGeometry *QgsGeos::symDifference( const QgsAbstractGeometry *geom, QString *errorMsg, const QgsGeometryParameters &parameters ) const
540{
541 return overlay( geom, OverlaySymDifference, errorMsg, parameters ).release();
542}
543
544static bool isZVerticalLine( const QgsAbstractGeometry *geom, double tolerance = 4 * std::numeric_limits<double>::epsilon() )
545{
546 // checks if the Geometry if a purely vertical 3D line LineString Z((X Y Z1, X Y Z2, ..., X Y Zn))
547 // This is needed because QgsGeos is not able to handle this type of geometry on distance computation.
548
550 {
551 return false;
552 }
553
554 bool isVertical = true;
555 if ( const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( geom ) )
556 {
557 const int nrPoints = line->numPoints();
558 if ( nrPoints == 1 )
559 {
560 return true;
561 }
562
563 // if the 2D part of two points of the line are different, this means
564 // that the line is not purely vertical
565 const double sqrTolerance = tolerance * tolerance;
566 const double *lineX = line->xData();
567 const double *lineY = line->yData();
568 for ( int iVert = nrPoints - 1, jVert = 0; jVert < nrPoints; iVert = jVert++ )
569 {
570 if ( QgsGeometryUtilsBase::sqrDistance2D( lineX[iVert], lineY[iVert], lineX[jVert], lineY[jVert] ) > sqrTolerance )
571 {
572 isVertical = false;
573 break;
574 }
575 }
576 }
577
578 return isVertical;
579}
580
581double QgsGeos::distance( const QgsAbstractGeometry *geom, QString *errorMsg ) const
582{
583 double distance = -1.0;
584 if ( !mGeos )
585 {
586 return distance;
587 }
588
589 geos::unique_ptr otherGeosGeom;
590
591 // GEOSPreparedDistance_r is not able to properly compute the distance if one
592 // of the geometries if a vertical line (LineString Z((X Y Z1, X Y Z2, ..., X Y Zn))).
593 // In that case, replace `geom` by a single point.
594 // However, GEOSDistance_r works.
595 if ( mGeosPrepared && isZVerticalLine( geom->simplifiedTypeRef() ) )
596 {
597 QgsPoint firstPoint = geom->vertexAt( QgsVertexId( 0, 0, 0 ) );
598 otherGeosGeom = asGeos( &firstPoint, mPrecision );
599 }
600 else
601 {
602 otherGeosGeom = asGeos( geom, mPrecision );
603 }
604
605 if ( !otherGeosGeom )
606 {
607 return distance;
608 }
609
610 GEOSContextHandle_t context = QgsGeosContext::get();
611 try
612 {
613 if ( mGeosPrepared && !isZVerticalLine( mGeometry->simplifiedTypeRef() ) )
614 {
615 GEOSPreparedDistance_r( context, mGeosPrepared.get(), otherGeosGeom.get(), &distance );
616 }
617 else
618 {
619 GEOSDistance_r( context, mGeos.get(), otherGeosGeom.get(), &distance );
620 }
621 }
623
624 return distance;
625}
626
627double QgsGeos::distance( double x, double y, QString *errorMsg ) const
628{
629 double distance = -1.0;
630 if ( !mGeos )
631 {
632 return distance;
633 }
634
635 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
636 if ( !point )
637 return distance;
638
639 GEOSContextHandle_t context = QgsGeosContext::get();
640 try
641 {
642 if ( mGeosPrepared )
643 {
644 GEOSPreparedDistance_r( context, mGeosPrepared.get(), point.get(), &distance );
645 }
646 else
647 {
648 GEOSDistance_r( context, mGeos.get(), point.get(), &distance );
649 }
650 }
652
653 return distance;
654}
655
656bool QgsGeos::distanceWithin( const QgsAbstractGeometry *geom, double maxdist, QString *errorMsg ) const
657{
658 if ( !mGeos )
659 {
660 return false;
661 }
662
663 geos::unique_ptr otherGeosGeom;
664
665 // GEOSPreparedDistanceWithin_r GEOSPreparedDistance_r are not able to properly compute the distance if one
666 // of the geometries if a vertical line (LineString Z((X Y Z1, X Y Z2, ..., X Y Zn))).
667 // In that case, replace `geom` by a single point.
668 // However, GEOSDistanceWithin_r and GEOSDistance_r work.
669 if ( mGeosPrepared && isZVerticalLine( geom->simplifiedTypeRef() ) )
670 {
671 QgsPoint firstPoint = geom->vertexAt( QgsVertexId( 0, 0, 0 ) );
672 otherGeosGeom = asGeos( &firstPoint );
673 }
674 else
675 {
676 otherGeosGeom = asGeos( geom, mPrecision );
677 }
678
679 if ( !otherGeosGeom )
680 {
681 return false;
682 }
683
684 // TODO: optimize implementation of this function to early-exit if
685 // any part of othergeosGeom is found to be within the given
686 // distance
687 double distance;
688
689 GEOSContextHandle_t context = QgsGeosContext::get();
690 try
691 {
692 if ( mGeosPrepared && !isZVerticalLine( mGeometry->simplifiedTypeRef() ) )
693 {
694#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
695 return GEOSPreparedDistanceWithin_r( context, mGeosPrepared.get(), otherGeosGeom.get(), maxdist );
696#else
697 GEOSPreparedDistance_r( context, mGeosPrepared.get(), otherGeosGeom.get(), &distance );
698#endif
699 }
700 else
701 {
702#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
703 return GEOSDistanceWithin_r( context, mGeos.get(), otherGeosGeom.get(), maxdist );
704#else
705 GEOSDistance_r( context, mGeos.get(), otherGeosGeom.get(), &distance );
706#endif
707 }
708 }
710
711 return distance <= maxdist;
712}
713
714bool QgsGeos::contains( double x, double y, QString *errorMsg ) const
715{
716 bool result = false;
717 GEOSContextHandle_t context = QgsGeosContext::get();
718 try
719 {
720#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
721 // defer point creation until after prepared geometry check, we may not need it
722#else
723 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
724 if ( !point )
725 return false;
726#endif
727 if ( mGeosPrepared ) //use faster version with prepared geometry
728 {
729#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
730 return GEOSPreparedContainsXY_r( context, mGeosPrepared.get(), x, y ) == 1;
731#else
732 return GEOSPreparedContains_r( context, mGeosPrepared.get(), point.get() ) == 1;
733#endif
734 }
735
736#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
737 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
738 if ( !point )
739 return false;
740#endif
741
742 result = ( GEOSContains_r( context, mGeos.get(), point.get() ) == 1 );
743 }
744 catch ( QgsGeosException &e )
745 {
746 logError( QStringLiteral( "GEOS" ), e.what() );
747 if ( errorMsg )
748 {
749 *errorMsg = e.what();
750 }
751 return false;
752 }
753
754 return result;
755}
756
757double QgsGeos::hausdorffDistance( const QgsAbstractGeometry *geom, QString *errorMsg ) const
758{
759 double distance = -1.0;
760 if ( !mGeos )
761 {
762 return distance;
763 }
764
765 geos::unique_ptr otherGeosGeom( asGeos( geom, mPrecision ) );
766 if ( !otherGeosGeom )
767 {
768 return distance;
769 }
770
771 try
772 {
773 GEOSHausdorffDistance_r( QgsGeosContext::get(), mGeos.get(), otherGeosGeom.get(), &distance );
774 }
776
777 return distance;
778}
779
780double QgsGeos::hausdorffDistanceDensify( const QgsAbstractGeometry *geom, double densifyFraction, QString *errorMsg ) const
781{
782 double distance = -1.0;
783 if ( !mGeos )
784 {
785 return distance;
786 }
787
788 geos::unique_ptr otherGeosGeom( asGeos( geom, mPrecision ) );
789 if ( !otherGeosGeom )
790 {
791 return distance;
792 }
793
794 try
795 {
796 GEOSHausdorffDistanceDensify_r( QgsGeosContext::get(), mGeos.get(), otherGeosGeom.get(), densifyFraction, &distance );
797 }
799
800 return distance;
801}
802
803double QgsGeos::frechetDistance( const QgsAbstractGeometry *geom, QString *errorMsg ) const
804{
805 double distance = -1.0;
806 if ( !mGeos )
807 {
808 return distance;
809 }
810
811 geos::unique_ptr otherGeosGeom( asGeos( geom, mPrecision ) );
812 if ( !otherGeosGeom )
813 {
814 return distance;
815 }
816
817 try
818 {
819 GEOSFrechetDistance_r( QgsGeosContext::get(), mGeos.get(), otherGeosGeom.get(), &distance );
820 }
822
823 return distance;
824}
825
826double QgsGeos::frechetDistanceDensify( const QgsAbstractGeometry *geom, double densifyFraction, QString *errorMsg ) const
827{
828 double distance = -1.0;
829 if ( !mGeos )
830 {
831 return distance;
832 }
833
834 geos::unique_ptr otherGeosGeom( asGeos( geom, mPrecision ) );
835 if ( !otherGeosGeom )
836 {
837 return distance;
838 }
839
840 try
841 {
842 GEOSFrechetDistanceDensify_r( QgsGeosContext::get(), mGeos.get(), otherGeosGeom.get(), densifyFraction, &distance );
843 }
845
846 return distance;
847}
848
849bool QgsGeos::intersects( const QgsAbstractGeometry *geom, QString *errorMsg ) const
850{
851 if ( !mGeos || !geom )
852 {
853 return false;
854 }
855
856#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
857 // special optimised case for point intersects
858 if ( const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( geom->simplifiedTypeRef() ) )
859 {
860 if ( mGeosPrepared )
861 {
862 try
863 {
864 return GEOSPreparedIntersectsXY_r( QgsGeosContext::get(), mGeosPrepared.get(), point->x(), point->y() ) == 1;
865 }
866 catch ( QgsGeosException &e )
867 {
868 logError( QStringLiteral( "GEOS" ), e.what() );
869 if ( errorMsg )
870 {
871 *errorMsg = e.what();
872 }
873 return false;
874 }
875 }
876 }
877#endif
878
879 return relation( geom, RelationIntersects, errorMsg );
880}
881
882bool QgsGeos::touches( const QgsAbstractGeometry *geom, QString *errorMsg ) const
883{
884 return relation( geom, RelationTouches, errorMsg );
885}
886
887bool QgsGeos::crosses( const QgsAbstractGeometry *geom, QString *errorMsg ) const
888{
889 return relation( geom, RelationCrosses, errorMsg );
890}
891
892bool QgsGeos::within( const QgsAbstractGeometry *geom, QString *errorMsg ) const
893{
894 return relation( geom, RelationWithin, errorMsg );
895}
896
897bool QgsGeos::overlaps( const QgsAbstractGeometry *geom, QString *errorMsg ) const
898{
899 return relation( geom, RelationOverlaps, errorMsg );
900}
901
902bool QgsGeos::contains( const QgsAbstractGeometry *geom, QString *errorMsg ) const
903{
904 if ( !mGeos || !geom )
905 {
906 return false;
907 }
908
909#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=12 )
910 // special optimised case for point containment
911 if ( const QgsPoint *point = qgsgeometry_cast< const QgsPoint * >( geom->simplifiedTypeRef() ) )
912 {
913 if ( mGeosPrepared )
914 {
915 try
916 {
917 return GEOSPreparedContainsXY_r( QgsGeosContext::get(), mGeosPrepared.get(), point->x(), point->y() ) == 1;
918 }
919 catch ( QgsGeosException &e )
920 {
921 logError( QStringLiteral( "GEOS" ), e.what() );
922 if ( errorMsg )
923 {
924 *errorMsg = e.what();
925 }
926 return false;
927 }
928 }
929 }
930#endif
931
932 return relation( geom, RelationContains, errorMsg );
933}
934
935bool QgsGeos::disjoint( const QgsAbstractGeometry *geom, QString *errorMsg ) const
936{
937 return relation( geom, RelationDisjoint, errorMsg );
938}
939
940QString QgsGeos::relate( const QgsAbstractGeometry *geom, QString *errorMsg ) const
941{
942 if ( !mGeos )
943 {
944 return QString();
945 }
946
947 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
948 if ( !geosGeom )
949 {
950 return QString();
951 }
952
953 QString result;
954 GEOSContextHandle_t context = QgsGeosContext::get();
955 try
956 {
957 char *r = GEOSRelate_r( context, mGeos.get(), geosGeom.get() );
958 if ( r )
959 {
960 result = QString( r );
961 GEOSFree_r( context, r );
962 }
963 }
964 catch ( QgsGeosException &e )
965 {
966 logError( QStringLiteral( "GEOS" ), e.what() );
967 if ( errorMsg )
968 {
969 *errorMsg = e.what();
970 }
971 }
972
973 return result;
974}
975
976bool QgsGeos::relatePattern( const QgsAbstractGeometry *geom, const QString &pattern, QString *errorMsg ) const
977{
978 if ( !mGeos || !geom )
979 {
980 return false;
981 }
982
983 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
984 if ( !geosGeom )
985 {
986 return false;
987 }
988
989 bool result = false;
990 GEOSContextHandle_t context = QgsGeosContext::get();
991 try
992 {
993 result = ( GEOSRelatePattern_r( context, mGeos.get(), geosGeom.get(), pattern.toLocal8Bit().constData() ) == 1 );
994 }
995 catch ( QgsGeosException &e )
996 {
997 logError( QStringLiteral( "GEOS" ), e.what() );
998 if ( errorMsg )
999 {
1000 *errorMsg = e.what();
1001 }
1002 }
1003
1004 return result;
1005}
1006
1007double QgsGeos::area( QString *errorMsg ) const
1008{
1009 double area = -1.0;
1010 if ( !mGeos )
1011 {
1012 return area;
1013 }
1014
1015 try
1016 {
1017 if ( GEOSArea_r( QgsGeosContext::get(), mGeos.get(), &area ) != 1 )
1018 return -1.0;
1019 }
1021 return area;
1022}
1023
1024double QgsGeos::length( QString *errorMsg ) const
1025{
1026 double length = -1.0;
1027 if ( !mGeos )
1028 {
1029 return length;
1030 }
1031 try
1032 {
1033 if ( GEOSLength_r( QgsGeosContext::get(), mGeos.get(), &length ) != 1 )
1034 return -1.0;
1035 }
1037 return length;
1038}
1039
1041 QVector<QgsGeometry> &newGeometries,
1042 bool topological,
1043 QgsPointSequence &topologyTestPoints,
1044 QString *errorMsg, bool skipIntersectionCheck ) const
1045{
1046
1047 EngineOperationResult returnCode = Success;
1048 if ( !mGeos || !mGeometry )
1049 {
1050 return InvalidBaseGeometry;
1051 }
1052
1053 //return if this type is point/multipoint
1054 if ( mGeometry->dimension() == 0 )
1055 {
1056 return SplitCannotSplitPoint; //cannot split points
1057 }
1058
1059 GEOSContextHandle_t context = QgsGeosContext::get();
1060 if ( !GEOSisValid_r( context, mGeos.get() ) )
1061 return InvalidBaseGeometry;
1062
1063 //make sure splitLine is valid
1064 if ( ( mGeometry->dimension() == 1 && splitLine.numPoints() < 1 ) ||
1065 ( mGeometry->dimension() == 2 && splitLine.numPoints() < 2 ) )
1066 return InvalidInput;
1067
1068 newGeometries.clear();
1069 geos::unique_ptr splitLineGeos;
1070
1071 try
1072 {
1073 if ( splitLine.numPoints() > 1 )
1074 {
1075 splitLineGeos = createGeosLinestring( &splitLine, mPrecision );
1076 }
1077 else if ( splitLine.numPoints() == 1 )
1078 {
1079 splitLineGeos = createGeosPointXY( splitLine.xAt( 0 ), splitLine.yAt( 0 ), false, 0, false, 0, 2, mPrecision );
1080 }
1081 else
1082 {
1083 return InvalidInput;
1084 }
1085
1086 if ( !GEOSisValid_r( context, splitLineGeos.get() ) || !GEOSisSimple_r( context, splitLineGeos.get() ) )
1087 {
1088 return InvalidInput;
1089 }
1090
1091 if ( topological )
1092 {
1093 //find out candidate points for topological corrections
1094 if ( !topologicalTestPointsSplit( splitLineGeos.get(), topologyTestPoints ) )
1095 {
1096 return InvalidInput; // TODO: is it really an invalid input?
1097 }
1098 }
1099
1100 //call split function depending on geometry type
1101 if ( mGeometry->dimension() == 1 )
1102 {
1103 returnCode = splitLinearGeometry( splitLineGeos.get(), newGeometries, skipIntersectionCheck );
1104 }
1105 else if ( mGeometry->dimension() == 2 )
1106 {
1107 returnCode = splitPolygonGeometry( splitLineGeos.get(), newGeometries, skipIntersectionCheck );
1108 }
1109 else
1110 {
1111 return InvalidInput;
1112 }
1113 }
1115
1116 return returnCode;
1117}
1118
1119
1120
1121bool QgsGeos::topologicalTestPointsSplit( const GEOSGeometry *splitLine, QgsPointSequence &testPoints, QString *errorMsg ) const
1122{
1123 //Find out the intersection points between splitLineGeos and this geometry.
1124 //These points need to be tested for topological correctness by the calling function
1125 //if topological editing is enabled
1126
1127 if ( !mGeos )
1128 {
1129 return false;
1130 }
1131
1132 GEOSContextHandle_t context = QgsGeosContext::get();
1133 try
1134 {
1135 testPoints.clear();
1136 geos::unique_ptr intersectionGeom( GEOSIntersection_r( context, mGeos.get(), splitLine ) );
1137 if ( !intersectionGeom )
1138 return false;
1139
1140 bool simple = false;
1141 int nIntersectGeoms = 1;
1142 if ( GEOSGeomTypeId_r( context, intersectionGeom.get() ) == GEOS_LINESTRING
1143 || GEOSGeomTypeId_r( context, intersectionGeom.get() ) == GEOS_POINT )
1144 simple = true;
1145
1146 if ( !simple )
1147 nIntersectGeoms = GEOSGetNumGeometries_r( context, intersectionGeom.get() );
1148
1149 for ( int i = 0; i < nIntersectGeoms; ++i )
1150 {
1151 const GEOSGeometry *currentIntersectGeom = nullptr;
1152 if ( simple )
1153 currentIntersectGeom = intersectionGeom.get();
1154 else
1155 currentIntersectGeom = GEOSGetGeometryN_r( context, intersectionGeom.get(), i );
1156
1157 const GEOSCoordSequence *lineSequence = GEOSGeom_getCoordSeq_r( context, currentIntersectGeom );
1158 unsigned int sequenceSize = 0;
1159 double x, y, z;
1160 if ( GEOSCoordSeq_getSize_r( context, lineSequence, &sequenceSize ) != 0 )
1161 {
1162 for ( unsigned int i = 0; i < sequenceSize; ++i )
1163 {
1164 if ( GEOSCoordSeq_getXYZ_r( context, lineSequence, i, &x, &y, &z ) )
1165 {
1166 testPoints.push_back( QgsPoint( x, y, z ) );
1167 }
1168 }
1169 }
1170 }
1171 }
1173
1174 return true;
1175}
1176
1177geos::unique_ptr QgsGeos::linePointDifference( GEOSGeometry *GEOSsplitPoint ) const
1178{
1179 GEOSContextHandle_t context = QgsGeosContext::get();
1180 int type = GEOSGeomTypeId_r( context, mGeos.get() );
1181
1182 std::unique_ptr< QgsMultiCurve > multiCurve;
1183 if ( type == GEOS_MULTILINESTRING )
1184 {
1185 multiCurve.reset( qgsgeometry_cast<QgsMultiCurve *>( mGeometry->clone() ) );
1186 }
1187 else if ( type == GEOS_LINESTRING )
1188 {
1189 multiCurve.reset( new QgsMultiCurve() );
1190 multiCurve->addGeometry( mGeometry->clone() );
1191 }
1192 else
1193 {
1194 return nullptr;
1195 }
1196
1197 if ( !multiCurve )
1198 {
1199 return nullptr;
1200 }
1201
1202
1203 // we might have a point or a multipoint, depending on number of
1204 // intersections between the geometry and the split geometry
1205 std::unique_ptr< QgsMultiPoint > splitPoints;
1206 {
1207 std::unique_ptr< QgsAbstractGeometry > splitGeom( fromGeos( GEOSsplitPoint ) );
1208
1209 if ( qgsgeometry_cast<QgsMultiPoint *>( splitGeom.get() ) )
1210 {
1211 splitPoints.reset( qgsgeometry_cast<QgsMultiPoint *>( splitGeom.release() ) );
1212 }
1213 else if ( qgsgeometry_cast<QgsPoint *>( splitGeom.get() ) )
1214 {
1215 splitPoints = std::make_unique< QgsMultiPoint >();
1216 if ( qgsgeometry_cast<QgsPoint *>( splitGeom.get() ) )
1217 {
1218 splitPoints->addGeometry( qgsgeometry_cast<QgsPoint *>( splitGeom.release() ) );
1219 }
1220 }
1221 }
1222
1223 QgsMultiCurve lines;
1224
1225 //For each part
1226 for ( int geometryIndex = 0; geometryIndex < multiCurve->numGeometries(); ++geometryIndex )
1227 {
1228 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( multiCurve->geometryN( geometryIndex ) );
1229 if ( !line )
1230 {
1231 const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( multiCurve->geometryN( geometryIndex ) );
1232 line = curve->curveToLine();
1233 }
1234 if ( !line )
1235 {
1236 return nullptr;
1237 }
1238 // we gather the intersection points and their distance from previous node grouped by segment
1239 QMap< int, QVector< QPair< double, QgsPoint > > >pointMap;
1240 for ( int splitPointIndex = 0; splitPointIndex < splitPoints->numGeometries(); ++splitPointIndex )
1241 {
1242 const QgsPoint *intersectionPoint = splitPoints->pointN( splitPointIndex );
1243
1244 QgsPoint segmentPoint2D;
1245 QgsVertexId nextVertex;
1246 // With closestSegment we only get a 2D point so we need to interpolate if we
1247 // don't want to lose Z data
1248 line->closestSegment( *intersectionPoint, segmentPoint2D, nextVertex );
1249
1250 // The intersection might belong to another part, skip it
1251 // Note: cannot test for equality because of Z
1252 if ( !qgsDoubleNear( intersectionPoint->x(), segmentPoint2D.x() ) || !qgsDoubleNear( intersectionPoint->y(), segmentPoint2D.y() ) )
1253 {
1254 continue;
1255 }
1256
1257 const QgsLineString segment = QgsLineString( line->pointN( nextVertex.vertex - 1 ), line->pointN( nextVertex.vertex ) );
1258 const double distance = segmentPoint2D.distance( line->pointN( nextVertex.vertex - 1 ) );
1259
1260 // Due to precision issues, distance can be a tad larger than the actual segment length, making interpolatePoint() return nullptr
1261 // In that case we'll use the segment's endpoint instead of interpolating
1262 std::unique_ptr< QgsPoint > correctSegmentPoint( distance > segment.length() ? segment.endPoint().clone() : segment.interpolatePoint( distance ) );
1263
1264 const QPair< double, QgsPoint > pair = qMakePair( distance, *correctSegmentPoint.get() );
1265 if ( pointMap.contains( nextVertex.vertex - 1 ) )
1266 pointMap[ nextVertex.vertex - 1 ].append( pair );
1267 else
1268 pointMap[ nextVertex.vertex - 1 ] = QVector< QPair< double, QgsPoint > >() << pair;
1269 }
1270
1271 // When we have more than one intersection point on a segment we need those points
1272 // to be sorted by their distance from the previous geometry vertex
1273 for ( auto &p : pointMap )
1274 {
1275 std::sort( p.begin(), p.end(), []( const QPair< double, QgsPoint > &a, const QPair< double, QgsPoint > &b ) { return a.first < b.first; } );
1276 }
1277
1278 //For each segment
1279 QgsLineString newLine;
1280 int nVertices = line->numPoints();
1281 QgsPoint splitPoint;
1282 for ( int vertexIndex = 0; vertexIndex < nVertices; ++vertexIndex )
1283 {
1284 QgsPoint currentPoint = line->pointN( vertexIndex );
1285 newLine.addVertex( currentPoint );
1286 if ( pointMap.contains( vertexIndex ) )
1287 {
1288 // For each intersecting point
1289 for ( int k = 0; k < pointMap[ vertexIndex ].size(); ++k )
1290 {
1291 splitPoint = pointMap[ vertexIndex ][k].second;
1292 if ( splitPoint == currentPoint )
1293 {
1294 lines.addGeometry( newLine.clone() );
1295 newLine = QgsLineString();
1296 newLine.addVertex( currentPoint );
1297 }
1298 else if ( splitPoint == line->pointN( vertexIndex + 1 ) )
1299 {
1300 newLine.addVertex( line->pointN( vertexIndex + 1 ) );
1301 lines.addGeometry( newLine.clone() );
1302 newLine = QgsLineString();
1303 }
1304 else
1305 {
1306 newLine.addVertex( splitPoint );
1307 lines.addGeometry( newLine.clone() );
1308 newLine = QgsLineString();
1309 newLine.addVertex( splitPoint );
1310 }
1311 }
1312 }
1313 }
1314 lines.addGeometry( newLine.clone() );
1315 }
1316
1317 return asGeos( &lines, mPrecision );
1318}
1319
1320QgsGeometryEngine::EngineOperationResult QgsGeos::splitLinearGeometry( const GEOSGeometry *splitLine, QVector<QgsGeometry> &newGeometries, bool skipIntersectionCheck ) const
1321{
1322 Q_UNUSED( skipIntersectionCheck )
1323 if ( !splitLine )
1324 return InvalidInput;
1325
1326 if ( !mGeos )
1327 return InvalidBaseGeometry;
1328
1329 GEOSContextHandle_t context = QgsGeosContext::get();
1330
1331 geos::unique_ptr intersectGeom( GEOSIntersection_r( context, splitLine, mGeos.get() ) );
1332 if ( !intersectGeom || GEOSisEmpty_r( context, intersectGeom.get() ) )
1333 return NothingHappened;
1334
1335 //check that split line has no linear intersection
1336 const int linearIntersect = GEOSRelatePattern_r( context, mGeos.get(), splitLine, "1********" );
1337 if ( linearIntersect > 0 )
1338 return InvalidInput;
1339
1340 geos::unique_ptr splitGeom = linePointDifference( intersectGeom.get() );
1341
1342 if ( !splitGeom )
1343 return InvalidBaseGeometry;
1344
1345 std::vector<geos::unique_ptr> lineGeoms;
1346
1347 const int splitType = GEOSGeomTypeId_r( context, splitGeom.get() );
1348 if ( splitType == GEOS_MULTILINESTRING )
1349 {
1350 const int nGeoms = GEOSGetNumGeometries_r( context, splitGeom.get() );
1351 lineGeoms.reserve( nGeoms );
1352 for ( int i = 0; i < nGeoms; ++i )
1353 lineGeoms.emplace_back( GEOSGeom_clone_r( context, GEOSGetGeometryN_r( context, splitGeom.get(), i ) ) );
1354
1355 }
1356 else
1357 {
1358 lineGeoms.emplace_back( GEOSGeom_clone_r( context, splitGeom.get() ) );
1359 }
1360
1361 mergeGeometriesMultiTypeSplit( lineGeoms );
1362
1363 for ( geos::unique_ptr &lineGeom : lineGeoms )
1364 {
1365 newGeometries << QgsGeometry( fromGeos( lineGeom.get() ) );
1366 }
1367
1368 return Success;
1369}
1370
1371QgsGeometryEngine::EngineOperationResult QgsGeos::splitPolygonGeometry( const GEOSGeometry *splitLine, QVector<QgsGeometry> &newGeometries, bool skipIntersectionCheck ) const
1372{
1373 if ( !splitLine )
1374 return InvalidInput;
1375
1376 if ( !mGeos )
1377 return InvalidBaseGeometry;
1378
1379 // we will need prepared geometry for intersection tests
1380 const_cast<QgsGeos *>( this )->prepareGeometry();
1381 if ( !mGeosPrepared )
1382 return EngineError;
1383
1384 GEOSContextHandle_t context = QgsGeosContext::get();
1385
1386 //first test if linestring intersects geometry. If not, return straight away
1387 if ( !skipIntersectionCheck && !GEOSPreparedIntersects_r( context, mGeosPrepared.get(), splitLine ) )
1388 return NothingHappened;
1389
1390 //first union all the polygon rings together (to get them noded, see JTS developer guide)
1391 geos::unique_ptr nodedGeometry = nodeGeometries( splitLine, mGeos.get() );
1392 if ( !nodedGeometry )
1393 return NodedGeometryError; //an error occurred during noding
1394
1395 const GEOSGeometry *noded = nodedGeometry.get();
1396 geos::unique_ptr polygons( GEOSPolygonize_r( context, &noded, 1 ) );
1397 if ( !polygons )
1398 {
1399 return InvalidBaseGeometry;
1400 }
1401 const int numberOfGeometriesPolygon = numberOfGeometries( polygons.get() );
1402 if ( numberOfGeometriesPolygon == 0 )
1403 {
1404 return InvalidBaseGeometry;
1405 }
1406
1407 //test every polygon is contained in original geometry
1408 //include in result if yes
1409 std::vector<geos::unique_ptr> testedGeometries;
1410
1411 // test whether the polygon parts returned from polygonize algorithm actually
1412 // belong to the source polygon geometry (if the source polygon contains some holes,
1413 // those would be also returned by polygonize and we need to skip them)
1414 for ( int i = 0; i < numberOfGeometriesPolygon; i++ )
1415 {
1416 const GEOSGeometry *polygon = GEOSGetGeometryN_r( context, polygons.get(), i );
1417
1418 geos::unique_ptr pointOnSurface( GEOSPointOnSurface_r( context, polygon ) );
1419 if ( pointOnSurface && GEOSPreparedIntersects_r( context, mGeosPrepared.get(), pointOnSurface.get() ) )
1420 testedGeometries.emplace_back( GEOSGeom_clone_r( context, polygon ) );
1421 }
1422
1423 const size_t nGeometriesThis = numberOfGeometries( mGeos.get() ); //original number of geometries
1424 if ( testedGeometries.empty() || testedGeometries.size() == nGeometriesThis )
1425 {
1426 //no split done, preserve original geometry
1427 return NothingHappened;
1428 }
1429
1430 // For multi-part geometries, try to identify parts that have been unchanged and try to merge them back
1431 // to a single multi-part geometry. For example, if there's a multi-polygon with three parts, but only
1432 // one part is being split, this function makes sure that the other two parts will be kept in a multi-part
1433 // geometry rather than being separated into two single-part geometries.
1434 mergeGeometriesMultiTypeSplit( testedGeometries );
1435
1436 size_t i;
1437 for ( i = 0; i < testedGeometries.size() && GEOSisValid_r( context, testedGeometries[i].get() ); ++i )
1438 ;
1439
1440 if ( i < testedGeometries.size() )
1441 {
1442 return InvalidBaseGeometry;
1443 }
1444
1445 for ( geos::unique_ptr &testedGeometry : testedGeometries )
1446 {
1447 newGeometries << QgsGeometry( fromGeos( testedGeometry.get() ) );
1448 }
1449
1450 return Success;
1451}
1452
1453geos::unique_ptr QgsGeos::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
1454{
1455 if ( !splitLine || !geom )
1456 return nullptr;
1457
1458 geos::unique_ptr geometryBoundary;
1459 GEOSContextHandle_t context = QgsGeosContext::get();
1460 if ( GEOSGeom_getDimensions_r( context, geom ) == 2 )
1461 geometryBoundary.reset( GEOSBoundary_r( context, geom ) );
1462 else
1463 geometryBoundary.reset( GEOSGeom_clone_r( context, geom ) );
1464
1465 geos::unique_ptr splitLineClone( GEOSGeom_clone_r( context, splitLine ) );
1466 geos::unique_ptr unionGeometry( GEOSUnion_r( context, splitLineClone.get(), geometryBoundary.get() ) );
1467
1468 return unionGeometry;
1469}
1470
1471int QgsGeos::mergeGeometriesMultiTypeSplit( std::vector<geos::unique_ptr> &splitResult ) const
1472{
1473 if ( !mGeos )
1474 return 1;
1475
1476 //convert mGeos to geometry collection
1477 GEOSContextHandle_t context = QgsGeosContext::get();
1478 int type = GEOSGeomTypeId_r( context, mGeos.get() );
1479 if ( type != GEOS_GEOMETRYCOLLECTION &&
1480 type != GEOS_MULTILINESTRING &&
1481 type != GEOS_MULTIPOLYGON &&
1482 type != GEOS_MULTIPOINT )
1483 return 0;
1484
1485 //collect all the geometries that belong to the initial multifeature
1486 std::vector<geos::unique_ptr> unionGeom;
1487
1488 std::vector<geos::unique_ptr> newSplitResult;
1489
1490 for ( size_t i = 0; i < splitResult.size(); ++i )
1491 {
1492 //is this geometry a part of the original multitype?
1493 bool isPart = false;
1494 for ( int j = 0; j < GEOSGetNumGeometries_r( context, mGeos.get() ); j++ )
1495 {
1496 if ( GEOSEquals_r( context, splitResult[i].get(), GEOSGetGeometryN_r( context, mGeos.get(), j ) ) )
1497 {
1498 isPart = true;
1499 break;
1500 }
1501 }
1502
1503 if ( isPart )
1504 {
1505 unionGeom.emplace_back( std::move( splitResult[i] ) );
1506 }
1507 else
1508 {
1509 std::vector<geos::unique_ptr> geomVector;
1510 geomVector.emplace_back( std::move( splitResult[i] ) );
1511
1512 if ( type == GEOS_MULTILINESTRING )
1513 newSplitResult.emplace_back( createGeosCollection( GEOS_MULTILINESTRING, geomVector ) );
1514 else if ( type == GEOS_MULTIPOLYGON )
1515 newSplitResult.emplace_back( createGeosCollection( GEOS_MULTIPOLYGON, geomVector ) );
1516 }
1517 }
1518
1519 splitResult = std::move( newSplitResult );
1520
1521 //make multifeature out of unionGeom
1522 if ( !unionGeom.empty() )
1523 {
1524 if ( type == GEOS_MULTILINESTRING )
1525 splitResult.emplace_back( createGeosCollection( GEOS_MULTILINESTRING, unionGeom ) );
1526 else if ( type == GEOS_MULTIPOLYGON )
1527 splitResult.emplace_back( createGeosCollection( GEOS_MULTIPOLYGON, unionGeom ) );
1528 }
1529
1530 return 0;
1531}
1532
1533geos::unique_ptr QgsGeos::createGeosCollection( int typeId, std::vector<geos::unique_ptr> &geoms )
1534{
1535 std::vector<GEOSGeometry *> geomarr;
1536 geomarr.reserve( geoms.size() );
1537
1538 GEOSContextHandle_t context = QgsGeosContext::get();
1539 for ( geos::unique_ptr &geomUniquePtr : geoms )
1540 {
1541 if ( geomUniquePtr )
1542 {
1543 if ( !GEOSisEmpty_r( context, geomUniquePtr.get() ) )
1544 {
1545 // don't add empty parts to a geos collection, it can cause crashes in GEOS
1546 // transfer ownership of the geometries to GEOSGeom_createCollection_r()
1547 geomarr.emplace_back( geomUniquePtr.release() );
1548 }
1549 }
1550 }
1551 geos::unique_ptr geomRes;
1552
1553 try
1554 {
1555 geomRes.reset( GEOSGeom_createCollection_r( context, typeId, geomarr.data(), geomarr.size() ) );
1556 }
1557 catch ( QgsGeosException & )
1558 {
1559 for ( GEOSGeometry *geom : geomarr )
1560 {
1561 GEOSGeom_destroy_r( context, geom );
1562 }
1563 }
1564
1565 return geomRes;
1566}
1567
1568std::unique_ptr<QgsAbstractGeometry> QgsGeos::fromGeos( const GEOSGeometry *geos )
1569{
1570 if ( !geos )
1571 {
1572 return nullptr;
1573 }
1574
1575 GEOSContextHandle_t context = QgsGeosContext::get();
1576 int nCoordDims = GEOSGeom_getCoordinateDimension_r( context, geos );
1577 int nDims = GEOSGeom_getDimensions_r( context, geos );
1578 bool hasZ = ( nCoordDims == 3 );
1579 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1580
1581 switch ( GEOSGeomTypeId_r( context, geos ) )
1582 {
1583 case GEOS_POINT: // a point
1584 {
1585 if ( GEOSisEmpty_r( context, geos ) )
1586 return nullptr;
1587
1588 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, geos );
1589 unsigned int nPoints = 0;
1590 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1591 if ( nPoints == 0 )
1592 {
1593 return nullptr;
1594 }
1595 // Since GEOS 3.13, Points with NAN coordinates are not considered empty anymore
1596 // See: https://github.com/libgeos/geos/pull/927
1597 // Handle this change by checking if QgsPoint is empty
1598 const QgsPoint point = coordSeqPoint( cs, 0, hasZ, hasM );
1599 return !point.isEmpty() ? std::unique_ptr<QgsAbstractGeometry>( point.clone() ) : nullptr;
1600 }
1601 case GEOS_LINESTRING:
1602 {
1603 return sequenceToLinestring( geos, hasZ, hasM );
1604 }
1605 case GEOS_POLYGON:
1606 {
1607 return fromGeosPolygon( geos );
1608 }
1609 case GEOS_MULTIPOINT:
1610 {
1611 auto multiPoint = std::make_unique<QgsMultiPoint>();
1612 int nParts = GEOSGetNumGeometries_r( context, geos );
1613 multiPoint->reserve( nParts );
1614 for ( int i = 0; i < nParts; ++i )
1615 {
1616 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, GEOSGetGeometryN_r( context, geos, i ) );
1617 if ( cs )
1618 {
1619 unsigned int nPoints = 0;
1620 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1621 if ( nPoints > 0 )
1622 multiPoint->addGeometry( coordSeqPoint( cs, 0, hasZ, hasM ).clone() );
1623 }
1624 }
1625 return std::move( multiPoint );
1626 }
1627 case GEOS_MULTILINESTRING:
1628 {
1629 auto multiLineString = std::make_unique<QgsMultiLineString>();
1630 int nParts = GEOSGetNumGeometries_r( context, geos );
1631 multiLineString->reserve( nParts );
1632 for ( int i = 0; i < nParts; ++i )
1633 {
1634 std::unique_ptr< QgsLineString >line( sequenceToLinestring( GEOSGetGeometryN_r( context, geos, i ), hasZ, hasM ) );
1635 if ( line )
1636 {
1637 multiLineString->addGeometry( line.release() );
1638 }
1639 }
1640 return std::move( multiLineString );
1641 }
1642 case GEOS_MULTIPOLYGON:
1643 {
1644 auto multiPolygon = std::make_unique<QgsMultiPolygon>();
1645
1646 int nParts = GEOSGetNumGeometries_r( context, geos );
1647 multiPolygon->reserve( nParts );
1648 for ( int i = 0; i < nParts; ++i )
1649 {
1650 std::unique_ptr< QgsPolygon > poly = fromGeosPolygon( GEOSGetGeometryN_r( context, geos, i ) );
1651 if ( poly )
1652 {
1653 multiPolygon->addGeometry( poly.release() );
1654 }
1655 }
1656 return std::move( multiPolygon );
1657 }
1658 case GEOS_GEOMETRYCOLLECTION:
1659 {
1660 auto geomCollection = std::make_unique<QgsGeometryCollection>();
1661 int nParts = GEOSGetNumGeometries_r( context, geos );
1662 geomCollection->reserve( nParts );
1663 for ( int i = 0; i < nParts; ++i )
1664 {
1665 std::unique_ptr< QgsAbstractGeometry > geom( fromGeos( GEOSGetGeometryN_r( context, geos, i ) ) );
1666 if ( geom )
1667 {
1668 geomCollection->addGeometry( geom.release() );
1669 }
1670 }
1671 return std::move( geomCollection );
1672 }
1673 }
1674 return nullptr;
1675}
1676
1677std::unique_ptr<QgsPolygon> QgsGeos::fromGeosPolygon( const GEOSGeometry *geos )
1678{
1679 GEOSContextHandle_t context = QgsGeosContext::get();
1680 if ( GEOSGeomTypeId_r( context, geos ) != GEOS_POLYGON )
1681 {
1682 return nullptr;
1683 }
1684
1685 int nCoordDims = GEOSGeom_getCoordinateDimension_r( context, geos );
1686 int nDims = GEOSGeom_getDimensions_r( context, geos );
1687 bool hasZ = ( nCoordDims == 3 );
1688 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1689
1690 auto polygon = std::make_unique<QgsPolygon>();
1691
1692 const GEOSGeometry *ring = GEOSGetExteriorRing_r( context, geos );
1693 if ( ring )
1694 {
1695 polygon->setExteriorRing( sequenceToLinestring( ring, hasZ, hasM ).release() );
1696 }
1697
1698 QVector<QgsCurve *> interiorRings;
1699 const int ringCount = GEOSGetNumInteriorRings_r( context, geos );
1700 interiorRings.reserve( ringCount );
1701 for ( int i = 0; i < ringCount; ++i )
1702 {
1703 ring = GEOSGetInteriorRingN_r( context, geos, i );
1704 if ( ring )
1705 {
1706 interiorRings.push_back( sequenceToLinestring( ring, hasZ, hasM ).release() );
1707 }
1708 }
1709 polygon->setInteriorRings( interiorRings );
1710
1711 return polygon;
1712}
1713
1714std::unique_ptr<QgsLineString> QgsGeos::sequenceToLinestring( const GEOSGeometry *geos, bool hasZ, bool hasM )
1715{
1716 GEOSContextHandle_t context = QgsGeosContext::get();
1717 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, geos );
1718
1719 unsigned int nPoints;
1720 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1721
1722 QVector< double > xOut( nPoints );
1723 QVector< double > yOut( nPoints );
1724 QVector< double > zOut;
1725 if ( hasZ )
1726 zOut.resize( nPoints );
1727 QVector< double > mOut;
1728 if ( hasM )
1729 mOut.resize( nPoints );
1730
1731 double *x = xOut.data();
1732 double *y = yOut.data();
1733 double *z = zOut.data();
1734 double *m = mOut.data();
1735
1736#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
1737 GEOSCoordSeq_copyToArrays_r( context, cs, x, y, hasZ ? z : nullptr, hasM ? m : nullptr );
1738#else
1739 for ( unsigned int i = 0; i < nPoints; ++i )
1740 {
1741 if ( hasZ )
1742 GEOSCoordSeq_getXYZ_r( context, cs, i, x++, y++, z++ );
1743 else
1744 GEOSCoordSeq_getXY_r( context, cs, i, x++, y++ );
1745 if ( hasM )
1746 {
1747 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, m++ );
1748 }
1749 }
1750#endif
1751 auto line = std::make_unique<QgsLineString>( xOut, yOut, zOut, mOut );
1752 return line;
1753}
1754
1755int QgsGeos::numberOfGeometries( GEOSGeometry *g )
1756{
1757 if ( !g )
1758 return 0;
1759
1760 GEOSContextHandle_t context = QgsGeosContext::get();
1761 int geometryType = GEOSGeomTypeId_r( context, g );
1762 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
1763 || geometryType == GEOS_POLYGON )
1764 return 1;
1765
1766 //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
1767 return GEOSGetNumGeometries_r( context, g );
1768}
1769
1770QgsPoint QgsGeos::coordSeqPoint( const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM )
1771{
1772 if ( !cs )
1773 {
1774 return QgsPoint();
1775 }
1776
1777 GEOSContextHandle_t context = QgsGeosContext::get();
1778
1779 double x, y;
1780 double z = 0;
1781 double m = 0;
1782 if ( hasZ )
1783 GEOSCoordSeq_getXYZ_r( context, cs, i, &x, &y, &z );
1784 else
1785 GEOSCoordSeq_getXY_r( context, cs, i, &x, &y );
1786 if ( hasM )
1787 {
1788 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, &m );
1789 }
1790
1792 if ( hasZ && hasM )
1793 {
1795 }
1796 else if ( hasZ )
1797 {
1799 }
1800 else if ( hasM )
1801 {
1803 }
1804 return QgsPoint( t, x, y, z, m );
1805}
1806
1807geos::unique_ptr QgsGeos::asGeos( const QgsAbstractGeometry *geom, double precision, Qgis::GeosCreationFlags flags )
1808{
1809 if ( !geom )
1810 return nullptr;
1811
1812 int coordDims = 2;
1813 if ( geom->is3D() )
1814 {
1815 ++coordDims;
1816 }
1817 if ( geom->isMeasure() )
1818 {
1819 ++coordDims;
1820 }
1821
1823 {
1824 int geosType = GEOS_GEOMETRYCOLLECTION;
1825
1827 {
1828 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1829 {
1831 geosType = GEOS_MULTIPOINT;
1832 break;
1833
1835 geosType = GEOS_MULTILINESTRING;
1836 break;
1837
1839 geosType = GEOS_MULTIPOLYGON;
1840 break;
1841
1844 return nullptr;
1845 }
1846 }
1847
1848
1849 const QgsGeometryCollection *c = qgsgeometry_cast<const QgsGeometryCollection *>( geom );
1850
1851 if ( !c )
1852 return nullptr;
1853
1854 std::vector<geos::unique_ptr> geomVector;
1855 geomVector.reserve( c->numGeometries() );
1856 for ( int i = 0; i < c->numGeometries(); ++i )
1857 {
1858 geos::unique_ptr geosGeom = asGeos( c->geometryN( i ), precision, flags );
1859 if ( flags & Qgis::GeosCreationFlag::RejectOnInvalidSubGeometry && !geosGeom )
1860 {
1861 return nullptr;
1862 }
1863 geomVector.emplace_back( std::move( geosGeom ) );
1864 }
1865 return createGeosCollection( geosType, geomVector );
1866 }
1869 {
1870 // PolyhedralSurface and TIN support
1871 // convert it to a geos MultiPolygon
1872 const QgsPolyhedralSurface *polyhedralSurface = qgsgeometry_cast<const QgsPolyhedralSurface *>( geom );
1873 if ( !polyhedralSurface )
1874 return nullptr;
1875
1876 std::vector<geos::unique_ptr> geomVector;
1877 geomVector.reserve( polyhedralSurface->numPatches() );
1878 for ( int i = 0; i < polyhedralSurface->numPatches(); ++i )
1879 {
1880 geos::unique_ptr geosPolygon = createGeosPolygon( polyhedralSurface->patchN( i ), precision );
1881 if ( flags & Qgis::GeosCreationFlag::RejectOnInvalidSubGeometry && !geosPolygon )
1882 {
1883 return nullptr;
1884 }
1885 geomVector.emplace_back( std::move( geosPolygon ) );
1886 }
1887
1888 return createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
1889 }
1890 else
1891 {
1892 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1893 {
1895 return createGeosPoint( static_cast<const QgsPoint *>( geom ), coordDims, precision, flags );
1896
1898 return createGeosLinestring( static_cast<const QgsLineString *>( geom ), precision, flags );
1899
1901 return createGeosPolygon( static_cast<const QgsPolygon *>( geom ), precision, flags );
1902
1905 return nullptr;
1906 }
1907 }
1908 return nullptr;
1909}
1910
1911std::unique_ptr<QgsAbstractGeometry> QgsGeos::overlay( const QgsAbstractGeometry *geom, Overlay op, QString *errorMsg, const QgsGeometryParameters &parameters ) const
1912{
1913 if ( !mGeos || !geom )
1914 {
1915 return nullptr;
1916 }
1917
1918 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
1919 if ( !geosGeom )
1920 {
1921 return nullptr;
1922 }
1923
1924 const double gridSize = parameters.gridSize();
1925
1926 GEOSContextHandle_t context = QgsGeosContext::get();
1927 try
1928 {
1929 geos::unique_ptr opGeom;
1930 switch ( op )
1931 {
1932 case OverlayIntersection:
1933 if ( gridSize > 0 )
1934 {
1935 opGeom.reset( GEOSIntersectionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1936 }
1937 else
1938 {
1939 opGeom.reset( GEOSIntersection_r( context, mGeos.get(), geosGeom.get() ) );
1940 }
1941 break;
1942
1943 case OverlayDifference:
1944 if ( gridSize > 0 )
1945 {
1946 opGeom.reset( GEOSDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1947 }
1948 else
1949 {
1950 opGeom.reset( GEOSDifference_r( context, mGeos.get(), geosGeom.get() ) );
1951 }
1952 break;
1953
1954 case OverlayUnion:
1955 {
1956 geos::unique_ptr unionGeometry;
1957 if ( gridSize > 0 )
1958 {
1959 unionGeometry.reset( GEOSUnionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1960 }
1961 else
1962 {
1963 unionGeometry.reset( GEOSUnion_r( context, mGeos.get(), geosGeom.get() ) );
1964 }
1965
1966 if ( unionGeometry && GEOSGeomTypeId_r( context, unionGeometry.get() ) == GEOS_MULTILINESTRING )
1967 {
1968 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, unionGeometry.get() ) );
1969 if ( mergedLines )
1970 {
1971 unionGeometry = std::move( mergedLines );
1972 }
1973 }
1974
1975 opGeom = std::move( unionGeometry );
1976 }
1977 break;
1978
1979 case OverlaySymDifference:
1980 if ( gridSize > 0 )
1981 {
1982 opGeom.reset( GEOSSymDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1983 }
1984 else
1985 {
1986 opGeom.reset( GEOSSymDifference_r( context, mGeos.get(), geosGeom.get() ) );
1987 }
1988 break;
1989 }
1990 return fromGeos( opGeom.get() );
1991 }
1992 catch ( QgsGeosException &e )
1993 {
1994 logError( QStringLiteral( "GEOS" ), e.what() );
1995 if ( errorMsg )
1996 {
1997 *errorMsg = e.what();
1998 }
1999 return nullptr;
2000 }
2001}
2002
2003bool QgsGeos::relation( const QgsAbstractGeometry *geom, Relation r, QString *errorMsg ) const
2004{
2005 if ( !mGeos || !geom )
2006 {
2007 return false;
2008 }
2009
2010 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2011 if ( !geosGeom )
2012 {
2013 return false;
2014 }
2015
2016 GEOSContextHandle_t context = QgsGeosContext::get();
2017 bool result = false;
2018 try
2019 {
2020 if ( mGeosPrepared ) //use faster version with prepared geometry
2021 {
2022 switch ( r )
2023 {
2024 case RelationIntersects:
2025 result = ( GEOSPreparedIntersects_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2026 break;
2027 case RelationTouches:
2028 result = ( GEOSPreparedTouches_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2029 break;
2030 case RelationCrosses:
2031 result = ( GEOSPreparedCrosses_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2032 break;
2033 case RelationWithin:
2034 result = ( GEOSPreparedWithin_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2035 break;
2036 case RelationContains:
2037 result = ( GEOSPreparedContains_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2038 break;
2039 case RelationDisjoint:
2040 result = ( GEOSPreparedDisjoint_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2041 break;
2042 case RelationOverlaps:
2043 result = ( GEOSPreparedOverlaps_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2044 break;
2045 }
2046 return result;
2047 }
2048
2049 switch ( r )
2050 {
2051 case RelationIntersects:
2052 result = ( GEOSIntersects_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2053 break;
2054 case RelationTouches:
2055 result = ( GEOSTouches_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2056 break;
2057 case RelationCrosses:
2058 result = ( GEOSCrosses_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2059 break;
2060 case RelationWithin:
2061 result = ( GEOSWithin_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2062 break;
2063 case RelationContains:
2064 result = ( GEOSContains_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2065 break;
2066 case RelationDisjoint:
2067 result = ( GEOSDisjoint_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2068 break;
2069 case RelationOverlaps:
2070 result = ( GEOSOverlaps_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2071 break;
2072 }
2073 }
2074 catch ( QgsGeosException &e )
2075 {
2076 logError( QStringLiteral( "GEOS" ), e.what() );
2077 if ( errorMsg )
2078 {
2079 *errorMsg = e.what();
2080 }
2081 return false;
2082 }
2083
2084 return result;
2085}
2086
2087QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, QString *errorMsg ) const
2088{
2089 if ( !mGeos )
2090 {
2091 return nullptr;
2092 }
2093
2094 geos::unique_ptr geos;
2095 try
2096 {
2097 geos.reset( GEOSBuffer_r( QgsGeosContext::get(), mGeos.get(), distance, segments ) );
2098 }
2099 CATCH_GEOS_WITH_ERRMSG( nullptr )
2100 return fromGeos( geos.get() ).release();
2101}
2102
2103QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, Qgis::EndCapStyle endCapStyle, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2104{
2105 geos::unique_ptr geos = buffer( mGeos.get(), distance, segments, endCapStyle, joinStyle, miterLimit, errorMsg );
2106 return fromGeos( geos.get() ).release();
2107}
2108
2109geos::unique_ptr QgsGeos::buffer( const GEOSGeometry *geometry, double distance, int segments, Qgis::EndCapStyle endCapStyle, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg )
2110{
2111 if ( !geometry )
2112 {
2113 return nullptr;
2114 }
2115
2116 geos::unique_ptr geos;
2117 try
2118 {
2119 geos.reset( GEOSBufferWithStyle_r( QgsGeosContext::get(), geometry, distance, segments, static_cast< int >( endCapStyle ), static_cast< int >( joinStyle ), miterLimit ) );
2120 }
2121 CATCH_GEOS_WITH_ERRMSG( nullptr )
2122 return geos;
2123}
2124
2125QgsAbstractGeometry *QgsGeos::simplify( double tolerance, QString *errorMsg ) const
2126{
2127 if ( !mGeos )
2128 {
2129 return nullptr;
2130 }
2131 geos::unique_ptr geos;
2132 try
2133 {
2134 geos.reset( GEOSTopologyPreserveSimplify_r( QgsGeosContext::get(), mGeos.get(), tolerance ) );
2135 }
2136 CATCH_GEOS_WITH_ERRMSG( nullptr )
2137 return fromGeos( geos.get() ).release();
2138}
2139
2140QgsAbstractGeometry *QgsGeos::interpolate( double distance, QString *errorMsg ) const
2141{
2142 if ( !mGeos )
2143 {
2144 return nullptr;
2145 }
2146 geos::unique_ptr geos;
2147 try
2148 {
2149 geos.reset( GEOSInterpolate_r( QgsGeosContext::get(), mGeos.get(), distance ) );
2150 }
2151 CATCH_GEOS_WITH_ERRMSG( nullptr )
2152 return fromGeos( geos.get() ).release();
2153}
2154
2155QgsPoint *QgsGeos::centroid( QString *errorMsg ) const
2156{
2157 if ( !mGeos )
2158 {
2159 return nullptr;
2160 }
2161
2162 geos::unique_ptr geos;
2163 double x;
2164 double y;
2165
2166 GEOSContextHandle_t context = QgsGeosContext::get();
2167 try
2168 {
2169 geos.reset( GEOSGetCentroid_r( context, mGeos.get() ) );
2170
2171 if ( !geos )
2172 return nullptr;
2173
2174 GEOSGeomGetX_r( context, geos.get(), &x );
2175 GEOSGeomGetY_r( context, geos.get(), &y );
2176 }
2177 CATCH_GEOS_WITH_ERRMSG( nullptr )
2178
2179 return new QgsPoint( x, y );
2180}
2181
2182QgsAbstractGeometry *QgsGeos::envelope( QString *errorMsg ) const
2183{
2184 if ( !mGeos )
2185 {
2186 return nullptr;
2187 }
2188 geos::unique_ptr geos;
2189 try
2190 {
2191 geos.reset( GEOSEnvelope_r( QgsGeosContext::get(), mGeos.get() ) );
2192 }
2193 CATCH_GEOS_WITH_ERRMSG( nullptr )
2194 return fromGeos( geos.get() ).release();
2195}
2196
2197QgsPoint *QgsGeos::pointOnSurface( QString *errorMsg ) const
2198{
2199 if ( !mGeos )
2200 {
2201 return nullptr;
2202 }
2203
2204 double x;
2205 double y;
2206
2207 GEOSContextHandle_t context = QgsGeosContext::get();
2208 geos::unique_ptr geos;
2209 try
2210 {
2211 geos.reset( GEOSPointOnSurface_r( context, mGeos.get() ) );
2212
2213 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
2214 {
2215 return nullptr;
2216 }
2217
2218 GEOSGeomGetX_r( context, geos.get(), &x );
2219 GEOSGeomGetY_r( context, geos.get(), &y );
2220 }
2221 CATCH_GEOS_WITH_ERRMSG( nullptr )
2222
2223 return new QgsPoint( x, y );
2224}
2225
2226QgsAbstractGeometry *QgsGeos::convexHull( QString *errorMsg ) const
2227{
2228 if ( !mGeos )
2229 {
2230 return nullptr;
2231 }
2232
2233 try
2234 {
2235 geos::unique_ptr cHull( GEOSConvexHull_r( QgsGeosContext::get(), mGeos.get() ) );
2236 std::unique_ptr< QgsAbstractGeometry > cHullGeom = fromGeos( cHull.get() );
2237 return cHullGeom.release();
2238 }
2239 CATCH_GEOS_WITH_ERRMSG( nullptr )
2240}
2241
2242std::unique_ptr< QgsAbstractGeometry > QgsGeos::concaveHull( double targetPercent, bool allowHoles, QString *errorMsg ) const
2243{
2244#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
2245 ( void )allowHoles;
2246 ( void )targetPercent;
2247 ( void )errorMsg;
2248 throw QgsNotSupportedException( QObject::tr( "Calculating concaveHull requires a QGIS build based on GEOS 3.11 or later" ) );
2249#else
2250 if ( !mGeos )
2251 {
2252 return nullptr;
2253 }
2254
2255 try
2256 {
2257 geos::unique_ptr concaveHull( GEOSConcaveHull_r( QgsGeosContext::get(), mGeos.get(), targetPercent, allowHoles ) );
2258 std::unique_ptr< QgsAbstractGeometry > concaveHullGeom = fromGeos( concaveHull.get() );
2259 return concaveHullGeom;
2260 }
2261 CATCH_GEOS_WITH_ERRMSG( nullptr )
2262#endif
2263}
2264
2265Qgis::CoverageValidityResult QgsGeos::validateCoverage( double gapWidth, std::unique_ptr<QgsAbstractGeometry> *invalidEdges, QString *errorMsg ) const
2266{
2267#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2268 ( void )gapWidth;
2269 ( void )invalidEdges;
2270 ( void )errorMsg;
2271 throw QgsNotSupportedException( QObject::tr( "Validating coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2272#else
2273 if ( !mGeos )
2274 {
2275 if ( errorMsg )
2276 *errorMsg = QStringLiteral( "Input geometry was not set" );
2278 }
2279
2280 GEOSContextHandle_t context = QgsGeosContext::get();
2281 try
2282 {
2283 GEOSGeometry *invalidEdgesGeos = nullptr;
2284 const int result = GEOSCoverageIsValid_r( context, mGeos.get(), gapWidth, invalidEdges ? &invalidEdgesGeos : nullptr );
2285 if ( invalidEdges && invalidEdgesGeos )
2286 {
2287 *invalidEdges = fromGeos( invalidEdgesGeos );
2288 }
2289 if ( invalidEdgesGeos )
2290 {
2291 GEOSGeom_destroy_r( context, invalidEdgesGeos );
2292 invalidEdgesGeos = nullptr;
2293 }
2294
2295 switch ( result )
2296 {
2297 case 0:
2299 case 1:
2301 case 2:
2302 break;
2303 }
2305 }
2307#endif
2308}
2309
2310std::unique_ptr<QgsAbstractGeometry> QgsGeos::simplifyCoverageVW( double tolerance, bool preserveBoundary, QString *errorMsg ) const
2311{
2312#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2313 ( void )tolerance;
2314 ( void )preserveBoundary;
2315 ( void )errorMsg;
2316 throw QgsNotSupportedException( QObject::tr( "Simplifying coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2317#else
2318 if ( !mGeos )
2319 {
2320 if ( errorMsg )
2321 *errorMsg = QStringLiteral( "Input geometry was not set" );
2322 return nullptr;
2323 }
2324
2325 try
2326 {
2327 geos::unique_ptr simplified( GEOSCoverageSimplifyVW_r( QgsGeosContext::get(), mGeos.get(), tolerance, preserveBoundary ? 1 : 0 ) );
2328 std::unique_ptr< QgsAbstractGeometry > simplifiedGeom = fromGeos( simplified.get() );
2329 return simplifiedGeom;
2330 }
2331 CATCH_GEOS_WITH_ERRMSG( nullptr )
2332#endif
2333}
2334
2335std::unique_ptr<QgsAbstractGeometry> QgsGeos::unionCoverage( QString *errorMsg ) const
2336{
2337 if ( !mGeos )
2338 {
2339 if ( errorMsg )
2340 *errorMsg = QStringLiteral( "Input geometry was not set" );
2341 return nullptr;
2342 }
2343
2344 try
2345 {
2346 geos::unique_ptr unioned( GEOSCoverageUnion_r( QgsGeosContext::get(), mGeos.get() ) );
2347 std::unique_ptr< QgsAbstractGeometry > result = fromGeos( unioned.get() );
2348 return result;
2349 }
2350 CATCH_GEOS_WITH_ERRMSG( nullptr )
2351}
2352
2353bool QgsGeos::isValid( QString *errorMsg, const bool allowSelfTouchingHoles, QgsGeometry *errorLoc ) const
2354{
2355 if ( !mGeos )
2356 {
2357 if ( errorMsg )
2358 *errorMsg = QObject::tr( "QGIS geometry cannot be converted to a GEOS geometry", "GEOS Error" );
2359 return false;
2360 }
2361
2362 GEOSContextHandle_t context = QgsGeosContext::get();
2363 try
2364 {
2365 GEOSGeometry *g1 = nullptr;
2366 char *r = nullptr;
2367 char res = GEOSisValidDetail_r( context, mGeos.get(), allowSelfTouchingHoles ? GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE : 0, &r, &g1 );
2368 const bool invalid = res != 1;
2369
2370 QString error;
2371 if ( r )
2372 {
2373 error = QString( r );
2374 GEOSFree_r( context, r );
2375 }
2376
2377 if ( invalid && errorMsg )
2378 {
2379 // Copied from https://github.com/libgeos/geos/blob/main/src/operation/valid/TopologyValidationError.cpp
2380 static const std::map< QString, QString > sTranslatedErrors
2381 {
2382 { QStringLiteral( "topology validation error" ), QObject::tr( "Topology validation error", "GEOS Error" ) },
2383 { QStringLiteral( "repeated point" ), QObject::tr( "Repeated point", "GEOS Error" ) },
2384 { QStringLiteral( "hole lies outside shell" ), QObject::tr( "Hole lies outside shell", "GEOS Error" ) },
2385 { QStringLiteral( "holes are nested" ), QObject::tr( "Holes are nested", "GEOS Error" ) },
2386 { QStringLiteral( "interior is disconnected" ), QObject::tr( "Interior is disconnected", "GEOS Error" ) },
2387 { QStringLiteral( "self-intersection" ), QObject::tr( "Self-intersection", "GEOS Error" ) },
2388 { QStringLiteral( "ring self-intersection" ), QObject::tr( "Ring self-intersection", "GEOS Error" ) },
2389 { QStringLiteral( "nested shells" ), QObject::tr( "Nested shells", "GEOS Error" ) },
2390 { QStringLiteral( "duplicate rings" ), QObject::tr( "Duplicate rings", "GEOS Error" ) },
2391 { QStringLiteral( "too few points in geometry component" ), QObject::tr( "Too few points in geometry component", "GEOS Error" ) },
2392 { QStringLiteral( "invalid coordinate" ), QObject::tr( "Invalid coordinate", "GEOS Error" ) },
2393 { QStringLiteral( "ring is not closed" ), QObject::tr( "Ring is not closed", "GEOS Error" ) },
2394 };
2395
2396 const auto translatedError = sTranslatedErrors.find( error.toLower() );
2397 if ( translatedError != sTranslatedErrors.end() )
2398 *errorMsg = translatedError->second;
2399 else
2400 *errorMsg = error;
2401
2402 if ( g1 && errorLoc )
2403 {
2404 *errorLoc = geometryFromGeos( g1 );
2405 }
2406 else if ( g1 )
2407 {
2408 GEOSGeom_destroy_r( context, g1 );
2409 }
2410 }
2411 return !invalid;
2412 }
2413 CATCH_GEOS_WITH_ERRMSG( false )
2414}
2415
2416bool QgsGeos::isEqual( const QgsAbstractGeometry *geom, QString *errorMsg ) const
2417{
2418 if ( !mGeos || !geom )
2419 {
2420 return false;
2421 }
2422
2423 try
2424 {
2425 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2426 if ( !geosGeom )
2427 {
2428 return false;
2429 }
2430 bool equal = GEOSEquals_r( QgsGeosContext::get(), mGeos.get(), geosGeom.get() );
2431 return equal;
2432 }
2433 CATCH_GEOS_WITH_ERRMSG( false )
2434}
2435
2436bool QgsGeos::isEmpty( QString *errorMsg ) const
2437{
2438 if ( !mGeos )
2439 {
2440 return false;
2441 }
2442
2443 try
2444 {
2445 return GEOSisEmpty_r( QgsGeosContext::get(), mGeos.get() );
2446 }
2447 CATCH_GEOS_WITH_ERRMSG( false )
2448}
2449
2450bool QgsGeos::isSimple( QString *errorMsg ) const
2451{
2452 if ( !mGeos )
2453 {
2454 return false;
2455 }
2456
2457 try
2458 {
2459 return GEOSisSimple_r( QgsGeosContext::get(), mGeos.get() );
2460 }
2461 CATCH_GEOS_WITH_ERRMSG( false )
2462}
2463
2464GEOSCoordSequence *QgsGeos::createCoordinateSequence( const QgsCurve *curve, double precision, bool forceClose )
2465{
2466 GEOSContextHandle_t context = QgsGeosContext::get();
2467
2468 std::unique_ptr< QgsLineString > segmentized;
2469 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( curve );
2470
2471 if ( !line )
2472 {
2473 segmentized.reset( curve->curveToLine() );
2474 line = segmentized.get();
2475 }
2476
2477 if ( !line )
2478 {
2479 return nullptr;
2480 }
2481 GEOSCoordSequence *coordSeq = nullptr;
2482
2483 const int numPoints = line->numPoints();
2484
2485 const bool hasZ = line->is3D();
2486
2487#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
2488 if ( qgsDoubleNear( precision, 0 ) )
2489 {
2490 if ( !forceClose || ( line->pointN( 0 ) == line->pointN( numPoints - 1 ) ) )
2491 {
2492 // use optimised method if we don't have to force close an open ring
2493 try
2494 {
2495 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, line->xData(), line->yData(), line->zData(), nullptr, numPoints );
2496 if ( !coordSeq )
2497 {
2498 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points" ).arg( numPoints ) );
2499 return nullptr;
2500 }
2501 }
2502 CATCH_GEOS( nullptr )
2503 }
2504 else
2505 {
2506 QVector< double > x = line->xVector();
2507 if ( numPoints > 0 )
2508 x.append( x.at( 0 ) );
2509 QVector< double > y = line->yVector();
2510 if ( numPoints > 0 )
2511 y.append( y.at( 0 ) );
2512 QVector< double > z = line->zVector();
2513 if ( hasZ && numPoints > 0 )
2514 z.append( z.at( 0 ) );
2515 try
2516 {
2517 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, x.constData(), y.constData(), !hasZ ? nullptr : z.constData(), nullptr, numPoints + 1 );
2518 if ( !coordSeq )
2519 {
2520 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create closed coordinate sequence for %1 points" ).arg( numPoints + 1 ) );
2521 return nullptr;
2522 }
2523 }
2524 CATCH_GEOS( nullptr )
2525 }
2526 return coordSeq;
2527 }
2528#endif
2529
2530 int coordDims = 2;
2531 const bool hasM = false; //line->isMeasure(); //disabled until geos supports m-coordinates
2532
2533 if ( hasZ )
2534 {
2535 ++coordDims;
2536 }
2537 if ( hasM )
2538 {
2539 ++coordDims;
2540 }
2541
2542 int numOutPoints = numPoints;
2543 if ( forceClose && ( line->pointN( 0 ) != line->pointN( numPoints - 1 ) ) )
2544 {
2545 ++numOutPoints;
2546 }
2547
2548 try
2549 {
2550 coordSeq = GEOSCoordSeq_create_r( context, numOutPoints, coordDims );
2551 if ( !coordSeq )
2552 {
2553 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ) );
2554 return nullptr;
2555 }
2556
2557 const double *xData = line->xData();
2558 const double *yData = line->yData();
2559 const double *zData = hasZ ? line->zData() : nullptr;
2560 const double *mData = hasM ? line->mData() : nullptr;
2561
2562 if ( precision > 0. )
2563 {
2564 for ( int i = 0; i < numOutPoints; ++i )
2565 {
2566 if ( i >= numPoints )
2567 {
2568 // start reading back from start of line
2569 xData = line->xData();
2570 yData = line->yData();
2571 zData = hasZ ? line->zData() : nullptr;
2572 mData = hasM ? line->mData() : nullptr;
2573 }
2574 if ( hasZ )
2575 {
2576 GEOSCoordSeq_setXYZ_r( context, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision, std::round( *zData++ / precision ) * precision );
2577 }
2578 else
2579 {
2580 GEOSCoordSeq_setXY_r( context, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision );
2581 }
2582 if ( hasM )
2583 {
2584 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2585 }
2586 }
2587 }
2588 else
2589 {
2590 for ( int i = 0; i < numOutPoints; ++i )
2591 {
2592 if ( i >= numPoints )
2593 {
2594 // start reading back from start of line
2595 xData = line->xData();
2596 yData = line->yData();
2597 zData = hasZ ? line->zData() : nullptr;
2598 mData = hasM ? line->mData() : nullptr;
2599 }
2600 if ( hasZ )
2601 {
2602 GEOSCoordSeq_setXYZ_r( context, coordSeq, i, *xData++, *yData++, *zData++ );
2603 }
2604 else
2605 {
2606 GEOSCoordSeq_setXY_r( context, coordSeq, i, *xData++, *yData++ );
2607 }
2608 if ( hasM )
2609 {
2610 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2611 }
2612 }
2613 }
2614 }
2615 CATCH_GEOS( nullptr )
2616
2617 return coordSeq;
2618}
2619
2620geos::unique_ptr QgsGeos::createGeosPoint( const QgsAbstractGeometry *point, int coordDims, double precision, Qgis::GeosCreationFlags )
2621{
2622 const QgsPoint *pt = qgsgeometry_cast<const QgsPoint *>( point );
2623 if ( !pt )
2624 return nullptr;
2625
2626 return createGeosPointXY( pt->x(), pt->y(), pt->is3D(), pt->z(), pt->isMeasure(), pt->m(), coordDims, precision );
2627}
2628
2629geos::unique_ptr QgsGeos::createGeosPointXY( double x, double y, bool hasZ, double z, bool hasM, double m, int coordDims, double precision, Qgis::GeosCreationFlags )
2630{
2631 Q_UNUSED( hasM )
2632 Q_UNUSED( m )
2633
2634 geos::unique_ptr geosPoint;
2635 GEOSContextHandle_t context = QgsGeosContext::get();
2636 try
2637 {
2638 if ( coordDims == 2 )
2639 {
2640 // optimised constructor
2641 if ( precision > 0. )
2642 geosPoint.reset( GEOSGeom_createPointFromXY_r( context, std::round( x / precision ) * precision, std::round( y / precision ) * precision ) );
2643 else
2644 geosPoint.reset( GEOSGeom_createPointFromXY_r( context, x, y ) );
2645 return geosPoint;
2646 }
2647
2648 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( context, 1, coordDims );
2649 if ( !coordSeq )
2650 {
2651 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ) );
2652 return nullptr;
2653 }
2654 if ( precision > 0. )
2655 {
2656 GEOSCoordSeq_setX_r( context, coordSeq, 0, std::round( x / precision ) * precision );
2657 GEOSCoordSeq_setY_r( context, coordSeq, 0, std::round( y / precision ) * precision );
2658 if ( hasZ )
2659 {
2660 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, std::round( z / precision ) * precision );
2661 }
2662 }
2663 else
2664 {
2665 GEOSCoordSeq_setX_r( context, coordSeq, 0, x );
2666 GEOSCoordSeq_setY_r( context, coordSeq, 0, y );
2667 if ( hasZ )
2668 {
2669 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, z );
2670 }
2671 }
2672#if 0 //disabled until geos supports m-coordinates
2673 if ( hasM )
2674 {
2675 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 3, m );
2676 }
2677#endif
2678 geosPoint.reset( GEOSGeom_createPoint_r( context, coordSeq ) );
2679 }
2680 CATCH_GEOS( nullptr )
2681 return geosPoint;
2682}
2683
2684geos::unique_ptr QgsGeos::createGeosLinestring( const QgsAbstractGeometry *curve, double precision, Qgis::GeosCreationFlags )
2685{
2686 const QgsCurve *c = qgsgeometry_cast<const QgsCurve *>( curve );
2687 if ( !c )
2688 return nullptr;
2689
2690 GEOSCoordSequence *coordSeq = createCoordinateSequence( c, precision );
2691 if ( !coordSeq )
2692 return nullptr;
2693
2694 geos::unique_ptr geosGeom;
2695 try
2696 {
2697 geosGeom.reset( GEOSGeom_createLineString_r( QgsGeosContext::get(), coordSeq ) );
2698 }
2699 CATCH_GEOS( nullptr )
2700 return geosGeom;
2701}
2702
2703geos::unique_ptr QgsGeos::createGeosPolygon( const QgsAbstractGeometry *poly, double precision, Qgis::GeosCreationFlags flags )
2704{
2705 const QgsCurvePolygon *polygon = qgsgeometry_cast<const QgsCurvePolygon *>( poly );
2706 if ( !polygon )
2707 return nullptr;
2708
2709 const QgsCurve *exteriorRing = polygon->exteriorRing();
2710 if ( !exteriorRing )
2711 {
2712 return nullptr;
2713 }
2714
2715 GEOSContextHandle_t context = QgsGeosContext::get();
2716 geos::unique_ptr geosPolygon;
2717 try
2718 {
2719 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( context, createCoordinateSequence( exteriorRing, precision, true ) ) );
2720
2721 int nHoles = 0;
2722 int nInteriorRings = polygon->numInteriorRings();
2724 {
2725 for ( int i = 0; i < nInteriorRings; ++i )
2726 {
2727 const QgsCurve *interiorRing = polygon->interiorRing( i );
2728 if ( !interiorRing->isEmpty() )
2729 {
2730 nHoles++;
2731 }
2732 }
2733 }
2734 else
2735 {
2736 nHoles = nInteriorRings;
2737 }
2738 GEOSGeometry **holes = nullptr;
2739 if ( nHoles > 0 )
2740 {
2741 holes = new GEOSGeometry*[ nHoles ];
2742 }
2743
2744 for ( int i = 0; i < nInteriorRings; ++i )
2745 {
2746 const QgsCurve *interiorRing = polygon->interiorRing( i );
2747 if ( !( flags & Qgis::GeosCreationFlag::SkipEmptyInteriorRings ) || !interiorRing->isEmpty() )
2748 {
2749 holes[i] = GEOSGeom_createLinearRing_r( context, createCoordinateSequence( interiorRing, precision, true ) );
2750 }
2751 }
2752 geosPolygon.reset( GEOSGeom_createPolygon_r( context, exteriorRingGeos.release(), holes, nHoles ) );
2753 delete[] holes;
2754 }
2755 CATCH_GEOS( nullptr )
2756
2757 return geosPolygon;
2758}
2759
2760geos::unique_ptr QgsGeos::offsetCurve( const GEOSGeometry *geometry, double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg )
2761{
2762 if ( !geometry )
2763 return nullptr;
2764
2765 geos::unique_ptr offset;
2766 try
2767 {
2768 // Force quadrant segments to be at least 8, see
2769 // https://github.com/qgis/QGIS/issues/53165#issuecomment-1563470832
2770 if ( segments < 8 )
2771 segments = 8;
2772 offset.reset( GEOSOffsetCurve_r( QgsGeosContext::get(), geometry, distance, segments, static_cast< int >( joinStyle ), miterLimit ) );
2773 }
2774 CATCH_GEOS_WITH_ERRMSG( nullptr )
2775 return offset;
2776}
2777
2778QgsAbstractGeometry *QgsGeos::offsetCurve( double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2779{
2780 geos::unique_ptr res = offsetCurve( mGeos.get(), distance, segments, joinStyle, miterLimit, errorMsg );
2781 if ( !res )
2782 return nullptr;
2783
2784 return fromGeos( res.get() ).release();
2785}
2786
2787std::unique_ptr<QgsAbstractGeometry> QgsGeos::singleSidedBuffer( double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2788{
2789 if ( !mGeos )
2790 {
2791 return nullptr;
2792 }
2793
2794 geos::unique_ptr geos;
2795 GEOSContextHandle_t context = QgsGeosContext::get();
2796 try
2797 {
2798 geos::buffer_params_unique_ptr bp( GEOSBufferParams_create_r( context ) );
2799 GEOSBufferParams_setSingleSided_r( context, bp.get(), 1 );
2800 GEOSBufferParams_setQuadrantSegments_r( context, bp.get(), segments );
2801 GEOSBufferParams_setJoinStyle_r( context, bp.get(), static_cast< int >( joinStyle ) );
2802 GEOSBufferParams_setMitreLimit_r( context, bp.get(), miterLimit ); //#spellok
2803
2804 if ( side == Qgis::BufferSide::Right )
2805 {
2806 distance = -distance;
2807 }
2808 geos.reset( GEOSBufferWithParams_r( context, mGeos.get(), bp.get(), distance ) );
2809 }
2810 CATCH_GEOS_WITH_ERRMSG( nullptr )
2811 return fromGeos( geos.get() );
2812}
2813
2814std::unique_ptr<QgsAbstractGeometry> QgsGeos::maximumInscribedCircle( double tolerance, QString *errorMsg ) const
2815{
2816 if ( !mGeos )
2817 {
2818 return nullptr;
2819 }
2820
2821 geos::unique_ptr geos;
2822 try
2823 {
2824 geos.reset( GEOSMaximumInscribedCircle_r( QgsGeosContext::get(), mGeos.get(), tolerance ) );
2825 }
2826 CATCH_GEOS_WITH_ERRMSG( nullptr )
2827 return fromGeos( geos.get() );
2828}
2829
2830std::unique_ptr<QgsAbstractGeometry> QgsGeos::largestEmptyCircle( double tolerance, const QgsAbstractGeometry *boundary, QString *errorMsg ) const
2831{
2832 if ( !mGeos )
2833 {
2834 return nullptr;
2835 }
2836
2837 geos::unique_ptr geos;
2838 try
2839 {
2840 geos::unique_ptr boundaryGeos;
2841 if ( boundary )
2842 boundaryGeos = asGeos( boundary );
2843
2844 geos.reset( GEOSLargestEmptyCircle_r( QgsGeosContext::get(), mGeos.get(), boundaryGeos.get(), tolerance ) );
2845 }
2846 CATCH_GEOS_WITH_ERRMSG( nullptr )
2847 return fromGeos( geos.get() );
2848}
2849
2850std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumWidth( QString *errorMsg ) const
2851{
2852 if ( !mGeos )
2853 {
2854 return nullptr;
2855 }
2856
2857 geos::unique_ptr geos;
2858 try
2859 {
2860 geos.reset( GEOSMinimumWidth_r( QgsGeosContext::get(), mGeos.get() ) );
2861 }
2862 CATCH_GEOS_WITH_ERRMSG( nullptr )
2863 return fromGeos( geos.get() );
2864}
2865
2866double QgsGeos::minimumClearance( QString *errorMsg ) const
2867{
2868 if ( !mGeos )
2869 {
2870 return std::numeric_limits< double >::quiet_NaN();
2871 }
2872
2873 geos::unique_ptr geos;
2874 double res = 0;
2875 try
2876 {
2877 if ( GEOSMinimumClearance_r( QgsGeosContext::get(), mGeos.get(), &res ) != 0 )
2878 return std::numeric_limits< double >::quiet_NaN();
2879 }
2880 CATCH_GEOS_WITH_ERRMSG( std::numeric_limits< double >::quiet_NaN() )
2881 return res;
2882}
2883
2884std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumClearanceLine( QString *errorMsg ) const
2885{
2886 if ( !mGeos )
2887 {
2888 return nullptr;
2889 }
2890
2891 geos::unique_ptr geos;
2892 try
2893 {
2894 geos.reset( GEOSMinimumClearanceLine_r( QgsGeosContext::get(), mGeos.get() ) );
2895 }
2896 CATCH_GEOS_WITH_ERRMSG( nullptr )
2897 return fromGeos( geos.get() );
2898}
2899
2900std::unique_ptr<QgsAbstractGeometry> QgsGeos::node( QString *errorMsg ) const
2901{
2902 if ( !mGeos )
2903 {
2904 return nullptr;
2905 }
2906
2907 geos::unique_ptr geos;
2908 try
2909 {
2910 geos.reset( GEOSNode_r( QgsGeosContext::get(), mGeos.get() ) );
2911 }
2912 CATCH_GEOS_WITH_ERRMSG( nullptr )
2913 return fromGeos( geos.get() );
2914}
2915
2916std::unique_ptr<QgsAbstractGeometry> QgsGeos::sharedPaths( const QgsAbstractGeometry *other, QString *errorMsg ) const
2917{
2918 if ( !mGeos || !other )
2919 {
2920 return nullptr;
2921 }
2922
2923 geos::unique_ptr geos;
2924 try
2925 {
2926 geos::unique_ptr otherGeos = asGeos( other );
2927 if ( !otherGeos )
2928 return nullptr;
2929
2930 geos.reset( GEOSSharedPaths_r( QgsGeosContext::get(), mGeos.get(), otherGeos.get() ) );
2931 }
2932 CATCH_GEOS_WITH_ERRMSG( nullptr )
2933 return fromGeos( geos.get() );
2934}
2935
2936std::unique_ptr<QgsAbstractGeometry> QgsGeos::reshapeGeometry( const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg ) const
2937{
2938 if ( !mGeos || mGeometry->dimension() == 0 )
2939 {
2940 if ( errorCode ) { *errorCode = InvalidBaseGeometry; }
2941 return nullptr;
2942 }
2943
2944 if ( reshapeWithLine.numPoints() < 2 )
2945 {
2946 if ( errorCode ) { *errorCode = InvalidInput; }
2947 return nullptr;
2948 }
2949
2950 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
2951
2952 GEOSContextHandle_t context = QgsGeosContext::get();
2953 //single or multi?
2954 int numGeoms = GEOSGetNumGeometries_r( context, mGeos.get() );
2955 if ( numGeoms == -1 )
2956 {
2957 if ( errorCode )
2958 {
2959 *errorCode = InvalidBaseGeometry;
2960 }
2961 return nullptr;
2962 }
2963
2964 bool isMultiGeom = false;
2965 int geosTypeId = GEOSGeomTypeId_r( context, mGeos.get() );
2966 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
2967 isMultiGeom = true;
2968
2969 bool isLine = ( mGeometry->dimension() == 1 );
2970
2971 if ( !isMultiGeom )
2972 {
2973 geos::unique_ptr reshapedGeometry;
2974 if ( isLine )
2975 {
2976 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2977 }
2978 else
2979 {
2980 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2981 }
2982
2983 if ( errorCode )
2984 {
2985 if ( reshapedGeometry )
2986 *errorCode = Success;
2987 else
2988 *errorCode = NothingHappened;
2989 }
2990
2991 std::unique_ptr< QgsAbstractGeometry > reshapeResult = fromGeos( reshapedGeometry.get() );
2992 return reshapeResult;
2993 }
2994 else
2995 {
2996 try
2997 {
2998 //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
2999 bool reshapeTookPlace = false;
3000
3001 geos::unique_ptr currentReshapeGeometry;
3002 GEOSGeometry **newGeoms = new GEOSGeometry*[numGeoms];
3003
3004 for ( int i = 0; i < numGeoms; ++i )
3005 {
3006 if ( isLine )
3007 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
3008 else
3009 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
3010
3011 if ( currentReshapeGeometry )
3012 {
3013 newGeoms[i] = currentReshapeGeometry.release();
3014 reshapeTookPlace = true;
3015 }
3016 else
3017 {
3018 newGeoms[i] = GEOSGeom_clone_r( context, GEOSGetGeometryN_r( context, mGeos.get(), i ) );
3019 }
3020 }
3021
3022 geos::unique_ptr newMultiGeom;
3023 if ( isLine )
3024 {
3025 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
3026 }
3027 else //multipolygon
3028 {
3029 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
3030 }
3031
3032 delete[] newGeoms;
3033 if ( !newMultiGeom )
3034 {
3035 if ( errorCode ) { *errorCode = EngineError; }
3036 return nullptr;
3037 }
3038
3039 if ( reshapeTookPlace )
3040 {
3041 if ( errorCode )
3042 *errorCode = Success;
3043 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom = fromGeos( newMultiGeom.get() );
3044 return reshapedMultiGeom;
3045 }
3046 else
3047 {
3048 if ( errorCode )
3049 {
3050 *errorCode = NothingHappened;
3051 }
3052 return nullptr;
3053 }
3054 }
3055 CATCH_GEOS_WITH_ERRMSG( nullptr )
3056 }
3057}
3058
3059std::unique_ptr< QgsAbstractGeometry > QgsGeos::mergeLines( QString *errorMsg, const QgsGeometryParameters &parameters ) const
3060{
3061 if ( !mGeos )
3062 {
3063 return nullptr;
3064 }
3065
3066 GEOSContextHandle_t context = QgsGeosContext::get();
3067 if ( GEOSGeomTypeId_r( context, mGeos.get() ) != GEOS_MULTILINESTRING )
3068 return nullptr;
3069
3070 geos::unique_ptr geos;
3071 try
3072 {
3073 double gridSize = parameters.gridSize();
3074 if ( gridSize > 0 )
3075 {
3076 geos::unique_ptr geosFixedSize( GEOSGeom_setPrecision_r( context, mGeos.get(), gridSize, 0 ) );
3077 geos.reset( GEOSLineMerge_r( context, geosFixedSize.get() ) );
3078 }
3079 else
3080 geos.reset( GEOSLineMerge_r( context, mGeos.get() ) );
3081 }
3082 CATCH_GEOS_WITH_ERRMSG( nullptr )
3083 return fromGeos( geos.get() );
3084}
3085
3086std::unique_ptr<QgsAbstractGeometry> QgsGeos::closestPoint( const QgsGeometry &other, QString *errorMsg ) const
3087{
3088 if ( !mGeos || isEmpty() || other.isEmpty() )
3089 {
3090 return nullptr;
3091 }
3092
3093 geos::unique_ptr otherGeom( asGeos( other.constGet(), mPrecision ) );
3094 if ( !otherGeom )
3095 {
3096 return nullptr;
3097 }
3098
3099 GEOSContextHandle_t context = QgsGeosContext::get();
3100 double nx = 0.0;
3101 double ny = 0.0;
3102 try
3103 {
3104 geos::coord_sequence_unique_ptr nearestCoord;
3105 if ( mGeosPrepared ) // use faster version with prepared geometry
3106 {
3107 nearestCoord.reset( GEOSPreparedNearestPoints_r( context, mGeosPrepared.get(), otherGeom.get() ) );
3108 }
3109 else
3110 {
3111 nearestCoord.reset( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3112 }
3113
3114 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx );
3115 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny );
3116 }
3117 catch ( QgsGeosException &e )
3118 {
3119 logError( QStringLiteral( "GEOS" ), e.what() );
3120 if ( errorMsg )
3121 {
3122 *errorMsg = e.what();
3123 }
3124 return nullptr;
3125 }
3126
3127 return std::make_unique< QgsPoint >( nx, ny );
3128}
3129
3130std::unique_ptr<QgsAbstractGeometry> QgsGeos::shortestLine( const QgsGeometry &other, QString *errorMsg ) const
3131{
3132 if ( !mGeos || other.isEmpty() )
3133 {
3134 return nullptr;
3135 }
3136
3137 return shortestLine( other.constGet(), errorMsg );
3138}
3139
3140std::unique_ptr< QgsAbstractGeometry > QgsGeos::shortestLine( const QgsAbstractGeometry *other, QString *errorMsg ) const
3141{
3142 if ( !other || other->isEmpty() )
3143 return nullptr;
3144
3145 geos::unique_ptr otherGeom( asGeos( other, mPrecision ) );
3146 if ( !otherGeom )
3147 {
3148 return nullptr;
3149 }
3150
3151 GEOSContextHandle_t context = QgsGeosContext::get();
3152 double nx1 = 0.0;
3153 double ny1 = 0.0;
3154 double nx2 = 0.0;
3155 double ny2 = 0.0;
3156 try
3157 {
3158 geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3159
3160 if ( !nearestCoord )
3161 {
3162 if ( errorMsg )
3163 *errorMsg = QStringLiteral( "GEOS returned no nearest points" );
3164 return nullptr;
3165 }
3166
3167 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx1 );
3168 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny1 );
3169 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 1, &nx2 );
3170 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 1, &ny2 );
3171 }
3172 catch ( QgsGeosException &e )
3173 {
3174 logError( QStringLiteral( "GEOS" ), e.what() );
3175 if ( errorMsg )
3176 {
3177 *errorMsg = e.what();
3178 }
3179 return nullptr;
3180 }
3181
3182 auto line = std::make_unique< QgsLineString >();
3183 line->addVertex( QgsPoint( nx1, ny1 ) );
3184 line->addVertex( QgsPoint( nx2, ny2 ) );
3185 return line;
3186}
3187
3188double QgsGeos::lineLocatePoint( const QgsPoint &point, QString *errorMsg ) const
3189{
3190 if ( !mGeos )
3191 {
3192 return -1;
3193 }
3194
3195 geos::unique_ptr otherGeom( asGeos( &point, mPrecision ) );
3196 if ( !otherGeom )
3197 {
3198 return -1;
3199 }
3200
3201 double distance = -1;
3202 try
3203 {
3204 distance = GEOSProject_r( QgsGeosContext::get(), mGeos.get(), otherGeom.get() );
3205 }
3206 catch ( QgsGeosException &e )
3207 {
3208 logError( QStringLiteral( "GEOS" ), e.what() );
3209 if ( errorMsg )
3210 {
3211 *errorMsg = e.what();
3212 }
3213 return -1;
3214 }
3215
3216 return distance;
3217}
3218
3219double QgsGeos::lineLocatePoint( double x, double y, QString *errorMsg ) const
3220{
3221 if ( !mGeos )
3222 {
3223 return -1;
3224 }
3225
3226 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
3227 if ( !point )
3228 return false;
3229
3230 double distance = -1;
3231 try
3232 {
3233 distance = GEOSProject_r( QgsGeosContext::get(), mGeos.get(), point.get() );
3234 }
3235 catch ( QgsGeosException &e )
3236 {
3237 logError( QStringLiteral( "GEOS" ), e.what() );
3238 if ( errorMsg )
3239 {
3240 *errorMsg = e.what();
3241 }
3242 return -1;
3243 }
3244
3245 return distance;
3246}
3247
3248QgsGeometry QgsGeos::polygonize( const QVector<const QgsAbstractGeometry *> &geometries, QString *errorMsg )
3249{
3250 GEOSGeometry **const lineGeosGeometries = new GEOSGeometry*[ geometries.size()];
3251 int validLines = 0;
3252 for ( const QgsAbstractGeometry *g : geometries )
3253 {
3254 geos::unique_ptr l = asGeos( g );
3255 if ( l )
3256 {
3257 lineGeosGeometries[validLines] = l.release();
3258 validLines++;
3259 }
3260 }
3261
3262 GEOSContextHandle_t context = QgsGeosContext::get();
3263 try
3264 {
3265 geos::unique_ptr result( GEOSPolygonize_r( context, lineGeosGeometries, validLines ) );
3266 for ( int i = 0; i < validLines; ++i )
3267 {
3268 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3269 }
3270 delete[] lineGeosGeometries;
3271 return QgsGeometry( fromGeos( result.get() ) );
3272 }
3273 catch ( QgsGeosException &e )
3274 {
3275 if ( errorMsg )
3276 {
3277 *errorMsg = e.what();
3278 }
3279 for ( int i = 0; i < validLines; ++i )
3280 {
3281 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3282 }
3283 delete[] lineGeosGeometries;
3284 return QgsGeometry();
3285 }
3286}
3287
3288std::unique_ptr<QgsAbstractGeometry> QgsGeos::voronoiDiagram( const QgsAbstractGeometry *extent, double tolerance, bool edgesOnly, QString *errorMsg ) const
3289{
3290 if ( !mGeos )
3291 {
3292 return nullptr;
3293 }
3294
3295 geos::unique_ptr extentGeosGeom;
3296 if ( extent )
3297 {
3298 extentGeosGeom = asGeos( extent, mPrecision );
3299 if ( !extentGeosGeom )
3300 {
3301 return nullptr;
3302 }
3303 }
3304
3305 geos::unique_ptr geos;
3306 GEOSContextHandle_t context = QgsGeosContext::get();
3307 try
3308 {
3309 geos.reset( GEOSVoronoiDiagram_r( context, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
3310
3311 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3312 {
3313 return nullptr;
3314 }
3315
3316 return fromGeos( geos.get() );
3317 }
3318 CATCH_GEOS_WITH_ERRMSG( nullptr )
3319}
3320
3321std::unique_ptr<QgsAbstractGeometry> QgsGeos::delaunayTriangulation( double tolerance, bool edgesOnly, QString *errorMsg ) const
3322{
3323 if ( !mGeos )
3324 {
3325 return nullptr;
3326 }
3327
3328 GEOSContextHandle_t context = QgsGeosContext::get();
3329 geos::unique_ptr geos;
3330 try
3331 {
3332 geos.reset( GEOSDelaunayTriangulation_r( context, mGeos.get(), tolerance, edgesOnly ) );
3333
3334 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3335 {
3336 return nullptr;
3337 }
3338
3339 return fromGeos( geos.get() );
3340 }
3341 CATCH_GEOS_WITH_ERRMSG( nullptr )
3342}
3343
3344std::unique_ptr<QgsAbstractGeometry> QgsGeos::constrainedDelaunayTriangulation( QString *errorMsg ) const
3345{
3346#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
3347 ( void )errorMsg;
3348 throw QgsNotSupportedException( QObject::tr( "Calculating constrainedDelaunayTriangulation requires a QGIS build based on GEOS 3.11 or later" ) );
3349#else
3350 if ( !mGeos )
3351 {
3352 return nullptr;
3353 }
3354
3355 geos::unique_ptr geos;
3356 GEOSContextHandle_t context = QgsGeosContext::get();
3357 try
3358 {
3359 geos.reset( GEOSConstrainedDelaunayTriangulation_r( context, mGeos.get() ) );
3360
3361 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3362 {
3363 return nullptr;
3364 }
3365
3366 std::unique_ptr< QgsAbstractGeometry > res = fromGeos( geos.get() );
3367 if ( const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( res.get() ) )
3368 {
3369 return std::unique_ptr< QgsAbstractGeometry >( collection->extractPartsByType( Qgis::WkbType::Polygon, true ) );
3370 }
3371 else
3372 {
3373 return res;
3374 }
3375 }
3376 CATCH_GEOS_WITH_ERRMSG( nullptr )
3377#endif
3378}
3379
3381static bool _linestringEndpoints( const GEOSGeometry *linestring, double &x1, double &y1, double &x2, double &y2 )
3382{
3383 GEOSContextHandle_t context = QgsGeosContext::get();
3384 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( context, linestring );
3385 if ( !coordSeq )
3386 return false;
3387
3388 unsigned int coordSeqSize;
3389 if ( GEOSCoordSeq_getSize_r( context, coordSeq, &coordSeqSize ) == 0 )
3390 return false;
3391
3392 if ( coordSeqSize < 2 )
3393 return false;
3394
3395 GEOSCoordSeq_getX_r( context, coordSeq, 0, &x1 );
3396 GEOSCoordSeq_getY_r( context, coordSeq, 0, &y1 );
3397 GEOSCoordSeq_getX_r( context, coordSeq, coordSeqSize - 1, &x2 );
3398 GEOSCoordSeq_getY_r( context, coordSeq, coordSeqSize - 1, &y2 );
3399 return true;
3400}
3401
3402
3404static geos::unique_ptr _mergeLinestrings( const GEOSGeometry *line1, const GEOSGeometry *line2, const QgsPointXY &intersectionPoint )
3405{
3406 double x1, y1, x2, y2;
3407 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
3408 return nullptr;
3409
3410 double rx1, ry1, rx2, ry2;
3411 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
3412 return nullptr;
3413
3414 bool intersectionAtOrigLineEndpoint =
3415 ( intersectionPoint.x() == x1 && intersectionPoint.y() == y1 ) !=
3416 ( intersectionPoint.x() == x2 && intersectionPoint.y() == y2 );
3417 bool intersectionAtReshapeLineEndpoint =
3418 ( intersectionPoint.x() == rx1 && intersectionPoint.y() == ry1 ) ||
3419 ( intersectionPoint.x() == rx2 && intersectionPoint.y() == ry2 );
3420
3421 GEOSContextHandle_t context = QgsGeosContext::get();
3422 // the intersection must be at the begin/end of both lines
3423 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
3424 {
3425 geos::unique_ptr g1( GEOSGeom_clone_r( context, line1 ) );
3426 geos::unique_ptr g2( GEOSGeom_clone_r( context, line2 ) );
3427 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
3428 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, geoms, 2 ) );
3429 geos::unique_ptr res( GEOSLineMerge_r( context, multiGeom.get() ) );
3430
3431 //keep the original orientation if the result has a start or end point in common with the original line
3432 //and this point is not the start or the end point for both lines
3433 double x1res, y1res, x2res, y2res;
3434 if ( !_linestringEndpoints( res.get(), x1res, y1res, x2res, y2res ) )
3435 return nullptr;
3436 if ( ( x1res == x2 && y1res == y2 ) || ( x2res == x1 && y2res == y1 ) )
3437 res.reset( GEOSReverse_r( context, res.get() ) );
3438
3439 return res;
3440 }
3441 else
3442 return nullptr;
3443}
3444
3445
3446geos::unique_ptr QgsGeos::reshapeLine( const GEOSGeometry *line, const GEOSGeometry *reshapeLineGeos, double precision )
3447{
3448 if ( !line || !reshapeLineGeos )
3449 return nullptr;
3450
3451 bool atLeastTwoIntersections = false;
3452 bool oneIntersection = false;
3453 QgsPointXY oneIntersectionPoint;
3454
3455 GEOSContextHandle_t context = QgsGeosContext::get();
3456 try
3457 {
3458 //make sure there are at least two intersection between line and reshape geometry
3459 geos::unique_ptr intersectGeom( GEOSIntersection_r( context, line, reshapeLineGeos ) );
3460 if ( intersectGeom )
3461 {
3462 const int geomType = GEOSGeomTypeId_r( context, intersectGeom.get() );
3463 atLeastTwoIntersections = ( geomType == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 1 )
3464 || ( geomType == GEOS_GEOMETRYCOLLECTION && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 ) // a collection implies at least two points!
3465 || ( geomType == GEOS_MULTILINESTRING && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 );
3466 // one point is enough when extending line at its endpoint
3467 if ( GEOSGeomTypeId_r( context, intersectGeom.get() ) == GEOS_POINT )
3468 {
3469 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( context, intersectGeom.get() );
3470 double xi, yi;
3471 GEOSCoordSeq_getX_r( context, intersectionCoordSeq, 0, &xi );
3472 GEOSCoordSeq_getY_r( context, intersectionCoordSeq, 0, &yi );
3473 oneIntersection = true;
3474 oneIntersectionPoint = QgsPointXY( xi, yi );
3475 }
3476 }
3477 }
3478 catch ( QgsGeosException & )
3479 {
3480 atLeastTwoIntersections = false;
3481 }
3482
3483 // special case when extending line at its endpoint
3484 if ( oneIntersection )
3485 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
3486
3487 if ( !atLeastTwoIntersections )
3488 return nullptr;
3489
3490 //begin and end point of original line
3491 double x1, y1, x2, y2;
3492 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
3493 return nullptr;
3494
3495 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1, false, 0, false, 0, 2, precision );
3496 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2, false, 0, false, 0, 2, precision );
3497
3498 bool isRing = false;
3499 if ( GEOSGeomTypeId_r( context, line ) == GEOS_LINEARRING
3500 || GEOSEquals_r( context, beginLineVertex.get(), endLineVertex.get() ) == 1 )
3501 isRing = true;
3502
3503 //node line and reshape line
3504 geos::unique_ptr nodedGeometry = nodeGeometries( reshapeLineGeos, line );
3505 if ( !nodedGeometry )
3506 {
3507 return nullptr;
3508 }
3509
3510 //and merge them together
3511 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, nodedGeometry.get() ) );
3512 if ( !mergedLines )
3513 {
3514 return nullptr;
3515 }
3516
3517 int numMergedLines = GEOSGetNumGeometries_r( context, mergedLines.get() );
3518 if ( numMergedLines < 2 ) //some special cases. Normally it is >2
3519 {
3520 if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
3521 {
3522 geos::unique_ptr result( GEOSGeom_clone_r( context, reshapeLineGeos ) );
3523 return result;
3524 }
3525 else
3526 return nullptr;
3527 }
3528
3529 QVector<GEOSGeometry *> resultLineParts; //collection with the line segments that will be contained in result
3530 QVector<GEOSGeometry *> probableParts; //parts where we can decide on inclusion only after going through all the candidates
3531
3532 for ( int i = 0; i < numMergedLines; ++i )
3533 {
3534 const GEOSGeometry *currentGeom = GEOSGetGeometryN_r( context, mergedLines.get(), i );
3535
3536 // have we already added this part?
3537 bool alreadyAdded = false;
3538 double distance = 0;
3539 double bufferDistance = std::pow( 10.0L, geomDigits( currentGeom ) - 11 );
3540 for ( const GEOSGeometry *other : std::as_const( resultLineParts ) )
3541 {
3542 GEOSHausdorffDistance_r( context, currentGeom, other, &distance );
3543 if ( distance < bufferDistance )
3544 {
3545 alreadyAdded = true;
3546 break;
3547 }
3548 }
3549 if ( alreadyAdded )
3550 continue;
3551
3552 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( context, currentGeom );
3553 unsigned int currentCoordSeqSize;
3554 GEOSCoordSeq_getSize_r( context, currentCoordSeq, &currentCoordSeqSize );
3555 if ( currentCoordSeqSize < 2 )
3556 continue;
3557
3558 //get the two endpoints of the current line merge result
3559 double xBegin, xEnd, yBegin, yEnd;
3560 GEOSCoordSeq_getX_r( context, currentCoordSeq, 0, &xBegin );
3561 GEOSCoordSeq_getY_r( context, currentCoordSeq, 0, &yBegin );
3562 GEOSCoordSeq_getX_r( context, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
3563 GEOSCoordSeq_getY_r( context, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
3564 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin, false, 0, false, 0, 2, precision );
3565 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd, false, 0, false, 0, 2, precision );
3566
3567 //check how many endpoints of the line merge result are on the (original) line
3568 int nEndpointsOnOriginalLine = 0;
3569 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
3570 nEndpointsOnOriginalLine += 1;
3571
3572 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
3573 nEndpointsOnOriginalLine += 1;
3574
3575 //check how many endpoints equal the endpoints of the original line
3576 int nEndpointsSameAsOriginalLine = 0;
3577 if ( GEOSEquals_r( context, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3578 || GEOSEquals_r( context, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3579 nEndpointsSameAsOriginalLine += 1;
3580
3581 if ( GEOSEquals_r( context, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3582 || GEOSEquals_r( context, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3583 nEndpointsSameAsOriginalLine += 1;
3584
3585 //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
3586 bool currentGeomOverlapsOriginalGeom = false;
3587 bool currentGeomOverlapsReshapeLine = false;
3588 if ( lineContainedInLine( currentGeom, line ) == 1 )
3589 currentGeomOverlapsOriginalGeom = true;
3590
3591 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
3592 currentGeomOverlapsReshapeLine = true;
3593
3594 //logic to decide if this part belongs to the result
3595 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3596 {
3597 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3598 }
3599 //for closed rings, we take one segment from the candidate list
3600 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3601 {
3602 probableParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3603 }
3604 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3605 {
3606 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3607 }
3608 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3609 {
3610 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3611 }
3612 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
3613 {
3614 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3615 }
3616 }
3617
3618 //add the longest segment from the probable list for rings (only used for polygon rings)
3619 if ( isRing && !probableParts.isEmpty() )
3620 {
3621 geos::unique_ptr maxGeom; //the longest geometry in the probabla list
3622 GEOSGeometry *currentGeom = nullptr;
3623 double maxLength = -std::numeric_limits<double>::max();
3624 double currentLength = 0;
3625 for ( int i = 0; i < probableParts.size(); ++i )
3626 {
3627 currentGeom = probableParts.at( i );
3628 GEOSLength_r( context, currentGeom, &currentLength );
3629 if ( currentLength > maxLength )
3630 {
3631 maxLength = currentLength;
3632 maxGeom.reset( currentGeom );
3633 }
3634 else
3635 {
3636 GEOSGeom_destroy_r( context, currentGeom );
3637 }
3638 }
3639 resultLineParts.push_back( maxGeom.release() );
3640 }
3641
3642 geos::unique_ptr result;
3643 if ( resultLineParts.empty() )
3644 return nullptr;
3645
3646 if ( resultLineParts.size() == 1 ) //the whole result was reshaped
3647 {
3648 result.reset( resultLineParts[0] );
3649 }
3650 else //>1
3651 {
3652 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
3653 for ( int i = 0; i < resultLineParts.size(); ++i )
3654 {
3655 lineArray[i] = resultLineParts[i];
3656 }
3657
3658 //create multiline from resultLineParts
3659 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
3660 delete [] lineArray;
3661
3662 //then do a linemerge with the newly combined partstrings
3663 result.reset( GEOSLineMerge_r( context, multiLineGeom.get() ) );
3664 }
3665
3666 //now test if the result is a linestring. Otherwise something went wrong
3667 if ( GEOSGeomTypeId_r( context, result.get() ) != GEOS_LINESTRING )
3668 {
3669 return nullptr;
3670 }
3671
3672 //keep the original orientation
3673 bool reverseLine = false;
3674 if ( isRing )
3675 {
3676 //for closed linestring check clockwise/counter-clockwise
3677 char isResultCCW = 0, isOriginCCW = 0;
3678 if ( GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, result.get() ), &isResultCCW ) &&
3679 GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, line ), &isOriginCCW )
3680 )
3681 {
3682 //reverse line if orientations are different
3683 reverseLine = ( isOriginCCW == 1 && isResultCCW == 0 ) || ( isOriginCCW == 0 && isResultCCW == 1 );
3684 }
3685 }
3686 else
3687 {
3688 //for linestring, check if the result has a start or end point in common with the original line
3689 double x1res, y1res, x2res, y2res;
3690 if ( !_linestringEndpoints( result.get(), x1res, y1res, x2res, y2res ) )
3691 return nullptr;
3692 geos::unique_ptr beginResultLineVertex = createGeosPointXY( x1res, y1res, false, 0, false, 0, 2, precision );
3693 geos::unique_ptr endResultLineVertex = createGeosPointXY( x2res, y2res, false, 0, false, 0, 2, precision );
3694 reverseLine = GEOSEquals_r( context, beginLineVertex.get(), endResultLineVertex.get() ) == 1 || GEOSEquals_r( context, endLineVertex.get(), beginResultLineVertex.get() ) == 1;
3695 }
3696 if ( reverseLine )
3697 result.reset( GEOSReverse_r( context, result.get() ) );
3698
3699 return result;
3700}
3701
3702geos::unique_ptr QgsGeos::reshapePolygon( const GEOSGeometry *polygon, const GEOSGeometry *reshapeLineGeos, double precision )
3703{
3704 //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
3705 int nIntersections = 0;
3706 int lastIntersectingRing = -2;
3707 const GEOSGeometry *lastIntersectingGeom = nullptr;
3708
3709 GEOSContextHandle_t context = QgsGeosContext::get();
3710 int nRings = GEOSGetNumInteriorRings_r( context, polygon );
3711 if ( nRings < 0 )
3712 return nullptr;
3713
3714 //does outer ring intersect?
3715 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( context, polygon );
3716 if ( GEOSIntersects_r( context, outerRing, reshapeLineGeos ) == 1 )
3717 {
3718 ++nIntersections;
3719 lastIntersectingRing = -1;
3720 lastIntersectingGeom = outerRing;
3721 }
3722
3723 //do inner rings intersect?
3724 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
3725
3726 try
3727 {
3728 for ( int i = 0; i < nRings; ++i )
3729 {
3730 innerRings[i] = GEOSGetInteriorRingN_r( context, polygon, i );
3731 if ( GEOSIntersects_r( context, innerRings[i], reshapeLineGeos ) == 1 )
3732 {
3733 ++nIntersections;
3734 lastIntersectingRing = i;
3735 lastIntersectingGeom = innerRings[i];
3736 }
3737 }
3738 }
3739 catch ( QgsGeosException & )
3740 {
3741 nIntersections = 0;
3742 }
3743
3744 if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
3745 {
3746 delete [] innerRings;
3747 return nullptr;
3748 }
3749
3750 //we have one intersecting ring, let's try to reshape it
3751 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos, precision );
3752 if ( !reshapeResult )
3753 {
3754 delete [] innerRings;
3755 return nullptr;
3756 }
3757
3758 //if reshaping took place, we need to reassemble the polygon and its rings
3759 GEOSGeometry *newRing = nullptr;
3760 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( context, reshapeResult.get() );
3761 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( context, reshapeSequence );
3762
3763 reshapeResult.reset();
3764
3765 try
3766 {
3767 newRing = GEOSGeom_createLinearRing_r( context, newCoordSequence );
3768 }
3769 catch ( QgsGeosException & )
3770 {
3771 // nothing to do: on exception newRing will be null
3772 }
3773
3774 if ( !newRing )
3775 {
3776 delete [] innerRings;
3777 return nullptr;
3778 }
3779
3780 GEOSGeometry *newOuterRing = nullptr;
3781 if ( lastIntersectingRing == -1 )
3782 newOuterRing = newRing;
3783 else
3784 newOuterRing = GEOSGeom_clone_r( context, outerRing );
3785
3786 //check if all the rings are still inside the outer boundary
3787 QVector<GEOSGeometry *> ringList;
3788 if ( nRings > 0 )
3789 {
3790 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( context, GEOSGeom_clone_r( context, newOuterRing ), nullptr, 0 );
3791 if ( outerRingPoly )
3792 {
3793 ringList.reserve( nRings );
3794 GEOSGeometry *currentRing = nullptr;
3795 for ( int i = 0; i < nRings; ++i )
3796 {
3797 if ( lastIntersectingRing == i )
3798 currentRing = newRing;
3799 else
3800 currentRing = GEOSGeom_clone_r( context, innerRings[i] );
3801
3802 //possibly a ring is no longer contained in the result polygon after reshape
3803 if ( GEOSContains_r( context, outerRingPoly, currentRing ) == 1 )
3804 ringList.push_back( currentRing );
3805 else
3806 GEOSGeom_destroy_r( context, currentRing );
3807 }
3808 }
3809 GEOSGeom_destroy_r( context, outerRingPoly );
3810 }
3811
3812 GEOSGeometry **newInnerRings = new GEOSGeometry*[ringList.size()];
3813 for ( int i = 0; i < ringList.size(); ++i )
3814 newInnerRings[i] = ringList.at( i );
3815
3816 delete [] innerRings;
3817
3818 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( context, newOuterRing, newInnerRings, ringList.size() ) );
3819 delete[] newInnerRings;
3820
3821 return reshapedPolygon;
3822}
3823
3824int QgsGeos::lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 )
3825{
3826 if ( !line1 || !line2 )
3827 {
3828 return -1;
3829 }
3830
3831 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
3832
3833 GEOSContextHandle_t context = QgsGeosContext::get();
3834 geos::unique_ptr bufferGeom( GEOSBuffer_r( context, line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ) );
3835 if ( !bufferGeom )
3836 return -2;
3837
3838 geos::unique_ptr intersectionGeom( GEOSIntersection_r( context, bufferGeom.get(), line1 ) );
3839
3840 //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
3841 double intersectGeomLength;
3842 double line1Length;
3843
3844 GEOSLength_r( context, intersectionGeom.get(), &intersectGeomLength );
3845 GEOSLength_r( context, line1, &line1Length );
3846
3847 double intersectRatio = line1Length / intersectGeomLength;
3848 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
3849 return 1;
3850
3851 return 0;
3852}
3853
3854int QgsGeos::pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line )
3855{
3856 if ( !point || !line )
3857 return -1;
3858
3859 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
3860
3861 GEOSContextHandle_t context = QgsGeosContext::get();
3862 geos::unique_ptr lineBuffer( GEOSBuffer_r( context, line, bufferDistance, 8 ) );
3863 if ( !lineBuffer )
3864 return -2;
3865
3866 bool contained = false;
3867 if ( GEOSContains_r( context, lineBuffer.get(), point ) == 1 )
3868 contained = true;
3869
3870 return contained;
3871}
3872
3873int QgsGeos::geomDigits( const GEOSGeometry *geom )
3874{
3875 GEOSContextHandle_t context = QgsGeosContext::get();
3876 geos::unique_ptr bbox( GEOSEnvelope_r( context, geom ) );
3877 if ( !bbox )
3878 return -1;
3879
3880 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( context, bbox.get() );
3881 if ( !bBoxRing )
3882 return -1;
3883
3884 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( context, bBoxRing );
3885
3886 if ( !bBoxCoordSeq )
3887 return -1;
3888
3889 unsigned int nCoords = 0;
3890 if ( !GEOSCoordSeq_getSize_r( context, bBoxCoordSeq, &nCoords ) )
3891 return -1;
3892
3893 int maxDigits = -1;
3894 for ( unsigned int i = 0; i < nCoords - 1; ++i )
3895 {
3896 double t;
3897 GEOSCoordSeq_getX_r( context, bBoxCoordSeq, i, &t );
3898
3899 int digits;
3900 digits = std::ceil( std::log10( std::fabs( t ) ) );
3901 if ( digits > maxDigits )
3902 maxDigits = digits;
3903
3904 GEOSCoordSeq_getY_r( context, bBoxCoordSeq, i, &t );
3905 digits = std::ceil( std::log10( std::fabs( t ) ) );
3906 if ( digits > maxDigits )
3907 maxDigits = digits;
3908 }
3909
3910 return maxDigits;
3911}
Provides global constants and enumerations for use throughout the application.
Definition qgis.h:54
BufferSide
Side of line to buffer.
Definition qgis.h:2059
@ Right
Buffer to right of line.
GeometryOperationResult
Success or failure of a geometry operation.
Definition qgis.h:2005
@ AddPartNotMultiGeometry
The source geometry is not multi.
@ InvalidBaseGeometry
The base geometry on which the operation is done is invalid or empty.
@ RejectOnInvalidSubGeometry
Don't allow geometries with invalid sub-geometries to be created.
@ SkipEmptyInteriorRings
Skip any empty polygon interior ring.
QFlags< GeosCreationFlag > GeosCreationFlags
Geos geometry creation behavior flags.
Definition qgis.h:2108
@ Polygon
Polygons.
@ Unknown
Unknown types.
@ Null
No geometry.
JoinStyle
Join styles for buffers.
Definition qgis.h:2084
EndCapStyle
End cap styles for buffers.
Definition qgis.h:2071
CoverageValidityResult
Coverage validity results.
Definition qgis.h:2117
@ Valid
Coverage is valid.
@ Invalid
Coverage is invalid. Invalidity includes polygons that overlap, that have gaps smaller than the gap w...
@ Error
An exception occurred while determining validity.
MakeValidMethod
Algorithms to use when repairing invalid geometries.
Definition qgis.h:2130
@ Linework
Combines all rings into a set of noded lines and then extracts valid polygons from that linework.
@ Structure
Structured method, first makes all rings valid and then merges shells and subtracts holes from shells...
WkbType
The WKB type describes the number of dimensions a geometry has.
Definition qgis.h:256
@ LineStringZM
LineStringZM.
@ Polygon
Polygon.
@ PointM
PointM.
@ PointZ
PointZ.
@ GeometryCollection
GeometryCollection.
@ PointZM
PointZM.
@ LineStringZ
LineStringZ.
@ PolyhedralSurface
PolyhedralSurface.
Abstract base class for all geometries.
virtual const QgsAbstractGeometry * simplifiedTypeRef() const
Returns a reference to the simplest lossless representation of this geometry, e.g.
bool isMeasure() const
Returns true if the geometry contains m values.
virtual QgsRectangle boundingBox() const
Returns the minimal bounding box for the geometry.
bool is3D() const
Returns true if the geometry is 3D and contains a z-value.
virtual QgsPoint vertexAt(QgsVertexId id) const =0
Returns the point corresponding to a specified vertex id.
Qgis::WkbType wkbType() const
Returns the WKB type of the geometry.
virtual bool isEmpty() const
Returns true if the geometry is empty.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
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.
Abstract base class for curved geometry type.
Definition qgscurve.h:35
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.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
static Qgis::GeometryOperationResult addPart(QgsAbstractGeometry *geometry, std::unique_ptr< QgsAbstractGeometry > part)
Add a part to multi type geometry.
A geometry engine is a low-level representation of a QgsAbstractGeometry object, optimised for use wi...
const QgsAbstractGeometry * mGeometry
EngineOperationResult
Success or failure of a geometry operation.
@ NothingHappened
Nothing happened, without any error.
@ InvalidBaseGeometry
The geometry on which the operation occurs is not valid.
@ InvalidInput
The input is not valid.
@ NodedGeometryError
Error occurred while creating a noded geometry.
@ EngineError
Error occurred in the geometry engine.
@ SplitCannotSplitPoint
Points cannot be split.
@ Success
Operation succeeded.
void logError(const QString &engineName, const QString &message) const
Logs an error message encountered during an operation.
static std::unique_ptr< QgsGeometryCollection > createCollectionOfType(Qgis::WkbType type)
Returns a new geometry collection matching a specified WKB type.
Encapsulates parameters under which a geometry operation is performed.
double gridSize() const
Returns the grid size which will be used to snap vertices of a geometry.
static double sqrDistance2D(double x1, double y1, double x2, double y2)
Returns the squared 2D distance between (x1, y1) and (x2, y2).
A geometry is the spatial representation of a feature.
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
const QgsAbstractGeometry * constGet() const
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Used to create and store a proj context object, correctly freeing the context upon destruction.
Definition qgsgeos.h:44
static GEOSContextHandle_t get()
Returns a thread local instance of a GEOS context, safe for use in the current thread.
Does vector analysis using the GEOS library and handles import, export, and exception handling.
Definition qgsgeos.h:139
double frechetDistanceDensify(const QgsAbstractGeometry *geometry, double densifyFraction, QString *errorMsg=nullptr) const
Returns the Fréchet distance between this geometry and another geometry, restricted to discrete point...
Definition qgsgeos.cpp:826
double hausdorffDistanceDensify(const QgsAbstractGeometry *geometry, double densifyFraction, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and another geometry.
Definition qgsgeos.cpp:780
bool touches(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom touches this.
Definition qgsgeos.cpp:882
static geos::unique_ptr offsetCurve(const GEOSGeometry *geometry, double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg=nullptr)
Directly calculates the offset curve for a GEOS geometry object and returns a GEOS geometry result.
Definition qgsgeos.cpp:2760
QgsAbstractGeometry * symDifference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const override
Calculate the symmetric difference of this and geom.
Definition qgsgeos.cpp:539
bool isValid(QString *errorMsg=nullptr, bool allowSelfTouchingHoles=false, QgsGeometry *errorLoc=nullptr) const override
Returns true if the geometry is valid.
Definition qgsgeos.cpp:2353
std::unique_ptr< QgsAbstractGeometry > subdivide(int maxNodes, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const
Subdivides the geometry.
Definition qgsgeos.cpp:448
QgsAbstractGeometry * intersection(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const override
Calculate the intersection of this and geom.
Definition qgsgeos.cpp:317
std::unique_ptr< QgsAbstractGeometry > constrainedDelaunayTriangulation(QString *errorMsg=nullptr) const
Returns a constrained Delaunay triangulation for the vertices of the geometry.
Definition qgsgeos.cpp:3344
QgsAbstractGeometry * buffer(double distance, int segments, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2087
std::unique_ptr< QgsAbstractGeometry > closestPoint(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the closest point on the geometry to the other geometry.
Definition qgsgeos.cpp:3086
std::unique_ptr< QgsAbstractGeometry > reshapeGeometry(const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg=nullptr) const
Reshapes the geometry using a line.
Definition qgsgeos.cpp:2936
std::unique_ptr< QgsAbstractGeometry > maximumInscribedCircle(double tolerance, QString *errorMsg=nullptr) const
Returns the maximum inscribed circle.
Definition qgsgeos.cpp:2814
static geos::unique_ptr asGeos(const QgsGeometry &geometry, double precision=0, Qgis::GeosCreationFlags flags=Qgis::GeosCreationFlags())
Returns a geos geometry - caller takes ownership of the object (should be deleted with GEOSGeom_destr...
Definition qgsgeos.cpp:257
QgsAbstractGeometry * difference(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const override
Calculate the difference of this and geom.
Definition qgsgeos.cpp:322
EngineOperationResult splitGeometry(const QgsLineString &splitLine, QVector< QgsGeometry > &newGeometries, bool topological, QgsPointSequence &topologyTestPoints, QString *errorMsg=nullptr, bool skipIntersectionCheck=false) const override
Splits this geometry according to a given line.
Definition qgsgeos.cpp:1040
std::unique_ptr< QgsAbstractGeometry > sharedPaths(const QgsAbstractGeometry *other, QString *errorMsg=nullptr) const
Find paths shared between the two given lineal geometries (this and other).
Definition qgsgeos.cpp:2916
std::unique_ptr< QgsAbstractGeometry > clip(const QgsRectangle &rectangle, QString *errorMsg=nullptr) const
Performs a fast, non-robust intersection between the geometry and a rectangle.
Definition qgsgeos.cpp:327
std::unique_ptr< QgsAbstractGeometry > node(QString *errorMsg=nullptr) const
Returns a (Multi)LineString representing the fully noded version of a collection of linestrings.
Definition qgsgeos.cpp:2900
double minimumClearance(QString *errorMsg=nullptr) const
Computes the minimum clearance of a geometry.
Definition qgsgeos.cpp:2866
std::unique_ptr< QgsAbstractGeometry > singleSidedBuffer(double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg=nullptr) const
Returns a single sided buffer for a geometry.
Definition qgsgeos.cpp:2787
bool disjoint(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is disjoint from this.
Definition qgsgeos.cpp:935
std::unique_ptr< QgsAbstractGeometry > mergeLines(QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
Definition qgsgeos.cpp:3059
std::unique_ptr< QgsAbstractGeometry > makeValid(Qgis::MakeValidMethod method=Qgis::MakeValidMethod::Linework, bool keepCollapsed=false, QString *errorMsg=nullptr) const
Repairs the geometry using GEOS make valid routine.
Definition qgsgeos.cpp:198
std::unique_ptr< QgsAbstractGeometry > voronoiDiagram(const QgsAbstractGeometry *extent=nullptr, double tolerance=0.0, bool edgesOnly=false, QString *errorMsg=nullptr) const
Creates a Voronoi diagram for the nodes contained within the geometry.
Definition qgsgeos.cpp:3288
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2125
std::unique_ptr< QgsAbstractGeometry > unionCoverage(QString *errorMsg=nullptr) const
Optimized union algorithm for polygonal inputs that are correctly noded and do not overlap.
Definition qgsgeos.cpp:2335
static std::unique_ptr< QgsPolygon > fromGeosPolygon(const GEOSGeometry *geos)
Definition qgsgeos.cpp:1677
double hausdorffDistance(const QgsAbstractGeometry *geometry, QString *errorMsg=nullptr) const
Returns the Hausdorff distance between this geometry and another geometry.
Definition qgsgeos.cpp:757
std::unique_ptr< QgsAbstractGeometry > minimumClearanceLine(QString *errorMsg=nullptr) const
Returns a LineString whose endpoints define the minimum clearance of a geometry.
Definition qgsgeos.cpp:2884
QgsAbstractGeometry * envelope(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2182
QString relate(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Returns the Dimensional Extended 9 Intersection Model (DE-9IM) representation of the relationship bet...
Definition qgsgeos.cpp:940
bool crosses(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom crosses this.
Definition qgsgeos.cpp:887
QgsAbstractGeometry * convexHull(QString *errorMsg=nullptr) const override
Calculate the convex hull of this.
Definition qgsgeos.cpp:2226
double distance(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Calculates the distance between this and geom.
Definition qgsgeos.cpp:581
std::unique_ptr< QgsAbstractGeometry > minimumWidth(QString *errorMsg=nullptr) const
Returns a linestring geometry which represents the minimum diameter of the geometry.
Definition qgsgeos.cpp:2850
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
Definition qgsgeos.cpp:2450
QgsGeos(const QgsAbstractGeometry *geometry, double precision=0, Qgis::GeosCreationFlags flags=Qgis::GeosCreationFlag::SkipEmptyInteriorRings)
GEOS geometry engine constructor.
Definition qgsgeos.cpp:177
QgsAbstractGeometry * combine(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const override
Calculate the combination of this and geom.
Definition qgsgeos.cpp:468
bool within(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom is within this.
Definition qgsgeos.cpp:892
std::unique_ptr< QgsAbstractGeometry > shortestLine(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the shortest line joining this geometry to the other geometry.
Definition qgsgeos.cpp:3130
std::unique_ptr< QgsAbstractGeometry > concaveHull(double targetPercent, bool allowHoles=false, QString *errorMsg=nullptr) const
Returns a possibly concave geometry that encloses the input geometry.
Definition qgsgeos.cpp:2242
void prepareGeometry() override
Prepares the geometry, so that subsequent calls to spatial relation methods are much faster.
Definition qgsgeos.cpp:289
static std::unique_ptr< QgsAbstractGeometry > fromGeos(const GEOSGeometry *geos)
Create a geometry from a GEOSGeometry.
Definition qgsgeos.cpp:1568
QgsAbstractGeometry * interpolate(double distance, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2140
bool distanceWithin(const QgsAbstractGeometry *geom, double maxdistance, QString *errorMsg=nullptr) const override
Checks if geom is within maxdistance distance from this geometry.
Definition qgsgeos.cpp:656
bool isEqual(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if this is equal to geom.
Definition qgsgeos.cpp:2416
bool contains(double x, double y, QString *errorMsg=nullptr) const
Returns true if the geometry contains the point at (x, y).
Definition qgsgeos.cpp:714
static QgsGeometry polygonize(const QVector< const QgsAbstractGeometry * > &geometries, QString *errorMsg=nullptr)
Creates a GeometryCollection geometry containing possible polygons formed from the constituent linewo...
Definition qgsgeos.cpp:3248
double frechetDistance(const QgsAbstractGeometry *geometry, QString *errorMsg=nullptr) const
Returns the Fréchet distance between this geometry and another geometry, restricted to discrete point...
Definition qgsgeos.cpp:803
bool relatePattern(const QgsAbstractGeometry *geom, const QString &pattern, QString *errorMsg=nullptr) const override
Tests whether two geometries are related by a specified Dimensional Extended 9 Intersection Model (DE...
Definition qgsgeos.cpp:976
QgsPoint * centroid(QString *errorMsg=nullptr) const override
Calculates the centroid of this.
Definition qgsgeos.cpp:2155
double lineLocatePoint(const QgsPoint &point, QString *errorMsg=nullptr) const
Returns a distance representing the location along this linestring of the closest point on this lines...
Definition qgsgeos.cpp:3188
std::unique_ptr< QgsAbstractGeometry > simplifyCoverageVW(double tolerance, bool preserveBoundary, QString *errorMsg=nullptr) const
Operates on a coverage (represented as a list of polygonal geometry with exactly matching edge geomet...
Definition qgsgeos.cpp:2310
bool isEmpty(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2436
static Qgis::GeometryOperationResult addPart(QgsGeometry &geometry, GEOSGeometry *newPart)
Adds a new island polygon to a multipolygon feature.
Definition qgsgeos.cpp:267
bool intersects(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom intersects this.
Definition qgsgeos.cpp:849
void geometryChanged() override
Should be called whenever the geometry associated with the engine has been modified and the engine mu...
Definition qgsgeos.cpp:282
bool overlaps(const QgsAbstractGeometry *geom, QString *errorMsg=nullptr) const override
Checks if geom overlaps this.
Definition qgsgeos.cpp:897
double area(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:1007
std::unique_ptr< QgsAbstractGeometry > delaunayTriangulation(double tolerance=0.0, bool edgesOnly=false, QString *errorMsg=nullptr) const
Returns the Delaunay triangulation for the vertices of the geometry.
Definition qgsgeos.cpp:3321
double length(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:1024
std::unique_ptr< QgsAbstractGeometry > largestEmptyCircle(double tolerance, const QgsAbstractGeometry *boundary=nullptr, QString *errorMsg=nullptr) const
Constructs the Largest Empty Circle for a set of obstacle geometries, up to a specified tolerance.
Definition qgsgeos.cpp:2830
Qgis::CoverageValidityResult validateCoverage(double gapWidth, std::unique_ptr< QgsAbstractGeometry > *invalidEdges, QString *errorMsg=nullptr) const
Analyze a coverage (represented as a collection of polygonal geometry with exactly matching edge geom...
Definition qgsgeos.cpp:2265
static QgsPoint coordSeqPoint(const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM)
Definition qgsgeos.cpp:1770
QgsPoint * pointOnSurface(QString *errorMsg=nullptr) const override
Calculate a point that is guaranteed to be on the surface of this.
Definition qgsgeos.cpp:2197
static QgsGeometry geometryFromGeos(GEOSGeometry *geos)
Creates a new QgsGeometry object, feeding in a geometry in GEOS format.
Definition qgsgeos.cpp:185
Line string geometry type, with support for z-dimension and m-values.
const double * yData() const
Returns a const pointer to the y vertex data.
const double * xData() const
Returns a const pointer to the x vertex data.
QVector< double > xVector() const
Returns the x vertex values as a vector.
const double * zData() const
Returns a const pointer to the z vertex data, or nullptr if the linestring does not have z values.
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.
QVector< double > yVector() const
Returns the y vertex values as a vector.
double yAt(int index) const override
Returns the y-coordinate of the specified node in the line string.
QVector< double > zVector() const
Returns the z vertex values as a vector.
double closestSegment(const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf=nullptr, double epsilon=4 *std::numeric_limits< double >::epsilon()) const override
Searches for the closest segment of the geometry to a given point.
const double * mData() const
Returns a const pointer to the m vertex data, or nullptr if the linestring does not have m values.
void addVertex(const QgsPoint &pt)
Adds a new vertex to the end of the line string.
QgsLineString * clone() const override
Clones the geometry by performing a deep copy.
double xAt(int index) const override
Returns the x-coordinate of the specified node in the line string.
Multi curve geometry collection.
bool addGeometry(QgsAbstractGeometry *g) override
Adds a geometry and takes ownership. Returns true in case of success.
Custom exception class which is raised when an operation is not supported.
Represents a 2D point.
Definition qgspointxy.h:60
double y
Definition qgspointxy.h:64
double x
Definition qgspointxy.h:63
Point geometry type, with support for z-dimension and m-values.
Definition qgspoint.h:49
QgsPoint * clone() const override
Clones the geometry by performing a deep copy.
Definition qgspoint.cpp:105
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
bool isEmpty() const override
Returns true if the geometry is empty.
Definition qgspoint.cpp:736
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition qgspoint.h:387
double m
Definition qgspoint.h:55
double y
Definition qgspoint.h:53
Polygon geometry type.
Definition qgspolygon.h:33
Polyhedral surface geometry type.
int numPatches() const
Returns the number of patches contained with the polyhedral surface.
const QgsPolygon * patchN(int i) const
Retrieves a patch from the polyhedral surface.
A rectangle specified with double values.
double xMinimum
double yMinimum
double xMaximum
void setYMinimum(double y)
Set the minimum y value.
void setXMinimum(double x)
Set the minimum x value.
void setYMaximum(double y)
Set the maximum y value.
void setXMaximum(double x)
Set the maximum x value.
double yMaximum
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 Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB type.
static Q_INVOKABLE bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
Contains geos related utilities and functions.
Definition qgsgeos.h:75
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
bool qgsDoubleNear(double a, double b, double epsilon=4 *std::numeric_limits< double >::epsilon())
Compare two doubles (but allow some difference)
Definition qgis.h:6392
QVector< QgsPoint > QgsPointSequence
#define CATCH_GEOS(r)
Definition qgsgeos.cpp:35
#define DEFAULT_QUADRANT_SEGMENTS
Definition qgsgeos.cpp:33
#define CATCH_GEOS_WITH_ERRMSG(r)
Definition qgsgeos.cpp:41
#define QgsDebugError(str)
Definition qgslogger.h:40
QLineF segment(int index, QRectF rect, double radius)
int precision
Utility class for identifying a unique vertex within a geometry.
Definition qgsvertexid.h:30
int vertex
Vertex number.
Definition qgsvertexid.h:94
struct GEOSGeom_t GEOSGeometry
Definition util.h:41