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