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