QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
qgsgeos.cpp
Go to the documentation of this file.
1/***************************************************************************
2 qgsgeos.cpp
3 -------------------------------------------------------------------
4Date : 22 Sept 2014
5Copyright : (C) 2014 by Marco Hugentobler
6email : marco.hugentobler at sourcepole dot com
7 ***************************************************************************
8 * *
9 * This program is free software; you can redistribute it and/or modify *
10 * it under the terms of the GNU General Public License as published by *
11 * the Free Software Foundation; either version 2 of the License, or *
12 * (at your option) any later version. *
13 * *
14 ***************************************************************************/
15
16#include "qgsgeos.h"
17#include "qgsabstractgeometry.h"
19#include "qgsgeometryfactory.h"
20#include "qgslinestring.h"
21#include "qgsmulticurve.h"
22#include "qgsmultilinestring.h"
23#include "qgsmultipoint.h"
24#include "qgsmultipolygon.h"
25#include "qgslogger.h"
26#include "qgspolygon.h"
28#include <limits>
29#include <cstdio>
30
31#define DEFAULT_QUADRANT_SEGMENTS 8
32
33#define CATCH_GEOS(r) \
34 catch (GEOSException &) \
35 { \
36 return r; \
37 }
38
39#define CATCH_GEOS_WITH_ERRMSG(r) \
40 catch (GEOSException &e) \
41 { \
42 if ( errorMsg ) \
43 { \
44 *errorMsg = e.what(); \
45 } \
46 return r; \
47 }
48
50
51static void throwGEOSException( const char *fmt, ... )
52{
53 va_list ap;
54 char buffer[1024];
55
56 va_start( ap, fmt );
57 vsnprintf( buffer, sizeof buffer, fmt, ap );
58 va_end( ap );
59
60 QString message = QString::fromUtf8( buffer );
61
62#ifdef _MSC_VER
63 // stupid stupid MSVC, *SOMETIMES* raises it's own exception if we throw GEOSException, resulting in a crash!
64 // see https://github.com/qgis/QGIS/issues/22709
65 // if you want to test alternative fixes for this, run the testqgsexpression.cpp test suite - that will crash
66 // and burn on the "line_interpolate_point point" test if a GEOSException is thrown.
67 // TODO - find a real fix for the underlying issue
68 try
69 {
70 throw GEOSException( message );
71 }
72 catch ( ... )
73 {
74 // oops, msvc threw an exception when we tried to throw the exception!
75 // just throw nothing instead (except your mouse at your monitor)
76 }
77#else
78 throw GEOSException( message );
79#endif
80}
81
82
83static void printGEOSNotice( const char *fmt, ... )
84{
85#if defined(QGISDEBUG)
86 va_list ap;
87 char buffer[1024];
88
89 va_start( ap, fmt );
90 vsnprintf( buffer, sizeof buffer, fmt, ap );
91 va_end( ap );
92#else
93 Q_UNUSED( fmt )
94#endif
95}
96
97class GEOSInit
98{
99 public:
100 GEOSContextHandle_t ctxt;
101
102 GEOSInit()
103 {
104 ctxt = initGEOS_r( printGEOSNotice, throwGEOSException );
105 }
106
107 ~GEOSInit()
108 {
109 finishGEOS_r( ctxt );
110 }
111
112 GEOSInit( const GEOSInit &rh ) = delete;
113 GEOSInit &operator=( const GEOSInit &rh ) = delete;
114};
115
116Q_GLOBAL_STATIC( GEOSInit, geosinit )
117
118void geos::GeosDeleter::operator()( GEOSGeometry *geom ) const
119{
120 GEOSGeom_destroy_r( geosinit()->ctxt, geom );
121}
122
123void geos::GeosDeleter::operator()( const GEOSPreparedGeometry *geom ) const
124{
125 GEOSPreparedGeom_destroy_r( geosinit()->ctxt, geom );
126}
127
128void geos::GeosDeleter::operator()( GEOSBufferParams *params ) const
129{
130 GEOSBufferParams_destroy_r( geosinit()->ctxt, params );
131}
132
133void geos::GeosDeleter::operator()( GEOSCoordSequence *sequence ) const
134{
135 GEOSCoordSeq_destroy_r( geosinit()->ctxt, sequence );
136}
137
138
140
141
143 : QgsGeometryEngine( geometry )
144 , mGeos( nullptr )
145 , mPrecision( precision )
146{
147 cacheGeos();
148}
149
151{
153 GEOSGeom_destroy_r( QgsGeos::getGEOSHandler(), geos );
154 return g;
155}
156
158{
159 QgsGeometry g( QgsGeos::fromGeos( geos.get() ) );
160 return g;
161}
162
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;
980 if ( GEOSCoordSeq_getSize_r( geosinit()->ctxt, lineSequence, &sequenceSize ) != 0 )
981 {
982 for ( unsigned int i = 0; i < sequenceSize; ++i )
983 {
984 if ( GEOSCoordSeq_getX_r( geosinit()->ctxt, lineSequence, i, &x ) != 0 )
985 {
986 if ( GEOSCoordSeq_getY_r( geosinit()->ctxt, lineSequence, i, &y ) != 0 )
987 {
988 testPoints.push_back( QgsPoint( x, y ) );
989 }
990 }
991 }
992 }
993 }
994 }
996
997 return true;
998}
999
1000geos::unique_ptr QgsGeos::linePointDifference( GEOSGeometry *GEOSsplitPoint ) const
1001{
1002 int type = GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() );
1003
1004 std::unique_ptr< QgsMultiCurve > multiCurve;
1005 if ( type == GEOS_MULTILINESTRING )
1006 {
1007 multiCurve.reset( qgsgeometry_cast<QgsMultiCurve *>( mGeometry->clone() ) );
1008 }
1009 else if ( type == GEOS_LINESTRING )
1010 {
1011 multiCurve.reset( new QgsMultiCurve() );
1012 multiCurve->addGeometry( mGeometry->clone() );
1013 }
1014 else
1015 {
1016 return nullptr;
1017 }
1018
1019 if ( !multiCurve )
1020 {
1021 return nullptr;
1022 }
1023
1024
1025 // we might have a point or a multipoint, depending on number of
1026 // intersections between the geometry and the split geometry
1027 std::unique_ptr< QgsAbstractGeometry > splitGeom( fromGeos( GEOSsplitPoint ) );
1028 std::unique_ptr< QgsMultiPoint > splitPoints( qgsgeometry_cast<QgsMultiPoint *>( splitGeom->clone() ) );
1029 if ( !splitPoints )
1030 {
1031 QgsPoint *splitPoint = qgsgeometry_cast<QgsPoint *>( splitGeom->clone() );
1032 if ( !splitPoint )
1033 {
1034 return nullptr;
1035 }
1036 splitPoints.reset( new QgsMultiPoint() );
1037 splitPoints->addGeometry( splitPoint );
1038 }
1039
1040 QgsMultiCurve lines;
1041
1042 //For each part
1043 for ( int i = 0; i < multiCurve->numGeometries(); ++i )
1044 {
1045 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( multiCurve->geometryN( i ) );
1046 if ( !line )
1047 {
1048 const QgsCurve *curve = qgsgeometry_cast<const QgsCurve *>( multiCurve->geometryN( i ) );
1049 line = curve->curveToLine();
1050 }
1051 if ( !line )
1052 {
1053 return nullptr;
1054 }
1055 // we gather the intersection points and their distance from previous node grouped by segment
1056 QMap< int, QVector< QPair< double, QgsPoint > > >pointMap;
1057 for ( int p = 0; p < splitPoints->numGeometries(); ++p )
1058 {
1059 const QgsPoint *intersectionPoint = splitPoints->pointN( p );
1060 QgsPoint segmentPoint2D;
1061 QgsVertexId nextVertex;
1062 // With closestSegment we only get a 2D point so we need to interpolate if we
1063 // don't want to lose Z data
1064 line->closestSegment( *intersectionPoint, segmentPoint2D, nextVertex );
1065 const QgsLineString segment = QgsLineString( line->pointN( nextVertex.vertex - 1 ), line->pointN( nextVertex.vertex ) );
1066 const double distance = segmentPoint2D.distance( line->pointN( nextVertex.vertex - 1 ) );
1067 std::unique_ptr< QgsPoint > correctSegmentPoint( segment.interpolatePoint( distance ) );
1068 const QPair< double, QgsPoint > pair = qMakePair( distance, *correctSegmentPoint.get() );
1069 if ( pointMap.contains( nextVertex.vertex - 1 ) )
1070 pointMap[ nextVertex.vertex - 1 ].append( pair );
1071 else
1072 pointMap[ nextVertex.vertex - 1 ] = QVector< QPair< double, QgsPoint > >() << pair;
1073 }
1074
1075 // When we have more than one intersection point on a segment we need those points
1076 // to be sorted by their distance from the previous geometry vertex
1077 for ( auto &p : pointMap )
1078 {
1079 std::sort( p.begin(), p.end(), []( const QPair< double, QgsPoint > &a, const QPair< double, QgsPoint > &b ) { return a.first < b.first; } );
1080 }
1081
1082 //For each segment
1083 QgsLineString newLine;
1084 int nVertices = line->numPoints();
1085 QgsPoint splitPoint;
1086 for ( int j = 0; j < nVertices; ++j )
1087 {
1088 QgsPoint currentPoint = line->pointN( j );
1089 newLine.addVertex( currentPoint );
1090 if ( pointMap.contains( j ) )
1091 {
1092 // For each intersecting point
1093 for ( int k = 0; k < pointMap[ j ].size(); ++k )
1094 {
1095 splitPoint = pointMap[ j ][k].second;
1096 if ( splitPoint == currentPoint )
1097 {
1098 lines.addGeometry( newLine.clone() );
1099 newLine = QgsLineString();
1100 newLine.addVertex( currentPoint );
1101 }
1102 else if ( splitPoint == line->pointN( j + 1 ) )
1103 {
1104 newLine.addVertex( line->pointN( j + 1 ) );
1105 lines.addGeometry( newLine.clone() );
1106 newLine = QgsLineString();
1107 }
1108 else
1109 {
1110 newLine.addVertex( splitPoint );
1111 lines.addGeometry( newLine.clone() );
1112 newLine = QgsLineString();
1113 newLine.addVertex( splitPoint );
1114 }
1115 }
1116 }
1117 }
1118 lines.addGeometry( newLine.clone() );
1119 }
1120
1121 return asGeos( &lines, mPrecision );
1122}
1123
1124QgsGeometryEngine::EngineOperationResult QgsGeos::splitLinearGeometry( GEOSGeometry *splitLine, QVector<QgsGeometry> &newGeometries, bool skipIntersectionCheck ) const
1125{
1126 Q_UNUSED( skipIntersectionCheck )
1127 if ( !splitLine )
1128 return InvalidInput;
1129
1130 if ( !mGeos )
1131 return InvalidBaseGeometry;
1132
1133 geos::unique_ptr intersectGeom( GEOSIntersection_r( geosinit()->ctxt, splitLine, mGeos.get() ) );
1134 if ( !intersectGeom )
1135 return NothingHappened;
1136
1137 //check that split line has no linear intersection
1138 int linearIntersect = GEOSRelatePattern_r( geosinit()->ctxt, mGeos.get(), splitLine, "1********" );
1139 if ( linearIntersect > 0 )
1140 return InvalidInput;
1141
1142 geos::unique_ptr splitGeom;
1143 splitGeom = linePointDifference( intersectGeom.get() );
1144
1145 if ( !splitGeom )
1146 return InvalidBaseGeometry;
1147
1148 QVector<GEOSGeometry *> lineGeoms;
1149
1150 int splitType = GEOSGeomTypeId_r( geosinit()->ctxt, splitGeom.get() );
1151 if ( splitType == GEOS_MULTILINESTRING )
1152 {
1153 int nGeoms = GEOSGetNumGeometries_r( geosinit()->ctxt, splitGeom.get() );
1154 lineGeoms.reserve( nGeoms );
1155 for ( int i = 0; i < nGeoms; ++i )
1156 lineGeoms << GEOSGeom_clone_r( geosinit()->ctxt, GEOSGetGeometryN_r( geosinit()->ctxt, splitGeom.get(), i ) );
1157
1158 }
1159 else
1160 {
1161 lineGeoms << GEOSGeom_clone_r( geosinit()->ctxt, splitGeom.get() );
1162 }
1163
1164 mergeGeometriesMultiTypeSplit( lineGeoms );
1165
1166 for ( int i = 0; i < lineGeoms.size(); ++i )
1167 {
1168 newGeometries << QgsGeometry( fromGeos( lineGeoms[i] ) );
1169 GEOSGeom_destroy_r( geosinit()->ctxt, lineGeoms[i] );
1170 }
1171
1172 return Success;
1173}
1174
1175QgsGeometryEngine::EngineOperationResult QgsGeos::splitPolygonGeometry( GEOSGeometry *splitLine, QVector<QgsGeometry> &newGeometries, bool skipIntersectionCheck ) const
1176{
1177 if ( !splitLine )
1178 return InvalidInput;
1179
1180 if ( !mGeos )
1181 return InvalidBaseGeometry;
1182
1183 // we will need prepared geometry for intersection tests
1184 const_cast<QgsGeos *>( this )->prepareGeometry();
1185 if ( !mGeosPrepared )
1186 return EngineError;
1187
1188 //first test if linestring intersects geometry. If not, return straight away
1189 if ( !skipIntersectionCheck && !GEOSPreparedIntersects_r( geosinit()->ctxt, mGeosPrepared.get(), splitLine ) )
1190 return NothingHappened;
1191
1192 //first union all the polygon rings together (to get them noded, see JTS developer guide)
1193 geos::unique_ptr nodedGeometry = nodeGeometries( splitLine, mGeos.get() );
1194 if ( !nodedGeometry )
1195 return NodedGeometryError; //an error occurred during noding
1196
1197 const GEOSGeometry *noded = nodedGeometry.get();
1198 geos::unique_ptr polygons( GEOSPolygonize_r( geosinit()->ctxt, &noded, 1 ) );
1199 if ( !polygons || numberOfGeometries( polygons.get() ) == 0 )
1200 {
1201 return InvalidBaseGeometry;
1202 }
1203
1204 //test every polygon is contained in original geometry
1205 //include in result if yes
1206 QVector<GEOSGeometry *> testedGeometries;
1207
1208 // test whether the polygon parts returned from polygonize algorithm actually
1209 // belong to the source polygon geometry (if the source polygon contains some holes,
1210 // those would be also returned by polygonize and we need to skip them)
1211 for ( int i = 0; i < numberOfGeometries( polygons.get() ); i++ )
1212 {
1213 const GEOSGeometry *polygon = GEOSGetGeometryN_r( geosinit()->ctxt, polygons.get(), i );
1214
1215 geos::unique_ptr pointOnSurface( GEOSPointOnSurface_r( geosinit()->ctxt, polygon ) );
1216 if ( pointOnSurface && GEOSPreparedIntersects_r( geosinit()->ctxt, mGeosPrepared.get(), pointOnSurface.get() ) )
1217 testedGeometries << GEOSGeom_clone_r( geosinit()->ctxt, polygon );
1218 }
1219
1220 int nGeometriesThis = numberOfGeometries( mGeos.get() ); //original number of geometries
1221 if ( testedGeometries.empty() || testedGeometries.size() == nGeometriesThis )
1222 {
1223 //no split done, preserve original geometry
1224 for ( int i = 0; i < testedGeometries.size(); ++i )
1225 {
1226 GEOSGeom_destroy_r( geosinit()->ctxt, testedGeometries[i] );
1227 }
1228 return NothingHappened;
1229 }
1230
1231 // For multi-part geometries, try to identify parts that have been unchanged and try to merge them back
1232 // to a single multi-part geometry. For example, if there's a multi-polygon with three parts, but only
1233 // one part is being split, this function makes sure that the other two parts will be kept in a multi-part
1234 // geometry rather than being separated into two single-part geometries.
1235 mergeGeometriesMultiTypeSplit( testedGeometries );
1236
1237 int i;
1238 for ( i = 0; i < testedGeometries.size() && GEOSisValid_r( geosinit()->ctxt, testedGeometries[i] ); ++i )
1239 ;
1240
1241 if ( i < testedGeometries.size() )
1242 {
1243 for ( i = 0; i < testedGeometries.size(); ++i )
1244 GEOSGeom_destroy_r( geosinit()->ctxt, testedGeometries[i] );
1245
1246 return InvalidBaseGeometry;
1247 }
1248
1249 for ( i = 0; i < testedGeometries.size(); ++i )
1250 {
1251 newGeometries << QgsGeometry( fromGeos( testedGeometries[i] ) );
1252 GEOSGeom_destroy_r( geosinit()->ctxt, testedGeometries[i] );
1253 }
1254
1255 return Success;
1256}
1257
1258geos::unique_ptr QgsGeos::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
1259{
1260 if ( !splitLine || !geom )
1261 return nullptr;
1262
1263 geos::unique_ptr geometryBoundary;
1264 if ( GEOSGeomTypeId_r( geosinit()->ctxt, geom ) == GEOS_POLYGON || GEOSGeomTypeId_r( geosinit()->ctxt, geom ) == GEOS_MULTIPOLYGON )
1265 geometryBoundary.reset( GEOSBoundary_r( geosinit()->ctxt, geom ) );
1266 else
1267 geometryBoundary.reset( GEOSGeom_clone_r( geosinit()->ctxt, geom ) );
1268
1269 geos::unique_ptr splitLineClone( GEOSGeom_clone_r( geosinit()->ctxt, splitLine ) );
1270 geos::unique_ptr unionGeometry( GEOSUnion_r( geosinit()->ctxt, splitLineClone.get(), geometryBoundary.get() ) );
1271
1272 return unionGeometry;
1273}
1274
1275int QgsGeos::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry *> &splitResult ) const
1276{
1277 if ( !mGeos )
1278 return 1;
1279
1280 //convert mGeos to geometry collection
1281 int type = GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() );
1282 if ( type != GEOS_GEOMETRYCOLLECTION &&
1283 type != GEOS_MULTILINESTRING &&
1284 type != GEOS_MULTIPOLYGON &&
1285 type != GEOS_MULTIPOINT )
1286 return 0;
1287
1288 QVector<GEOSGeometry *> copyList = splitResult;
1289 splitResult.clear();
1290
1291 //collect all the geometries that belong to the initial multifeature
1292 QVector<GEOSGeometry *> unionGeom;
1293
1294 for ( int i = 0; i < copyList.size(); ++i )
1295 {
1296 //is this geometry a part of the original multitype?
1297 bool isPart = false;
1298 for ( int j = 0; j < GEOSGetNumGeometries_r( geosinit()->ctxt, mGeos.get() ); j++ )
1299 {
1300 if ( GEOSEquals_r( geosinit()->ctxt, copyList[i], GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), j ) ) )
1301 {
1302 isPart = true;
1303 break;
1304 }
1305 }
1306
1307 if ( isPart )
1308 {
1309 unionGeom << copyList[i];
1310 }
1311 else
1312 {
1313 QVector<GEOSGeometry *> geomVector;
1314 geomVector << copyList[i];
1315
1316 if ( type == GEOS_MULTILINESTRING )
1317 splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector ).release();
1318 else if ( type == GEOS_MULTIPOLYGON )
1319 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector ).release();
1320 else
1321 GEOSGeom_destroy_r( geosinit()->ctxt, copyList[i] );
1322 }
1323 }
1324
1325 //make multifeature out of unionGeom
1326 if ( !unionGeom.isEmpty() )
1327 {
1328 if ( type == GEOS_MULTILINESTRING )
1329 splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom ).release();
1330 else if ( type == GEOS_MULTIPOLYGON )
1331 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom ).release();
1332 }
1333 else
1334 {
1335 unionGeom.clear();
1336 }
1337
1338 return 0;
1339}
1340
1341geos::unique_ptr QgsGeos::createGeosCollection( int typeId, const QVector<GEOSGeometry *> &geoms )
1342{
1343 int nNullGeoms = geoms.count( nullptr );
1344 int nNotNullGeoms = geoms.size() - nNullGeoms;
1345
1346 GEOSGeometry **geomarr = new GEOSGeometry*[ nNotNullGeoms ];
1347 if ( !geomarr )
1348 {
1349 return nullptr;
1350 }
1351
1352 int i = 0;
1353 QVector<GEOSGeometry *>::const_iterator geomIt = geoms.constBegin();
1354 for ( ; geomIt != geoms.constEnd(); ++geomIt )
1355 {
1356 if ( *geomIt )
1357 {
1358 if ( GEOSisEmpty_r( geosinit()->ctxt, *geomIt ) )
1359 {
1360 // don't add empty parts to a geos collection, it can cause crashes in GEOS
1361 nNullGeoms++;
1362 nNotNullGeoms--;
1363 GEOSGeom_destroy_r( geosinit()->ctxt, *geomIt );
1364 }
1365 else
1366 {
1367 geomarr[i] = *geomIt;
1368 ++i;
1369 }
1370 }
1371 }
1372 geos::unique_ptr geom;
1373
1374 try
1375 {
1376 geom.reset( GEOSGeom_createCollection_r( geosinit()->ctxt, typeId, geomarr, nNotNullGeoms ) );
1377 }
1378 catch ( GEOSException & )
1379 {
1380 }
1381
1382 delete [] geomarr;
1383
1384 return geom;
1385}
1386
1387std::unique_ptr<QgsAbstractGeometry> QgsGeos::fromGeos( const GEOSGeometry *geos )
1388{
1389 if ( !geos )
1390 {
1391 return nullptr;
1392 }
1393
1394 int nCoordDims = GEOSGeom_getCoordinateDimension_r( geosinit()->ctxt, geos );
1395 int nDims = GEOSGeom_getDimensions_r( geosinit()->ctxt, geos );
1396 bool hasZ = ( nCoordDims == 3 );
1397 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1398
1399 switch ( GEOSGeomTypeId_r( geosinit()->ctxt, geos ) )
1400 {
1401 case GEOS_POINT: // a point
1402 {
1403 if ( GEOSisEmpty_r( geosinit()->ctxt, geos ) )
1404 return nullptr;
1405
1406 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, geos );
1407 unsigned int nPoints = 0;
1408 GEOSCoordSeq_getSize_r( geosinit()->ctxt, cs, &nPoints );
1409 return nPoints > 0 ? std::unique_ptr<QgsAbstractGeometry>( coordSeqPoint( cs, 0, hasZ, hasM ).clone() ) : nullptr;
1410 }
1411 case GEOS_LINESTRING:
1412 {
1413 return sequenceToLinestring( geos, hasZ, hasM );
1414 }
1415 case GEOS_POLYGON:
1416 {
1417 return fromGeosPolygon( geos );
1418 }
1419 case GEOS_MULTIPOINT:
1420 {
1421 std::unique_ptr< QgsMultiPoint > multiPoint( new QgsMultiPoint() );
1422 int nParts = GEOSGetNumGeometries_r( geosinit()->ctxt, geos );
1423 multiPoint->reserve( nParts );
1424 for ( int i = 0; i < nParts; ++i )
1425 {
1426 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, GEOSGetGeometryN_r( geosinit()->ctxt, geos, i ) );
1427 if ( cs )
1428 {
1429 unsigned int nPoints = 0;
1430 GEOSCoordSeq_getSize_r( geosinit()->ctxt, cs, &nPoints );
1431 if ( nPoints > 0 )
1432 multiPoint->addGeometry( coordSeqPoint( cs, 0, hasZ, hasM ).clone() );
1433 }
1434 }
1435 return std::move( multiPoint );
1436 }
1437 case GEOS_MULTILINESTRING:
1438 {
1439 std::unique_ptr< QgsMultiLineString > multiLineString( new QgsMultiLineString() );
1440 int nParts = GEOSGetNumGeometries_r( geosinit()->ctxt, geos );
1441 multiLineString->reserve( nParts );
1442 for ( int i = 0; i < nParts; ++i )
1443 {
1444 std::unique_ptr< QgsLineString >line( sequenceToLinestring( GEOSGetGeometryN_r( geosinit()->ctxt, geos, i ), hasZ, hasM ) );
1445 if ( line )
1446 {
1447 multiLineString->addGeometry( line.release() );
1448 }
1449 }
1450 return std::move( multiLineString );
1451 }
1452 case GEOS_MULTIPOLYGON:
1453 {
1454 std::unique_ptr< QgsMultiPolygon > multiPolygon( new QgsMultiPolygon() );
1455
1456 int nParts = GEOSGetNumGeometries_r( geosinit()->ctxt, geos );
1457 multiPolygon->reserve( nParts );
1458 for ( int i = 0; i < nParts; ++i )
1459 {
1460 std::unique_ptr< QgsPolygon > poly = fromGeosPolygon( GEOSGetGeometryN_r( geosinit()->ctxt, geos, i ) );
1461 if ( poly )
1462 {
1463 multiPolygon->addGeometry( poly.release() );
1464 }
1465 }
1466 return std::move( multiPolygon );
1467 }
1468 case GEOS_GEOMETRYCOLLECTION:
1469 {
1470 std::unique_ptr< QgsGeometryCollection > geomCollection( new QgsGeometryCollection() );
1471 int nParts = GEOSGetNumGeometries_r( geosinit()->ctxt, geos );
1472 geomCollection->reserve( nParts );
1473 for ( int i = 0; i < nParts; ++i )
1474 {
1475 std::unique_ptr< QgsAbstractGeometry > geom( fromGeos( GEOSGetGeometryN_r( geosinit()->ctxt, geos, i ) ) );
1476 if ( geom )
1477 {
1478 geomCollection->addGeometry( geom.release() );
1479 }
1480 }
1481 return std::move( geomCollection );
1482 }
1483 }
1484 return nullptr;
1485}
1486
1487std::unique_ptr<QgsPolygon> QgsGeos::fromGeosPolygon( const GEOSGeometry *geos )
1488{
1489 if ( GEOSGeomTypeId_r( geosinit()->ctxt, geos ) != GEOS_POLYGON )
1490 {
1491 return nullptr;
1492 }
1493
1494 int nCoordDims = GEOSGeom_getCoordinateDimension_r( geosinit()->ctxt, geos );
1495 int nDims = GEOSGeom_getDimensions_r( geosinit()->ctxt, geos );
1496 bool hasZ = ( nCoordDims == 3 );
1497 bool hasM = ( ( nDims - nCoordDims ) == 1 );
1498
1499 std::unique_ptr< QgsPolygon > polygon( new QgsPolygon() );
1500
1501 const GEOSGeometry *ring = GEOSGetExteriorRing_r( geosinit()->ctxt, geos );
1502 if ( ring )
1503 {
1504 polygon->setExteriorRing( sequenceToLinestring( ring, hasZ, hasM ).release() );
1505 }
1506
1507 QVector<QgsCurve *> interiorRings;
1508 const int ringCount = GEOSGetNumInteriorRings_r( geosinit()->ctxt, geos );
1509 interiorRings.reserve( ringCount );
1510 for ( int i = 0; i < ringCount; ++i )
1511 {
1512 ring = GEOSGetInteriorRingN_r( geosinit()->ctxt, geos, i );
1513 if ( ring )
1514 {
1515 interiorRings.push_back( sequenceToLinestring( ring, hasZ, hasM ).release() );
1516 }
1517 }
1518 polygon->setInteriorRings( interiorRings );
1519
1520 return polygon;
1521}
1522
1523std::unique_ptr<QgsLineString> QgsGeos::sequenceToLinestring( const GEOSGeometry *geos, bool hasZ, bool hasM )
1524{
1525 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, geos );
1526
1527 unsigned int nPoints;
1528 GEOSCoordSeq_getSize_r( geosinit()->ctxt, cs, &nPoints );
1529
1530 QVector< double > xOut( nPoints );
1531 QVector< double > yOut( nPoints );
1532 QVector< double > zOut;
1533 if ( hasZ )
1534 zOut.resize( nPoints );
1535 QVector< double > mOut;
1536 if ( hasM )
1537 mOut.resize( nPoints );
1538
1539 double *x = xOut.data();
1540 double *y = yOut.data();
1541 double *z = zOut.data();
1542 double *m = mOut.data();
1543
1544#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
1545 GEOSCoordSeq_copyToArrays_r( geosinit()->ctxt, cs, x, y, hasZ ? z : nullptr, hasM ? m : nullptr );
1546#else
1547 for ( unsigned int i = 0; i < nPoints; ++i )
1548 {
1549 if ( hasZ )
1550 GEOSCoordSeq_getXYZ_r( geosinit()->ctxt, cs, i, x++, y++, z++ );
1551 else
1552 GEOSCoordSeq_getXY_r( geosinit()->ctxt, cs, i, x++, y++ );
1553 if ( hasM )
1554 {
1555 GEOSCoordSeq_getOrdinate_r( geosinit()->ctxt, cs, i, 3, m++ );
1556 }
1557 }
1558#endif
1559 std::unique_ptr< QgsLineString > line( new QgsLineString( xOut, yOut, zOut, mOut ) );
1560 return line;
1561}
1562
1563int QgsGeos::numberOfGeometries( GEOSGeometry *g )
1564{
1565 if ( !g )
1566 return 0;
1567
1568 int geometryType = GEOSGeomTypeId_r( geosinit()->ctxt, g );
1569 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
1570 || geometryType == GEOS_POLYGON )
1571 return 1;
1572
1573 //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
1574 return GEOSGetNumGeometries_r( geosinit()->ctxt, g );
1575}
1576
1577QgsPoint QgsGeos::coordSeqPoint( const GEOSCoordSequence *cs, int i, bool hasZ, bool hasM )
1578{
1579 if ( !cs )
1580 {
1581 return QgsPoint();
1582 }
1583
1584 double x, y;
1585 double z = 0;
1586 double m = 0;
1587 if ( hasZ )
1588 GEOSCoordSeq_getXYZ_r( geosinit()->ctxt, cs, i, &x, &y, &z );
1589 else
1590 GEOSCoordSeq_getXY_r( geosinit()->ctxt, cs, i, &x, &y );
1591 if ( hasM )
1592 {
1593 GEOSCoordSeq_getOrdinate_r( geosinit()->ctxt, cs, i, 3, &m );
1594 }
1595
1597 if ( hasZ && hasM )
1598 {
1600 }
1601 else if ( hasZ )
1602 {
1604 }
1605 else if ( hasM )
1606 {
1608 }
1609 return QgsPoint( t, x, y, z, m );
1610}
1611
1613{
1614 if ( !geom )
1615 return nullptr;
1616
1617 int coordDims = 2;
1618 if ( geom->is3D() )
1619 {
1620 ++coordDims;
1621 }
1622 if ( geom->isMeasure() )
1623 {
1624 ++coordDims;
1625 }
1626
1628 {
1629 int geosType = GEOS_GEOMETRYCOLLECTION;
1630
1632 {
1633 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1634 {
1636 geosType = GEOS_MULTIPOINT;
1637 break;
1638
1640 geosType = GEOS_MULTILINESTRING;
1641 break;
1642
1644 geosType = GEOS_MULTIPOLYGON;
1645 break;
1646
1649 return nullptr;
1650 }
1651 }
1652
1653
1654 const QgsGeometryCollection *c = qgsgeometry_cast<const QgsGeometryCollection *>( geom );
1655
1656 if ( !c )
1657 return nullptr;
1658
1659 QVector< GEOSGeometry * > geomVector( c->numGeometries() );
1660 for ( int i = 0; i < c->numGeometries(); ++i )
1661 {
1662 geomVector[i] = asGeos( c->geometryN( i ), precision ).release();
1663 }
1664 return createGeosCollection( geosType, geomVector );
1665 }
1666 else
1667 {
1668 switch ( QgsWkbTypes::geometryType( geom->wkbType() ) )
1669 {
1671 return createGeosPoint( static_cast<const QgsPoint *>( geom ), coordDims, precision );
1672
1674 return createGeosLinestring( static_cast<const QgsLineString *>( geom ), precision );
1675
1677 return createGeosPolygon( static_cast<const QgsPolygon *>( geom ), precision );
1678
1681 return nullptr;
1682 }
1683 }
1684 return nullptr;
1685}
1686
1687std::unique_ptr<QgsAbstractGeometry> QgsGeos::overlay( const QgsAbstractGeometry *geom, Overlay op, QString *errorMsg, const QgsGeometryParameters &parameters ) const
1688{
1689 if ( !mGeos || !geom )
1690 {
1691 return nullptr;
1692 }
1693
1694 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
1695 if ( !geosGeom )
1696 {
1697 return nullptr;
1698 }
1699
1700 const double gridSize = parameters.gridSize();
1701
1702 try
1703 {
1704 geos::unique_ptr opGeom;
1705 switch ( op )
1706 {
1707 case OverlayIntersection:
1708 if ( gridSize > 0 )
1709 {
1710 opGeom.reset( GEOSIntersectionPrec_r( geosinit()->ctxt, mGeos.get(), geosGeom.get(), gridSize ) );
1711 }
1712 else
1713 {
1714 opGeom.reset( GEOSIntersection_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) );
1715 }
1716 break;
1717
1718 case OverlayDifference:
1719 if ( gridSize > 0 )
1720 {
1721 opGeom.reset( GEOSDifferencePrec_r( geosinit()->ctxt, mGeos.get(), geosGeom.get(), gridSize ) );
1722 }
1723 else
1724 {
1725 opGeom.reset( GEOSDifference_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) );
1726 }
1727 break;
1728
1729 case OverlayUnion:
1730 {
1731 geos::unique_ptr unionGeometry;
1732 if ( gridSize > 0 )
1733 {
1734 unionGeometry.reset( GEOSUnionPrec_r( geosinit()->ctxt, mGeos.get(), geosGeom.get(), gridSize ) );
1735 }
1736 else
1737 {
1738 unionGeometry.reset( GEOSUnion_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) );
1739 }
1740
1741 if ( unionGeometry && GEOSGeomTypeId_r( geosinit()->ctxt, unionGeometry.get() ) == GEOS_MULTILINESTRING )
1742 {
1743 geos::unique_ptr mergedLines( GEOSLineMerge_r( geosinit()->ctxt, unionGeometry.get() ) );
1744 if ( mergedLines )
1745 {
1746 unionGeometry = std::move( mergedLines );
1747 }
1748 }
1749
1750 opGeom = std::move( unionGeometry );
1751 }
1752 break;
1753
1754 case OverlaySymDifference:
1755 if ( gridSize > 0 )
1756 {
1757 opGeom.reset( GEOSSymDifferencePrec_r( geosinit()->ctxt, mGeos.get(), geosGeom.get(), gridSize ) );
1758 }
1759 else
1760 {
1761 opGeom.reset( GEOSSymDifference_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) );
1762 }
1763 break;
1764 }
1765 return fromGeos( opGeom.get() );
1766 }
1767 catch ( GEOSException &e )
1768 {
1769 logError( QStringLiteral( "GEOS" ), e.what() );
1770 if ( errorMsg )
1771 {
1772 *errorMsg = e.what();
1773 }
1774 return nullptr;
1775 }
1776}
1777
1778bool QgsGeos::relation( const QgsAbstractGeometry *geom, Relation r, QString *errorMsg ) const
1779{
1780 if ( !mGeos || !geom )
1781 {
1782 return false;
1783 }
1784
1785 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
1786 if ( !geosGeom )
1787 {
1788 return false;
1789 }
1790
1791 bool result = false;
1792 try
1793 {
1794 if ( mGeosPrepared ) //use faster version with prepared geometry
1795 {
1796 switch ( r )
1797 {
1798 case RelationIntersects:
1799 result = ( GEOSPreparedIntersects_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1800 break;
1801 case RelationTouches:
1802 result = ( GEOSPreparedTouches_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1803 break;
1804 case RelationCrosses:
1805 result = ( GEOSPreparedCrosses_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1806 break;
1807 case RelationWithin:
1808 result = ( GEOSPreparedWithin_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1809 break;
1810 case RelationContains:
1811 result = ( GEOSPreparedContains_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1812 break;
1813 case RelationDisjoint:
1814 result = ( GEOSPreparedDisjoint_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1815 break;
1816 case RelationOverlaps:
1817 result = ( GEOSPreparedOverlaps_r( geosinit()->ctxt, mGeosPrepared.get(), geosGeom.get() ) == 1 );
1818 break;
1819 }
1820 return result;
1821 }
1822
1823 switch ( r )
1824 {
1825 case RelationIntersects:
1826 result = ( GEOSIntersects_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1827 break;
1828 case RelationTouches:
1829 result = ( GEOSTouches_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1830 break;
1831 case RelationCrosses:
1832 result = ( GEOSCrosses_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1833 break;
1834 case RelationWithin:
1835 result = ( GEOSWithin_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1836 break;
1837 case RelationContains:
1838 result = ( GEOSContains_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1839 break;
1840 case RelationDisjoint:
1841 result = ( GEOSDisjoint_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1842 break;
1843 case RelationOverlaps:
1844 result = ( GEOSOverlaps_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() ) == 1 );
1845 break;
1846 }
1847 }
1848 catch ( GEOSException &e )
1849 {
1850 logError( QStringLiteral( "GEOS" ), e.what() );
1851 if ( errorMsg )
1852 {
1853 *errorMsg = e.what();
1854 }
1855 return false;
1856 }
1857
1858 return result;
1859}
1860
1861QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, QString *errorMsg ) const
1862{
1863 if ( !mGeos )
1864 {
1865 return nullptr;
1866 }
1867
1869 try
1870 {
1871 geos.reset( GEOSBuffer_r( geosinit()->ctxt, mGeos.get(), distance, segments ) );
1872 }
1873 CATCH_GEOS_WITH_ERRMSG( nullptr )
1874 return fromGeos( geos.get() ).release();
1875}
1876
1877QgsAbstractGeometry *QgsGeos::buffer( double distance, int segments, Qgis::EndCapStyle endCapStyle, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
1878{
1879 if ( !mGeos )
1880 {
1881 return nullptr;
1882 }
1883
1885 try
1886 {
1887 geos.reset( GEOSBufferWithStyle_r( geosinit()->ctxt, mGeos.get(), distance, segments, static_cast< int >( endCapStyle ), static_cast< int >( joinStyle ), miterLimit ) );
1888 }
1889 CATCH_GEOS_WITH_ERRMSG( nullptr )
1890 return fromGeos( geos.get() ).release();
1891}
1892
1893QgsAbstractGeometry *QgsGeos::simplify( double tolerance, QString *errorMsg ) const
1894{
1895 if ( !mGeos )
1896 {
1897 return nullptr;
1898 }
1900 try
1901 {
1902 geos.reset( GEOSTopologyPreserveSimplify_r( geosinit()->ctxt, mGeos.get(), tolerance ) );
1903 }
1904 CATCH_GEOS_WITH_ERRMSG( nullptr )
1905 return fromGeos( geos.get() ).release();
1906}
1907
1908QgsAbstractGeometry *QgsGeos::interpolate( double distance, QString *errorMsg ) const
1909{
1910 if ( !mGeos )
1911 {
1912 return nullptr;
1913 }
1915 try
1916 {
1917 geos.reset( GEOSInterpolate_r( geosinit()->ctxt, mGeos.get(), distance ) );
1918 }
1919 CATCH_GEOS_WITH_ERRMSG( nullptr )
1920 return fromGeos( geos.get() ).release();
1921}
1922
1923QgsPoint *QgsGeos::centroid( QString *errorMsg ) const
1924{
1925 if ( !mGeos )
1926 {
1927 return nullptr;
1928 }
1929
1931 double x;
1932 double y;
1933
1934 try
1935 {
1936 geos.reset( GEOSGetCentroid_r( geosinit()->ctxt, mGeos.get() ) );
1937
1938 if ( !geos )
1939 return nullptr;
1940
1941 GEOSGeomGetX_r( geosinit()->ctxt, geos.get(), &x );
1942 GEOSGeomGetY_r( geosinit()->ctxt, geos.get(), &y );
1943 }
1944 CATCH_GEOS_WITH_ERRMSG( nullptr )
1945
1946 return new QgsPoint( x, y );
1947}
1948
1949QgsAbstractGeometry *QgsGeos::envelope( QString *errorMsg ) const
1950{
1951 if ( !mGeos )
1952 {
1953 return nullptr;
1954 }
1956 try
1957 {
1958 geos.reset( GEOSEnvelope_r( geosinit()->ctxt, mGeos.get() ) );
1959 }
1960 CATCH_GEOS_WITH_ERRMSG( nullptr )
1961 return fromGeos( geos.get() ).release();
1962}
1963
1964QgsPoint *QgsGeos::pointOnSurface( QString *errorMsg ) const
1965{
1966 if ( !mGeos )
1967 {
1968 return nullptr;
1969 }
1970
1971 double x;
1972 double y;
1973
1975 try
1976 {
1977 geos.reset( GEOSPointOnSurface_r( geosinit()->ctxt, mGeos.get() ) );
1978
1979 if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
1980 {
1981 return nullptr;
1982 }
1983
1984 GEOSGeomGetX_r( geosinit()->ctxt, geos.get(), &x );
1985 GEOSGeomGetY_r( geosinit()->ctxt, geos.get(), &y );
1986 }
1987 CATCH_GEOS_WITH_ERRMSG( nullptr )
1988
1989 return new QgsPoint( x, y );
1990}
1991
1992QgsAbstractGeometry *QgsGeos::convexHull( QString *errorMsg ) const
1993{
1994 if ( !mGeos )
1995 {
1996 return nullptr;
1997 }
1998
1999 try
2000 {
2001 geos::unique_ptr cHull( GEOSConvexHull_r( geosinit()->ctxt, mGeos.get() ) );
2002 std::unique_ptr< QgsAbstractGeometry > cHullGeom = fromGeos( cHull.get() );
2003 return cHullGeom.release();
2004 }
2005 CATCH_GEOS_WITH_ERRMSG( nullptr )
2006}
2007
2008QgsAbstractGeometry *QgsGeos::concaveHull( double targetPercent, bool allowHoles, QString *errorMsg ) const
2009{
2010#if GEOS_VERSION_MAJOR==3 && GEOS_VERSION_MINOR<11
2011 ( void )allowHoles;
2012 ( void )targetPercent;
2013 ( void )errorMsg;
2014 throw QgsNotSupportedException( QObject::tr( "Calculating concaveHull requires a QGIS build based on GEOS 3.11 or later" ) );
2015#else
2016 if ( !mGeos )
2017 {
2018 return nullptr;
2019 }
2020
2021 try
2022 {
2023 geos::unique_ptr concaveHull( GEOSConcaveHull_r( geosinit()->ctxt, mGeos.get(), targetPercent, allowHoles ) );
2024 std::unique_ptr< QgsAbstractGeometry > concaveHullGeom = fromGeos( concaveHull.get() );
2025 return concaveHullGeom.release();
2026 }
2027 CATCH_GEOS_WITH_ERRMSG( nullptr )
2028#endif
2029}
2030
2031bool QgsGeos::isValid( QString *errorMsg, const bool allowSelfTouchingHoles, QgsGeometry *errorLoc ) const
2032{
2033 if ( !mGeos )
2034 {
2035 return false;
2036 }
2037
2038 try
2039 {
2040 GEOSGeometry *g1 = nullptr;
2041 char *r = nullptr;
2042 char res = GEOSisValidDetail_r( geosinit()->ctxt, mGeos.get(), allowSelfTouchingHoles ? GEOSVALID_ALLOW_SELFTOUCHING_RING_FORMING_HOLE : 0, &r, &g1 );
2043 const bool invalid = res != 1;
2044
2045 QString error;
2046 if ( r )
2047 {
2048 error = QString( r );
2049 GEOSFree_r( geosinit()->ctxt, r );
2050 }
2051
2052 if ( invalid && errorMsg )
2053 {
2054 static QgsStringMap translatedErrors;
2055
2056 if ( translatedErrors.empty() )
2057 {
2058 // Copied from https://git.osgeo.org/gitea/geos/geos/src/branch/master/src/operation/valid/TopologyValidationError.cpp
2059 translatedErrors.insert( QStringLiteral( "topology validation error" ), QObject::tr( "Topology validation error", "GEOS Error" ) );
2060 translatedErrors.insert( QStringLiteral( "repeated point" ), QObject::tr( "Repeated point", "GEOS Error" ) );
2061 translatedErrors.insert( QStringLiteral( "hole lies outside shell" ), QObject::tr( "Hole lies outside shell", "GEOS Error" ) );
2062 translatedErrors.insert( QStringLiteral( "holes are nested" ), QObject::tr( "Holes are nested", "GEOS Error" ) );
2063 translatedErrors.insert( QStringLiteral( "interior is disconnected" ), QObject::tr( "Interior is disconnected", "GEOS Error" ) );
2064 translatedErrors.insert( QStringLiteral( "self-intersection" ), QObject::tr( "Self-intersection", "GEOS Error" ) );
2065 translatedErrors.insert( QStringLiteral( "ring self-intersection" ), QObject::tr( "Ring self-intersection", "GEOS Error" ) );
2066 translatedErrors.insert( QStringLiteral( "nested shells" ), QObject::tr( "Nested shells", "GEOS Error" ) );
2067 translatedErrors.insert( QStringLiteral( "duplicate rings" ), QObject::tr( "Duplicate rings", "GEOS Error" ) );
2068 translatedErrors.insert( QStringLiteral( "too few points in geometry component" ), QObject::tr( "Too few points in geometry component", "GEOS Error" ) );
2069 translatedErrors.insert( QStringLiteral( "invalid coordinate" ), QObject::tr( "Invalid coordinate", "GEOS Error" ) );
2070 translatedErrors.insert( QStringLiteral( "ring is not closed" ), QObject::tr( "Ring is not closed", "GEOS Error" ) );
2071 }
2072
2073 *errorMsg = translatedErrors.value( error.toLower(), error );
2074
2075 if ( g1 && errorLoc )
2076 {
2077 *errorLoc = geometryFromGeos( g1 );
2078 }
2079 else if ( g1 )
2080 {
2081 GEOSGeom_destroy_r( geosinit()->ctxt, g1 );
2082 }
2083 }
2084 return !invalid;
2085 }
2086 CATCH_GEOS_WITH_ERRMSG( false )
2087}
2088
2089bool QgsGeos::isEqual( const QgsAbstractGeometry *geom, QString *errorMsg ) const
2090{
2091 if ( !mGeos || !geom )
2092 {
2093 return false;
2094 }
2095
2096 try
2097 {
2098 geos::unique_ptr geosGeom( asGeos( geom, mPrecision ) );
2099 if ( !geosGeom )
2100 {
2101 return false;
2102 }
2103 bool equal = GEOSEquals_r( geosinit()->ctxt, mGeos.get(), geosGeom.get() );
2104 return equal;
2105 }
2106 CATCH_GEOS_WITH_ERRMSG( false )
2107}
2108
2109bool QgsGeos::isEmpty( QString *errorMsg ) const
2110{
2111 if ( !mGeos )
2112 {
2113 return false;
2114 }
2115
2116 try
2117 {
2118 return GEOSisEmpty_r( geosinit()->ctxt, mGeos.get() );
2119 }
2120 CATCH_GEOS_WITH_ERRMSG( false )
2121}
2122
2123bool QgsGeos::isSimple( QString *errorMsg ) const
2124{
2125 if ( !mGeos )
2126 {
2127 return false;
2128 }
2129
2130 try
2131 {
2132 return GEOSisSimple_r( geosinit()->ctxt, mGeos.get() );
2133 }
2134 CATCH_GEOS_WITH_ERRMSG( false )
2135}
2136
2137GEOSCoordSequence *QgsGeos::createCoordinateSequence( const QgsCurve *curve, double precision, bool forceClose )
2138{
2139 GEOSContextHandle_t ctxt = geosinit()->ctxt;
2140
2141 std::unique_ptr< QgsLineString > segmentized;
2142 const QgsLineString *line = qgsgeometry_cast<const QgsLineString *>( curve );
2143
2144 if ( !line )
2145 {
2146 segmentized.reset( curve->curveToLine() );
2147 line = segmentized.get();
2148 }
2149
2150 if ( !line )
2151 {
2152 return nullptr;
2153 }
2154 GEOSCoordSequence *coordSeq = nullptr;
2155
2156 const int numPoints = line->numPoints();
2157
2158 const bool hasZ = line->is3D();
2159
2160#if GEOS_VERSION_MAJOR>3 || ( GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR>=10 )
2161 if ( qgsDoubleNear( precision, 0 ) )
2162 {
2163 if ( !forceClose || ( line->pointN( 0 ) == line->pointN( numPoints - 1 ) ) )
2164 {
2165 // use optimised method if we don't have to force close an open ring
2166 try
2167 {
2168 coordSeq = GEOSCoordSeq_copyFromArrays_r( ctxt, line->xData(), line->yData(), line->zData(), nullptr, numPoints );
2169 if ( !coordSeq )
2170 {
2171 QgsDebugMsg( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points" ).arg( numPoints ) );
2172 return nullptr;
2173 }
2174 }
2175 CATCH_GEOS( nullptr )
2176 }
2177 else
2178 {
2179 QVector< double > x = line->xVector();
2180 if ( numPoints > 0 )
2181 x.append( x.at( 0 ) );
2182 QVector< double > y = line->yVector();
2183 if ( numPoints > 0 )
2184 y.append( y.at( 0 ) );
2185 QVector< double > z = line->zVector();
2186 if ( hasZ && numPoints > 0 )
2187 z.append( z.at( 0 ) );
2188 try
2189 {
2190 coordSeq = GEOSCoordSeq_copyFromArrays_r( ctxt, x.constData(), y.constData(), !hasZ ? nullptr : z.constData(), nullptr, numPoints + 1 );
2191 if ( !coordSeq )
2192 {
2193 QgsDebugMsg( QStringLiteral( "GEOS Exception: Could not create closed coordinate sequence for %1 points" ).arg( numPoints + 1 ) );
2194 return nullptr;
2195 }
2196 }
2197 CATCH_GEOS( nullptr )
2198 }
2199 return coordSeq;
2200 }
2201#endif
2202
2203 int coordDims = 2;
2204 const bool hasM = false; //line->isMeasure(); //disabled until geos supports m-coordinates
2205
2206 if ( hasZ )
2207 {
2208 ++coordDims;
2209 }
2210 if ( hasM )
2211 {
2212 ++coordDims;
2213 }
2214
2215 int numOutPoints = numPoints;
2216 if ( forceClose && ( line->pointN( 0 ) != line->pointN( numPoints - 1 ) ) )
2217 {
2218 ++numOutPoints;
2219 }
2220
2221 try
2222 {
2223 coordSeq = GEOSCoordSeq_create_r( ctxt, numOutPoints, coordDims );
2224 if ( !coordSeq )
2225 {
2226 QgsDebugMsg( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for %1 points in %2 dimensions" ).arg( numPoints ).arg( coordDims ) );
2227 return nullptr;
2228 }
2229
2230 const double *xData = line->xData();
2231 const double *yData = line->yData();
2232 const double *zData = hasZ ? line->zData() : nullptr;
2233 const double *mData = hasM ? line->mData() : nullptr;
2234
2235 if ( precision > 0. )
2236 {
2237 for ( int i = 0; i < numOutPoints; ++i )
2238 {
2239 if ( i >= numPoints )
2240 {
2241 // start reading back from start of line
2242 xData = line->xData();
2243 yData = line->yData();
2244 zData = hasZ ? line->zData() : nullptr;
2245 mData = hasM ? line->mData() : nullptr;
2246 }
2247 if ( hasZ )
2248 {
2249 GEOSCoordSeq_setXYZ_r( ctxt, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision, std::round( *zData++ / precision ) * precision );
2250 }
2251 else
2252 {
2253 GEOSCoordSeq_setXY_r( ctxt, coordSeq, i, std::round( *xData++ / precision ) * precision, std::round( *yData++ / precision ) * precision );
2254 }
2255 if ( hasM )
2256 {
2257 GEOSCoordSeq_setOrdinate_r( ctxt, coordSeq, i, 3, *mData++ );
2258 }
2259 }
2260 }
2261 else
2262 {
2263 for ( int i = 0; i < numOutPoints; ++i )
2264 {
2265 if ( i >= numPoints )
2266 {
2267 // start reading back from start of line
2268 xData = line->xData();
2269 yData = line->yData();
2270 zData = hasZ ? line->zData() : nullptr;
2271 mData = hasM ? line->mData() : nullptr;
2272 }
2273 if ( hasZ )
2274 {
2275 GEOSCoordSeq_setXYZ_r( ctxt, coordSeq, i, *xData++, *yData++, *zData++ );
2276 }
2277 else
2278 {
2279 GEOSCoordSeq_setXY_r( ctxt, coordSeq, i, *xData++, *yData++ );
2280 }
2281 if ( hasM )
2282 {
2283 GEOSCoordSeq_setOrdinate_r( ctxt, coordSeq, i, 3, *mData++ );
2284 }
2285 }
2286 }
2287 }
2288 CATCH_GEOS( nullptr )
2289
2290 return coordSeq;
2291}
2292
2293geos::unique_ptr QgsGeos::createGeosPoint( const QgsAbstractGeometry *point, int coordDims, double precision )
2294{
2295 const QgsPoint *pt = qgsgeometry_cast<const QgsPoint *>( point );
2296 if ( !pt )
2297 return nullptr;
2298
2299 return createGeosPointXY( pt->x(), pt->y(), pt->is3D(), pt->z(), pt->isMeasure(), pt->m(), coordDims, precision );
2300}
2301
2302geos::unique_ptr QgsGeos::createGeosPointXY( double x, double y, bool hasZ, double z, bool hasM, double m, int coordDims, double precision )
2303{
2304 Q_UNUSED( hasM )
2305 Q_UNUSED( m )
2306
2307 geos::unique_ptr geosPoint;
2308 try
2309 {
2310 if ( coordDims == 2 )
2311 {
2312 // optimised constructor
2313 if ( precision > 0. )
2314 geosPoint.reset( GEOSGeom_createPointFromXY_r( geosinit()->ctxt, std::round( x / precision ) * precision, std::round( y / precision ) * precision ) );
2315 else
2316 geosPoint.reset( GEOSGeom_createPointFromXY_r( geosinit()->ctxt, x, y ) );
2317 return geosPoint;
2318 }
2319
2320 GEOSCoordSequence *coordSeq = GEOSCoordSeq_create_r( geosinit()->ctxt, 1, coordDims );
2321 if ( !coordSeq )
2322 {
2323 QgsDebugMsg( QStringLiteral( "GEOS Exception: Could not create coordinate sequence for point with %1 dimensions" ).arg( coordDims ) );
2324 return nullptr;
2325 }
2326 if ( precision > 0. )
2327 {
2328 GEOSCoordSeq_setX_r( geosinit()->ctxt, coordSeq, 0, std::round( x / precision ) * precision );
2329 GEOSCoordSeq_setY_r( geosinit()->ctxt, coordSeq, 0, std::round( y / precision ) * precision );
2330 if ( hasZ )
2331 {
2332 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 2, std::round( z / precision ) * precision );
2333 }
2334 }
2335 else
2336 {
2337 GEOSCoordSeq_setX_r( geosinit()->ctxt, coordSeq, 0, x );
2338 GEOSCoordSeq_setY_r( geosinit()->ctxt, coordSeq, 0, y );
2339 if ( hasZ )
2340 {
2341 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 2, z );
2342 }
2343 }
2344#if 0 //disabled until geos supports m-coordinates
2345 if ( hasM )
2346 {
2347 GEOSCoordSeq_setOrdinate_r( geosinit()->ctxt, coordSeq, 0, 3, m );
2348 }
2349#endif
2350 geosPoint.reset( GEOSGeom_createPoint_r( geosinit()->ctxt, coordSeq ) );
2351 }
2352 CATCH_GEOS( nullptr )
2353 return geosPoint;
2354}
2355
2356geos::unique_ptr QgsGeos::createGeosLinestring( const QgsAbstractGeometry *curve, double precision )
2357{
2358 const QgsCurve *c = qgsgeometry_cast<const QgsCurve *>( curve );
2359 if ( !c )
2360 return nullptr;
2361
2362 GEOSCoordSequence *coordSeq = createCoordinateSequence( c, precision );
2363 if ( !coordSeq )
2364 return nullptr;
2365
2366 geos::unique_ptr geosGeom;
2367 try
2368 {
2369 geosGeom.reset( GEOSGeom_createLineString_r( geosinit()->ctxt, coordSeq ) );
2370 }
2371 CATCH_GEOS( nullptr )
2372 return geosGeom;
2373}
2374
2375geos::unique_ptr QgsGeos::createGeosPolygon( const QgsAbstractGeometry *poly, double precision )
2376{
2377 const QgsCurvePolygon *polygon = qgsgeometry_cast<const QgsCurvePolygon *>( poly );
2378 if ( !polygon )
2379 return nullptr;
2380
2381 const QgsCurve *exteriorRing = polygon->exteriorRing();
2382 if ( !exteriorRing )
2383 {
2384 return nullptr;
2385 }
2386
2387 geos::unique_ptr geosPolygon;
2388 try
2389 {
2390 geos::unique_ptr exteriorRingGeos( GEOSGeom_createLinearRing_r( geosinit()->ctxt, createCoordinateSequence( exteriorRing, precision, true ) ) );
2391
2392 int nHoles = polygon->numInteriorRings();
2393 GEOSGeometry **holes = nullptr;
2394 if ( nHoles > 0 )
2395 {
2396 holes = new GEOSGeometry*[ nHoles ];
2397 }
2398
2399 for ( int i = 0; i < nHoles; ++i )
2400 {
2401 const QgsCurve *interiorRing = polygon->interiorRing( i );
2402 holes[i] = GEOSGeom_createLinearRing_r( geosinit()->ctxt, createCoordinateSequence( interiorRing, precision, true ) );
2403 }
2404 geosPolygon.reset( GEOSGeom_createPolygon_r( geosinit()->ctxt, exteriorRingGeos.release(), holes, nHoles ) );
2405 delete[] holes;
2406 }
2407 CATCH_GEOS( nullptr )
2408
2409 return geosPolygon;
2410}
2411
2412QgsAbstractGeometry *QgsGeos::offsetCurve( double distance, int segments, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2413{
2414 if ( !mGeos )
2415 return nullptr;
2416
2417 geos::unique_ptr offset;
2418 try
2419 {
2420 offset.reset( GEOSOffsetCurve_r( geosinit()->ctxt, mGeos.get(), distance, segments, static_cast< int >( joinStyle ), miterLimit ) );
2421 }
2422 CATCH_GEOS_WITH_ERRMSG( nullptr )
2423 std::unique_ptr< QgsAbstractGeometry > offsetGeom = fromGeos( offset.get() );
2424 return offsetGeom.release();
2425}
2426
2427std::unique_ptr<QgsAbstractGeometry> QgsGeos::singleSidedBuffer( double distance, int segments, Qgis::BufferSide side, Qgis::JoinStyle joinStyle, double miterLimit, QString *errorMsg ) const
2428{
2429 if ( !mGeos )
2430 {
2431 return nullptr;
2432 }
2433
2435 try
2436 {
2437 geos::buffer_params_unique_ptr bp( GEOSBufferParams_create_r( geosinit()->ctxt ) );
2438 GEOSBufferParams_setSingleSided_r( geosinit()->ctxt, bp.get(), 1 );
2439 GEOSBufferParams_setQuadrantSegments_r( geosinit()->ctxt, bp.get(), segments );
2440 GEOSBufferParams_setJoinStyle_r( geosinit()->ctxt, bp.get(), static_cast< int >( joinStyle ) );
2441 GEOSBufferParams_setMitreLimit_r( geosinit()->ctxt, bp.get(), miterLimit ); //#spellok
2442
2443 if ( side == Qgis::BufferSide::Right )
2444 {
2445 distance = -distance;
2446 }
2447 geos.reset( GEOSBufferWithParams_r( geosinit()->ctxt, mGeos.get(), bp.get(), distance ) );
2448 }
2449 CATCH_GEOS_WITH_ERRMSG( nullptr )
2450 return fromGeos( geos.get() );
2451}
2452
2453std::unique_ptr<QgsAbstractGeometry> QgsGeos::maximumInscribedCircle( double tolerance, QString *errorMsg ) const
2454{
2455 if ( !mGeos )
2456 {
2457 return nullptr;
2458 }
2459
2461 try
2462 {
2463 geos.reset( GEOSMaximumInscribedCircle_r( geosinit()->ctxt, mGeos.get(), tolerance ) );
2464 }
2465 CATCH_GEOS_WITH_ERRMSG( nullptr )
2466 return fromGeos( geos.get() );
2467}
2468
2469std::unique_ptr<QgsAbstractGeometry> QgsGeos::largestEmptyCircle( double tolerance, const QgsAbstractGeometry *boundary, QString *errorMsg ) const
2470{
2471 if ( !mGeos )
2472 {
2473 return nullptr;
2474 }
2475
2477 try
2478 {
2479 geos::unique_ptr boundaryGeos;
2480 if ( boundary )
2481 boundaryGeos = asGeos( boundary );
2482
2483 geos.reset( GEOSLargestEmptyCircle_r( geosinit()->ctxt, mGeos.get(), boundaryGeos.get(), tolerance ) );
2484 }
2485 CATCH_GEOS_WITH_ERRMSG( nullptr )
2486 return fromGeos( geos.get() );
2487}
2488
2489std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumWidth( QString *errorMsg ) const
2490{
2491 if ( !mGeos )
2492 {
2493 return nullptr;
2494 }
2495
2497 try
2498 {
2499 geos.reset( GEOSMinimumWidth_r( geosinit()->ctxt, mGeos.get() ) );
2500 }
2501 CATCH_GEOS_WITH_ERRMSG( nullptr )
2502 return fromGeos( geos.get() );
2503}
2504
2505double QgsGeos::minimumClearance( QString *errorMsg ) const
2506{
2507 if ( !mGeos )
2508 {
2509 return std::numeric_limits< double >::quiet_NaN();
2510 }
2511
2513 double res = 0;
2514 try
2515 {
2516 if ( GEOSMinimumClearance_r( geosinit()->ctxt, mGeos.get(), &res ) != 0 )
2517 return std::numeric_limits< double >::quiet_NaN();
2518 }
2519 CATCH_GEOS_WITH_ERRMSG( std::numeric_limits< double >::quiet_NaN() )
2520 return res;
2521}
2522
2523std::unique_ptr<QgsAbstractGeometry> QgsGeos::minimumClearanceLine( QString *errorMsg ) const
2524{
2525 if ( !mGeos )
2526 {
2527 return nullptr;
2528 }
2529
2531 try
2532 {
2533 geos.reset( GEOSMinimumClearanceLine_r( geosinit()->ctxt, mGeos.get() ) );
2534 }
2535 CATCH_GEOS_WITH_ERRMSG( nullptr )
2536 return fromGeos( geos.get() );
2537}
2538
2539std::unique_ptr<QgsAbstractGeometry> QgsGeos::node( QString *errorMsg ) const
2540{
2541 if ( !mGeos )
2542 {
2543 return nullptr;
2544 }
2545
2547 try
2548 {
2549 geos.reset( GEOSNode_r( geosinit()->ctxt, mGeos.get() ) );
2550 }
2551 CATCH_GEOS_WITH_ERRMSG( nullptr )
2552 return fromGeos( geos.get() );
2553}
2554
2555std::unique_ptr<QgsAbstractGeometry> QgsGeos::sharedPaths( const QgsAbstractGeometry *other, QString *errorMsg ) const
2556{
2557 if ( !mGeos || !other )
2558 {
2559 return nullptr;
2560 }
2561
2563 try
2564 {
2565 geos::unique_ptr otherGeos = asGeos( other );
2566 if ( !otherGeos )
2567 return nullptr;
2568
2569 geos.reset( GEOSSharedPaths_r( geosinit()->ctxt, mGeos.get(), otherGeos.get() ) );
2570 }
2571 CATCH_GEOS_WITH_ERRMSG( nullptr )
2572 return fromGeos( geos.get() );
2573}
2574
2575std::unique_ptr<QgsAbstractGeometry> QgsGeos::reshapeGeometry( const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg ) const
2576{
2577 if ( !mGeos || mGeometry->dimension() == 0 )
2578 {
2579 if ( errorCode ) { *errorCode = InvalidBaseGeometry; }
2580 return nullptr;
2581 }
2582
2583 if ( reshapeWithLine.numPoints() < 2 )
2584 {
2585 if ( errorCode ) { *errorCode = InvalidInput; }
2586 return nullptr;
2587 }
2588
2589 geos::unique_ptr reshapeLineGeos = createGeosLinestring( &reshapeWithLine, mPrecision );
2590
2591 //single or multi?
2592 int numGeoms = GEOSGetNumGeometries_r( geosinit()->ctxt, mGeos.get() );
2593 if ( numGeoms == -1 )
2594 {
2595 if ( errorCode )
2596 {
2597 *errorCode = InvalidBaseGeometry;
2598 }
2599 return nullptr;
2600 }
2601
2602 bool isMultiGeom = false;
2603 int geosTypeId = GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() );
2604 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
2605 isMultiGeom = true;
2606
2607 bool isLine = ( mGeometry->dimension() == 1 );
2608
2609 if ( !isMultiGeom )
2610 {
2611 geos::unique_ptr reshapedGeometry;
2612 if ( isLine )
2613 {
2614 reshapedGeometry = reshapeLine( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2615 }
2616 else
2617 {
2618 reshapedGeometry = reshapePolygon( mGeos.get(), reshapeLineGeos.get(), mPrecision );
2619 }
2620
2621 if ( errorCode )
2622 *errorCode = Success;
2623 std::unique_ptr< QgsAbstractGeometry > reshapeResult = fromGeos( reshapedGeometry.get() );
2624 return reshapeResult;
2625 }
2626 else
2627 {
2628 try
2629 {
2630 //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
2631 bool reshapeTookPlace = false;
2632
2633 geos::unique_ptr currentReshapeGeometry;
2634 GEOSGeometry **newGeoms = new GEOSGeometry*[numGeoms];
2635
2636 for ( int i = 0; i < numGeoms; ++i )
2637 {
2638 if ( isLine )
2639 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2640 else
2641 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ), reshapeLineGeos.get(), mPrecision );
2642
2643 if ( currentReshapeGeometry )
2644 {
2645 newGeoms[i] = currentReshapeGeometry.release();
2646 reshapeTookPlace = true;
2647 }
2648 else
2649 {
2650 newGeoms[i] = GEOSGeom_clone_r( geosinit()->ctxt, GEOSGetGeometryN_r( geosinit()->ctxt, mGeos.get(), i ) );
2651 }
2652 }
2653
2654 geos::unique_ptr newMultiGeom;
2655 if ( isLine )
2656 {
2657 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, newGeoms, numGeoms ) );
2658 }
2659 else //multipolygon
2660 {
2661 newMultiGeom.reset( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTIPOLYGON, newGeoms, numGeoms ) );
2662 }
2663
2664 delete[] newGeoms;
2665 if ( !newMultiGeom )
2666 {
2667 if ( errorCode ) { *errorCode = EngineError; }
2668 return nullptr;
2669 }
2670
2671 if ( reshapeTookPlace )
2672 {
2673 if ( errorCode )
2674 *errorCode = Success;
2675 std::unique_ptr< QgsAbstractGeometry > reshapedMultiGeom = fromGeos( newMultiGeom.get() );
2676 return reshapedMultiGeom;
2677 }
2678 else
2679 {
2680 if ( errorCode )
2681 {
2682 *errorCode = NothingHappened;
2683 }
2684 return nullptr;
2685 }
2686 }
2687 CATCH_GEOS_WITH_ERRMSG( nullptr )
2688 }
2689}
2690
2691QgsGeometry QgsGeos::mergeLines( QString *errorMsg ) const
2692{
2693 if ( !mGeos )
2694 {
2695 return QgsGeometry();
2696 }
2697
2698 if ( GEOSGeomTypeId_r( geosinit()->ctxt, mGeos.get() ) != GEOS_MULTILINESTRING )
2699 return QgsGeometry();
2700
2702 try
2703 {
2704 geos.reset( GEOSLineMerge_r( geosinit()->ctxt, mGeos.get() ) );
2705 }
2707 return QgsGeometry( fromGeos( geos.get() ) );
2708}
2709
2710QgsGeometry QgsGeos::closestPoint( const QgsGeometry &other, QString *errorMsg ) const
2711{
2712 if ( !mGeos || isEmpty() || other.isEmpty() )
2713 {
2714 return QgsGeometry();
2715 }
2716
2717 geos::unique_ptr otherGeom( asGeos( other.constGet(), mPrecision ) );
2718 if ( !otherGeom )
2719 {
2720 return QgsGeometry();
2721 }
2722
2723 double nx = 0.0;
2724 double ny = 0.0;
2725 try
2726 {
2728 if ( mGeosPrepared ) // use faster version with prepared geometry
2729 {
2730 nearestCoord.reset( GEOSPreparedNearestPoints_r( geosinit()->ctxt, mGeosPrepared.get(), otherGeom.get() ) );
2731 }
2732 else
2733 {
2734 nearestCoord.reset( GEOSNearestPoints_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() ) );
2735 }
2736
2737 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 0, &nx );
2738 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 0, &ny );
2739 }
2740 catch ( GEOSException &e )
2741 {
2742 logError( QStringLiteral( "GEOS" ), e.what() );
2743 if ( errorMsg )
2744 {
2745 *errorMsg = e.what();
2746 }
2747 return QgsGeometry();
2748 }
2749
2750 return QgsGeometry( new QgsPoint( nx, ny ) );
2751}
2752
2753QgsGeometry QgsGeos::shortestLine( const QgsGeometry &other, QString *errorMsg ) const
2754{
2755 if ( !mGeos || other.isEmpty() )
2756 {
2757 return QgsGeometry();
2758 }
2759
2760 return shortestLine( other.constGet(), errorMsg );
2761}
2762
2763QgsGeometry QgsGeos::shortestLine( const QgsAbstractGeometry *other, QString *errorMsg ) const
2764{
2765 if ( !other || other->isEmpty() )
2766 return QgsGeometry();
2767
2768 geos::unique_ptr otherGeom( asGeos( other, mPrecision ) );
2769 if ( !otherGeom )
2770 {
2771 return QgsGeometry();
2772 }
2773
2774 double nx1 = 0.0;
2775 double ny1 = 0.0;
2776 double nx2 = 0.0;
2777 double ny2 = 0.0;
2778 try
2779 {
2780 geos::coord_sequence_unique_ptr nearestCoord( GEOSNearestPoints_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() ) );
2781
2782 if ( !nearestCoord )
2783 {
2784 if ( errorMsg )
2785 *errorMsg = QStringLiteral( "GEOS returned no nearest points" );
2786 return QgsGeometry();
2787 }
2788
2789 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 0, &nx1 );
2790 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 0, &ny1 );
2791 ( void )GEOSCoordSeq_getX_r( geosinit()->ctxt, nearestCoord.get(), 1, &nx2 );
2792 ( void )GEOSCoordSeq_getY_r( geosinit()->ctxt, nearestCoord.get(), 1, &ny2 );
2793 }
2794 catch ( GEOSException &e )
2795 {
2796 logError( QStringLiteral( "GEOS" ), e.what() );
2797 if ( errorMsg )
2798 {
2799 *errorMsg = e.what();
2800 }
2801 return QgsGeometry();
2802 }
2803
2804 QgsLineString *line = new QgsLineString();
2805 line->addVertex( QgsPoint( nx1, ny1 ) );
2806 line->addVertex( QgsPoint( nx2, ny2 ) );
2807 return QgsGeometry( line );
2808}
2809
2810double QgsGeos::lineLocatePoint( const QgsPoint &point, QString *errorMsg ) const
2811{
2812 if ( !mGeos )
2813 {
2814 return -1;
2815 }
2816
2817 geos::unique_ptr otherGeom( asGeos( &point, mPrecision ) );
2818 if ( !otherGeom )
2819 {
2820 return -1;
2821 }
2822
2823 double distance = -1;
2824 try
2825 {
2826 distance = GEOSProject_r( geosinit()->ctxt, mGeos.get(), otherGeom.get() );
2827 }
2828 catch ( GEOSException &e )
2829 {
2830 logError( QStringLiteral( "GEOS" ), e.what() );
2831 if ( errorMsg )
2832 {
2833 *errorMsg = e.what();
2834 }
2835 return -1;
2836 }
2837
2838 return distance;
2839}
2840
2841double QgsGeos::lineLocatePoint( double x, double y, QString *errorMsg ) const
2842{
2843 if ( !mGeos )
2844 {
2845 return -1;
2846 }
2847
2848 geos::unique_ptr point = createGeosPointXY( x, y, false, 0, false, 0, 2, 0 );
2849 if ( !point )
2850 return false;
2851
2852 double distance = -1;
2853 try
2854 {
2855 distance = GEOSProject_r( geosinit()->ctxt, mGeos.get(), point.get() );
2856 }
2857 catch ( GEOSException &e )
2858 {
2859 logError( QStringLiteral( "GEOS" ), e.what() );
2860 if ( errorMsg )
2861 {
2862 *errorMsg = e.what();
2863 }
2864 return -1;
2865 }
2866
2867 return distance;
2868}
2869
2870QgsGeometry QgsGeos::polygonize( const QVector<const QgsAbstractGeometry *> &geometries, QString *errorMsg )
2871{
2872 GEOSGeometry **const lineGeosGeometries = new GEOSGeometry*[ geometries.size()];
2873 int validLines = 0;
2874 for ( const QgsAbstractGeometry *g : geometries )
2875 {
2876 geos::unique_ptr l = asGeos( g );
2877 if ( l )
2878 {
2879 lineGeosGeometries[validLines] = l.release();
2880 validLines++;
2881 }
2882 }
2883
2884 try
2885 {
2886 geos::unique_ptr result( GEOSPolygonize_r( geosinit()->ctxt, lineGeosGeometries, validLines ) );
2887 for ( int i = 0; i < validLines; ++i )
2888 {
2889 GEOSGeom_destroy_r( geosinit()->ctxt, lineGeosGeometries[i] );
2890 }
2891 delete[] lineGeosGeometries;
2892 return QgsGeometry( fromGeos( result.get() ) );
2893 }
2894 catch ( GEOSException &e )
2895 {
2896 if ( errorMsg )
2897 {
2898 *errorMsg = e.what();
2899 }
2900 for ( int i = 0; i < validLines; ++i )
2901 {
2902 GEOSGeom_destroy_r( geosinit()->ctxt, lineGeosGeometries[i] );
2903 }
2904 delete[] lineGeosGeometries;
2905 return QgsGeometry();
2906 }
2907}
2908
2909QgsGeometry QgsGeos::voronoiDiagram( const QgsAbstractGeometry *extent, double tolerance, bool edgesOnly, QString *errorMsg ) const
2910{
2911 if ( !mGeos )
2912 {
2913 return QgsGeometry();
2914 }
2915
2916 geos::unique_ptr extentGeosGeom;
2917 if ( extent )
2918 {
2919 extentGeosGeom = asGeos( extent, mPrecision );
2920 if ( !extentGeosGeom )
2921 {
2922 return QgsGeometry();
2923 }
2924 }
2925
2927 try
2928 {
2929 geos.reset( GEOSVoronoiDiagram_r( geosinit()->ctxt, mGeos.get(), extentGeosGeom.get(), tolerance, edgesOnly ) );
2930
2931 if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
2932 {
2933 return QgsGeometry();
2934 }
2935
2936 return QgsGeometry( fromGeos( geos.get() ) );
2937 }
2939}
2940
2941QgsGeometry QgsGeos::delaunayTriangulation( double tolerance, bool edgesOnly, QString *errorMsg ) const
2942{
2943 if ( !mGeos )
2944 {
2945 return QgsGeometry();
2946 }
2947
2949 try
2950 {
2951 geos.reset( GEOSDelaunayTriangulation_r( geosinit()->ctxt, mGeos.get(), tolerance, edgesOnly ) );
2952
2953 if ( !geos || GEOSisEmpty_r( geosinit()->ctxt, geos.get() ) != 0 )
2954 {
2955 return QgsGeometry();
2956 }
2957
2958 return QgsGeometry( fromGeos( geos.get() ) );
2959 }
2961}
2962
2963
2965static bool _linestringEndpoints( const GEOSGeometry *linestring, double &x1, double &y1, double &x2, double &y2 )
2966{
2967 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, linestring );
2968 if ( !coordSeq )
2969 return false;
2970
2971 unsigned int coordSeqSize;
2972 if ( GEOSCoordSeq_getSize_r( geosinit()->ctxt, coordSeq, &coordSeqSize ) == 0 )
2973 return false;
2974
2975 if ( coordSeqSize < 2 )
2976 return false;
2977
2978 GEOSCoordSeq_getX_r( geosinit()->ctxt, coordSeq, 0, &x1 );
2979 GEOSCoordSeq_getY_r( geosinit()->ctxt, coordSeq, 0, &y1 );
2980 GEOSCoordSeq_getX_r( geosinit()->ctxt, coordSeq, coordSeqSize - 1, &x2 );
2981 GEOSCoordSeq_getY_r( geosinit()->ctxt, coordSeq, coordSeqSize - 1, &y2 );
2982 return true;
2983}
2984
2985
2987static geos::unique_ptr _mergeLinestrings( const GEOSGeometry *line1, const GEOSGeometry *line2, const QgsPointXY &intersectionPoint )
2988{
2989 double x1, y1, x2, y2;
2990 if ( !_linestringEndpoints( line1, x1, y1, x2, y2 ) )
2991 return nullptr;
2992
2993 double rx1, ry1, rx2, ry2;
2994 if ( !_linestringEndpoints( line2, rx1, ry1, rx2, ry2 ) )
2995 return nullptr;
2996
2997 bool intersectionAtOrigLineEndpoint =
2998 ( intersectionPoint.x() == x1 && intersectionPoint.y() == y1 ) !=
2999 ( intersectionPoint.x() == x2 && intersectionPoint.y() == y2 );
3000 bool intersectionAtReshapeLineEndpoint =
3001 ( intersectionPoint.x() == rx1 && intersectionPoint.y() == ry1 ) ||
3002 ( intersectionPoint.x() == rx2 && intersectionPoint.y() == ry2 );
3003
3004 // the intersection must be at the begin/end of both lines
3005 if ( intersectionAtOrigLineEndpoint && intersectionAtReshapeLineEndpoint )
3006 {
3007 geos::unique_ptr g1( GEOSGeom_clone_r( geosinit()->ctxt, line1 ) );
3008 geos::unique_ptr g2( GEOSGeom_clone_r( geosinit()->ctxt, line2 ) );
3009 GEOSGeometry *geoms[2] = { g1.release(), g2.release() };
3010 geos::unique_ptr multiGeom( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, geoms, 2 ) );
3011 geos::unique_ptr res( GEOSLineMerge_r( geosinit()->ctxt, multiGeom.get() ) );
3012 return res;
3013 }
3014 else
3015 return nullptr;
3016}
3017
3018
3019geos::unique_ptr QgsGeos::reshapeLine( const GEOSGeometry *line, const GEOSGeometry *reshapeLineGeos, double precision )
3020{
3021 if ( !line || !reshapeLineGeos )
3022 return nullptr;
3023
3024 bool atLeastTwoIntersections = false;
3025 bool oneIntersection = false;
3026 QgsPointXY oneIntersectionPoint;
3027
3028 try
3029 {
3030 //make sure there are at least two intersection between line and reshape geometry
3031 geos::unique_ptr intersectGeom( GEOSIntersection_r( geosinit()->ctxt, line, reshapeLineGeos ) );
3032 if ( intersectGeom )
3033 {
3034 const int geomType = GEOSGeomTypeId_r( geosinit()->ctxt, intersectGeom.get() );
3035 atLeastTwoIntersections = ( geomType == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 1 )
3036 || ( geomType == GEOS_GEOMETRYCOLLECTION && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 0 ) // a collection implies at least two points!
3037 || ( geomType == GEOS_MULTILINESTRING && GEOSGetNumGeometries_r( geosinit()->ctxt, intersectGeom.get() ) > 0 );
3038 // one point is enough when extending line at its endpoint
3039 if ( GEOSGeomTypeId_r( geosinit()->ctxt, intersectGeom.get() ) == GEOS_POINT )
3040 {
3041 const GEOSCoordSequence *intersectionCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, intersectGeom.get() );
3042 double xi, yi;
3043 GEOSCoordSeq_getX_r( geosinit()->ctxt, intersectionCoordSeq, 0, &xi );
3044 GEOSCoordSeq_getY_r( geosinit()->ctxt, intersectionCoordSeq, 0, &yi );
3045 oneIntersection = true;
3046 oneIntersectionPoint = QgsPointXY( xi, yi );
3047 }
3048 }
3049 }
3050 catch ( GEOSException & )
3051 {
3052 atLeastTwoIntersections = false;
3053 }
3054
3055 // special case when extending line at its endpoint
3056 if ( oneIntersection )
3057 return _mergeLinestrings( line, reshapeLineGeos, oneIntersectionPoint );
3058
3059 if ( !atLeastTwoIntersections )
3060 return nullptr;
3061
3062 //begin and end point of original line
3063 double x1, y1, x2, y2;
3064 if ( !_linestringEndpoints( line, x1, y1, x2, y2 ) )
3065 return nullptr;
3066
3067 geos::unique_ptr beginLineVertex = createGeosPointXY( x1, y1, false, 0, false, 0, 2, precision );
3068 geos::unique_ptr endLineVertex = createGeosPointXY( x2, y2, false, 0, false, 0, 2, precision );
3069
3070 bool isRing = false;
3071 if ( GEOSGeomTypeId_r( geosinit()->ctxt, line ) == GEOS_LINEARRING
3072 || GEOSEquals_r( geosinit()->ctxt, beginLineVertex.get(), endLineVertex.get() ) == 1 )
3073 isRing = true;
3074
3075 //node line and reshape line
3076 geos::unique_ptr nodedGeometry = nodeGeometries( reshapeLineGeos, line );
3077 if ( !nodedGeometry )
3078 {
3079 return nullptr;
3080 }
3081
3082 //and merge them together
3083 geos::unique_ptr mergedLines( GEOSLineMerge_r( geosinit()->ctxt, nodedGeometry.get() ) );
3084 if ( !mergedLines )
3085 {
3086 return nullptr;
3087 }
3088
3089 int numMergedLines = GEOSGetNumGeometries_r( geosinit()->ctxt, mergedLines.get() );
3090 if ( numMergedLines < 2 ) //some special cases. Normally it is >2
3091 {
3092 if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
3093 {
3094 geos::unique_ptr result( GEOSGeom_clone_r( geosinit()->ctxt, reshapeLineGeos ) );
3095 return result;
3096 }
3097 else
3098 return nullptr;
3099 }
3100
3101 QVector<GEOSGeometry *> resultLineParts; //collection with the line segments that will be contained in result
3102 QVector<GEOSGeometry *> probableParts; //parts where we can decide on inclusion only after going through all the candidates
3103
3104 for ( int i = 0; i < numMergedLines; ++i )
3105 {
3106 const GEOSGeometry *currentGeom = GEOSGetGeometryN_r( geosinit()->ctxt, mergedLines.get(), i );
3107
3108 // have we already added this part?
3109 bool alreadyAdded = false;
3110 double distance = 0;
3111 double bufferDistance = std::pow( 10.0L, geomDigits( currentGeom ) - 11 );
3112 for ( const GEOSGeometry *other : std::as_const( resultLineParts ) )
3113 {
3114 GEOSHausdorffDistance_r( geosinit()->ctxt, currentGeom, other, &distance );
3115 if ( distance < bufferDistance )
3116 {
3117 alreadyAdded = true;
3118 break;
3119 }
3120 }
3121 if ( alreadyAdded )
3122 continue;
3123
3124 const GEOSCoordSequence *currentCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, currentGeom );
3125 unsigned int currentCoordSeqSize;
3126 GEOSCoordSeq_getSize_r( geosinit()->ctxt, currentCoordSeq, &currentCoordSeqSize );
3127 if ( currentCoordSeqSize < 2 )
3128 continue;
3129
3130 //get the two endpoints of the current line merge result
3131 double xBegin, xEnd, yBegin, yEnd;
3132 GEOSCoordSeq_getX_r( geosinit()->ctxt, currentCoordSeq, 0, &xBegin );
3133 GEOSCoordSeq_getY_r( geosinit()->ctxt, currentCoordSeq, 0, &yBegin );
3134 GEOSCoordSeq_getX_r( geosinit()->ctxt, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
3135 GEOSCoordSeq_getY_r( geosinit()->ctxt, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
3136 geos::unique_ptr beginCurrentGeomVertex = createGeosPointXY( xBegin, yBegin, false, 0, false, 0, 2, precision );
3137 geos::unique_ptr endCurrentGeomVertex = createGeosPointXY( xEnd, yEnd, false, 0, false, 0, 2, precision );
3138
3139 //check how many endpoints of the line merge result are on the (original) line
3140 int nEndpointsOnOriginalLine = 0;
3141 if ( pointContainedInLine( beginCurrentGeomVertex.get(), line ) == 1 )
3142 nEndpointsOnOriginalLine += 1;
3143
3144 if ( pointContainedInLine( endCurrentGeomVertex.get(), line ) == 1 )
3145 nEndpointsOnOriginalLine += 1;
3146
3147 //check how many endpoints equal the endpoints of the original line
3148 int nEndpointsSameAsOriginalLine = 0;
3149 if ( GEOSEquals_r( geosinit()->ctxt, beginCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3150 || GEOSEquals_r( geosinit()->ctxt, beginCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3151 nEndpointsSameAsOriginalLine += 1;
3152
3153 if ( GEOSEquals_r( geosinit()->ctxt, endCurrentGeomVertex.get(), beginLineVertex.get() ) == 1
3154 || GEOSEquals_r( geosinit()->ctxt, endCurrentGeomVertex.get(), endLineVertex.get() ) == 1 )
3155 nEndpointsSameAsOriginalLine += 1;
3156
3157 //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
3158 bool currentGeomOverlapsOriginalGeom = false;
3159 bool currentGeomOverlapsReshapeLine = false;
3160 if ( lineContainedInLine( currentGeom, line ) == 1 )
3161 currentGeomOverlapsOriginalGeom = true;
3162
3163 if ( lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
3164 currentGeomOverlapsReshapeLine = true;
3165
3166 //logic to decide if this part belongs to the result
3167 if ( !isRing && nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3168 {
3169 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3170 }
3171 //for closed rings, we take one segment from the candidate list
3172 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
3173 {
3174 probableParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3175 }
3176 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3177 {
3178 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3179 }
3180 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
3181 {
3182 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3183 }
3184 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
3185 {
3186 resultLineParts.push_back( GEOSGeom_clone_r( geosinit()->ctxt, currentGeom ) );
3187 }
3188 }
3189
3190 //add the longest segment from the probable list for rings (only used for polygon rings)
3191 if ( isRing && !probableParts.isEmpty() )
3192 {
3193 geos::unique_ptr maxGeom; //the longest geometry in the probabla list
3194 GEOSGeometry *currentGeom = nullptr;
3195 double maxLength = -std::numeric_limits<double>::max();
3196 double currentLength = 0;
3197 for ( int i = 0; i < probableParts.size(); ++i )
3198 {
3199 currentGeom = probableParts.at( i );
3200 GEOSLength_r( geosinit()->ctxt, currentGeom, &currentLength );
3201 if ( currentLength > maxLength )
3202 {
3203 maxLength = currentLength;
3204 maxGeom.reset( currentGeom );
3205 }
3206 else
3207 {
3208 GEOSGeom_destroy_r( geosinit()->ctxt, currentGeom );
3209 }
3210 }
3211 resultLineParts.push_back( maxGeom.release() );
3212 }
3213
3214 geos::unique_ptr result;
3215 if ( resultLineParts.empty() )
3216 return nullptr;
3217
3218 if ( resultLineParts.size() == 1 ) //the whole result was reshaped
3219 {
3220 result.reset( resultLineParts[0] );
3221 }
3222 else //>1
3223 {
3224 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
3225 for ( int i = 0; i < resultLineParts.size(); ++i )
3226 {
3227 lineArray[i] = resultLineParts[i];
3228 }
3229
3230 //create multiline from resultLineParts
3231 geos::unique_ptr multiLineGeom( GEOSGeom_createCollection_r( geosinit()->ctxt, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ) );
3232 delete [] lineArray;
3233
3234 //then do a linemerge with the newly combined partstrings
3235 result.reset( GEOSLineMerge_r( geosinit()->ctxt, multiLineGeom.get() ) );
3236 }
3237
3238 //now test if the result is a linestring. Otherwise something went wrong
3239 if ( GEOSGeomTypeId_r( geosinit()->ctxt, result.get() ) != GEOS_LINESTRING )
3240 {
3241 return nullptr;
3242 }
3243
3244 return result;
3245}
3246
3247geos::unique_ptr QgsGeos::reshapePolygon( const GEOSGeometry *polygon, const GEOSGeometry *reshapeLineGeos, double precision )
3248{
3249 //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
3250 int nIntersections = 0;
3251 int lastIntersectingRing = -2;
3252 const GEOSGeometry *lastIntersectingGeom = nullptr;
3253
3254 int nRings = GEOSGetNumInteriorRings_r( geosinit()->ctxt, polygon );
3255 if ( nRings < 0 )
3256 return nullptr;
3257
3258 //does outer ring intersect?
3259 const GEOSGeometry *outerRing = GEOSGetExteriorRing_r( geosinit()->ctxt, polygon );
3260 if ( GEOSIntersects_r( geosinit()->ctxt, outerRing, reshapeLineGeos ) == 1 )
3261 {
3262 ++nIntersections;
3263 lastIntersectingRing = -1;
3264 lastIntersectingGeom = outerRing;
3265 }
3266
3267 //do inner rings intersect?
3268 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
3269
3270 try
3271 {
3272 for ( int i = 0; i < nRings; ++i )
3273 {
3274 innerRings[i] = GEOSGetInteriorRingN_r( geosinit()->ctxt, polygon, i );
3275 if ( GEOSIntersects_r( geosinit()->ctxt, innerRings[i], reshapeLineGeos ) == 1 )
3276 {
3277 ++nIntersections;
3278 lastIntersectingRing = i;
3279 lastIntersectingGeom = innerRings[i];
3280 }
3281 }
3282 }
3283 catch ( GEOSException & )
3284 {
3285 nIntersections = 0;
3286 }
3287
3288 if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
3289 {
3290 delete [] innerRings;
3291 return nullptr;
3292 }
3293
3294 //we have one intersecting ring, let's try to reshape it
3295 geos::unique_ptr reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos, precision );
3296 if ( !reshapeResult )
3297 {
3298 delete [] innerRings;
3299 return nullptr;
3300 }
3301
3302 //if reshaping took place, we need to reassemble the polygon and its rings
3303 GEOSGeometry *newRing = nullptr;
3304 const GEOSCoordSequence *reshapeSequence = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, reshapeResult.get() );
3305 GEOSCoordSequence *newCoordSequence = GEOSCoordSeq_clone_r( geosinit()->ctxt, reshapeSequence );
3306
3307 reshapeResult.reset();
3308
3309 newRing = GEOSGeom_createLinearRing_r( geosinit()->ctxt, newCoordSequence );
3310 if ( !newRing )
3311 {
3312 delete [] innerRings;
3313 return nullptr;
3314 }
3315
3316 GEOSGeometry *newOuterRing = nullptr;
3317 if ( lastIntersectingRing == -1 )
3318 newOuterRing = newRing;
3319 else
3320 newOuterRing = GEOSGeom_clone_r( geosinit()->ctxt, outerRing );
3321
3322 //check if all the rings are still inside the outer boundary
3323 QVector<GEOSGeometry *> ringList;
3324 if ( nRings > 0 )
3325 {
3326 GEOSGeometry *outerRingPoly = GEOSGeom_createPolygon_r( geosinit()->ctxt, GEOSGeom_clone_r( geosinit()->ctxt, newOuterRing ), nullptr, 0 );
3327 if ( outerRingPoly )
3328 {
3329 ringList.reserve( nRings );
3330 GEOSGeometry *currentRing = nullptr;
3331 for ( int i = 0; i < nRings; ++i )
3332 {
3333 if ( lastIntersectingRing == i )
3334 currentRing = newRing;
3335 else
3336 currentRing = GEOSGeom_clone_r( geosinit()->ctxt, innerRings[i] );
3337
3338 //possibly a ring is no longer contained in the result polygon after reshape
3339 if ( GEOSContains_r( geosinit()->ctxt, outerRingPoly, currentRing ) == 1 )
3340 ringList.push_back( currentRing );
3341 else
3342 GEOSGeom_destroy_r( geosinit()->ctxt, currentRing );
3343 }
3344 }
3345 GEOSGeom_destroy_r( geosinit()->ctxt, outerRingPoly );
3346 }
3347
3348 GEOSGeometry **newInnerRings = new GEOSGeometry*[ringList.size()];
3349 for ( int i = 0; i < ringList.size(); ++i )
3350 newInnerRings[i] = ringList.at( i );
3351
3352 delete [] innerRings;
3353
3354 geos::unique_ptr reshapedPolygon( GEOSGeom_createPolygon_r( geosinit()->ctxt, newOuterRing, newInnerRings, ringList.size() ) );
3355 delete[] newInnerRings;
3356
3357 return reshapedPolygon;
3358}
3359
3360int QgsGeos::lineContainedInLine( const GEOSGeometry *line1, const GEOSGeometry *line2 )
3361{
3362 if ( !line1 || !line2 )
3363 {
3364 return -1;
3365 }
3366
3367 double bufferDistance = std::pow( 10.0L, geomDigits( line2 ) - 11 );
3368
3369 geos::unique_ptr bufferGeom( GEOSBuffer_r( geosinit()->ctxt, line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ) );
3370 if ( !bufferGeom )
3371 return -2;
3372
3373 geos::unique_ptr intersectionGeom( GEOSIntersection_r( geosinit()->ctxt, bufferGeom.get(), line1 ) );
3374
3375 //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
3376 double intersectGeomLength;
3377 double line1Length;
3378
3379 GEOSLength_r( geosinit()->ctxt, intersectionGeom.get(), &intersectGeomLength );
3380 GEOSLength_r( geosinit()->ctxt, line1, &line1Length );
3381
3382 double intersectRatio = line1Length / intersectGeomLength;
3383 if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
3384 return 1;
3385
3386 return 0;
3387}
3388
3389int QgsGeos::pointContainedInLine( const GEOSGeometry *point, const GEOSGeometry *line )
3390{
3391 if ( !point || !line )
3392 return -1;
3393
3394 double bufferDistance = std::pow( 10.0L, geomDigits( line ) - 11 );
3395
3396 geos::unique_ptr lineBuffer( GEOSBuffer_r( geosinit()->ctxt, line, bufferDistance, 8 ) );
3397 if ( !lineBuffer )
3398 return -2;
3399
3400 bool contained = false;
3401 if ( GEOSContains_r( geosinit()->ctxt, lineBuffer.get(), point ) == 1 )
3402 contained = true;
3403
3404 return contained;
3405}
3406
3407int QgsGeos::geomDigits( const GEOSGeometry *geom )
3408{
3409 geos::unique_ptr bbox( GEOSEnvelope_r( geosinit()->ctxt, geom ) );
3410 if ( !bbox.get() )
3411 return -1;
3412
3413 const GEOSGeometry *bBoxRing = GEOSGetExteriorRing_r( geosinit()->ctxt, bbox.get() );
3414 if ( !bBoxRing )
3415 return -1;
3416
3417 const GEOSCoordSequence *bBoxCoordSeq = GEOSGeom_getCoordSeq_r( geosinit()->ctxt, bBoxRing );
3418
3419 if ( !bBoxCoordSeq )
3420 return -1;
3421
3422 unsigned int nCoords = 0;
3423 if ( !GEOSCoordSeq_getSize_r( geosinit()->ctxt, bBoxCoordSeq, &nCoords ) )
3424 return -1;
3425
3426 int maxDigits = -1;
3427 for ( unsigned int i = 0; i < nCoords - 1; ++i )
3428 {
3429 double t;
3430 GEOSCoordSeq_getX_r( geosinit()->ctxt, bBoxCoordSeq, i, &t );
3431
3432 int digits;
3433 digits = std::ceil( std::log10( std::fabs( t ) ) );
3434 if ( digits > maxDigits )
3435 maxDigits = digits;
3436
3437 GEOSCoordSeq_getY_r( geosinit()->ctxt, bBoxCoordSeq, i, &t );
3438 digits = std::ceil( std::log10( std::fabs( t ) ) );
3439 if ( digits > maxDigits )
3440 maxDigits = digits;
3441 }
3442
3443 return maxDigits;
3444}
3445
3446GEOSContextHandle_t QgsGeos::getGEOSHandler()
3447{
3448 return geosinit()->ctxt;
3449}
The Qgis class provides global constants for use throughout the application.
Definition: qgis.h:72
BufferSide
Side of line to buffer.
Definition: qgis.h:1008
GeometryOperationResult
Success or failure of a geometry operation.
Definition: qgis.h:955
@ AddPartNotMultiGeometry
The source geometry is not multi.
@ InvalidBaseGeometry
The base geometry on which the operation is done is invalid or empty.
JoinStyle
Join styles for buffers.
Definition: qgis.h:1033
EndCapStyle
End cap styles for buffers.
Definition: qgis.h:1020
MakeValidMethod
Algorithms to use when repairing invalid geometries.
Definition: qgis.h:1046
@ 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...
Abstract base class for all geometries.
bool is3D() const SIP_HOLDGIL
Returns true if the geometry is 3D and contains a z-value.
virtual QgsRectangle boundingBox() const =0
Returns the minimal bounding box for the geometry.
virtual bool isEmpty() const
Returns true if the geometry is empty.
QgsWkbTypes::Type wkbType() const SIP_HOLDGIL
Returns the WKB type of the geometry.
virtual int dimension() const =0
Returns the inherent dimension of the geometry.
bool isMeasure() const SIP_HOLDGIL
Returns true if the geometry contains m values.
virtual QgsAbstractGeometry * clone() const =0
Clones the geometry by performing a deep copy.
Curve polygon geometry type.
const QgsCurve * interiorRing(int i) const SIP_HOLDGIL
Retrieves an interior ring from the curve polygon.
const QgsCurve * exteriorRing() const SIP_HOLDGIL
Returns the curve polygon's exterior ring.
int numInteriorRings() const SIP_HOLDGIL
Returns the number of interior rings contained with the curve polygon.
Abstract base class for curved geometry type.
Definition: qgscurve.h:36
virtual QgsLineString * curveToLine(double tolerance=M_PI_2/90, SegmentationToleranceType toleranceType=MaximumAngle) const =0
Returns a new line string geometry corresponding to a segmentized approximation of the curve.
Geometry collection.
virtual bool addGeometry(QgsAbstractGeometry *g)
Adds a geometry and takes ownership. Returns true in case of success.
static Qgis::GeometryOperationResult addPart(QgsAbstractGeometry *geometry, std::unique_ptr< QgsAbstractGeometry > part)
Add a part to multi type geometry.
A geometry engine is a low-level representation of a QgsAbstractGeometry object, optimised for use wi...
const QgsAbstractGeometry * mGeometry
EngineOperationResult
Success or failure of a geometry operation.
@ NothingHappened
Nothing happened, without any error.
@ InvalidBaseGeometry
The geometry on which the operation occurs is not valid.
@ InvalidInput
The input is not valid.
@ NodedGeometryError
Error occurred while creating a noded geometry.
@ EngineError
Error occurred in the geometry engine.
@ SplitCannotSplitPoint
Points cannot be split.
@ Success
Operation succeeded.
void logError(const QString &engineName, const QString &message) const
Logs an error message encountered during an operation.
static std::unique_ptr< QgsGeometryCollection > createCollectionOfType(QgsWkbTypes::Type type)
Returns a new geometry collection matching a specified WKB type.
Encapsulates parameters under which a geometry operation is performed.
Definition: qgsgeometry.h:111
double gridSize() const
Returns the grid size which will be used to snap vertices of a geometry.
Definition: qgsgeometry.h:124
A geometry is the spatial representation of a feature.
Definition: qgsgeometry.h:164
const QgsAbstractGeometry * constGet() const SIP_HOLDGIL
Returns a non-modifiable (const) reference to the underlying abstract geometry primitive.
Q_GADGET bool isNull
Definition: qgsgeometry.h:166
QgsAbstractGeometry * get()
Returns a modifiable (non-const) reference to the underlying abstract geometry primitive.
bool isEmpty() const
Returns true if the geometry is empty (eg a linestring with no vertices, or a collection with no geom...
Does vector analysis using the geos library and handles import, export, exception handling*.
Definition: qgsgeos.h: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:2941
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:2031
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:1861
QgsAbstractGeometry * concaveHull(double targetPercent, bool allowHoles=false, QString *errorMsg=nullptr) const
Returns a possibly concave geometry that encloses the input geometry.
Definition: qgsgeos.cpp:2008
std::unique_ptr< QgsAbstractGeometry > reshapeGeometry(const QgsLineString &reshapeWithLine, EngineOperationResult *errorCode, QString *errorMsg=nullptr) const
Reshapes the geometry using a line.
Definition: qgsgeos.cpp:2575
std::unique_ptr< QgsAbstractGeometry > maximumInscribedCircle(double tolerance, QString *errorMsg=nullptr) const
Returns the maximum inscribed circle.
Definition: qgsgeos.cpp:2453
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:2555
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:2539
double minimumClearance(QString *errorMsg=nullptr) const
Computes the minimum clearance of a geometry.
Definition: qgsgeos.cpp:2505
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:2427
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:2753
QgsAbstractGeometry * simplify(double tolerance, QString *errorMsg=nullptr) const override
Definition: qgsgeos.cpp:1893
QgsGeometry closestPoint(const QgsGeometry &other, QString *errorMsg=nullptr) const
Returns the closest point on the geometry to the other geometry.
Definition: qgsgeos.cpp:2710
static std::unique_ptr< QgsPolygon > fromGeosPolygon(const GEOSGeometry *geos)
Definition: qgsgeos.cpp:1487
std::unique_ptr< QgsAbstractGeometry > minimumClearanceLine(QString *errorMsg=nullptr) const
Returns a LineString whose endpoints define the minimum clearance of a geometry.
Definition: qgsgeos.cpp:2523
QgsAbstractGeometry * envelope(QString *errorMsg=nullptr) const override
Definition: qgsgeos.cpp:1949
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:1992
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:2489
bool isSimple(QString *errorMsg=nullptr) const override
Determines whether the geometry is simple (according to OGC definition).
Definition: qgsgeos.cpp:2123
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:1387
QgsAbstractGeometry * interpolate(double distance, QString *errorMsg=nullptr) const override
Definition: qgsgeos.cpp:1908
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:2089
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:2691
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:2870
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:2412
static GEOSContextHandle_t getGEOSHandler()
Definition: qgsgeos.cpp:3446
QgsPoint * centroid(QString *errorMsg=nullptr) const override
Calculates the centroid of this.
Definition: qgsgeos.cpp:1923
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:2810
bool isEmpty(QString *errorMsg=nullptr) const override
Definition: qgsgeos.cpp:2109
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:2909
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:2469
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:1577
QgsPoint * pointOnSurface(QString *errorMsg=nullptr) const override
Calculate a point that is guaranteed to be on the surface of this.
Definition: qgsgeos.cpp:1964
static QgsGeometry geometryFromGeos(GEOSGeometry *geos)
Creates a new QgsGeometry object, feeding in a geometry in GEOS format.
Definition: qgsgeos.cpp:150
Line string geometry type, with support for z-dimension and m-values.
Definition: qgslinestring.h:45
const double * yData() const
Returns a const pointer to the y vertex data.
const double * xData() const
Returns a const pointer to the x vertex data.
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 SIP_HOLDGIL
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.
Definition: qgsmulticurve.h:30
bool addGeometry(QgsAbstractGeometry *g) override
Adds a geometry and takes ownership. Returns true in case of success.
Multi line string geometry collection.
Multi point geometry collection.
Definition: qgsmultipoint.h:30
Multi polygon geometry collection.
Custom exception class which is raised when an operation is not supported.
Definition: qgsexception.h:118
A class to represent a 2D point.
Definition: qgspointxy.h:59
double y
Definition: qgspointxy.h:63
Q_GADGET double x
Definition: qgspointxy.h:62
Point geometry type, with support for z-dimension and m-values.
Definition: qgspoint.h:49
double distance(double x, double y) const SIP_HOLDGIL
Returns the Cartesian 2D distance between this point and a specified x, y coordinate.
Definition: qgspoint.h:343
Q_GADGET double x
Definition: qgspoint.h:52
double z
Definition: qgspoint.h:54
double m
Definition: qgspoint.h:55
double y
Definition: qgspoint.h:53
Polygon geometry type.
Definition: qgspolygon.h:34
A rectangle specified with double values.
Definition: qgsrectangle.h:42
double yMaximum() const SIP_HOLDGIL
Returns the y maximum value (top side of rectangle).
Definition: qgsrectangle.h:193
double xMaximum() const SIP_HOLDGIL
Returns the x maximum value (right side of rectangle).
Definition: qgsrectangle.h:183
double xMinimum() const SIP_HOLDGIL
Returns the x minimum value (left side of rectangle).
Definition: qgsrectangle.h:188
double yMinimum() const SIP_HOLDGIL
Returns the y minimum value (bottom side of rectangle).
Definition: qgsrectangle.h:198
void setYMinimum(double y) SIP_HOLDGIL
Set the minimum y value.
Definition: qgsrectangle.h:161
bool isNull() const
Test if the rectangle is null (all coordinates zero or after call to setMinimal()).
Definition: qgsrectangle.h:479
void setXMaximum(double x) SIP_HOLDGIL
Set the maximum x value.
Definition: qgsrectangle.h:156
void setXMinimum(double x) SIP_HOLDGIL
Set the minimum x value.
Definition: qgsrectangle.h:151
double height() const SIP_HOLDGIL
Returns the height of the rectangle.
Definition: qgsrectangle.h:230
void setYMaximum(double y) SIP_HOLDGIL
Set the maximum y value.
Definition: qgsrectangle.h:166
double width() const SIP_HOLDGIL
Returns the width of the rectangle.
Definition: qgsrectangle.h:223
bool isEmpty() const
Returns true if the rectangle is empty.
Definition: qgsrectangle.h:469
static GeometryType geometryType(Type type) SIP_HOLDGIL
Returns the geometry type for a WKB type, e.g., both MultiPolygon and CurvePolygon would have a Polyg...
Definition: qgswkbtypes.h:968
static bool isMultiType(Type type) SIP_HOLDGIL
Returns true if the WKB type is a multi type.
Definition: qgswkbtypes.h:862
Type
The WKB type describes the number of dimensions a geometry has.
Definition: qgswkbtypes.h:70
@ GeometryCollection
Definition: qgswkbtypes.h:79
static Type flatType(Type type) SIP_HOLDGIL
Returns the flat type for a WKB type.
Definition: qgswkbtypes.h:732
Contains geos related utilities and functions.
Definition: qgsgeos.h: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:2527
QMap< QString, QString > QgsStringMap
Definition: qgis.h:3022
QVector< QgsPoint > QgsPointSequence
Q_GLOBAL_STATIC(QReadWriteLock, sDefinitionCacheLock)
#define CATCH_GEOS(r)
Definition: qgsgeos.cpp:33
#define DEFAULT_QUADRANT_SEGMENTS
Definition: qgsgeos.cpp:31
#define CATCH_GEOS_WITH_ERRMSG(r)
Definition: qgsgeos.cpp:39
#define QgsDebugMsg(str)
Definition: qgslogger.h:38
QLineF segment(int index, QRectF rect, double radius)
int precision
Utility class for identifying a unique vertex within a geometry.
Definition: qgsvertexid.h:31
int vertex
Vertex number.
Definition: qgsvertexid.h:95
void CORE_EXPORT operator()(GEOSGeometry *geom) const
Destroys the GEOS geometry geom, using the static QGIS geos context.