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