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