QGIS API Documentation 3.34.0-Prizren (ffbdd678812)
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
2046bool QgsGeos::isValid( QString *errorMsg, const bool allowSelfTouchingHoles, QgsGeometry *errorLoc ) const
2047{
2048 if ( !mGeos )
2049 {
2050 return false;
2051 }
2052
2053 try
2054 {
2055 GEOSGeometry *g1 = nullptr;
2056 char *r = nullptr;
2057 char res = GEOSisValidDetail_r( geosinit()->ctxt, mGeos.get(), allowSelfTouchingHoles ? GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE : 0, &r, &g1 );
2058 const bool invalid = res != 1;
2059
2060 QString error;
2061 if ( r )
2062 {
2063 error = QString( r );
2064 GEOSFree_r( geosinit()->ctxt, r );
2065 }
2066
2067 if ( invalid && errorMsg )
2068 {
2069 static QgsStringMap translatedErrors;
2070
2071 if ( translatedErrors.empty() )
2072 {
2073 // Copied from https://git.osgeo.org/gitea/geos/geos/src/branch/master/src/operation/valid/TopologyValidationError.cpp
2074 translatedErrors.insert( QStringLiteral( "topology validation error" ), QObject::tr( "Topology validation error", "GEOS Error" ) );
2075 translatedErrors.insert( QStringLiteral( "repeated point" ), QObject::tr( "Repeated point", "GEOS Error" ) );
2076 translatedErrors.insert( QStringLiteral( "hole lies outside shell" ), QObject::tr( "Hole lies outside shell", "GEOS Error" ) );
2077 translatedErrors.insert( QStringLiteral( "holes are nested" ), QObject::tr( "Holes are nested", "GEOS Error" ) );
2078 translatedErrors.insert( QStringLiteral( "interior is disconnected" ), QObject::tr( "Interior is disconnected", "GEOS Error" ) );
2079 translatedErrors.insert( QStringLiteral( "self-intersection" ), QObject::tr( "Self-intersection", "GEOS Error" ) );
2080 translatedErrors.insert( QStringLiteral( "ring self-intersection" ), QObject::tr( "Ring self-intersection", "GEOS Error" ) );
2081 translatedErrors.insert( QStringLiteral( "nested shells" ), QObject::tr( "Nested shells", "GEOS Error" ) );
2082 translatedErrors.insert( QStringLiteral( "duplicate rings" ), QObject::tr( "Duplicate rings", "GEOS Error" ) );
2083 translatedErrors.insert( QStringLiteral( "too few points in geometry component" ), QObject::tr( "Too few points in geometry component", "GEOS Error" ) );
2084 translatedErrors.insert( QStringLiteral( "invalid coordinate" ), QObject::tr( "Invalid coordinate", "GEOS Error" ) );
2085 translatedErrors.insert( QStringLiteral( "ring is not closed" ), QObject::tr( "Ring is not closed", "GEOS Error" ) );
2086 }
2087
2088 *errorMsg = translatedErrors.value( error.toLower(), error );
2089
2090 if ( g1 && errorLoc )
2091 {
2092 *errorLoc = geometryFromGeos( g1 );
2093 }
2094 else if ( g1 )
2095 {
2096 GEOSGeom_destroy_r( geosinit()->ctxt, g1 );
2097 }
2098 }
2099 return !invalid;
2100 }
2101 CATCH_GEOS_WITH_ERRMSG( false )
2102}
2103
2104bool QgsGeos::isEqual( const QgsAbstractGeometry *geom, QString *errorMsg ) const
2105{
2106 if ( !mGeos || !geom )
2107 {
2108 return false;
2109 }
2110
2111 try
2112 {
2113 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2114 if ( !geosGeom )
2115 {
2116 return false;
2117 }
2118 bool equal = GEOSEquals_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() );
2119 return equal;
2120 }
2121 CATCH_GEOS_WITH_ERRMSG( false )
2122}
2123
2124bool QgsGeos::isEmpty( QString *errorMsg ) const
2125{
2126 if ( !mGeos )
2127 {
2128 return false;
2129 }
2130
2131 try
2132 {
2133 return GEOSisEmpty_r( geosinit()->ctxt, mGeos.get() );
2134 }
2135 CATCH_GEOS_WITH_ERRMSG( false )
2136}
2137
2138bool QgsGeos::isSimple( QString *errorMsg ) const
2139{
2140 if ( !mGeos )
2141 {
2142 return false;
2143 }
2144
2145 try
2146 {
2147 return GEOSisSimple_r( geosinit()->ctxt, mGeos.get() );
2148 }
2149 CATCH_GEOS_WITH_ERRMSG( false )
2150}
2151
2152GEOSCoordSequence *QgsGeos::createCoordinateSequence( const QgsCurve *curve, double precision, bool forceClose )
2153{
2154 GEOSContextHandle_t ctxt = geosinit()->ctxt;
2155
2156 std::unique_ptr< QgsLineString > segmentized;
2157 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( curve );
2158
2159 if ( !line )
2160 {
2161 segmentized.reset( curve->curveToLine() );
2162 line = segmentized.get();
2163 }
2164
2165 if ( !line )
2166 {
2167 return nullptr;
2168 }
2169 GEOSCoordSequence *coordSeq = nullptr;
2170
2171 const int numPoints = line->numPoints();
2172
2173 const bool hasZ = line->is3D();
2174
2175#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
2176 if ( qgsDoubleNear( precision, 0 ) )
2177 {
2178 if ( !forceClose || ( line->pointN( 0 ) == line->pointN( numPoints - 1 ) ) )
2179 {
2180 // use optimised method if we don't have to force close an open ring
2181 try
2182 {
2183 coordSeq = GEOSCoordSeq_copyFromArrays_r( ctxt, line->xData(), line->yData(), line->zData(), nullptr, numPoints );
2184 if ( !coordSeq )
2185 {
2186 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points" ).arg( numPoints ) );
2187 return nullptr;
2188 }
2189 }
2190 CATCH_GEOS( nullptr )
2191 }
2192 else
2193 {
2194 QVector< double > x = line->xVector();
2195 if ( numPoints > 0 )
2196 x.append( x.at( 0 ) );
2197 QVector< double > y = line->yVector();
2198 if ( numPoints > 0 )
2199 y.append( y.at( 0 ) );
2200 QVector< double > z = line->zVector();
2201 if ( hasZ && numPoints > 0 )
2202 z.append( z.at( 0 ) );
2203 try
2204 {
2205 coordSeq = GEOSCoordSeq_copyFromArrays_r( ctxt, x.constData(), y.constData(), !hasZ ? nullptr : z.constData(), nullptr, numPoints + 1 );
2206 if ( !coordSeq )
2207 {
2208 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create closed coordinate sequence for %1 points" ).arg( numPoints + 1 ) );
2209 return nullptr;
2210 }
2211 }
2212 CATCH_GEOS( nullptr )
2213 }
2214 return coordSeq;
2215 }
2216#endif
2217
2218 int coordDims = 2;
2219 const bool hasM = false; //line->isMeasure(); //disabled until geos supports m-coordinates
2220
2221 if ( hasZ )
2222 {
2223 ++coordDims;
2224 }
2225 if ( hasM )
2226 {
2227 ++coordDims;
2228 }
2229
2230 int numOutPoints = numPoints;
2231 if ( forceClose && ( line->pointN( 0 ) != line->pointN( numPoints - 1 ) ) )
2232 {
2233 ++numOutPoints;
2234 }
2235
2236 try
2237 {
2238 coordSeq = GEOSCoordSeq_create_r( ctxt, numOutPoints, coordDims );
2239 if ( !coordSeq )
2240 {
2241 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ) );
2242 return nullptr;
2243 }
2244
2245 const double *xData = line->xData();
2246 const double *yData = line->yData();
2247 const double *zData = hasZ ? line->zData() : nullptr;
2248 const double *mData = hasM ? line->mData() : nullptr;
2249
2250 if ( precision > 0. )
2251 {
2252 for ( int i = 0; i < numOutPoints; ++i )
2253 {
2254 if ( i >= numPoints )
2255 {
2256 // start reading back from start of line
2257 xData = line->xData();
2258 yData = line->yData();
2259 zData = hasZ ? line->zData() : nullptr;
2260 mData = hasM ? line->mData() : nullptr;
2261 }
2262 if ( hasZ )
2263 {
2264 GEOSCoordSeq_setXYZ_r( ctxt, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision, std::round( *zData++ / precision ) * precision );
2265 }
2266 else
2267 {
2268 GEOSCoordSeq_setXY_r( ctxt, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision );
2269 }
2270 if ( hasM )
2271 {
2272 GEOSCoordSeq_setOrdinate_r( ctxt, coordSeq, i, 3, *mData++ );
2273 }
2274 }
2275 }
2276 else
2277 {
2278 for ( int i = 0; i < numOutPoints; ++i )
2279 {
2280 if ( i >= numPoints )
2281 {
2282 // start reading back from start of line
2283 xData = line->xData();
2284 yData = line->yData();
2285 zData = hasZ ? line->zData() : nullptr;
2286 mData = hasM ? line->mData() : nullptr;
2287 }
2288 if ( hasZ )
2289 {
2290 GEOSCoordSeq_setXYZ_r( ctxt, coordSeq, i, *xData++, *yData++, *zData++ );
2291 }
2292 else
2293 {
2294 GEOSCoordSeq_setXY_r( ctxt, coordSeq, i, *xData++, *yData++ );
2295 }
2296 if ( hasM )
2297 {
2298 GEOSCoordSeq_setOrdinate_r( ctxt, coordSeq, i, 3, *mData++ );
2299 }
2300 }
2301 }
2302 }
2303 CATCH_GEOS( nullptr )
2304
2305 return coordSeq;
2306}
2307
2308geos::unique_ptr QgsGeos::createGeosPoint( const QgsAbstractGeometry *point, int coordDims, double precision )
2309{
2310 const QgsPoint *pt = qgsgeometry_cast<const QgsPoint *>( point );
2311 if ( !pt )
2312 return nullptr;
2313
2314 return createGeosPointXY( pt->x(), pt->y(), pt->is3D(), pt->z(), pt->isMeasure(), pt->m(), coordDims, precision );
2315}
2316
2317geos::unique_ptr QgsGeos::createGeosPointXY( double x, double y, bool hasZ, double z, bool hasM, double m, int coordDims, double precision )
2318{
2319 Q_UNUSED( hasM )
2320 Q_UNUSED( m )
2321
2322 geos::unique_ptr geosPoint;
2323 try
2324 {
2325 if ( coordDims == 2 )
2326 {
2327 // optimised constructor
2328 if ( precision > 0. )
2329 geosPoint.reset( GEOSGeom_createPointFromXY_r( geosinit()->ctxt, std::round( x / precision ) * precision, std::round( y / precision ) * precision ) );
2330 else
2331 geosPoint.reset( GEOSGeom_createPointFromXY_r( geosinit()->ctxt, x, y ) );
2332 return geosPoint;
2333 }
2334
2335 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( geosinit()->ctxt, 1, coordDims );
2336 if ( !coordSeq )
2337 {
2338 QgsDebugError( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ) );
2339 return nullptr;
2340 }
2341 if ( precision > 0. )
2342 {
2343 GEOSCoordSeq_setX_r( geosinit()->ctxt, coordSeq, 0, std::round( x / precision ) * precision );
2344 GEOSCoordSeq_setY_r( geosinit()->ctxt, coordSeq, 0, std::round( y / precision ) * precision );
2345 if ( hasZ )
2346 {
2347 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 2, std::round( z / precision ) * precision );
2348 }
2349 }
2350 else
2351 {
2352 GEOSCoordSeq_setX_r( geosinit()->ctxt, coordSeq, 0, x );
2353 GEOSCoordSeq_setY_r( geosinit()->ctxt, coordSeq, 0, y );
2354 if ( hasZ )
2355 {
2356 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 2, z );
2357 }
2358 }
2359#if 0 //disabled until geos supports m-coordinates
2360 if ( hasM )
2361 {
2362 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 3, m );
2363 }
2364#endif
2365 geosPoint.reset( GEOSGeom_createPoint_r( geosinit()->ctxt, coordSeq ) );
2366 }
2367 CATCH_GEOS( nullptr )
2368 return geosPoint;
2369}
2370
2371geos::unique_ptr QgsGeos::createGeosLinestring( const QgsAbstractGeometry *curve, double precision )
2372{
2373 const QgsCurve *c = qgsgeometry_cast<const QgsCurve *>( curve );
2374 if ( !c )
2375 return nullptr;
2376
2377 GEOSCoordSequence *coordSeq = createCoordinateSequence( c, precision );
2378 if ( !coordSeq )
2379 return nullptr;
2380
2381 geos::unique_ptr geosGeom;
2382 try
2383 {
2384 geosGeom.reset( GEOSGeom_createLineString_r( geosinit()->ctxt, coordSeq ) );
2385 }
2386 CATCH_GEOS( nullptr )
2387 return geosGeom;
2388}
2389
2390geos::unique_ptr QgsGeos::createGeosPolygon( const QgsAbstractGeometry *poly, double precision )
2391{
2392 const QgsCurvePolygon *polygon = qgsgeometry_cast<const QgsCurvePolygon *>( poly );
2393 if ( !polygon )
2394 return nullptr;
2395
2396 const QgsCurve *exteriorRing = polygon->exteriorRing();
2397 if ( !exteriorRing )
2398 {
2399 return nullptr;
2400 }
2401
2402 geos::unique_ptr geosPolygon;
2403 try
2404 {
2405 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( geosinit()->ctxt, createCoordinateSequence( exteriorRing, precision, true ) ) );
2406
2407 int nHoles = polygon->numInteriorRings();
2408 GEOSGeometry **holes = nullptr;
2409 if ( nHoles > 0 )
2410 {
2411 holes = new GEOSGeometry*[ nHoles ];
2412 }
2413
2414 for ( int i = 0; i < nHoles; ++i )
2415 {
2416 const QgsCurve *interiorRing = polygon->interiorRing( i );
2417 holes[i] = GEOSGeom_createLinearRing_r( geosinit()->ctxt, createCoordinateSequence( interiorRing, precision, true ) );
2418 }
2419 geosPolygon.reset( GEOSGeom_createPolygon_r( geosinit()->ctxt, exteriorRingGeos.release(), holes, nHoles ) );
2420 delete[] holes;
2421 }
2422 CATCH_GEOS( nullptr )
2423
2424 return geosPolygon;
2425}
2426
2427QgsAbstractGeometry *QgsGeos::offsetCurve( double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2428{
2429 if ( !mGeos )
2430 return nullptr;
2431
2432 geos::unique_ptr offset;
2433 try
2434 {
2435 // Force quadrant segments to be at least 8, see
2436 // https://github.com/qgis/QGIS/issues/53165#issuecomment-1563470832
2437 if ( segments < 8 )
2438 segments = 8;
2439 offset.reset( GEOSOffsetCurve_r( geosinit()->ctxt, mGeos.get(), distance, segments, static_cast< int >( joinStyle ), miterLimit ) );
2440 }
2441 CATCH_GEOS_WITH_ERRMSG( nullptr )
2442 std::unique_ptr< QgsAbstractGeometry > offsetGeom = fromGeos( offset.get() );
2443 return offsetGeom.release();
2444}
2445
2446std::unique_ptr<QgsAbstractGeometry> QgsGeos::singleSidedBuffer( double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2447{
2448 if ( !mGeos )
2449 {
2450 return nullptr;
2451 }
2452
2454 try
2455 {
2456 geos::buffer_params_unique_ptr bp( GEOSBufferParams_create_r( geosinit()->ctxt ) );
2457 GEOSBufferParams_setSingleSided_r( geosinit()->ctxt, bp.get(), 1 );
2458 GEOSBufferParams_setQuadrantSegments_r( geosinit()->ctxt, bp.get(), segments );
2459 GEOSBufferParams_setJoinStyle_r( geosinit()->ctxt, bp.get(), static_cast< int >( joinStyle ) );
2460 GEOSBufferParams_setMitreLimit_r( geosinit()->ctxt, bp.get(), miterLimit ); //#spellok
2461
2462 if ( side == Qgis::BufferSide::Right )
2463 {
2464 distance = -distance;
2465 }
2466 geos.reset( GEOSBufferWithParams_r( geosinit()->ctxt, mGeos.get(), bp.get(), distance ) );
2467 }
2468 CATCH_GEOS_WITH_ERRMSG( nullptr )
2469 return fromGeos( geos.get() );
2470}
2471
2472std::unique_ptr<QgsAbstractGeometry> QgsGeos::maximumInscribedCircle( double tolerance, QString *errorMsg ) const
2473{
2474 if ( !mGeos )
2475 {
2476 return nullptr;
2477 }
2478
2480 try
2481 {
2482 geos.reset( GEOSMaximumInscribedCircle_r( geosinit()->ctxt, mGeos.get(), tolerance ) );
2483 }
2484 CATCH_GEOS_WITH_ERRMSG( nullptr )
2485 return fromGeos( geos.get() );
2486}
2487
2488std::unique_ptr<QgsAbstractGeometry> QgsGeos::largestEmptyCircle( double tolerance, const QgsAbstractGeometry *boundary, QString *errorMsg ) const
2489{
2490 if ( !mGeos )
2491 {
2492 return nullptr;
2493 }
2494
2496 try
2497 {
2498 geos::unique_ptr boundaryGeos;
2499 if ( boundary )
2500 boundaryGeos = asGeos( boundary );
2501
2502 geos.reset( GEOSLargestEmptyCircle_r( geosinit()->ctxt, mGeos.get(), boundaryGeos.get(), tolerance ) );
2503 }
2504 CATCH_GEOS_WITH_ERRMSG( nullptr )
2505 return fromGeos( geos.get() );
2506}
2507
2508std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumWidth( QString *errorMsg ) const
2509{
2510 if ( !mGeos )
2511 {
2512 return nullptr;
2513 }
2514
2516 try
2517 {
2518 geos.reset( GEOSMinimumWidth_r( geosinit()->ctxt, mGeos.get() ) );
2519 }
2520 CATCH_GEOS_WITH_ERRMSG( nullptr )
2521 return fromGeos( geos.get() );
2522}
2523
2524double QgsGeos::minimumClearance( QString *errorMsg ) const
2525{
2526 if ( !mGeos )
2527 {
2528 return std::numeric_limits< double >::quiet_NaN();
2529 }
2530
2532 double res = 0;
2533 try
2534 {
2535 if ( GEOSMinimumClearance_r( geosinit()->ctxt, mGeos.get(), &res ) != 0 )
2536 return std::numeric_limits< double >::quiet_NaN();
2537 }
2538 CATCH_GEOS_WITH_ERRMSG( std::numeric_limits< double >::quiet_NaN() )
2539 return res;
2540}
2541
2542std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumClearanceLine( QString *errorMsg ) const
2543{
2544 if ( !mGeos )
2545 {
2546 return nullptr;
2547 }
2548
2550 try
2551 {
2552 geos.reset( GEOSMinimumClearanceLine_r( geosinit()->ctxt, mGeos.get() ) );
2553 }
2554 CATCH_GEOS_WITH_ERRMSG( nullptr )
2555 return fromGeos( geos.get() );
2556}
2557
2558std::unique_ptr<QgsAbstractGeometry> QgsGeos::node( QString *errorMsg ) const
2559{
2560 if ( !mGeos )
2561 {
2562 return nullptr;
2563 }
2564
2566 try
2567 {
2568 geos.reset( GEOSNode_r( geosinit()->ctxt, mGeos.get() ) );
2569 }
2570 CATCH_GEOS_WITH_ERRMSG( nullptr )
2571 return fromGeos( geos.get() );
2572}
2573
2574std::unique_ptr<QgsAbstractGeometry> QgsGeos::sharedPaths( const QgsAbstractGeometry *other, QString *errorMsg ) const
2575{
2576 if ( !mGeos || !other )
2577 {
2578 return nullptr;
2579 }
2580
2582 try
2583 {
2584 geos::unique_ptr otherGeos = asGeos( other );
2585 if ( !otherGeos )
2586 return nullptr;
2587
2588 geos.reset( GEOSSharedPaths_r( geosinit()->ctxt, mGeos.get(), otherGeos.get() ) );
2589 }
2590 CATCH_GEOS_WITH_ERRMSG( nullptr )
2591 return fromGeos( geos.get() );
2592}
2593
2594std::unique_ptr<QgsAbstractGeometry> QgsGeos::reshapeGeometry( const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg ) const
2595{
2596 if ( !mGeos || mGeometry->dimension() == 0 )
2597 {
2598 if ( errorCode ) { *errorCode = InvalidBaseGeometry; }
2599 return nullptr;
2600 }
2601
2602 if ( reshapeWithLine.numPoints() < 2 )
2603 {
2604 if ( errorCode ) { *errorCode = InvalidInput; }
2605 return nullptr;
2606 }
2607
2608 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
2609
2610 //single or multi?
2611 int numGeoms = GEOSGetNumGeometries_r( geosinit()->ctxt, mGeos.get() );
2612 if ( numGeoms == -1 )
2613 {
2614 if ( errorCode )
2615 {
2616 *errorCode = InvalidBaseGeometry;
2617 }
2618 return nullptr;
2619 }
2620
2621 bool isMultiGeom = false;
2622 int geosTypeId = GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() );
2623 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
2624 isMultiGeom = true;
2625
2626 bool isLine = ( mGeometry->dimension() == 1 );
2627
2628 if ( !isMultiGeom )
2629 {
2630 geos::unique_ptr reshapedGeometry;
2631 if ( isLine )
2632 {
2633 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2634 }
2635 else
2636 {
2637 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2638 }
2639
2640 if ( errorCode )
2641 *errorCode = Success;
2642 std::unique_ptr< QgsAbstractGeometry > reshapeResult = fromGeos( reshapedGeometry.get() );
2643 return reshapeResult;
2644 }
2645 else
2646 {
2647 try
2648 {
2649 //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
2650 bool reshapeTookPlace = false;
2651
2652 geos::unique_ptr currentReshapeGeometry;
2653 GEOSGeometry **newGeoms = new GEOSGeometry*[numGeoms];
2654
2655 for ( int i = 0; i < numGeoms; ++i )
2656 {
2657 if ( isLine )
2658 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2659 else
2660 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2661
2662 if ( currentReshapeGeometry )
2663 {
2664 newGeoms[i] = currentReshapeGeometry.release();
2665 reshapeTookPlace = true;
2666 }
2667 else
2668 {
2669 newGeoms[i] = GEOSGeom_clone_r( geosinit()->ctxt, GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ) );
2670 }
2671 }
2672
2673 geos::unique_ptr newMultiGeom;
2674 if ( isLine )
2675 {
2676 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
2677 }
2678 else //multipolygon
2679 {
2680 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
2681 }
2682
2683 delete[] newGeoms;
2684 if ( !newMultiGeom )
2685 {
2686 if ( errorCode ) { *errorCode = EngineError; }
2687 return nullptr;
2688 }
2689
2690 if ( reshapeTookPlace )
2691 {
2692 if ( errorCode )
2693 *errorCode = Success;
2694 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom = fromGeos( newMultiGeom.get() );
2695 return reshapedMultiGeom;
2696 }
2697 else
2698 {
2699 if ( errorCode )
2700 {
2701 *errorCode = NothingHappened;
2702 }
2703 return nullptr;
2704 }
2705 }
2706 CATCH_GEOS_WITH_ERRMSG( nullptr )
2707 }
2708}
2709
2710QgsGeometry QgsGeos::mergeLines( QString *errorMsg ) const
2711{
2712 if ( !mGeos )
2713 {
2714 return QgsGeometry();
2715 }
2716
2717 if ( GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() ) != GEOS_MULTILINESTRING )
2718 return QgsGeometry();
2719
2721 try
2722 {
2723 geos.reset( GEOSLineMerge_r( geosinit()->ctxt, mGeos.get() ) );
2724 }
2726 return QgsGeometry( fromGeos( geos.get() ) );
2727}
2728
2729QgsGeometry QgsGeos::closestPoint( const QgsGeometry &other, QString *errorMsg ) const
2730{
2731 if ( !mGeos || isEmpty() || other.isEmpty() )
2732 {
2733 return QgsGeometry();
2734 }
2735
2736 geos::unique_ptr otherGeom( asGeos( other.constGet(), mPrecision ) );
2737 if ( !otherGeom )
2738 {
2739 return QgsGeometry();
2740 }
2741
2742 double nx = 0.0;
2743 double ny = 0.0;
2744 try
2745 {
2747 if ( mGeosPrepared ) // use faster version with prepared geometry
2748 {
2749 nearestCoord.reset( GEOSPreparedNearestPoints_r( geosinit()->ctxt, mGeosPrepared.get(), otherGeom.get() ) );
2750 }
2751 else
2752 {
2753 nearestCoord.reset( GEOSNearestPoints_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() ) );
2754 }
2755
2756 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 0, &nx );
2757 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 0, &ny );
2758 }
2759 catch ( GEOSException &e )
2760 {
2761 logError( QStringLiteral( "GEOS" ), e.what() );
2762 if ( errorMsg )
2763 {
2764 *errorMsg = e.what();
2765 }
2766 return QgsGeometry();
2767 }
2768
2769 return QgsGeometry( new QgsPoint( nx, ny ) );
2770}
2771
2772QgsGeometry QgsGeos::shortestLine( const QgsGeometry &other, QString *errorMsg ) const
2773{
2774 if ( !mGeos || other.isEmpty() )
2775 {
2776 return QgsGeometry();
2777 }
2778
2779 return shortestLine( other.constGet(), errorMsg );
2780}
2781
2782QgsGeometry QgsGeos::shortestLine( const QgsAbstractGeometry *other, QString *errorMsg ) const
2783{
2784 if ( !other || other->isEmpty() )
2785 return QgsGeometry();
2786
2787 geos::unique_ptr otherGeom( asGeos( other, mPrecision ) );
2788 if ( !otherGeom )
2789 {
2790 return QgsGeometry();
2791 }
2792
2793 double nx1 = 0.0;
2794 double ny1 = 0.0;
2795 double nx2 = 0.0;
2796 double ny2 = 0.0;
2797 try
2798 {
2799 geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() ) );
2800
2801 if ( !nearestCoord )
2802 {
2803 if ( errorMsg )
2804 *errorMsg = QStringLiteral( "GEOS returned no nearest points" );
2805 return QgsGeometry();
2806 }
2807
2808 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 0, &nx1 );
2809 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 0, &ny1 );
2810 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 1, &nx2 );
2811 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 1, &ny2 );
2812 }
2813 catch ( GEOSException &e )
2814 {
2815 logError( QStringLiteral( "GEOS" ), e.what() );
2816 if ( errorMsg )
2817 {
2818 *errorMsg = e.what();
2819 }
2820 return QgsGeometry();
2821 }
2822
2823 QgsLineString *line = new QgsLineString();
2824 line->addVertex( QgsPoint( nx1, ny1 ) );
2825 line->addVertex( QgsPoint( nx2, ny2 ) );
2826 return QgsGeometry( line );
2827}
2828
2829double QgsGeos::lineLocatePoint( const QgsPoint &point, QString *errorMsg ) const
2830{
2831 if ( !mGeos )
2832 {
2833 return -1;
2834 }
2835
2836 geos::unique_ptr otherGeom( asGeos( &point, mPrecision ) );
2837 if ( !otherGeom )
2838 {
2839 return -1;
2840 }
2841
2842 double distance = -1;
2843 try
2844 {
2845 distance = GEOSProject_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() );
2846 }
2847 catch ( GEOSException &e )
2848 {
2849 logError( QStringLiteral( "GEOS" ), e.what() );
2850 if ( errorMsg )
2851 {
2852 *errorMsg = e.what();
2853 }
2854 return -1;
2855 }
2856
2857 return distance;
2858}
2859
2860double QgsGeos::lineLocatePoint( double x, double y, QString *errorMsg ) const
2861{
2862 if ( !mGeos )
2863 {
2864 return -1;
2865 }
2866
2867 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
2868 if ( !point )
2869 return false;
2870
2871 double distance = -1;
2872 try
2873 {
2874 distance = GEOSProject_r( geosinit()->ctxt, mGeos.get(), point.get() );
2875 }
2876 catch ( GEOSException &e )
2877 {
2878 logError( QStringLiteral( "GEOS" ), e.what() );
2879 if ( errorMsg )
2880 {
2881 *errorMsg = e.what();
2882 }
2883 return -1;
2884 }
2885
2886 return distance;
2887}
2888
2889QgsGeometry QgsGeos::polygonize( const QVector<const QgsAbstractGeometry *> &geometries, QString *errorMsg )
2890{
2891 GEOSGeometry **const lineGeosGeometries = new GEOSGeometry*[ geometries.size()];
2892 int validLines = 0;
2893 for ( const QgsAbstractGeometry *g : geometries )
2894 {
2895 geos::unique_ptr l = asGeos( g );
2896 if ( l )
2897 {
2898 lineGeosGeometries[validLines] = l.release();
2899 validLines++;
2900 }
2901 }
2902
2903 try
2904 {
2905 geos::unique_ptr result( GEOSPolygonize_r( geosinit()->ctxt, lineGeosGeometries, validLines ) );
2906 for ( int i = 0; i < validLines; ++i )
2907 {
2908 GEOSGeom_destroy_r( geosinit()->ctxt, lineGeosGeometries[i] );
2909 }
2910 delete[] lineGeosGeometries;
2911 return QgsGeometry( fromGeos( result.get() ) );
2912 }
2913 catch ( GEOSException &e )
2914 {
2915 if ( errorMsg )
2916 {
2917 *errorMsg = e.what();
2918 }
2919 for ( int i = 0; i < validLines; ++i )
2920 {
2921 GEOSGeom_destroy_r( geosinit()->ctxt, lineGeosGeometries[i] );
2922 }
2923 delete[] lineGeosGeometries;
2924 return QgsGeometry();
2925 }
2926}
2927
2928QgsGeometry QgsGeos::voronoiDiagram( const QgsAbstractGeometry *extent, double tolerance, bool edgesOnly, QString *errorMsg ) const
2929{
2930 if ( !mGeos )
2931 {
2932 return QgsGeometry();
2933 }
2934
2935 geos::unique_ptr extentGeosGeom;
2936 if ( extent )
2937 {
2938 extentGeosGeom = asGeos( extent, mPrecision );
2939 if ( !extentGeosGeom )
2940 {
2941 return QgsGeometry();
2942 }
2943 }
2944
2946 try
2947 {
2948 geos.reset( GEOSVoronoiDiagram_r( geosinit()->ctxt, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
2949
2950 if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
2951 {
2952 return QgsGeometry();
2953 }
2954
2955 return QgsGeometry( fromGeos( geos.get() ) );
2956 }
2958}
2959
2960QgsGeometry QgsGeos::delaunayTriangulation( double tolerance, bool edgesOnly, QString *errorMsg ) const
2961{
2962 if ( !mGeos )
2963 {
2964 return QgsGeometry();
2965 }
2966
2968 try
2969 {
2970 geos.reset( GEOSDelaunayTriangulation_r( geosinit()->ctxt, mGeos.get(), tolerance, edgesOnly ) );
2971
2972 if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
2973 {
2974 return QgsGeometry();
2975 }
2976
2977 return QgsGeometry( fromGeos( geos.get() ) );
2978 }
2980}
2981
2982
2984static bool _linestringEndpoints( const GEOSGeometry *linestring, double &x1, double &y1, double &x2, double &y2 )
2985{
2986 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, linestring );
2987 if ( !coordSeq )
2988 return false;
2989
2990 unsigned int coordSeqSize;
2991 if ( GEOSCoordSeq_getSize_r( geosinit()->ctxt, coordSeq, &coordSeqSize ) == 0 )
2992 return false;
2993
2994 if ( coordSeqSize < 2 )
2995 return false;
2996
2997 GEOSCoordSeq_getX_r( geosinit()->ctxt, coordSeq, 0, &x1 );
2998 GEOSCoordSeq_getY_r( geosinit()->ctxt, coordSeq, 0, &y1 );
2999 GEOSCoordSeq_getX_r( geosinit()->ctxt, coordSeq, coordSeqSize - 1, &x2 );
3000 GEOSCoordSeq_getY_r( geosinit()->ctxt, coordSeq, coordSeqSize - 1, &y2 );
3001 return true;
3002}
3003
3004
3006static geos::unique_ptr _mergeLinestrings( const GEOSGeometry *line1, const GEOSGeometry *line2, const QgsPointXY &intersectionPoint )
3007{
3008 double x1, y1, x2, y2;
3009 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
3010 return nullptr;
3011
3012 double rx1, ry1, rx2, ry2;
3013 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
3014 return nullptr;
3015
3016 bool intersectionAtOrigLineEndpoint =
3017 ( intersectionPoint.x() == x1 && intersectionPoint.y() == y1 ) !=
3018 ( intersectionPoint.x() == x2 && intersectionPoint.y() == y2 );
3019 bool intersectionAtReshapeLineEndpoint =
3020 ( intersectionPoint.x() == rx1 && intersectionPoint.y() == ry1 ) ||
3021 ( intersectionPoint.x() == rx2 && intersectionPoint.y() == ry2 );
3022
3023 // the intersection must be at the begin/end of both lines
3024 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
3025 {
3026 geos::unique_ptr g1( GEOSGeom_clone_r( geosinit()->ctxt, line1 ) );
3027 geos::unique_ptr g2( GEOSGeom_clone_r( geosinit()->ctxt, line2 ) );
3028 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
3029 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, geoms, 2 ) );
3030 geos::unique_ptr res( GEOSLineMerge_r( geosinit()->ctxt, multiGeom.get() ) );
3031 return res;
3032 }
3033 else
3034 return nullptr;
3035}
3036
3037
3038geos::unique_ptr QgsGeos::reshapeLine( const GEOSGeometry *line, const GEOSGeometry *reshapeLineGeos, double precision )
3039{
3040 if ( !line || !reshapeLineGeos )
3041 return nullptr;
3042
3043 bool atLeastTwoIntersections = false;
3044 bool oneIntersection = false;
3045 QgsPointXY oneIntersectionPoint;
3046
3047 try
3048 {
3049 //make sure there are at least two intersection between line and reshape geometry
3050 geos::unique_ptr intersectGeom( GEOSIntersection_r( geosinit()->ctxt, line, reshapeLineGeos ) );
3051 if ( intersectGeom )
3052 {
3053 const int geomType = GEOSGeomTypeId_r( geosinit()->ctxt, intersectGeom.get() );
3054 atLeastTwoIntersections = ( geomType == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 1 )
3055 || ( geomType == GEOS_GEOMETRYCOLLECTION && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 0 ) // a collection implies at least two points!
3056 || ( geomType == GEOS_MULTILINESTRING && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 0 );
3057 // one point is enough when extending line at its endpoint
3058 if ( GEOSGeomTypeId_r( geosinit()->ctxt, intersectGeom.get() ) == GEOS_POINT )
3059 {
3060 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, intersectGeom.get() );
3061 double xi, yi;
3062 GEOSCoordSeq_getX_r( geosinit()->ctxt, intersectionCoordSeq, 0, &xi );
3063 GEOSCoordSeq_getY_r( geosinit()->ctxt, intersectionCoordSeq, 0, &yi );
3064 oneIntersection = true;
3065 oneIntersectionPoint = QgsPointXY( xi, yi );
3066 }
3067 }
3068 }
3069 catch ( GEOSException & )
3070 {
3071 atLeastTwoIntersections = false;
3072 }
3073
3074 // special case when extending line at its endpoint
3075 if ( oneIntersection )
3076 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
3077
3078 if ( !atLeastTwoIntersections )
3079 return nullptr;
3080
3081 //begin and end point of original line
3082 double x1, y1, x2, y2;
3083 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
3084 return nullptr;
3085
3086 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1, false, 0, false, 0, 2, precision );
3087 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2, false, 0, false, 0, 2, precision );
3088
3089 bool isRing = false;
3090 if ( GEOSGeomTypeId_r( geosinit()->ctxt, line ) == GEOS_LINEARRING
3091 || GEOSEquals_r( geosinit()->ctxt, beginLineVertex.get(), endLineVertex.get() ) == 1 )
3092 isRing = true;
3093
3094 //node line and reshape line
3095 geos::unique_ptr nodedGeometry = nodeGeometries( reshapeLineGeos, line );
3096 if ( !nodedGeometry )
3097 {
3098 return nullptr;
3099 }
3100
3101 //and merge them together
3102 geos::unique_ptr mergedLines( GEOSLineMerge_r( geosinit()->ctxt, nodedGeometry.get() ) );
3103 if ( !mergedLines )
3104 {
3105 return nullptr;
3106 }
3107
3108 int numMergedLines = GEOSGetNumGeometries_r( geosinit()->ctxt, mergedLines.get() );
3109 if ( numMergedLines < 2 ) //some special cases. Normally it is >2
3110 {
3111 if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
3112 {
3113 geos::unique_ptr result( GEOSGeom_clone_r( geosinit()->ctxt, reshapeLineGeos ) );
3114 return result;
3115 }
3116 else
3117 return nullptr;
3118 }
3119
3120 QVector<GEOSGeometry *> resultLineParts; //collection with the line segments that will be contained in result
3121 QVector<GEOSGeometry *> probableParts; //parts where we can decide on inclusion only after going through all the candidates
3122
3123 for ( int i = 0; i < numMergedLines; ++i )
3124 {
3125 const GEOSGeometry *currentGeom = GEOSGetGeometryN_r( geosinit()->ctxt, mergedLines.get(), i );
3126
3127 // have we already added this part?
3128 bool alreadyAdded = false;
3129 double distance = 0;
3130 double bufferDistance = std::pow( 10.0L, geomDigits( currentGeom ) - 11 );
3131 for ( const GEOSGeometry *other : std::as_const( resultLineParts ) )
3132 {
3133 GEOSHausdorffDistance_r( geosinit()->ctxt, currentGeom, other, &distance );
3134 if ( distance < bufferDistance )
3135 {
3136 alreadyAdded = true;
3137 break;
3138 }
3139 }
3140 if ( alreadyAdded )
3141 continue;
3142
3143 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, currentGeom );
3144 unsigned int currentCoordSeqSize;
3145 GEOSCoordSeq_getSize_r( geosinit()->ctxt, currentCoordSeq, &currentCoordSeqSize );
3146 if ( currentCoordSeqSize < 2 )
3147 continue;
3148
3149 //get the two endpoints of the current line merge result
3150 double xBegin, xEnd, yBegin, yEnd;
3151 GEOSCoordSeq_getX_r( geosinit()->ctxt, currentCoordSeq, 0, &xBegin );
3152 GEOSCoordSeq_getY_r( geosinit()->ctxt, currentCoordSeq, 0, &yBegin );
3153 GEOSCoordSeq_getX_r( geosinit()->ctxt, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
3154 GEOSCoordSeq_getY_r( geosinit()->ctxt, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
3155 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin, false, 0, false, 0, 2, precision );
3156 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd, false, 0, false, 0, 2, precision );
3157
3158 //check how many endpoints of the line merge result are on the (original) line
3159 int nEndpointsOnOriginalLine = 0;
3160 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
3161 nEndpointsOnOriginalLine += 1;
3162
3163 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
3164 nEndpointsOnOriginalLine += 1;
3165
3166 //check how many endpoints equal the endpoints of the original line
3167 int nEndpointsSameAsOriginalLine = 0;
3168 if ( GEOSEquals_r( geosinit()->ctxt, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3169 || GEOSEquals_r( geosinit()->ctxt, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3170 nEndpointsSameAsOriginalLine += 1;
3171
3172 if ( GEOSEquals_r( geosinit()->ctxt, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3173 || GEOSEquals_r( geosinit()->ctxt, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3174 nEndpointsSameAsOriginalLine += 1;
3175
3176 //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
3177 bool currentGeomOverlapsOriginalGeom = false;
3178 bool currentGeomOverlapsReshapeLine = false;
3179 if ( lineContainedInLine( currentGeom, line ) == 1 )
3180 currentGeomOverlapsOriginalGeom = true;
3181
3182 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
3183 currentGeomOverlapsReshapeLine = true;
3184
3185 //logic to decide if this part belongs to the result
3186 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3187 {
3188 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3189 }
3190 //for closed rings, we take one segment from the candidate list
3191 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3192 {
3193 probableParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3194 }
3195 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3196 {
3197 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3198 }
3199 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3200 {
3201 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3202 }
3203 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
3204 {
3205 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3206 }
3207 }
3208
3209 //add the longest segment from the probable list for rings (only used for polygon rings)
3210 if ( isRing && !probableParts.isEmpty() )
3211 {
3212 geos::unique_ptr maxGeom; //the longest geometry in the probabla list
3213 GEOSGeometry *currentGeom = nullptr;
3214 double maxLength = -std::numeric_limits<double>::max();
3215 double currentLength = 0;
3216 for ( int i = 0; i < probableParts.size(); ++i )
3217 {
3218 currentGeom = probableParts.at( i );
3219 GEOSLength_r( geosinit()->ctxt, currentGeom, &currentLength );
3220 if ( currentLength > maxLength )
3221 {
3222 maxLength = currentLength;
3223 maxGeom.reset( currentGeom );
3224 }
3225 else
3226 {
3227 GEOSGeom_destroy_r( geosinit()->ctxt, currentGeom );
3228 }
3229 }
3230 resultLineParts.push_back( maxGeom.release() );
3231 }
3232
3233 geos::unique_ptr result;
3234 if ( resultLineParts.empty() )
3235 return nullptr;
3236
3237 if ( resultLineParts.size() == 1 ) //the whole result was reshaped
3238 {
3239 result.reset( resultLineParts[0] );
3240 }
3241 else //>1
3242 {
3243 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
3244 for ( int i = 0; i < resultLineParts.size(); ++i )
3245 {
3246 lineArray[i] = resultLineParts[i];
3247 }
3248
3249 //create multiline from resultLineParts
3250 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
3251 delete [] lineArray;
3252
3253 //then do a linemerge with the newly combined partstrings
3254 result.reset( GEOSLineMerge_r( geosinit()->ctxt, multiLineGeom.get() ) );
3255 }
3256
3257 //now test if the result is a linestring. Otherwise something went wrong
3258 if ( GEOSGeomTypeId_r( geosinit()->ctxt, result.get() ) != GEOS_LINESTRING )
3259 {
3260 return nullptr;
3261 }
3262
3263 return result;
3264}
3265
3266geos::unique_ptr QgsGeos::reshapePolygon( const GEOSGeometry *polygon, const GEOSGeometry *reshapeLineGeos, double precision )
3267{
3268 //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
3269 int nIntersections = 0;
3270 int lastIntersectingRing = -2;
3271 const GEOSGeometry *lastIntersectingGeom = nullptr;
3272
3273 int nRings = GEOSGetNumInteriorRings_r( geosinit()->ctxt, polygon );
3274 if ( nRings < 0 )
3275 return nullptr;
3276
3277 //does outer ring intersect?
3278 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( geosinit()->ctxt, polygon );
3279 if ( GEOSIntersects_r( geosinit()->ctxt, outerRing, reshapeLineGeos ) == 1 )
3280 {
3281 ++nIntersections;
3282 lastIntersectingRing = -1;
3283 lastIntersectingGeom = outerRing;
3284 }
3285
3286 //do inner rings intersect?
3287 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
3288
3289 try
3290 {
3291 for ( int i = 0; i < nRings; ++i )
3292 {
3293 innerRings[i] = GEOSGetInteriorRingN_r( geosinit()->ctxt, polygon, i );
3294 if ( GEOSIntersects_r( geosinit()->ctxt, innerRings[i], reshapeLineGeos ) == 1 )
3295 {
3296 ++nIntersections;
3297 lastIntersectingRing = i;
3298 lastIntersectingGeom = innerRings[i];
3299 }
3300 }
3301 }
3302 catch ( GEOSException & )
3303 {
3304 nIntersections = 0;
3305 }
3306
3307 if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
3308 {
3309 delete [] innerRings;
3310 return nullptr;
3311 }
3312
3313 //we have one intersecting ring, let's try to reshape it
3314 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos, precision );
3315 if ( !reshapeResult )
3316 {
3317 delete [] innerRings;
3318 return nullptr;
3319 }
3320
3321 //if reshaping took place, we need to reassemble the polygon and its rings
3322 GEOSGeometry *newRing = nullptr;
3323 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, reshapeResult.get() );
3324 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( geosinit()->ctxt, reshapeSequence );
3325
3326 reshapeResult.reset();
3327
3328 newRing = GEOSGeom_createLinearRing_r( geosinit()->ctxt, newCoordSequence );
3329 if ( !newRing )
3330 {
3331 delete [] innerRings;
3332 return nullptr;
3333 }
3334
3335 GEOSGeometry *newOuterRing = nullptr;
3336 if ( lastIntersectingRing == -1 )
3337 newOuterRing = newRing;
3338 else
3339 newOuterRing = GEOSGeom_clone_r( geosinit()->ctxt, outerRing );
3340
3341 //check if all the rings are still inside the outer boundary
3342 QVector<GEOSGeometry *> ringList;
3343 if ( nRings > 0 )
3344 {
3345 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( geosinit()->ctxt, GEOSGeom_clone_r( geosinit()->ctxt, newOuterRing ), nullptr, 0 );
3346 if ( outerRingPoly )
3347 {
3348 ringList.reserve( nRings );
3349 GEOSGeometry *currentRing = nullptr;
3350 for ( int i = 0; i < nRings; ++i )
3351 {
3352 if ( lastIntersectingRing == i )
3353 currentRing = newRing;
3354 else
3355 currentRing = GEOSGeom_clone_r( geosinit()->ctxt, innerRings[i] );
3356
3357 //possibly a ring is no longer contained in the result polygon after reshape
3358 if ( GEOSContains_r( geosinit()->ctxt, outerRingPoly, currentRing ) == 1 )
3359 ringList.push_back( currentRing );
3360 else
3361 GEOSGeom_destroy_r( geosinit()->ctxt, currentRing );
3362 }
3363 }
3364 GEOSGeom_destroy_r( geosinit()->ctxt, outerRingPoly );
3365 }
3366
3367 GEOSGeometry **newInnerRings = new GEOSGeometry*[ringList.size()];
3368 for ( int i = 0; i < ringList.size(); ++i )
3369 newInnerRings[i] = ringList.at( i );
3370
3371 delete [] innerRings;
3372
3373 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( geosinit()->ctxt, newOuterRing, newInnerRings, ringList.size() ) );
3374 delete[] newInnerRings;
3375
3376 return reshapedPolygon;
3377}
3378
3379int QgsGeos::lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 )
3380{
3381 if ( !line1 || !line2 )
3382 {
3383 return -1;
3384 }
3385
3386 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
3387
3388 geos::unique_ptr bufferGeom( GEOSBuffer_r( geosinit()->ctxt, line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ) );
3389 if ( !bufferGeom )
3390 return -2;
3391
3392 geos::unique_ptr intersectionGeom( GEOSIntersection_r( geosinit()->ctxt, bufferGeom.get(), line1 ) );
3393
3394 //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
3395 double intersectGeomLength;
3396 double line1Length;
3397
3398 GEOSLength_r( geosinit()->ctxt, intersectionGeom.get(), &intersectGeomLength );
3399 GEOSLength_r( geosinit()->ctxt, line1, &line1Length );
3400
3401 double intersectRatio = line1Length / intersectGeomLength;
3402 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
3403 return 1;
3404
3405 return 0;
3406}
3407
3408int QgsGeos::pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line )
3409{
3410 if ( !point || !line )
3411 return -1;
3412
3413 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
3414
3415 geos::unique_ptr lineBuffer( GEOSBuffer_r( geosinit()->ctxt, line, bufferDistance, 8 ) );
3416 if ( !lineBuffer )
3417 return -2;
3418
3419 bool contained = false;
3420 if ( GEOSContains_r( geosinit()->ctxt, lineBuffer.get(), point ) == 1 )
3421 contained = true;
3422
3423 return contained;
3424}
3425
3426int QgsGeos::geomDigits( const GEOSGeometry *geom )
3427{
3428 geos::unique_ptr bbox( GEOSEnvelope_r( geosinit()->ctxt, geom ) );
3429 if ( !bbox.get() )
3430 return -1;
3431
3432 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( geosinit()->ctxt, bbox.get() );
3433 if ( !bBoxRing )
3434 return -1;
3435
3436 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, bBoxRing );
3437
3438 if ( !bBoxCoordSeq )
3439 return -1;
3440
3441 unsigned int nCoords = 0;
3442 if ( !GEOSCoordSeq_getSize_r( geosinit()->ctxt, bBoxCoordSeq, &nCoords ) )
3443 return -1;
3444
3445 int maxDigits = -1;
3446 for ( unsigned int i = 0; i < nCoords - 1; ++i )
3447 {
3448 double t;
3449 GEOSCoordSeq_getX_r( geosinit()->ctxt, bBoxCoordSeq, i, &t );
3450
3451 int digits;
3452 digits = std::ceil( std::log10( std::fabs( t ) ) );
3453 if ( digits > maxDigits )
3454 maxDigits = digits;
3455
3456 GEOSCoordSeq_getY_r( geosinit()->ctxt, bBoxCoordSeq, i, &t );
3457 digits = std::ceil( std::log10( std::fabs( t ) ) );
3458 if ( digits > maxDigits )
3459 maxDigits = digits;
3460 }
3461
3462 return maxDigits;
3463}
3464
3465GEOSContextHandle_t QgsGeos::getGEOSHandler()
3466{
3467 return geosinit()->ctxt;
3468}
The Qgis class provides global constants for use throughout the application.
Definition qgis.h:54
BufferSide
Side of line to buffer.
Definition qgis.h:1574
@ Right
Buffer to right of line.
GeometryOperationResult
Success or failure of a geometry operation.
Definition qgis.h:1520
@ 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:1599
EndCapStyle
End cap styles for buffers.
Definition qgis.h:1586
MakeValidMethod
Algorithms to use when repairing invalid geometries.
Definition qgis.h:1612
@ 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
@ 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:36
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:99
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:2960
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:2046
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
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:2594
std::unique_ptr< QgsAbstractGeometry > maximumInscribedCircle(double tolerance, QString *errorMsg=nullptr) const
Returns the maximum inscribed circle.
Definition qgsgeos.cpp:2472
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:2574
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:2558
double minimumClearance(QString *errorMsg=nullptr) const
Computes the minimum clearance of a geometry.
Definition qgsgeos.cpp:2524
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:2446
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:2772
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:1908
QgsGeometry closestPoint(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the closest point on the geometry to the other geometry.
Definition qgsgeos.cpp:2729
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:2542
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:2508
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
Definition qgsgeos.cpp:2138
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:2104
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:2710
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:2889
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:2427
static GEOSContextHandle_t getGEOSHandler()
Definition qgsgeos.cpp:3465
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:2829
bool isEmpty(QString *errorMsg=nullptr) const override
Definition qgsgeos.cpp:2124
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:2928
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:2488
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:59
double y
Definition qgspointxy.h:63
double x
Definition qgspointxy.h:62
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:343
double m
Definition qgspoint.h:55
double y
Definition qgspoint.h:53
Polygon geometry type.
Definition qgspolygon.h:34
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:37
std::unique_ptr< GEOSGeometry, GeosDeleter > unique_ptr
Scoped GEOS pointer.
Definition qgsgeos.h:74
std::unique_ptr< GEOSCoordSequence, GeosDeleter > coord_sequence_unique_ptr
Scoped GEOS coordinate sequence pointer.
Definition qgsgeos.h:89
std::unique_ptr< GEOSBufferParams, GeosDeleter > buffer_params_unique_ptr
Scoped GEOS buffer params pointer.
Definition qgsgeos.h:84
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:4332
QMap< QString, QString > QgsStringMap
Definition qgis.h:4877
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:31
int vertex
Vertex number.
Definition qgsvertexid.h:95
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