QGIS API Documentation 3.41.0-Master (88383c3d16f)
Loading...
Searching...
No Matches
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 (GEOSException &) \
37 { \
38 return r; \
39 }
40
41#define CATCH_GEOS_WITH_ERRMSG(r) \
42 catch (GEOSException &e) \
43 { \
44 if ( errorMsg ) \
45 { \
46 *errorMsg = e.what(); \
47 } \
48 return r; \
49 }
50
52
53static void throwGEOSException( 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 GEOSException, 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 GEOSException is thrown.
69 // TODO - find a real fix for the underlying issue
70 try
71 {
72 throw GEOSException( 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 GEOSException( 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, throwGEOSException );
116#else
117 mContext = initGEOS_r( printGEOSNotice, throwGEOSException );
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 ( GEOSException &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 ( GEOSException &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 ( GEOSException &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 ( GEOSException &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 ( GEOSException &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 ( GEOSException &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 ( GEOSException &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 ( GEOSGeomTypeId_r( context, geom ) == GEOS_POLYGON || GEOSGeomTypeId_r( context, geom ) == GEOS_MULTIPOLYGON )
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 ( GEOSException & )
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 return nPoints > 0 ? std::unique_ptr<QgsAbstractGeometry>( coordSeqPoint( cs, 0, hasZ, hasM ).clone() ) : nullptr;
1592 }
1593 case GEOS_LINESTRING:
1594 {
1595 return sequenceToLinestring( geos, hasZ, hasM );
1596 }
1597 case GEOS_POLYGON:
1598 {
1599 return fromGeosPolygon( geos );
1600 }
1601 case GEOS_MULTIPOINT:
1602 {
1603 auto multiPoint = std::make_unique<QgsMultiPoint>();
1604 int nParts = GEOSGetNumGeometries_r( context, geos );
1605 multiPoint->reserve( nParts );
1606 for ( int i = 0; i < nParts; ++i )
1607 {
1608 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, GEOSGetGeometryN_r( context, geos, i ) );
1609 if ( cs )
1610 {
1611 unsigned int nPoints = 0;
1612 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1613 if ( nPoints > 0 )
1614 multiPoint->addGeometry( coordSeqPoint( cs, 0, hasZ, hasM ).clone() );
1615 }
1616 }
1617 return std::move( multiPoint );
1618 }
1619 case GEOS_MULTILINESTRING:
1620 {
1621 auto multiLineString = std::make_unique<QgsMultiLineString>();
1622 int nParts = GEOSGetNumGeometries_r( context, geos );
1623 multiLineString->reserve( nParts );
1624 for ( int i = 0; i < nParts; ++i )
1625 {
1626 std::unique_ptr< QgsLineString >line( sequenceToLinestring( GEOSGetGeometryN_r( context, geos, i ), hasZ, hasM ) );
1627 if ( line )
1628 {
1629 multiLineString->addGeometry( line.release() );
1630 }
1631 }
1632 return std::move( multiLineString );
1633 }
1634 case GEOS_MULTIPOLYGON:
1635 {
1636 auto multiPolygon = std::make_unique<QgsMultiPolygon>();
1637
1638 int nParts = GEOSGetNumGeometries_r( context, geos );
1639 multiPolygon->reserve( nParts );
1640 for ( int i = 0; i < nParts; ++i )
1641 {
1642 std::unique_ptr< QgsPolygon > poly = fromGeosPolygon( GEOSGetGeometryN_r( context, geos, i ) );
1643 if ( poly )
1644 {
1645 multiPolygon->addGeometry( poly.release() );
1646 }
1647 }
1648 return std::move( multiPolygon );
1649 }
1650 case GEOS_GEOMETRYCOLLECTION:
1651 {
1652 auto geomCollection = std::make_unique<QgsGeometryCollection>();
1653 int nParts = GEOSGetNumGeometries_r( context, geos );
1654 geomCollection->reserve( nParts );
1655 for ( int i = 0; i < nParts; ++i )
1656 {
1657 std::unique_ptr< QgsAbstractGeometry > geom( fromGeos( GEOSGetGeometryN_r( context, geos, i ) ) );
1658 if ( geom )
1659 {
1660 geomCollection->addGeometry( geom.release() );
1661 }
1662 }
1663 return std::move( geomCollection );
1664 }
1665 }
1666 return nullptr;
1667}
1668
1669std::unique_ptr<QgsPolygon> QgsGeos::fromGeosPolygon( const GEOSGeometry *geos )
1670{
1671 GEOSContextHandle_t context = QgsGeosContext::get();
1672 if ( GEOSGeomTypeId_r( context, geos ) != GEOS_POLYGON )
1673 {
1674 return nullptr;
1675 }
1676
1677 int nCoordDims = GEOSGeom_getCoordinateDimension_r( context, geos );
1678 int nDims = GEOSGeom_getDimensions_r( context, geos );
1679 bool hasZ = ( nCoordDims == 3 );
1680 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1681
1682 auto polygon = std::make_unique<QgsPolygon>();
1683
1684 const GEOSGeometry *ring = GEOSGetExteriorRing_r( context, geos );
1685 if ( ring )
1686 {
1687 polygon->setExteriorRing( sequenceToLinestring( ring, hasZ, hasM ).release() );
1688 }
1689
1690 QVector<QgsCurve *> interiorRings;
1691 const int ringCount = GEOSGetNumInteriorRings_r( context, geos );
1692 interiorRings.reserve( ringCount );
1693 for ( int i = 0; i < ringCount; ++i )
1694 {
1695 ring = GEOSGetInteriorRingN_r( context, geos, i );
1696 if ( ring )
1697 {
1698 interiorRings.push_back( sequenceToLinestring( ring, hasZ, hasM ).release() );
1699 }
1700 }
1701 polygon->setInteriorRings( interiorRings );
1702
1703 return polygon;
1704}
1705
1706std::unique_ptr<QgsLineString> QgsGeos::sequenceToLinestring( const GEOSGeometry *geos, bool hasZ, bool hasM )
1707{
1708 GEOSContextHandle_t context = QgsGeosContext::get();
1709 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( context, geos );
1710
1711 unsigned int nPoints;
1712 GEOSCoordSeq_getSize_r( context, cs, &nPoints );
1713
1714 QVector< double > xOut( nPoints );
1715 QVector< double > yOut( nPoints );
1716 QVector< double > zOut;
1717 if ( hasZ )
1718 zOut.resize( nPoints );
1719 QVector< double > mOut;
1720 if ( hasM )
1721 mOut.resize( nPoints );
1722
1723 double *x = xOut.data();
1724 double *y = yOut.data();
1725 double *z = zOut.data();
1726 double *m = mOut.data();
1727
1728#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
1729 GEOSCoordSeq_copyToArrays_r( context, cs, x, y, hasZ ? z : nullptr, hasM ? m : nullptr );
1730#else
1731 for ( unsigned int i = 0; i < nPoints; ++i )
1732 {
1733 if ( hasZ )
1734 GEOSCoordSeq_getXYZ_r( context, cs, i, x++, y++, z++ );
1735 else
1736 GEOSCoordSeq_getXY_r( context, cs, i, x++, y++ );
1737 if ( hasM )
1738 {
1739 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, m++ );
1740 }
1741 }
1742#endif
1743 auto line = std::make_unique<QgsLineString>( xOut, yOut, zOut, mOut );
1744 return line;
1745}
1746
1747int QgsGeos::numberOfGeometries( GEOSGeometry *g )
1748{
1749 if ( !g )
1750 return 0;
1751
1752 GEOSContextHandle_t context = QgsGeosContext::get();
1753 int geometryType = GEOSGeomTypeId_r( context, g );
1754 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
1755 || geometryType == GEOS_POLYGON )
1756 return 1;
1757
1758 //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
1759 return GEOSGetNumGeometries_r( context, g );
1760}
1761
1762QgsPoint QgsGeos::coordSeqPoint( const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM )
1763{
1764 if ( !cs )
1765 {
1766 return QgsPoint();
1767 }
1768
1769 GEOSContextHandle_t context = QgsGeosContext::get();
1770
1771 double x, y;
1772 double z = 0;
1773 double m = 0;
1774 if ( hasZ )
1775 GEOSCoordSeq_getXYZ_r( context, cs, i, &x, &y, &z );
1776 else
1777 GEOSCoordSeq_getXY_r( context, cs, i, &x, &y );
1778 if ( hasM )
1779 {
1780 GEOSCoordSeq_getOrdinate_r( context, cs, i, 3, &m );
1781 }
1782
1784 if ( hasZ && hasM )
1785 {
1787 }
1788 else if ( hasZ )
1789 {
1791 }
1792 else if ( hasM )
1793 {
1795 }
1796 return QgsPoint( t, x, y, z, m );
1797}
1798
1799geos::unique_ptr QgsGeos::asGeos( const QgsAbstractGeometry *geom, double precision, Qgis::GeosCreationFlags flags )
1800{
1801 if ( !geom )
1802 return nullptr;
1803
1804 int coordDims = 2;
1805 if ( geom->is3D() )
1806 {
1807 ++coordDims;
1808 }
1809 if ( geom->isMeasure() )
1810 {
1811 ++coordDims;
1812 }
1813
1815 {
1816 int geosType = GEOS_GEOMETRYCOLLECTION;
1817
1819 {
1820 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1821 {
1823 geosType = GEOS_MULTIPOINT;
1824 break;
1825
1827 geosType = GEOS_MULTILINESTRING;
1828 break;
1829
1831 geosType = GEOS_MULTIPOLYGON;
1832 break;
1833
1836 return nullptr;
1837 }
1838 }
1839
1840
1841 const QgsGeometryCollection *c = qgsgeometry_cast<const QgsGeometryCollection *>( geom );
1842
1843 if ( !c )
1844 return nullptr;
1845
1846 std::vector<geos::unique_ptr> geomVector;
1847 geomVector.reserve( c->numGeometries() );
1848 for ( int i = 0; i < c->numGeometries(); ++i )
1849 {
1850 geos::unique_ptr geosGeom = asGeos( c->geometryN( i ), precision, flags );
1851 if ( flags & Qgis::GeosCreationFlag::RejectOnInvalidSubGeometry && !geosGeom )
1852 {
1853 return nullptr;
1854 }
1855 geomVector.emplace_back( std::move( geosGeom ) );
1856 }
1857 return createGeosCollection( geosType, geomVector );
1858 }
1861 {
1862 // PolyhedralSurface and TIN support
1863 // convert it to a geos MultiPolygon
1864 const QgsPolyhedralSurface *polyhedralSurface = qgsgeometry_cast<const QgsPolyhedralSurface *>( geom );
1865 if ( !polyhedralSurface )
1866 return nullptr;
1867
1868 std::vector<geos::unique_ptr> geomVector;
1869 geomVector.reserve( polyhedralSurface->numPatches() );
1870 for ( int i = 0; i < polyhedralSurface->numPatches(); ++i )
1871 {
1872 geos::unique_ptr geosPolygon = createGeosPolygon( polyhedralSurface->patchN( i ), precision );
1873 if ( flags & Qgis::GeosCreationFlag::RejectOnInvalidSubGeometry && !geosPolygon )
1874 {
1875 return nullptr;
1876 }
1877 geomVector.emplace_back( std::move( geosPolygon ) );
1878 }
1879
1880 return createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
1881 }
1882 else
1883 {
1884 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1885 {
1887 return createGeosPoint( static_cast<const QgsPoint *>( geom ), coordDims, precision, flags );
1888
1890 return createGeosLinestring( static_cast<const QgsLineString *>( geom ), precision, flags );
1891
1893 return createGeosPolygon( static_cast<const QgsPolygon *>( geom ), precision, flags );
1894
1897 return nullptr;
1898 }
1899 }
1900 return nullptr;
1901}
1902
1903std::unique_ptr<QgsAbstractGeometry> QgsGeos::overlay( const QgsAbstractGeometry *geom, Overlay op, QString *errorMsg, const QgsGeometryParameters &parameters ) const
1904{
1905 if ( !mGeos || !geom )
1906 {
1907 return nullptr;
1908 }
1909
1910 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
1911 if ( !geosGeom )
1912 {
1913 return nullptr;
1914 }
1915
1916 const double gridSize = parameters.gridSize();
1917
1918 GEOSContextHandle_t context = QgsGeosContext::get();
1919 try
1920 {
1921 geos::unique_ptr opGeom;
1922 switch ( op )
1923 {
1924 case OverlayIntersection:
1925 if ( gridSize > 0 )
1926 {
1927 opGeom.reset( GEOSIntersectionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1928 }
1929 else
1930 {
1931 opGeom.reset( GEOSIntersection_r( context, mGeos.get(), geosGeom.get() ) );
1932 }
1933 break;
1934
1935 case OverlayDifference:
1936 if ( gridSize > 0 )
1937 {
1938 opGeom.reset( GEOSDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1939 }
1940 else
1941 {
1942 opGeom.reset( GEOSDifference_r( context, mGeos.get(), geosGeom.get() ) );
1943 }
1944 break;
1945
1946 case OverlayUnion:
1947 {
1948 geos::unique_ptr unionGeometry;
1949 if ( gridSize > 0 )
1950 {
1951 unionGeometry.reset( GEOSUnionPrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1952 }
1953 else
1954 {
1955 unionGeometry.reset( GEOSUnion_r( context, mGeos.get(), geosGeom.get() ) );
1956 }
1957
1958 if ( unionGeometry && GEOSGeomTypeId_r( context, unionGeometry.get() ) == GEOS_MULTILINESTRING )
1959 {
1960 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, unionGeometry.get() ) );
1961 if ( mergedLines )
1962 {
1963 unionGeometry = std::move( mergedLines );
1964 }
1965 }
1966
1967 opGeom = std::move( unionGeometry );
1968 }
1969 break;
1970
1971 case OverlaySymDifference:
1972 if ( gridSize > 0 )
1973 {
1974 opGeom.reset( GEOSSymDifferencePrec_r( context, mGeos.get(), geosGeom.get(), gridSize ) );
1975 }
1976 else
1977 {
1978 opGeom.reset( GEOSSymDifference_r( context, mGeos.get(), geosGeom.get() ) );
1979 }
1980 break;
1981 }
1982 return fromGeos( opGeom.get() );
1983 }
1984 catch ( GEOSException &e )
1985 {
1986 logError( QStringLiteral( "GEOS" ), e.what() );
1987 if ( errorMsg )
1988 {
1989 *errorMsg = e.what();
1990 }
1991 return nullptr;
1992 }
1993}
1994
1995bool QgsGeos::relation( const QgsAbstractGeometry *geom, Relation r, QString *errorMsg ) const
1996{
1997 if ( !mGeos || !geom )
1998 {
1999 return false;
2000 }
2001
2002 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2003 if ( !geosGeom )
2004 {
2005 return false;
2006 }
2007
2008 GEOSContextHandle_t context = QgsGeosContext::get();
2009 bool result = false;
2010 try
2011 {
2012 if ( mGeosPrepared ) //use faster version with prepared geometry
2013 {
2014 switch ( r )
2015 {
2016 case RelationIntersects:
2017 result = ( GEOSPreparedIntersects_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2018 break;
2019 case RelationTouches:
2020 result = ( GEOSPreparedTouches_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2021 break;
2022 case RelationCrosses:
2023 result = ( GEOSPreparedCrosses_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2024 break;
2025 case RelationWithin:
2026 result = ( GEOSPreparedWithin_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2027 break;
2028 case RelationContains:
2029 result = ( GEOSPreparedContains_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2030 break;
2031 case RelationDisjoint:
2032 result = ( GEOSPreparedDisjoint_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2033 break;
2034 case RelationOverlaps:
2035 result = ( GEOSPreparedOverlaps_r( context, mGeosPrepared.get(), geosGeom.get() ) == 1 );
2036 break;
2037 }
2038 return result;
2039 }
2040
2041 switch ( r )
2042 {
2043 case RelationIntersects:
2044 result = ( GEOSIntersects_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2045 break;
2046 case RelationTouches:
2047 result = ( GEOSTouches_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2048 break;
2049 case RelationCrosses:
2050 result = ( GEOSCrosses_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2051 break;
2052 case RelationWithin:
2053 result = ( GEOSWithin_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2054 break;
2055 case RelationContains:
2056 result = ( GEOSContains_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2057 break;
2058 case RelationDisjoint:
2059 result = ( GEOSDisjoint_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2060 break;
2061 case RelationOverlaps:
2062 result = ( GEOSOverlaps_r( context, mGeos.get(), geosGeom.get() ) == 1 );
2063 break;
2064 }
2065 }
2066 catch ( GEOSException &e )
2067 {
2068 logError( QStringLiteral( "GEOS" ), e.what() );
2069 if ( errorMsg )
2070 {
2071 *errorMsg = e.what();
2072 }
2073 return false;
2074 }
2075
2076 return result;
2077}
2078
2079QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, QString *errorMsg ) const
2080{
2081 if ( !mGeos )
2082 {
2083 return nullptr;
2084 }
2085
2086 geos::unique_ptr geos;
2087 try
2088 {
2089 geos.reset( GEOSBuffer_r( QgsGeosContext::get(), mGeos.get(), distance, segments ) );
2090 }
2091 CATCH_GEOS_WITH_ERRMSG( nullptr )
2092 return fromGeos( geos.get() ).release();
2093}
2094
2095QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, Qgis::EndCapStyle endCapStyle, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2096{
2097 geos::unique_ptr geos = buffer( mGeos.get(), distance, segments, endCapStyle, joinStyle, miterLimit, errorMsg );
2098 return fromGeos( geos.get() ).release();
2099}
2100
2101geos::unique_ptr QgsGeos::buffer( const GEOSGeometry *geometry, double distance, int segments, Qgis::EndCapStyle endCapStyle, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg )
2102{
2103 if ( !geometry )
2104 {
2105 return nullptr;
2106 }
2107
2108 geos::unique_ptr geos;
2109 try
2110 {
2111 geos.reset( GEOSBufferWithStyle_r( QgsGeosContext::get(), geometry, distance, segments, static_cast< int >( endCapStyle ), static_cast< int >( joinStyle ), miterLimit ) );
2112 }
2113 CATCH_GEOS_WITH_ERRMSG( nullptr )
2114 return geos;
2115}
2116
2117QgsAbstractGeometry *QgsGeos::simplify( double tolerance, QString *errorMsg ) const
2118{
2119 if ( !mGeos )
2120 {
2121 return nullptr;
2122 }
2123 geos::unique_ptr geos;
2124 try
2125 {
2126 geos.reset( GEOSTopologyPreserveSimplify_r( QgsGeosContext::get(), mGeos.get(), tolerance ) );
2127 }
2128 CATCH_GEOS_WITH_ERRMSG( nullptr )
2129 return fromGeos( geos.get() ).release();
2130}
2131
2132QgsAbstractGeometry *QgsGeos::interpolate( double distance, QString *errorMsg ) const
2133{
2134 if ( !mGeos )
2135 {
2136 return nullptr;
2137 }
2138 geos::unique_ptr geos;
2139 try
2140 {
2141 geos.reset( GEOSInterpolate_r( QgsGeosContext::get(), mGeos.get(), distance ) );
2142 }
2143 CATCH_GEOS_WITH_ERRMSG( nullptr )
2144 return fromGeos( geos.get() ).release();
2145}
2146
2147QgsPoint *QgsGeos::centroid( QString *errorMsg ) const
2148{
2149 if ( !mGeos )
2150 {
2151 return nullptr;
2152 }
2153
2154 geos::unique_ptr geos;
2155 double x;
2156 double y;
2157
2158 GEOSContextHandle_t context = QgsGeosContext::get();
2159 try
2160 {
2161 geos.reset( GEOSGetCentroid_r( context, mGeos.get() ) );
2162
2163 if ( !geos )
2164 return nullptr;
2165
2166 GEOSGeomGetX_r( context, geos.get(), &x );
2167 GEOSGeomGetY_r( context, geos.get(), &y );
2168 }
2169 CATCH_GEOS_WITH_ERRMSG( nullptr )
2170
2171 return new QgsPoint( x, y );
2172}
2173
2174QgsAbstractGeometry *QgsGeos::envelope( QString *errorMsg ) const
2175{
2176 if ( !mGeos )
2177 {
2178 return nullptr;
2179 }
2180 geos::unique_ptr geos;
2181 try
2182 {
2183 geos.reset( GEOSEnvelope_r( QgsGeosContext::get(), mGeos.get() ) );
2184 }
2185 CATCH_GEOS_WITH_ERRMSG( nullptr )
2186 return fromGeos( geos.get() ).release();
2187}
2188
2189QgsPoint *QgsGeos::pointOnSurface( QString *errorMsg ) const
2190{
2191 if ( !mGeos )
2192 {
2193 return nullptr;
2194 }
2195
2196 double x;
2197 double y;
2198
2199 GEOSContextHandle_t context = QgsGeosContext::get();
2200 geos::unique_ptr geos;
2201 try
2202 {
2203 geos.reset( GEOSPointOnSurface_r( context, mGeos.get() ) );
2204
2205 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
2206 {
2207 return nullptr;
2208 }
2209
2210 GEOSGeomGetX_r( context, geos.get(), &x );
2211 GEOSGeomGetY_r( context, geos.get(), &y );
2212 }
2213 CATCH_GEOS_WITH_ERRMSG( nullptr )
2214
2215 return new QgsPoint( x, y );
2216}
2217
2218QgsAbstractGeometry *QgsGeos::convexHull( QString *errorMsg ) const
2219{
2220 if ( !mGeos )
2221 {
2222 return nullptr;
2223 }
2224
2225 try
2226 {
2227 geos::unique_ptr cHull( GEOSConvexHull_r( QgsGeosContext::get(), mGeos.get() ) );
2228 std::unique_ptr< QgsAbstractGeometry > cHullGeom = fromGeos( cHull.get() );
2229 return cHullGeom.release();
2230 }
2231 CATCH_GEOS_WITH_ERRMSG( nullptr )
2232}
2233
2234std::unique_ptr< QgsAbstractGeometry > QgsGeos::concaveHull( double targetPercent, bool allowHoles, QString *errorMsg ) const
2235{
2236#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
2237 ( void )allowHoles;
2238 ( void )targetPercent;
2239 ( void )errorMsg;
2240 throw QgsNotSupportedException( QObject::tr( "Calculating concaveHull requires a QGIS build based on GEOS 3.11 or later" ) );
2241#else
2242 if ( !mGeos )
2243 {
2244 return nullptr;
2245 }
2246
2247 try
2248 {
2249 geos::unique_ptr concaveHull( GEOSConcaveHull_r( QgsGeosContext::get(), mGeos.get(), targetPercent, allowHoles ) );
2250 std::unique_ptr< QgsAbstractGeometry > concaveHullGeom = fromGeos( concaveHull.get() );
2251 return concaveHullGeom;
2252 }
2253 CATCH_GEOS_WITH_ERRMSG( nullptr )
2254#endif
2255}
2256
2257Qgis::CoverageValidityResult QgsGeos::validateCoverage( double gapWidth, std::unique_ptr<QgsAbstractGeometry> *invalidEdges, QString *errorMsg ) const
2258{
2259#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2260 ( void )gapWidth;
2261 ( void )invalidEdges;
2262 ( void )errorMsg;
2263 throw QgsNotSupportedException( QObject::tr( "Validating coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2264#else
2265 if ( !mGeos )
2266 {
2267 if ( errorMsg )
2268 *errorMsg = QStringLiteral( "Input geometry was not set" );
2270 }
2271
2272 GEOSContextHandle_t context = QgsGeosContext::get();
2273 try
2274 {
2275 GEOSGeometry *invalidEdgesGeos = nullptr;
2276 const int result = GEOSCoverageIsValid_r( context, mGeos.get(), gapWidth, invalidEdges ? &invalidEdgesGeos : nullptr );
2277 if ( invalidEdges && invalidEdgesGeos )
2278 {
2279 *invalidEdges = fromGeos( invalidEdgesGeos );
2280 }
2281 if ( invalidEdgesGeos )
2282 {
2283 GEOSGeom_destroy_r( context, invalidEdgesGeos );
2284 invalidEdgesGeos = nullptr;
2285 }
2286
2287 switch ( result )
2288 {
2289 case 0:
2291 case 1:
2293 case 2:
2294 break;
2295 }
2297 }
2299#endif
2300}
2301
2302std::unique_ptr<QgsAbstractGeometry> QgsGeos::simplifyCoverageVW( double tolerance, bool preserveBoundary, QString *errorMsg ) const
2303{
2304#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<12
2305 ( void )tolerance;
2306 ( void )preserveBoundary;
2307 ( void )errorMsg;
2308 throw QgsNotSupportedException( QObject::tr( "Simplifying coverages requires a QGIS build based on GEOS 3.12 or later" ) );
2309#else
2310 if ( !mGeos )
2311 {
2312 if ( errorMsg )
2313 *errorMsg = QStringLiteral( "Input geometry was not set" );
2314 return nullptr;
2315 }
2316
2317 try
2318 {
2319 geos::unique_ptr simplified( GEOSCoverageSimplifyVW_r( QgsGeosContext::get(), mGeos.get(), tolerance, preserveBoundary ? 1 : 0 ) );
2320 std::unique_ptr< QgsAbstractGeometry > simplifiedGeom = fromGeos( simplified.get() );
2321 return simplifiedGeom;
2322 }
2323 CATCH_GEOS_WITH_ERRMSG( nullptr )
2324#endif
2325}
2326
2327std::unique_ptr<QgsAbstractGeometry> QgsGeos::unionCoverage( QString *errorMsg ) const
2328{
2329 if ( !mGeos )
2330 {
2331 if ( errorMsg )
2332 *errorMsg = QStringLiteral( "Input geometry was not set" );
2333 return nullptr;
2334 }
2335
2336 try
2337 {
2338 geos::unique_ptr unioned( GEOSCoverageUnion_r( QgsGeosContext::get(), mGeos.get() ) );
2339 std::unique_ptr< QgsAbstractGeometry > result = fromGeos( unioned.get() );
2340 return result;
2341 }
2342 CATCH_GEOS_WITH_ERRMSG( nullptr )
2343}
2344
2345bool QgsGeos::isValid( QString *errorMsg, const bool allowSelfTouchingHoles, QgsGeometry *errorLoc ) const
2346{
2347 if ( !mGeos )
2348 {
2349 if ( errorMsg )
2350 *errorMsg = QObject::tr( "QGIS geometry cannot be converted to a GEOS geometry", "GEOS Error" );
2351 return false;
2352 }
2353
2354 GEOSContextHandle_t context = QgsGeosContext::get();
2355 try
2356 {
2357 GEOSGeometry *g1 = nullptr;
2358 char *r = nullptr;
2359 char res = GEOSisValidDetail_r( context, mGeos.get(), allowSelfTouchingHoles ? GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE : 0, &r, &g1 );
2360 const bool invalid = res != 1;
2361
2362 QString error;
2363 if ( r )
2364 {
2365 error = QString( r );
2366 GEOSFree_r( context, r );
2367 }
2368
2369 if ( invalid && errorMsg )
2370 {
2371 // Copied from https://github.com/libgeos/geos/blob/main/src/operation/valid/TopologyValidationError.cpp
2372 static const std::map< QString, QString > sTranslatedErrors
2373 {
2374 { QStringLiteral( "topology validation error" ), QObject::tr( "Topology validation error", "GEOS Error" ) },
2375 { QStringLiteral( "repeated point" ), QObject::tr( "Repeated point", "GEOS Error" ) },
2376 { QStringLiteral( "hole lies outside shell" ), QObject::tr( "Hole lies outside shell", "GEOS Error" ) },
2377 { QStringLiteral( "holes are nested" ), QObject::tr( "Holes are nested", "GEOS Error" ) },
2378 { QStringLiteral( "interior is disconnected" ), QObject::tr( "Interior is disconnected", "GEOS Error" ) },
2379 { QStringLiteral( "self-intersection" ), QObject::tr( "Self-intersection", "GEOS Error" ) },
2380 { QStringLiteral( "ring self-intersection" ), QObject::tr( "Ring self-intersection", "GEOS Error" ) },
2381 { QStringLiteral( "nested shells" ), QObject::tr( "Nested shells", "GEOS Error" ) },
2382 { QStringLiteral( "duplicate rings" ), QObject::tr( "Duplicate rings", "GEOS Error" ) },
2383 { QStringLiteral( "too few points in geometry component" ), QObject::tr( "Too few points in geometry component", "GEOS Error" ) },
2384 { QStringLiteral( "invalid coordinate" ), QObject::tr( "Invalid coordinate", "GEOS Error" ) },
2385 { QStringLiteral( "ring is not closed" ), QObject::tr( "Ring is not closed", "GEOS Error" ) },
2386 };
2387
2388 const auto translatedError = sTranslatedErrors.find( error.toLower() );
2389 if ( translatedError != sTranslatedErrors.end() )
2390 *errorMsg = translatedError->second;
2391 else
2392 *errorMsg = error;
2393
2394 if ( g1 && errorLoc )
2395 {
2396 *errorLoc = geometryFromGeos( g1 );
2397 }
2398 else if ( g1 )
2399 {
2400 GEOSGeom_destroy_r( context, g1 );
2401 }
2402 }
2403 return !invalid;
2404 }
2405 CATCH_GEOS_WITH_ERRMSG( false )
2406}
2407
2408bool QgsGeos::isEqual( const QgsAbstractGeometry *geom, QString *errorMsg ) const
2409{
2410 if ( !mGeos || !geom )
2411 {
2412 return false;
2413 }
2414
2415 try
2416 {
2417 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2418 if ( !geosGeom )
2419 {
2420 return false;
2421 }
2422 bool equal = GEOSEquals_r( QgsGeosContext::get(), mGeos.get(), geosGeom.get() );
2423 return equal;
2424 }
2425 CATCH_GEOS_WITH_ERRMSG( false )
2426}
2427
2428bool QgsGeos::isEmpty( QString *errorMsg ) const
2429{
2430 if ( !mGeos )
2431 {
2432 return false;
2433 }
2434
2435 try
2436 {
2437 return GEOSisEmpty_r( QgsGeosContext::get(), mGeos.get() );
2438 }
2439 CATCH_GEOS_WITH_ERRMSG( false )
2440}
2441
2442bool QgsGeos::isSimple( QString *errorMsg ) const
2443{
2444 if ( !mGeos )
2445 {
2446 return false;
2447 }
2448
2449 try
2450 {
2451 return GEOSisSimple_r( QgsGeosContext::get(), mGeos.get() );
2452 }
2453 CATCH_GEOS_WITH_ERRMSG( false )
2454}
2455
2456GEOSCoordSequence *QgsGeos::createCoordinateSequence( const QgsCurve *curve, double precision, bool forceClose )
2457{
2458 GEOSContextHandle_t context = QgsGeosContext::get();
2459
2460 std::unique_ptr< QgsLineString > segmentized;
2461 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( curve );
2462
2463 if ( !line )
2464 {
2465 segmentized.reset( curve->curveToLine() );
2466 line = segmentized.get();
2467 }
2468
2469 if ( !line )
2470 {
2471 return nullptr;
2472 }
2473 GEOSCoordSequence *coordSeq = nullptr;
2474
2475 const int numPoints = line->numPoints();
2476
2477 const bool hasZ = line->is3D();
2478
2479#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
2480 if ( qgsDoubleNear( precision, 0 ) )
2481 {
2482 if ( !forceClose || ( line->pointN( 0 ) == line->pointN( numPoints - 1 ) ) )
2483 {
2484 // use optimised method if we don't have to force close an open ring
2485 try
2486 {
2487 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, line->xData(), line->yData(), line->zData(), nullptr, numPoints );
2488 if ( !coordSeq )
2489 {
2490 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points" ).arg( numPoints ) );
2491 return nullptr;
2492 }
2493 }
2494 CATCH_GEOS( nullptr )
2495 }
2496 else
2497 {
2498 QVector< double > x = line->xVector();
2499 if ( numPoints > 0 )
2500 x.append( x.at( 0 ) );
2501 QVector< double > y = line->yVector();
2502 if ( numPoints > 0 )
2503 y.append( y.at( 0 ) );
2504 QVector< double > z = line->zVector();
2505 if ( hasZ && numPoints > 0 )
2506 z.append( z.at( 0 ) );
2507 try
2508 {
2509 coordSeq = GEOSCoordSeq_copyFromArrays_r( context, x.constData(), y.constData(), !hasZ ? nullptr : z.constData(), nullptr, numPoints + 1 );
2510 if ( !coordSeq )
2511 {
2512 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create closed coordinate sequence for %1 points" ).arg( numPoints + 1 ) );
2513 return nullptr;
2514 }
2515 }
2516 CATCH_GEOS( nullptr )
2517 }
2518 return coordSeq;
2519 }
2520#endif
2521
2522 int coordDims = 2;
2523 const bool hasM = false; //line->isMeasure(); //disabled until geos supports m-coordinates
2524
2525 if ( hasZ )
2526 {
2527 ++coordDims;
2528 }
2529 if ( hasM )
2530 {
2531 ++coordDims;
2532 }
2533
2534 int numOutPoints = numPoints;
2535 if ( forceClose && ( line->pointN( 0 ) != line->pointN( numPoints - 1 ) ) )
2536 {
2537 ++numOutPoints;
2538 }
2539
2540 try
2541 {
2542 coordSeq = GEOSCoordSeq_create_r( context, numOutPoints, coordDims );
2543 if ( !coordSeq )
2544 {
2545 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ) );
2546 return nullptr;
2547 }
2548
2549 const double *xData = line->xData();
2550 const double *yData = line->yData();
2551 const double *zData = hasZ ? line->zData() : nullptr;
2552 const double *mData = hasM ? line->mData() : nullptr;
2553
2554 if ( precision > 0. )
2555 {
2556 for ( int i = 0; i < numOutPoints; ++i )
2557 {
2558 if ( i >= numPoints )
2559 {
2560 // start reading back from start of line
2561 xData = line->xData();
2562 yData = line->yData();
2563 zData = hasZ ? line->zData() : nullptr;
2564 mData = hasM ? line->mData() : nullptr;
2565 }
2566 if ( hasZ )
2567 {
2568 GEOSCoordSeq_setXYZ_r( context, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision, std::round( *zData++ / precision ) * precision );
2569 }
2570 else
2571 {
2572 GEOSCoordSeq_setXY_r( context, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision );
2573 }
2574 if ( hasM )
2575 {
2576 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2577 }
2578 }
2579 }
2580 else
2581 {
2582 for ( int i = 0; i < numOutPoints; ++i )
2583 {
2584 if ( i >= numPoints )
2585 {
2586 // start reading back from start of line
2587 xData = line->xData();
2588 yData = line->yData();
2589 zData = hasZ ? line->zData() : nullptr;
2590 mData = hasM ? line->mData() : nullptr;
2591 }
2592 if ( hasZ )
2593 {
2594 GEOSCoordSeq_setXYZ_r( context, coordSeq, i, *xData++, *yData++, *zData++ );
2595 }
2596 else
2597 {
2598 GEOSCoordSeq_setXY_r( context, coordSeq, i, *xData++, *yData++ );
2599 }
2600 if ( hasM )
2601 {
2602 GEOSCoordSeq_setOrdinate_r( context, coordSeq, i, 3, *mData++ );
2603 }
2604 }
2605 }
2606 }
2607 CATCH_GEOS( nullptr )
2608
2609 return coordSeq;
2610}
2611
2612geos::unique_ptr QgsGeos::createGeosPoint( const QgsAbstractGeometry *point, int coordDims, double precision, Qgis::GeosCreationFlags )
2613{
2614 const QgsPoint *pt = qgsgeometry_cast<const QgsPoint *>( point );
2615 if ( !pt )
2616 return nullptr;
2617
2618 return createGeosPointXY( pt->x(), pt->y(), pt->is3D(), pt->z(), pt->isMeasure(), pt->m(), coordDims, precision );
2619}
2620
2621geos::unique_ptr QgsGeos::createGeosPointXY( double x, double y, bool hasZ, double z, bool hasM, double m, int coordDims, double precision, Qgis::GeosCreationFlags )
2622{
2623 Q_UNUSED( hasM )
2624 Q_UNUSED( m )
2625
2626 geos::unique_ptr geosPoint;
2627 GEOSContextHandle_t context = QgsGeosContext::get();
2628 try
2629 {
2630 if ( coordDims == 2 )
2631 {
2632 // optimised constructor
2633 if ( precision > 0. )
2634 geosPoint.reset( GEOSGeom_createPointFromXY_r( context, std::round( x / precision ) * precision, std::round( y / precision ) * precision ) );
2635 else
2636 geosPoint.reset( GEOSGeom_createPointFromXY_r( context, x, y ) );
2637 return geosPoint;
2638 }
2639
2640 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( context, 1, coordDims );
2641 if ( !coordSeq )
2642 {
2643 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ) );
2644 return nullptr;
2645 }
2646 if ( precision > 0. )
2647 {
2648 GEOSCoordSeq_setX_r( context, coordSeq, 0, std::round( x / precision ) * precision );
2649 GEOSCoordSeq_setY_r( context, coordSeq, 0, std::round( y / precision ) * precision );
2650 if ( hasZ )
2651 {
2652 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, std::round( z / precision ) * precision );
2653 }
2654 }
2655 else
2656 {
2657 GEOSCoordSeq_setX_r( context, coordSeq, 0, x );
2658 GEOSCoordSeq_setY_r( context, coordSeq, 0, y );
2659 if ( hasZ )
2660 {
2661 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 2, z );
2662 }
2663 }
2664#if 0 //disabled until geos supports m-coordinates
2665 if ( hasM )
2666 {
2667 GEOSCoordSeq_setOrdinate_r( context, coordSeq, 0, 3, m );
2668 }
2669#endif
2670 geosPoint.reset( GEOSGeom_createPoint_r( context, coordSeq ) );
2671 }
2672 CATCH_GEOS( nullptr )
2673 return geosPoint;
2674}
2675
2676geos::unique_ptr QgsGeos::createGeosLinestring( const QgsAbstractGeometry *curve, double precision, Qgis::GeosCreationFlags )
2677{
2678 const QgsCurve *c = qgsgeometry_cast<const QgsCurve *>( curve );
2679 if ( !c )
2680 return nullptr;
2681
2682 GEOSCoordSequence *coordSeq = createCoordinateSequence( c, precision );
2683 if ( !coordSeq )
2684 return nullptr;
2685
2686 geos::unique_ptr geosGeom;
2687 try
2688 {
2689 geosGeom.reset( GEOSGeom_createLineString_r( QgsGeosContext::get(), coordSeq ) );
2690 }
2691 CATCH_GEOS( nullptr )
2692 return geosGeom;
2693}
2694
2695geos::unique_ptr QgsGeos::createGeosPolygon( const QgsAbstractGeometry *poly, double precision, Qgis::GeosCreationFlags flags )
2696{
2697 const QgsCurvePolygon *polygon = qgsgeometry_cast<const QgsCurvePolygon *>( poly );
2698 if ( !polygon )
2699 return nullptr;
2700
2701 const QgsCurve *exteriorRing = polygon->exteriorRing();
2702 if ( !exteriorRing )
2703 {
2704 return nullptr;
2705 }
2706
2707 GEOSContextHandle_t context = QgsGeosContext::get();
2708 geos::unique_ptr geosPolygon;
2709 try
2710 {
2711 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( context, createCoordinateSequence( exteriorRing, precision, true ) ) );
2712
2713 int nHoles = 0;
2714 int nInteriorRings = polygon->numInteriorRings();
2716 {
2717 for ( int i = 0; i < nInteriorRings; ++i )
2718 {
2719 const QgsCurve *interiorRing = polygon->interiorRing( i );
2720 if ( !interiorRing->isEmpty() )
2721 {
2722 nHoles++;
2723 }
2724 }
2725 }
2726 else
2727 {
2728 nHoles = nInteriorRings;
2729 }
2730 GEOSGeometry **holes = nullptr;
2731 if ( nHoles > 0 )
2732 {
2733 holes = new GEOSGeometry*[ nHoles ];
2734 }
2735
2736 for ( int i = 0; i < nInteriorRings; ++i )
2737 {
2738 const QgsCurve *interiorRing = polygon->interiorRing( i );
2739 if ( !( flags & Qgis::GeosCreationFlag::SkipEmptyInteriorRings ) || !interiorRing->isEmpty() )
2740 {
2741 holes[i] = GEOSGeom_createLinearRing_r( context, createCoordinateSequence( interiorRing, precision, true ) );
2742 }
2743 }
2744 geosPolygon.reset( GEOSGeom_createPolygon_r( context, exteriorRingGeos.release(), holes, nHoles ) );
2745 delete[] holes;
2746 }
2747 CATCH_GEOS( nullptr )
2748
2749 return geosPolygon;
2750}
2751
2752geos::unique_ptr QgsGeos::offsetCurve( const GEOSGeometry *geometry, double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg )
2753{
2754 if ( !geometry )
2755 return nullptr;
2756
2757 geos::unique_ptr offset;
2758 try
2759 {
2760 // Force quadrant segments to be at least 8, see
2761 // https://github.com/qgis/QGIS/issues/53165#issuecomment-1563470832
2762 if ( segments < 8 )
2763 segments = 8;
2764 offset.reset( GEOSOffsetCurve_r( QgsGeosContext::get(), geometry, distance, segments, static_cast< int >( joinStyle ), miterLimit ) );
2765 }
2766 CATCH_GEOS_WITH_ERRMSG( nullptr )
2767 return offset;
2768}
2769
2770QgsAbstractGeometry *QgsGeos::offsetCurve( double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2771{
2772 geos::unique_ptr res = offsetCurve( mGeos.get(), distance, segments, joinStyle, miterLimit, errorMsg );
2773 if ( !res )
2774 return nullptr;
2775
2776 return fromGeos( res.get() ).release();
2777}
2778
2779std::unique_ptr<QgsAbstractGeometry> QgsGeos::singleSidedBuffer( double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2780{
2781 if ( !mGeos )
2782 {
2783 return nullptr;
2784 }
2785
2786 geos::unique_ptr geos;
2787 GEOSContextHandle_t context = QgsGeosContext::get();
2788 try
2789 {
2790 geos::buffer_params_unique_ptr bp( GEOSBufferParams_create_r( context ) );
2791 GEOSBufferParams_setSingleSided_r( context, bp.get(), 1 );
2792 GEOSBufferParams_setQuadrantSegments_r( context, bp.get(), segments );
2793 GEOSBufferParams_setJoinStyle_r( context, bp.get(), static_cast< int >( joinStyle ) );
2794 GEOSBufferParams_setMitreLimit_r( context, bp.get(), miterLimit ); //#spellok
2795
2796 if ( side == Qgis::BufferSide::Right )
2797 {
2798 distance = -distance;
2799 }
2800 geos.reset( GEOSBufferWithParams_r( context, mGeos.get(), bp.get(), distance ) );
2801 }
2802 CATCH_GEOS_WITH_ERRMSG( nullptr )
2803 return fromGeos( geos.get() );
2804}
2805
2806std::unique_ptr<QgsAbstractGeometry> QgsGeos::maximumInscribedCircle( double tolerance, QString *errorMsg ) const
2807{
2808 if ( !mGeos )
2809 {
2810 return nullptr;
2811 }
2812
2813 geos::unique_ptr geos;
2814 try
2815 {
2816 geos.reset( GEOSMaximumInscribedCircle_r( QgsGeosContext::get(), mGeos.get(), tolerance ) );
2817 }
2818 CATCH_GEOS_WITH_ERRMSG( nullptr )
2819 return fromGeos( geos.get() );
2820}
2821
2822std::unique_ptr<QgsAbstractGeometry> QgsGeos::largestEmptyCircle( double tolerance, const QgsAbstractGeometry *boundary, QString *errorMsg ) const
2823{
2824 if ( !mGeos )
2825 {
2826 return nullptr;
2827 }
2828
2829 geos::unique_ptr geos;
2830 try
2831 {
2832 geos::unique_ptr boundaryGeos;
2833 if ( boundary )
2834 boundaryGeos = asGeos( boundary );
2835
2836 geos.reset( GEOSLargestEmptyCircle_r( QgsGeosContext::get(), mGeos.get(), boundaryGeos.get(), tolerance ) );
2837 }
2838 CATCH_GEOS_WITH_ERRMSG( nullptr )
2839 return fromGeos( geos.get() );
2840}
2841
2842std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumWidth( QString *errorMsg ) const
2843{
2844 if ( !mGeos )
2845 {
2846 return nullptr;
2847 }
2848
2849 geos::unique_ptr geos;
2850 try
2851 {
2852 geos.reset( GEOSMinimumWidth_r( QgsGeosContext::get(), mGeos.get() ) );
2853 }
2854 CATCH_GEOS_WITH_ERRMSG( nullptr )
2855 return fromGeos( geos.get() );
2856}
2857
2858double QgsGeos::minimumClearance( QString *errorMsg ) const
2859{
2860 if ( !mGeos )
2861 {
2862 return std::numeric_limits< double >::quiet_NaN();
2863 }
2864
2865 geos::unique_ptr geos;
2866 double res = 0;
2867 try
2868 {
2869 if ( GEOSMinimumClearance_r( QgsGeosContext::get(), mGeos.get(), &res ) != 0 )
2870 return std::numeric_limits< double >::quiet_NaN();
2871 }
2872 CATCH_GEOS_WITH_ERRMSG( std::numeric_limits< double >::quiet_NaN() )
2873 return res;
2874}
2875
2876std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumClearanceLine( QString *errorMsg ) const
2877{
2878 if ( !mGeos )
2879 {
2880 return nullptr;
2881 }
2882
2883 geos::unique_ptr geos;
2884 try
2885 {
2886 geos.reset( GEOSMinimumClearanceLine_r( QgsGeosContext::get(), mGeos.get() ) );
2887 }
2888 CATCH_GEOS_WITH_ERRMSG( nullptr )
2889 return fromGeos( geos.get() );
2890}
2891
2892std::unique_ptr<QgsAbstractGeometry> QgsGeos::node( QString *errorMsg ) const
2893{
2894 if ( !mGeos )
2895 {
2896 return nullptr;
2897 }
2898
2899 geos::unique_ptr geos;
2900 try
2901 {
2902 geos.reset( GEOSNode_r( QgsGeosContext::get(), mGeos.get() ) );
2903 }
2904 CATCH_GEOS_WITH_ERRMSG( nullptr )
2905 return fromGeos( geos.get() );
2906}
2907
2908std::unique_ptr<QgsAbstractGeometry> QgsGeos::sharedPaths( const QgsAbstractGeometry *other, QString *errorMsg ) const
2909{
2910 if ( !mGeos || !other )
2911 {
2912 return nullptr;
2913 }
2914
2915 geos::unique_ptr geos;
2916 try
2917 {
2918 geos::unique_ptr otherGeos = asGeos( other );
2919 if ( !otherGeos )
2920 return nullptr;
2921
2922 geos.reset( GEOSSharedPaths_r( QgsGeosContext::get(), mGeos.get(), otherGeos.get() ) );
2923 }
2924 CATCH_GEOS_WITH_ERRMSG( nullptr )
2925 return fromGeos( geos.get() );
2926}
2927
2928std::unique_ptr<QgsAbstractGeometry> QgsGeos::reshapeGeometry( const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg ) const
2929{
2930 if ( !mGeos || mGeometry->dimension() == 0 )
2931 {
2932 if ( errorCode ) { *errorCode = InvalidBaseGeometry; }
2933 return nullptr;
2934 }
2935
2936 if ( reshapeWithLine.numPoints() < 2 )
2937 {
2938 if ( errorCode ) { *errorCode = InvalidInput; }
2939 return nullptr;
2940 }
2941
2942 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
2943
2944 GEOSContextHandle_t context = QgsGeosContext::get();
2945 //single or multi?
2946 int numGeoms = GEOSGetNumGeometries_r( context, mGeos.get() );
2947 if ( numGeoms == -1 )
2948 {
2949 if ( errorCode )
2950 {
2951 *errorCode = InvalidBaseGeometry;
2952 }
2953 return nullptr;
2954 }
2955
2956 bool isMultiGeom = false;
2957 int geosTypeId = GEOSGeomTypeId_r( context, mGeos.get() );
2958 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
2959 isMultiGeom = true;
2960
2961 bool isLine = ( mGeometry->dimension() == 1 );
2962
2963 if ( !isMultiGeom )
2964 {
2965 geos::unique_ptr reshapedGeometry;
2966 if ( isLine )
2967 {
2968 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2969 }
2970 else
2971 {
2972 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2973 }
2974
2975 if ( errorCode )
2976 *errorCode = Success;
2977 std::unique_ptr< QgsAbstractGeometry > reshapeResult = fromGeos( reshapedGeometry.get() );
2978 return reshapeResult;
2979 }
2980 else
2981 {
2982 try
2983 {
2984 //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
2985 bool reshapeTookPlace = false;
2986
2987 geos::unique_ptr currentReshapeGeometry;
2988 GEOSGeometry **newGeoms = new GEOSGeometry*[numGeoms];
2989
2990 for ( int i = 0; i < numGeoms; ++i )
2991 {
2992 if ( isLine )
2993 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2994 else
2995 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( context, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2996
2997 if ( currentReshapeGeometry )
2998 {
2999 newGeoms[i] = currentReshapeGeometry.release();
3000 reshapeTookPlace = true;
3001 }
3002 else
3003 {
3004 newGeoms[i] = GEOSGeom_clone_r( context, GEOSGetGeometryN_r( context, mGeos.get(), i ) );
3005 }
3006 }
3007
3008 geos::unique_ptr newMultiGeom;
3009 if ( isLine )
3010 {
3011 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
3012 }
3013 else //multipolygon
3014 {
3015 newMultiGeom.reset( GEOSGeom_createCollection_r( context, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
3016 }
3017
3018 delete[] newGeoms;
3019 if ( !newMultiGeom )
3020 {
3021 if ( errorCode ) { *errorCode = EngineError; }
3022 return nullptr;
3023 }
3024
3025 if ( reshapeTookPlace )
3026 {
3027 if ( errorCode )
3028 *errorCode = Success;
3029 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom = fromGeos( newMultiGeom.get() );
3030 return reshapedMultiGeom;
3031 }
3032 else
3033 {
3034 if ( errorCode )
3035 {
3036 *errorCode = NothingHappened;
3037 }
3038 return nullptr;
3039 }
3040 }
3041 CATCH_GEOS_WITH_ERRMSG( nullptr )
3042 }
3043}
3044
3045std::unique_ptr< QgsAbstractGeometry > QgsGeos::mergeLines( QString *errorMsg ) const
3046{
3047 if ( !mGeos )
3048 {
3049 return nullptr;
3050 }
3051
3052 GEOSContextHandle_t context = QgsGeosContext::get();
3053 if ( GEOSGeomTypeId_r( context, mGeos.get() ) != GEOS_MULTILINESTRING )
3054 return nullptr;
3055
3056 geos::unique_ptr geos;
3057 try
3058 {
3059 geos.reset( GEOSLineMerge_r( context, mGeos.get() ) );
3060 }
3061 CATCH_GEOS_WITH_ERRMSG( nullptr )
3062 return fromGeos( geos.get() );
3063}
3064
3065std::unique_ptr<QgsAbstractGeometry> QgsGeos::closestPoint( const QgsGeometry &other, QString *errorMsg ) const
3066{
3067 if ( !mGeos || isEmpty() || other.isEmpty() )
3068 {
3069 return nullptr;
3070 }
3071
3072 geos::unique_ptr otherGeom( asGeos( other.constGet(), mPrecision ) );
3073 if ( !otherGeom )
3074 {
3075 return nullptr;
3076 }
3077
3078 GEOSContextHandle_t context = QgsGeosContext::get();
3079 double nx = 0.0;
3080 double ny = 0.0;
3081 try
3082 {
3083 geos::coord_sequence_unique_ptr nearestCoord;
3084 if ( mGeosPrepared ) // use faster version with prepared geometry
3085 {
3086 nearestCoord.reset( GEOSPreparedNearestPoints_r( context, mGeosPrepared.get(), otherGeom.get() ) );
3087 }
3088 else
3089 {
3090 nearestCoord.reset( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3091 }
3092
3093 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx );
3094 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny );
3095 }
3096 catch ( GEOSException &e )
3097 {
3098 logError( QStringLiteral( "GEOS" ), e.what() );
3099 if ( errorMsg )
3100 {
3101 *errorMsg = e.what();
3102 }
3103 return nullptr;
3104 }
3105
3106 return std::make_unique< QgsPoint >( nx, ny );
3107}
3108
3109std::unique_ptr<QgsAbstractGeometry> QgsGeos::shortestLine( const QgsGeometry &other, QString *errorMsg ) const
3110{
3111 if ( !mGeos || other.isEmpty() )
3112 {
3113 return nullptr;
3114 }
3115
3116 return shortestLine( other.constGet(), errorMsg );
3117}
3118
3119std::unique_ptr< QgsAbstractGeometry > QgsGeos::shortestLine( const QgsAbstractGeometry *other, QString *errorMsg ) const
3120{
3121 if ( !other || other->isEmpty() )
3122 return nullptr;
3123
3124 geos::unique_ptr otherGeom( asGeos( other, mPrecision ) );
3125 if ( !otherGeom )
3126 {
3127 return nullptr;
3128 }
3129
3130 GEOSContextHandle_t context = QgsGeosContext::get();
3131 double nx1 = 0.0;
3132 double ny1 = 0.0;
3133 double nx2 = 0.0;
3134 double ny2 = 0.0;
3135 try
3136 {
3137 geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( context, mGeos.get(), otherGeom.get() ) );
3138
3139 if ( !nearestCoord )
3140 {
3141 if ( errorMsg )
3142 *errorMsg = QStringLiteral( "GEOS returned no nearest points" );
3143 return nullptr;
3144 }
3145
3146 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 0, &nx1 );
3147 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 0, &ny1 );
3148 ( void )GEOSCoordSeq_getX_r( context, nearestCoord.get(), 1, &nx2 );
3149 ( void )GEOSCoordSeq_getY_r( context, nearestCoord.get(), 1, &ny2 );
3150 }
3151 catch ( GEOSException &e )
3152 {
3153 logError( QStringLiteral( "GEOS" ), e.what() );
3154 if ( errorMsg )
3155 {
3156 *errorMsg = e.what();
3157 }
3158 return nullptr;
3159 }
3160
3161 auto line = std::make_unique< QgsLineString >();
3162 line->addVertex( QgsPoint( nx1, ny1 ) );
3163 line->addVertex( QgsPoint( nx2, ny2 ) );
3164 return line;
3165}
3166
3167double QgsGeos::lineLocatePoint( const QgsPoint &point, QString *errorMsg ) const
3168{
3169 if ( !mGeos )
3170 {
3171 return -1;
3172 }
3173
3174 geos::unique_ptr otherGeom( asGeos( &point, mPrecision ) );
3175 if ( !otherGeom )
3176 {
3177 return -1;
3178 }
3179
3180 double distance = -1;
3181 try
3182 {
3183 distance = GEOSProject_r( QgsGeosContext::get(), mGeos.get(), otherGeom.get() );
3184 }
3185 catch ( GEOSException &e )
3186 {
3187 logError( QStringLiteral( "GEOS" ), e.what() );
3188 if ( errorMsg )
3189 {
3190 *errorMsg = e.what();
3191 }
3192 return -1;
3193 }
3194
3195 return distance;
3196}
3197
3198double QgsGeos::lineLocatePoint( double x, double y, QString *errorMsg ) const
3199{
3200 if ( !mGeos )
3201 {
3202 return -1;
3203 }
3204
3205 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
3206 if ( !point )
3207 return false;
3208
3209 double distance = -1;
3210 try
3211 {
3212 distance = GEOSProject_r( QgsGeosContext::get(), mGeos.get(), point.get() );
3213 }
3214 catch ( GEOSException &e )
3215 {
3216 logError( QStringLiteral( "GEOS" ), e.what() );
3217 if ( errorMsg )
3218 {
3219 *errorMsg = e.what();
3220 }
3221 return -1;
3222 }
3223
3224 return distance;
3225}
3226
3227QgsGeometry QgsGeos::polygonize( const QVector<const QgsAbstractGeometry *> &geometries, QString *errorMsg )
3228{
3229 GEOSGeometry **const lineGeosGeometries = new GEOSGeometry*[ geometries.size()];
3230 int validLines = 0;
3231 for ( const QgsAbstractGeometry *g : geometries )
3232 {
3233 geos::unique_ptr l = asGeos( g );
3234 if ( l )
3235 {
3236 lineGeosGeometries[validLines] = l.release();
3237 validLines++;
3238 }
3239 }
3240
3241 GEOSContextHandle_t context = QgsGeosContext::get();
3242 try
3243 {
3244 geos::unique_ptr result( GEOSPolygonize_r( context, lineGeosGeometries, validLines ) );
3245 for ( int i = 0; i < validLines; ++i )
3246 {
3247 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3248 }
3249 delete[] lineGeosGeometries;
3250 return QgsGeometry( fromGeos( result.get() ) );
3251 }
3252 catch ( GEOSException &e )
3253 {
3254 if ( errorMsg )
3255 {
3256 *errorMsg = e.what();
3257 }
3258 for ( int i = 0; i < validLines; ++i )
3259 {
3260 GEOSGeom_destroy_r( context, lineGeosGeometries[i] );
3261 }
3262 delete[] lineGeosGeometries;
3263 return QgsGeometry();
3264 }
3265}
3266
3267std::unique_ptr<QgsAbstractGeometry> QgsGeos::voronoiDiagram( const QgsAbstractGeometry *extent, double tolerance, bool edgesOnly, QString *errorMsg ) const
3268{
3269 if ( !mGeos )
3270 {
3271 return nullptr;
3272 }
3273
3274 geos::unique_ptr extentGeosGeom;
3275 if ( extent )
3276 {
3277 extentGeosGeom = asGeos( extent, mPrecision );
3278 if ( !extentGeosGeom )
3279 {
3280 return nullptr;
3281 }
3282 }
3283
3284 geos::unique_ptr geos;
3285 GEOSContextHandle_t context = QgsGeosContext::get();
3286 try
3287 {
3288 geos.reset( GEOSVoronoiDiagram_r( context, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
3289
3290 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3291 {
3292 return nullptr;
3293 }
3294
3295 return fromGeos( geos.get() );
3296 }
3297 CATCH_GEOS_WITH_ERRMSG( nullptr )
3298}
3299
3300std::unique_ptr<QgsAbstractGeometry> QgsGeos::delaunayTriangulation( double tolerance, bool edgesOnly, QString *errorMsg ) const
3301{
3302 if ( !mGeos )
3303 {
3304 return nullptr;
3305 }
3306
3307 GEOSContextHandle_t context = QgsGeosContext::get();
3308 geos::unique_ptr geos;
3309 try
3310 {
3311 geos.reset( GEOSDelaunayTriangulation_r( context, mGeos.get(), tolerance, edgesOnly ) );
3312
3313 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3314 {
3315 return nullptr;
3316 }
3317
3318 return fromGeos( geos.get() );
3319 }
3320 CATCH_GEOS_WITH_ERRMSG( nullptr )
3321}
3322
3323std::unique_ptr<QgsAbstractGeometry> QgsGeos::constrainedDelaunayTriangulation( QString *errorMsg ) const
3324{
3325#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
3326 ( void )errorMsg;
3327 throw QgsNotSupportedException( QObject::tr( "Calculating constrainedDelaunayTriangulation requires a QGIS build based on GEOS 3.11 or later" ) );
3328#else
3329 if ( !mGeos )
3330 {
3331 return nullptr;
3332 }
3333
3334 geos::unique_ptr geos;
3335 GEOSContextHandle_t context = QgsGeosContext::get();
3336 try
3337 {
3338 geos.reset( GEOSConstrainedDelaunayTriangulation_r( context, mGeos.get() ) );
3339
3340 if ( !geos || GEOSisEmpty_r( context, geos.get() ) != 0 )
3341 {
3342 return nullptr;
3343 }
3344
3345 std::unique_ptr< QgsAbstractGeometry > res = fromGeos( geos.get() );
3346 if ( const QgsGeometryCollection *collection = qgsgeometry_cast< const QgsGeometryCollection * >( res.get() ) )
3347 {
3348 return std::unique_ptr< QgsAbstractGeometry >( collection->extractPartsByType( Qgis::WkbType::Polygon, true ) );
3349 }
3350 else
3351 {
3352 return res;
3353 }
3354 }
3355 CATCH_GEOS_WITH_ERRMSG( nullptr )
3356#endif
3357}
3358
3360static bool _linestringEndpoints( const GEOSGeometry *linestring, double &x1, double &y1, double &x2, double &y2 )
3361{
3362 GEOSContextHandle_t context = QgsGeosContext::get();
3363 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( context, linestring );
3364 if ( !coordSeq )
3365 return false;
3366
3367 unsigned int coordSeqSize;
3368 if ( GEOSCoordSeq_getSize_r( context, coordSeq, &coordSeqSize ) == 0 )
3369 return false;
3370
3371 if ( coordSeqSize < 2 )
3372 return false;
3373
3374 GEOSCoordSeq_getX_r( context, coordSeq, 0, &x1 );
3375 GEOSCoordSeq_getY_r( context, coordSeq, 0, &y1 );
3376 GEOSCoordSeq_getX_r( context, coordSeq, coordSeqSize - 1, &x2 );
3377 GEOSCoordSeq_getY_r( context, coordSeq, coordSeqSize - 1, &y2 );
3378 return true;
3379}
3380
3381
3383static geos::unique_ptr _mergeLinestrings( const GEOSGeometry *line1, const GEOSGeometry *line2, const QgsPointXY &intersectionPoint )
3384{
3385 double x1, y1, x2, y2;
3386 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
3387 return nullptr;
3388
3389 double rx1, ry1, rx2, ry2;
3390 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
3391 return nullptr;
3392
3393 bool intersectionAtOrigLineEndpoint =
3394 ( intersectionPoint.x() == x1 && intersectionPoint.y() == y1 ) !=
3395 ( intersectionPoint.x() == x2 && intersectionPoint.y() == y2 );
3396 bool intersectionAtReshapeLineEndpoint =
3397 ( intersectionPoint.x() == rx1 && intersectionPoint.y() == ry1 ) ||
3398 ( intersectionPoint.x() == rx2 && intersectionPoint.y() == ry2 );
3399
3400 GEOSContextHandle_t context = QgsGeosContext::get();
3401 // the intersection must be at the begin/end of both lines
3402 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
3403 {
3404 geos::unique_ptr g1( GEOSGeom_clone_r( context, line1 ) );
3405 geos::unique_ptr g2( GEOSGeom_clone_r( context, line2 ) );
3406 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
3407 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, geoms, 2 ) );
3408 geos::unique_ptr res( GEOSLineMerge_r( context, multiGeom.get() ) );
3409
3410 //keep the original orientation if the result has a start or end point in common with the original line
3411 //and this point is not the start or the end point for both lines
3412 double x1res, y1res, x2res, y2res;
3413 if ( !_linestringEndpoints( res.get(), x1res, y1res, x2res, y2res ) )
3414 return nullptr;
3415 if ( ( x1res == x2 && y1res == y2 ) || ( x2res == x1 && y2res == y1 ) )
3416 res.reset( GEOSReverse_r( context, res.get() ) );
3417
3418 return res;
3419 }
3420 else
3421 return nullptr;
3422}
3423
3424
3425geos::unique_ptr QgsGeos::reshapeLine( const GEOSGeometry *line, const GEOSGeometry *reshapeLineGeos, double precision )
3426{
3427 if ( !line || !reshapeLineGeos )
3428 return nullptr;
3429
3430 bool atLeastTwoIntersections = false;
3431 bool oneIntersection = false;
3432 QgsPointXY oneIntersectionPoint;
3433
3434 GEOSContextHandle_t context = QgsGeosContext::get();
3435 try
3436 {
3437 //make sure there are at least two intersection between line and reshape geometry
3438 geos::unique_ptr intersectGeom( GEOSIntersection_r( context, line, reshapeLineGeos ) );
3439 if ( intersectGeom )
3440 {
3441 const int geomType = GEOSGeomTypeId_r( context, intersectGeom.get() );
3442 atLeastTwoIntersections = ( geomType == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 1 )
3443 || ( geomType == GEOS_GEOMETRYCOLLECTION && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 ) // a collection implies at least two points!
3444 || ( geomType == GEOS_MULTILINESTRING && GEOSGetNumGeometries_r( context, intersectGeom.get() ) > 0 );
3445 // one point is enough when extending line at its endpoint
3446 if ( GEOSGeomTypeId_r( context, intersectGeom.get() ) == GEOS_POINT )
3447 {
3448 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( context, intersectGeom.get() );
3449 double xi, yi;
3450 GEOSCoordSeq_getX_r( context, intersectionCoordSeq, 0, &xi );
3451 GEOSCoordSeq_getY_r( context, intersectionCoordSeq, 0, &yi );
3452 oneIntersection = true;
3453 oneIntersectionPoint = QgsPointXY( xi, yi );
3454 }
3455 }
3456 }
3457 catch ( GEOSException & )
3458 {
3459 atLeastTwoIntersections = false;
3460 }
3461
3462 // special case when extending line at its endpoint
3463 if ( oneIntersection )
3464 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
3465
3466 if ( !atLeastTwoIntersections )
3467 return nullptr;
3468
3469 //begin and end point of original line
3470 double x1, y1, x2, y2;
3471 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
3472 return nullptr;
3473
3474 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1, false, 0, false, 0, 2, precision );
3475 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2, false, 0, false, 0, 2, precision );
3476
3477 bool isRing = false;
3478 if ( GEOSGeomTypeId_r( context, line ) == GEOS_LINEARRING
3479 || GEOSEquals_r( context, beginLineVertex.get(), endLineVertex.get() ) == 1 )
3480 isRing = true;
3481
3482 //node line and reshape line
3483 geos::unique_ptr nodedGeometry = nodeGeometries( reshapeLineGeos, line );
3484 if ( !nodedGeometry )
3485 {
3486 return nullptr;
3487 }
3488
3489 //and merge them together
3490 geos::unique_ptr mergedLines( GEOSLineMerge_r( context, nodedGeometry.get() ) );
3491 if ( !mergedLines )
3492 {
3493 return nullptr;
3494 }
3495
3496 int numMergedLines = GEOSGetNumGeometries_r( context, mergedLines.get() );
3497 if ( numMergedLines < 2 ) //some special cases. Normally it is >2
3498 {
3499 if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
3500 {
3501 geos::unique_ptr result( GEOSGeom_clone_r( context, reshapeLineGeos ) );
3502 return result;
3503 }
3504 else
3505 return nullptr;
3506 }
3507
3508 QVector<GEOSGeometry *> resultLineParts; //collection with the line segments that will be contained in result
3509 QVector<GEOSGeometry *> probableParts; //parts where we can decide on inclusion only after going through all the candidates
3510
3511 for ( int i = 0; i < numMergedLines; ++i )
3512 {
3513 const GEOSGeometry *currentGeom = GEOSGetGeometryN_r( context, mergedLines.get(), i );
3514
3515 // have we already added this part?
3516 bool alreadyAdded = false;
3517 double distance = 0;
3518 double bufferDistance = std::pow( 10.0L, geomDigits( currentGeom ) - 11 );
3519 for ( const GEOSGeometry *other : std::as_const( resultLineParts ) )
3520 {
3521 GEOSHausdorffDistance_r( context, currentGeom, other, &distance );
3522 if ( distance < bufferDistance )
3523 {
3524 alreadyAdded = true;
3525 break;
3526 }
3527 }
3528 if ( alreadyAdded )
3529 continue;
3530
3531 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( context, currentGeom );
3532 unsigned int currentCoordSeqSize;
3533 GEOSCoordSeq_getSize_r( context, currentCoordSeq, &currentCoordSeqSize );
3534 if ( currentCoordSeqSize < 2 )
3535 continue;
3536
3537 //get the two endpoints of the current line merge result
3538 double xBegin, xEnd, yBegin, yEnd;
3539 GEOSCoordSeq_getX_r( context, currentCoordSeq, 0, &xBegin );
3540 GEOSCoordSeq_getY_r( context, currentCoordSeq, 0, &yBegin );
3541 GEOSCoordSeq_getX_r( context, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
3542 GEOSCoordSeq_getY_r( context, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
3543 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin, false, 0, false, 0, 2, precision );
3544 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd, false, 0, false, 0, 2, precision );
3545
3546 //check how many endpoints of the line merge result are on the (original) line
3547 int nEndpointsOnOriginalLine = 0;
3548 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
3549 nEndpointsOnOriginalLine += 1;
3550
3551 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
3552 nEndpointsOnOriginalLine += 1;
3553
3554 //check how many endpoints equal the endpoints of the original line
3555 int nEndpointsSameAsOriginalLine = 0;
3556 if ( GEOSEquals_r( context, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3557 || GEOSEquals_r( context, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3558 nEndpointsSameAsOriginalLine += 1;
3559
3560 if ( GEOSEquals_r( context, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3561 || GEOSEquals_r( context, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3562 nEndpointsSameAsOriginalLine += 1;
3563
3564 //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
3565 bool currentGeomOverlapsOriginalGeom = false;
3566 bool currentGeomOverlapsReshapeLine = false;
3567 if ( lineContainedInLine( currentGeom, line ) == 1 )
3568 currentGeomOverlapsOriginalGeom = true;
3569
3570 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
3571 currentGeomOverlapsReshapeLine = true;
3572
3573 //logic to decide if this part belongs to the result
3574 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3575 {
3576 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3577 }
3578 //for closed rings, we take one segment from the candidate list
3579 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3580 {
3581 probableParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3582 }
3583 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3584 {
3585 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3586 }
3587 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3588 {
3589 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3590 }
3591 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
3592 {
3593 resultLineParts.push_back( GEOSGeom_clone_r( context, currentGeom ) );
3594 }
3595 }
3596
3597 //add the longest segment from the probable list for rings (only used for polygon rings)
3598 if ( isRing && !probableParts.isEmpty() )
3599 {
3600 geos::unique_ptr maxGeom; //the longest geometry in the probabla list
3601 GEOSGeometry *currentGeom = nullptr;
3602 double maxLength = -std::numeric_limits<double>::max();
3603 double currentLength = 0;
3604 for ( int i = 0; i < probableParts.size(); ++i )
3605 {
3606 currentGeom = probableParts.at( i );
3607 GEOSLength_r( context, currentGeom, &currentLength );
3608 if ( currentLength > maxLength )
3609 {
3610 maxLength = currentLength;
3611 maxGeom.reset( currentGeom );
3612 }
3613 else
3614 {
3615 GEOSGeom_destroy_r( context, currentGeom );
3616 }
3617 }
3618 resultLineParts.push_back( maxGeom.release() );
3619 }
3620
3621 geos::unique_ptr result;
3622 if ( resultLineParts.empty() )
3623 return nullptr;
3624
3625 if ( resultLineParts.size() == 1 ) //the whole result was reshaped
3626 {
3627 result.reset( resultLineParts[0] );
3628 }
3629 else //>1
3630 {
3631 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
3632 for ( int i = 0; i < resultLineParts.size(); ++i )
3633 {
3634 lineArray[i] = resultLineParts[i];
3635 }
3636
3637 //create multiline from resultLineParts
3638 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( context, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
3639 delete [] lineArray;
3640
3641 //then do a linemerge with the newly combined partstrings
3642 result.reset( GEOSLineMerge_r( context, multiLineGeom.get() ) );
3643 }
3644
3645 //now test if the result is a linestring. Otherwise something went wrong
3646 if ( GEOSGeomTypeId_r( context, result.get() ) != GEOS_LINESTRING )
3647 {
3648 return nullptr;
3649 }
3650
3651 //keep the original orientation
3652 bool reverseLine = false;
3653 if ( isRing )
3654 {
3655 //for closed linestring check clockwise/counter-clockwise
3656 char isResultCCW = 0, isOriginCCW = 0;
3657 if ( GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, result.get() ), &isResultCCW ) &&
3658 GEOSCoordSeq_isCCW_r( context, GEOSGeom_getCoordSeq_r( context, line ), &isOriginCCW )
3659 )
3660 {
3661 //reverse line if orientations are different
3662 reverseLine = ( isOriginCCW == 1 && isResultCCW == 0 ) || ( isOriginCCW == 0 && isResultCCW == 1 );
3663 }
3664 }
3665 else
3666 {
3667 //for linestring, check if the result has a start or end point in common with the original line
3668 double x1res, y1res, x2res, y2res;
3669 if ( !_linestringEndpoints( result.get(), x1res, y1res, x2res, y2res ) )
3670 return nullptr;
3671 geos::unique_ptr beginResultLineVertex = createGeosPointXY( x1res, y1res, false, 0, false, 0, 2, precision );
3672 geos::unique_ptr endResultLineVertex = createGeosPointXY( x2res, y2res, false, 0, false, 0, 2, precision );
3673 reverseLine = GEOSEquals_r( context, beginLineVertex.get(), endResultLineVertex.get() ) == 1 || GEOSEquals_r( context, endLineVertex.get(), beginResultLineVertex.get() ) == 1;
3674 }
3675 if ( reverseLine )
3676 result.reset( GEOSReverse_r( context, result.get() ) );
3677
3678 return result;
3679}
3680
3681geos::unique_ptr QgsGeos::reshapePolygon( const GEOSGeometry *polygon, const GEOSGeometry *reshapeLineGeos, double precision )
3682{
3683 //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
3684 int nIntersections = 0;
3685 int lastIntersectingRing = -2;
3686 const GEOSGeometry *lastIntersectingGeom = nullptr;
3687
3688 GEOSContextHandle_t context = QgsGeosContext::get();
3689 int nRings = GEOSGetNumInteriorRings_r( context, polygon );
3690 if ( nRings < 0 )
3691 return nullptr;
3692
3693 //does outer ring intersect?
3694 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( context, polygon );
3695 if ( GEOSIntersects_r( context, outerRing, reshapeLineGeos ) == 1 )
3696 {
3697 ++nIntersections;
3698 lastIntersectingRing = -1;
3699 lastIntersectingGeom = outerRing;
3700 }
3701
3702 //do inner rings intersect?
3703 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
3704
3705 try
3706 {
3707 for ( int i = 0; i < nRings; ++i )
3708 {
3709 innerRings[i] = GEOSGetInteriorRingN_r( context, polygon, i );
3710 if ( GEOSIntersects_r( context, innerRings[i], reshapeLineGeos ) == 1 )
3711 {
3712 ++nIntersections;
3713 lastIntersectingRing = i;
3714 lastIntersectingGeom = innerRings[i];
3715 }
3716 }
3717 }
3718 catch ( GEOSException & )
3719 {
3720 nIntersections = 0;
3721 }
3722
3723 if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
3724 {
3725 delete [] innerRings;
3726 return nullptr;
3727 }
3728
3729 //we have one intersecting ring, let's try to reshape it
3730 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos, precision );
3731 if ( !reshapeResult )
3732 {
3733 delete [] innerRings;
3734 return nullptr;
3735 }
3736
3737 //if reshaping took place, we need to reassemble the polygon and its rings
3738 GEOSGeometry *newRing = nullptr;
3739 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( context, reshapeResult.get() );
3740 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( context, reshapeSequence );
3741
3742 reshapeResult.reset();
3743
3744 newRing = GEOSGeom_createLinearRing_r( context, newCoordSequence );
3745 if ( !newRing )
3746 {
3747 delete [] innerRings;
3748 return nullptr;
3749 }
3750
3751 GEOSGeometry *newOuterRing = nullptr;
3752 if ( lastIntersectingRing == -1 )
3753 newOuterRing = newRing;
3754 else
3755 newOuterRing = GEOSGeom_clone_r( context, outerRing );
3756
3757 //check if all the rings are still inside the outer boundary
3758 QVector<GEOSGeometry *> ringList;
3759 if ( nRings > 0 )
3760 {
3761 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( context, GEOSGeom_clone_r( context, newOuterRing ), nullptr, 0 );
3762 if ( outerRingPoly )
3763 {
3764 ringList.reserve( nRings );
3765 GEOSGeometry *currentRing = nullptr;
3766 for ( int i = 0; i < nRings; ++i )
3767 {
3768 if ( lastIntersectingRing == i )
3769 currentRing = newRing;
3770 else
3771 currentRing = GEOSGeom_clone_r( context, innerRings[i] );
3772
3773 //possibly a ring is no longer contained in the result polygon after reshape
3774 if ( GEOSContains_r( context, outerRingPoly, currentRing ) == 1 )
3775 ringList.push_back( currentRing );
3776 else
3777 GEOSGeom_destroy_r( context, currentRing );
3778 }
3779 }
3780 GEOSGeom_destroy_r( context, outerRingPoly );
3781 }
3782
3783 GEOSGeometry **newInnerRings = new GEOSGeometry*[ringList.size()];
3784 for ( int i = 0; i < ringList.size(); ++i )
3785 newInnerRings[i] = ringList.at( i );
3786
3787 delete [] innerRings;
3788
3789 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( context, newOuterRing, newInnerRings, ringList.size() ) );
3790 delete[] newInnerRings;
3791
3792 return reshapedPolygon;
3793}
3794
3795int QgsGeos::lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 )
3796{
3797 if ( !line1 || !line2 )
3798 {
3799 return -1;
3800 }
3801
3802 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
3803
3804 GEOSContextHandle_t context = QgsGeosContext::get();
3805 geos::unique_ptr bufferGeom( GEOSBuffer_r( context, line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ) );
3806 if ( !bufferGeom )
3807 return -2;
3808
3809 geos::unique_ptr intersectionGeom( GEOSIntersection_r( context, bufferGeom.get(), line1 ) );
3810
3811 //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
3812 double intersectGeomLength;
3813 double line1Length;
3814
3815 GEOSLength_r( context, intersectionGeom.get(), &intersectGeomLength );
3816 GEOSLength_r( context, line1, &line1Length );
3817
3818 double intersectRatio = line1Length / intersectGeomLength;
3819 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
3820 return 1;
3821
3822 return 0;
3823}
3824
3825int QgsGeos::pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line )
3826{
3827 if ( !point || !line )
3828 return -1;
3829
3830 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
3831
3832 GEOSContextHandle_t context = QgsGeosContext::get();
3833 geos::unique_ptr lineBuffer( GEOSBuffer_r( context, line, bufferDistance, 8 ) );
3834 if ( !lineBuffer )
3835 return -2;
3836
3837 bool contained = false;
3838 if ( GEOSContains_r( context, lineBuffer.get(), point ) == 1 )
3839 contained = true;
3840
3841 return contained;
3842}
3843
3844int QgsGeos::geomDigits( const GEOSGeometry *geom )
3845{
3846 GEOSContextHandle_t context = QgsGeosContext::get();
3847 geos::unique_ptr bbox( GEOSEnvelope_r( context, geom ) );
3848 if ( !bbox.get() )
3849 return -1;
3850
3851 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( context, bbox.get() );
3852 if ( !bBoxRing )
3853 return -1;
3854
3855 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( context, bBoxRing );
3856
3857 if ( !bBoxCoordSeq )
3858 return -1;
3859
3860 unsigned int nCoords = 0;
3861 if ( !GEOSCoordSeq_getSize_r( context, bBoxCoordSeq, &nCoords ) )
3862 return -1;
3863
3864 int maxDigits = -1;
3865 for ( unsigned int i = 0; i < nCoords - 1; ++i )
3866 {
3867 double t;
3868 GEOSCoordSeq_getX_r( context, bBoxCoordSeq, i, &t );
3869
3870 int digits;
3871 digits = std::ceil( std::log10( std::fabs( t ) ) );
3872 if ( digits > maxDigits )
3873 maxDigits = digits;
3874
3875 GEOSCoordSeq_getY_r( context, bBoxCoordSeq, i, &t );
3876 digits = std::ceil( std::log10( std::fabs( t ) ) );
3877 if ( digits > maxDigits )
3878 maxDigits = digits;
3879 }
3880
3881 return maxDigits;
3882}
The Qgis class provides global constants for use throughout the application.
Definition qgis.h:54
BufferSide
Side of line to buffer.
Definition qgis.h:2026
@ Right
Buffer to right of line.
GeometryOperationResult
Success or failure of a geometry operation.
Definition qgis.h:1972
@ 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:2075
@ Polygon
Polygons.
@ Unknown
Unknown types.
@ Null
No geometry.
JoinStyle
Join styles for buffers.
Definition qgis.h:2051
EndCapStyle
End cap styles for buffers.
Definition qgis.h:2038
CoverageValidityResult
Coverage validity results.
Definition qgis.h:2084
@ 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:2097
@ 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:2752
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:2345
std::unique_ptr< QgsAbstractGeometry > subdivide(int maxNodes, QString *errorMsg=nullptr, const QgsGeometryParameters &parameters=QgsGeometryParameters()) const
Subdivides the geometry.
Definition qgsgeos.cpp:448
std::unique_ptr< QgsAbstractGeometry > mergeLines(QString *errorMsg=nullptr) const
Merges any connected lines in a LineString/MultiLineString geometry and converts them to single line ...
Definition qgsgeos.cpp:3045
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:3323
QgsAbstractGeometry * buffer(double distance, int segments, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2079
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:3065
std::unique_ptr< QgsAbstractGeometry > reshapeGeometry(const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg=nullptr) const
Reshapes the geometry using a line.
Definition qgsgeos.cpp:2928
std::unique_ptr< QgsAbstractGeometry > maximumInscribedCircle(double tolerance, QString *errorMsg=nullptr) const
Returns the maximum inscribed circle.
Definition qgsgeos.cpp:2806
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:2908
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:2892
double minimumClearance(QString *errorMsg=nullptr) const
Computes the minimum clearance of a geometry.
Definition qgsgeos.cpp:2858
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:2779
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 > 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:3267
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2117
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:2327
static std::unique_ptr< QgsPolygon > fromGeosPolygon(const GEOSGeometry *geos)
Definition qgsgeos.cpp:1669
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:2876
QgsAbstractGeometry * envelope(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2174
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:2218
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:2842
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
Definition qgsgeos.cpp:2442
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:3109
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:2234
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:2132
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:2408
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:3227
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:2147
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:3167
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:2302
bool isEmpty(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2428
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:3300
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:2822
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:2257
static QgsPoint coordSeqPoint(const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM)
Definition qgsgeos.cpp:1762
QgsPoint * pointOnSurface(QString *errorMsg=nullptr) const override
Calculate a point that is guaranteed to be on the surface of this.
Definition qgsgeos.cpp:2189
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.
A class to represent 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
double z
Definition qgspoint.h:54
double x
Definition qgspoint.h:52
double distance(double x, double y) const
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition qgspoint.h:393
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 bool isMultiType(Qgis::WkbType type)
Returns true if the WKB type is a multi type.
static Qgis::WkbType flatType(Qgis::WkbType type)
Returns the flat type for a WKB 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:6125
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