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