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