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