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