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