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