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