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