QGIS API Documentation  2.9.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
qgsgeometry.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometry.cpp - Geometry (stored as Open Geospatial Consortium WKB)
3  -------------------------------------------------------------------
4 Date : 02 May 2005
5 Copyright : (C) 2005 by Brendan Morley
6 email : morb at ozemail dot com dot au
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 <limits>
17 #include <cstdarg>
18 #include <cstdio>
19 #include <cmath>
20 
21 #include "qgis.h"
22 #include "qgsgeometry.h"
23 #include "qgsapplication.h"
24 #include "qgslogger.h"
25 #include "qgsmessagelog.h"
26 #include "qgspoint.h"
27 #include "qgsrectangle.h"
28 
29 #include "qgsmaplayerregistry.h"
30 #include "qgsvectorlayer.h"
31 #include "qgsproject.h"
32 #include "qgsmessagelog.h"
33 #include "qgsgeometryvalidator.h"
34 
35 #include <QDebug>
36 
37 #ifndef Q_OS_WIN
38 #include <netinet/in.h>
39 #else
40 #include <winsock.h>
41 #endif
42 
43 #define DEFAULT_QUADRANT_SEGMENTS 8
44 
45 #define CATCH_GEOS(r) \
46  catch (GEOSException &e) \
47  { \
48  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr("GEOS") ); \
49  return r; \
50  }
51 
53 {
54  public:
55  GEOSException( QString theMsg )
56  {
57  if ( theMsg == "Unknown exception thrown" && lastMsg.isNull() )
58  {
59  msg = theMsg;
60  }
61  else
62  {
63  msg = theMsg;
64  lastMsg = msg;
65  }
66  }
67 
68  // copy constructor
70  {
71  *this = rhs;
72  }
73 
75  {
76  if ( lastMsg == msg )
77  lastMsg = QString::null;
78  }
79 
80  QString what()
81  {
82  return msg;
83  }
84 
85  private:
86  QString msg;
87  static QString lastMsg;
88 };
89 
90 QString GEOSException::lastMsg;
91 
92 static void throwGEOSException( const char *fmt, ... )
93 {
94  va_list ap;
95  char buffer[1024];
96 
97  va_start( ap, fmt );
98  vsnprintf( buffer, sizeof buffer, fmt, ap );
99  va_end( ap );
100 
101  qWarning() << QString( "GEOS exception: %1" ).arg( buffer );
102 
103  throw GEOSException( QString::fromUtf8( buffer ) );
104 }
105 
106 static void printGEOSNotice( const char *fmt, ... )
107 {
108 #if defined(QGISDEBUG)
109  va_list ap;
110  char buffer[1024];
111 
112  va_start( ap, fmt );
113  vsnprintf( buffer, sizeof buffer, fmt, ap );
114  va_end( ap );
115 
116  QgsDebugMsg( QString( "GEOS notice: %1" ).arg( QString::fromUtf8( buffer ) ) );
117 #else
118  Q_UNUSED( fmt );
119 #endif
120 }
121 
122 class GEOSInit
123 {
124  public:
125  GEOSContextHandle_t ctxt;
126 
128  {
129  ctxt = initGEOS_r( printGEOSNotice, throwGEOSException );
130  }
131 
133  {
134  finishGEOS_r( ctxt );
135  }
136 };
137 
139 
140 GEOSContextHandle_t QgsGeometry::getGEOSHandler()
141 {
142  return geosinit.ctxt;
143 }
144 
146  : mGeometry( 0 )
147  , mGeometrySize( 0 )
148  , mGeos( 0 )
149  , mDirtyWkb( false )
150  , mDirtyGeos( false )
151 {
152 }
153 
155  : mGeometry( 0 )
156  , mGeometrySize( rhs.mGeometrySize )
157  , mDirtyWkb( rhs.mDirtyWkb )
158  , mDirtyGeos( rhs.mDirtyGeos )
159 {
160  if ( mGeometrySize && rhs.mGeometry )
161  {
162  mGeometry = new unsigned char[mGeometrySize];
163  memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
164  }
165 
166  // deep-copy the GEOS Geometry if appropriate
167  if ( rhs.mGeos )
168  mGeos = GEOSGeom_clone_r( geosinit.ctxt, rhs.mGeos );
169  else
170  mGeos = 0;
171 
172 }
173 
176 {
177  if ( mGeometry )
178  delete [] mGeometry;
179 
180  if ( mGeos )
181  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
182 
183 }
184 
185 static unsigned int getNumGeosPoints( const GEOSGeometry *geom )
186 {
187  unsigned int n;
188  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, geom );
189  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &n );
190  return n;
191 }
192 
193 static GEOSGeometry *createGeosPoint( const double x, const double y )
194 {
195  GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosinit.ctxt, 1, 2 );
196  GEOSCoordSeq_setX_r( geosinit.ctxt, coord, 0, x );
197  GEOSCoordSeq_setY_r( geosinit.ctxt, coord, 0, y );
198  return GEOSGeom_createPoint_r( geosinit.ctxt, coord );
199 }
200 
201 static GEOSGeometry *createGeosPoint( const QgsPoint &point )
202 {
203  return createGeosPoint( point.x(), point.y() );
204 }
205 
206 static GEOSCoordSequence *createGeosCoordSequence( const QgsPolyline& points )
207 {
208  GEOSCoordSequence *coord = 0;
209 
210  try
211  {
212  coord = GEOSCoordSeq_create_r( geosinit.ctxt, points.count(), 2 );
213  int i;
214  for ( i = 0; i < points.count(); i++ )
215  {
216  GEOSCoordSeq_setX_r( geosinit.ctxt, coord, i, points[i].x() );
217  GEOSCoordSeq_setY_r( geosinit.ctxt, coord, i, points[i].y() );
218  }
219  return coord;
220  }
221  catch ( GEOSException &e )
222  {
223  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
224  /*if ( coord )
225  GEOSCoordSeq_destroy( coord );*/
226  throw;
227  }
228 }
229 
230 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms )
231 {
232  GEOSGeometry **geomarr = new GEOSGeometry*[ geoms.size()];
233  if ( !geomarr )
234  return 0;
235 
236  for ( int i = 0; i < geoms.size(); i++ )
237  geomarr[i] = geoms[i];
238 
239  GEOSGeometry *geom = 0;
240 
241  try
242  {
243  geom = GEOSGeom_createCollection_r( geosinit.ctxt, typeId, geomarr, geoms.size() );
244  }
245  catch ( GEOSException &e )
246  {
247  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
248  }
249 
250  delete [] geomarr;
251 
252  return geom;
253 }
254 
255 static GEOSGeometry *createGeosLineString( const QgsPolyline& polyline )
256 {
257  GEOSCoordSequence *coord = 0;
258 
259  try
260  {
261  coord = createGeosCoordSequence( polyline );
262  return GEOSGeom_createLineString_r( geosinit.ctxt, coord );
263  }
264  catch ( GEOSException &e )
265  {
266  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
267  //MH: for strange reasons, geos3 crashes when removing the coordinate sequence
268  //if ( coord )
269  //GEOSCoordSeq_destroy( coord );
270  return 0;
271  }
272 }
273 
274 static GEOSGeometry *createGeosLinearRing( const QgsPolyline& polyline )
275 {
276  GEOSCoordSequence *coord = 0;
277 
278  if ( polyline.count() <= 2 )
279  return 0;
280 
281  try
282  {
283  if ( polyline[0] != polyline[polyline.size()-1] )
284  {
285  // Ring not closed
286  QgsPolyline closed( polyline );
287  closed << closed[0];
288  coord = createGeosCoordSequence( closed );
289  }
290  else
291  {
292  // XXX [MD] this exception should not be silenced!
293  // this is here just because maptopixel simplification can return invalid linear rings
294  if ( polyline.count() == 3 ) //-> Avoid 'GEOS::IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4'
295  return 0;
296 
297  coord = createGeosCoordSequence( polyline );
298  }
299 
300  return GEOSGeom_createLinearRing_r( geosinit.ctxt, coord );
301  }
302  catch ( GEOSException &e )
303  {
304  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
305  /* as MH has noticed ^, this crashes geos
306  if ( coord )
307  GEOSCoordSeq_destroy( coord );*/
308  return 0;
309  }
310 }
311 
312 static GEOSGeometry *createGeosPolygon( const QVector<GEOSGeometry*> &rings )
313 {
314  GEOSGeometry *shell;
315 
316  if ( rings.size() == 0 )
317  {
318 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
319  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)))
320  return GEOSGeom_createEmptyPolygon_r( geosinit.ctxt );
321 #else
322  shell = GEOSGeom_createLinearRing_r( geosinit.ctxt, GEOSCoordSeq_create_r( geosinit.ctxt, 0, 2 ) );
323 #endif
324  }
325  else
326  {
327  shell = rings[0];
328  }
329 
330  GEOSGeometry **holes = NULL;
331  int nHoles = 0;
332 
333  if ( rings.size() > 1 )
334  {
335  nHoles = rings.size() - 1;
336  holes = new GEOSGeometry*[ nHoles ];
337  if ( !holes )
338  return 0;
339 
340  for ( int i = 0; i < nHoles; i++ )
341  holes[i] = rings[i+1];
342  }
343 
344  GEOSGeometry *geom = GEOSGeom_createPolygon_r( geosinit.ctxt, shell, holes, nHoles );
345 
346  if ( holes )
347  delete [] holes;
348 
349  return geom;
350 }
351 
352 static GEOSGeometry *createGeosPolygon( GEOSGeometry *shell )
353 {
354  return createGeosPolygon( QVector<GEOSGeometry*>() << shell );
355 }
356 
357 static GEOSGeometry *createGeosPolygon( const QgsPolygon& polygon )
358 {
359  if ( polygon.count() == 0 )
360  return 0;
361 
362  QVector<GEOSGeometry *> geoms;
363 
364  try
365  {
366  for ( int i = 0; i < polygon.count(); i++ )
367  {
368  GEOSGeometry *ring = createGeosLinearRing( polygon[i] );
369  if ( !ring )
370  {
371  // something went really wrong - exit
372  for ( int j = 0; j < geoms.count(); j++ )
373  GEOSGeom_destroy_r( geosinit.ctxt, geoms[j] );
374  // XXX [MD] we just silently return here - but we shouldn't
375  // this is just because maptopixel simplification can return invalid linear rings
376  return 0;
377  }
378  geoms << ring;
379  }
380 
381  return createGeosPolygon( geoms );
382  }
383  catch ( GEOSException &e )
384  {
385  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
386  for ( int i = 0; i < geoms.count(); i++ )
387  GEOSGeom_destroy_r( geosinit.ctxt, geoms[i] );
388  return 0;
389  }
390 }
391 
392 static QgsGeometry *fromGeosGeom( GEOSGeometry *geom )
393 {
394  if ( !geom )
395  return 0;
396 
397  QgsGeometry *g = new QgsGeometry;
398  g->fromGeos( geom );
399  return g;
400 }
401 
403 {
404  try
405  {
406  GEOSWKTReader *reader = GEOSWKTReader_create_r( geosinit.ctxt );
407  QgsGeometry *g = fromGeosGeom( GEOSWKTReader_read_r( geosinit.ctxt, reader, wkt.toLocal8Bit().data() ) );
408  GEOSWKTReader_destroy_r( geosinit.ctxt, reader );
409  return g;
410  }
411  catch ( GEOSException &e )
412  {
413  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
414  return 0;
415  }
416 }
417 
419 {
420  return fromGeosGeom( createGeosPoint( point ) );
421 }
422 
424 {
425  return fromGeosGeom( createGeosLineString( polyline ) );
426 }
427 
429 {
430  return fromGeosGeom( createGeosPolygon( polygon ) );
431 }
432 
434 {
435  QVector<GEOSGeometry *> geoms;
436 
437  try
438  {
439  for ( int i = 0; i < multipoint.size(); ++i )
440  geoms << createGeosPoint( multipoint[i] );
441 
442  return fromGeosGeom( createGeosCollection( GEOS_MULTIPOINT, geoms ) );
443  }
444  catch ( GEOSException &e )
445  {
446  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
447 
448  for ( int i = 0; i < geoms.size(); ++i )
449  GEOSGeom_destroy_r( geosinit.ctxt, geoms[i] );
450 
451  return 0;
452  }
453 }
454 
456 {
457  QVector<GEOSGeometry *> geoms;
458 
459  try
460  {
461  for ( int i = 0; i < multiline.count(); i++ )
462  geoms << createGeosLineString( multiline[i] );
463 
464  return fromGeosGeom( createGeosCollection( GEOS_MULTILINESTRING, geoms ) );
465  }
466  catch ( GEOSException &e )
467  {
468  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
469 
470  for ( int i = 0; i < geoms.count(); i++ )
471  GEOSGeom_destroy_r( geosinit.ctxt, geoms[i] );
472 
473  return 0;
474  }
475 }
476 
478 {
479  if ( multipoly.count() == 0 )
480  return 0;
481 
482  QVector<GEOSGeometry *> geoms;
483 
484  try
485  {
486  for ( int i = 0; i < multipoly.count(); i++ )
487  geoms << createGeosPolygon( multipoly[i] );
488 
489  return fromGeosGeom( createGeosCollection( GEOS_MULTIPOLYGON, geoms ) );
490  }
491  catch ( GEOSException &e )
492  {
493  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
494 
495  for ( int i = 0; i < geoms.count(); i++ )
496  GEOSGeom_destroy_r( geosinit.ctxt, geoms[i] );
497 
498  return 0;
499  }
500 }
501 
503 {
504  QgsPolyline ring;
505  ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
506  ring.append( QgsPoint( rect.xMaximum(), rect.yMinimum() ) );
507  ring.append( QgsPoint( rect.xMaximum(), rect.yMaximum() ) );
508  ring.append( QgsPoint( rect.xMinimum(), rect.yMaximum() ) );
509  ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) );
510 
511  QgsPolygon polygon;
512  polygon.append( ring );
513 
514  return fromPolygon( polygon );
515 }
516 
517 QgsGeometry *QgsGeometry::fromQPointF( const QPointF &point )
518 {
519  return fromGeosGeom( createGeosPoint( point.x(), point.y() ) );
520 }
521 
522 QgsGeometry *QgsGeometry::fromQPolygonF( const QPolygonF &polygon )
523 {
524  if ( polygon.isClosed() )
525  {
526  return QgsGeometry::fromPolygon( createPolygonFromQPolygonF( polygon ) );
527  }
528  else
529  {
530  return QgsGeometry::fromPolyline( createPolylineFromQPolygonF( polygon ) );
531  }
532 }
533 
534 QgsPolygon QgsGeometry::createPolygonFromQPolygonF( const QPolygonF &polygon )
535 {
536  QgsPolygon result;
537  result << createPolylineFromQPolygonF( polygon );
538  return result;
539 }
540 
541 QgsPolyline QgsGeometry::createPolylineFromQPolygonF( const QPolygonF &polygon )
542 {
543  QgsPolyline result;
544  QPolygonF::const_iterator it = polygon.constBegin();
545  for ( ; it != polygon.constEnd(); ++it )
546  {
547  result.append( QgsPoint( *it ) );
548  }
549  return result;
550 }
551 
553 {
554  if ( &rhs == this )
555  return *this;
556 
557  // remove old geometry if it exists
558  if ( mGeometry )
559  {
560  delete [] mGeometry;
561  mGeometry = 0;
562  }
563 
564  mGeometrySize = rhs.mGeometrySize;
565 
566  // deep-copy the GEOS Geometry if appropriate
567  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
568  mGeos = rhs.mGeos ? GEOSGeom_clone_r( geosinit.ctxt, rhs.mGeos ) : 0;
569 
570  mDirtyGeos = rhs.mDirtyGeos;
571  mDirtyWkb = rhs.mDirtyWkb;
572 
573  if ( mGeometrySize && rhs.mGeometry )
574  {
575  mGeometry = new unsigned char[mGeometrySize];
576  memcpy( mGeometry, rhs.mGeometry, mGeometrySize );
577  }
578 
579  return *this;
580 } // QgsGeometry::operator=( QgsGeometry const & rhs )
581 
582 
583 void QgsGeometry::fromWkb( unsigned char *wkb, size_t length )
584 {
585  // delete any existing WKB geometry before assigning new one
586  if ( mGeometry )
587  {
588  delete [] mGeometry;
589  mGeometry = 0;
590  }
591 
592  if ( mGeos )
593  {
594  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
595  mGeos = 0;
596  }
597 
598  mGeometry = wkb;
599  mGeometrySize = length;
600 
601  mDirtyWkb = false;
602  mDirtyGeos = true;
603 }
604 
605 const unsigned char *QgsGeometry::asWkb() const
606 {
607  if ( mDirtyWkb )
608  exportGeosToWkb();
609 
610  return mGeometry;
611 }
612 
613 size_t QgsGeometry::wkbSize() const
614 {
615  if ( mDirtyWkb )
616  exportGeosToWkb();
617 
618  return mGeometrySize;
619 }
620 
621 const GEOSGeometry* QgsGeometry::asGeos() const
622 {
623  if ( mDirtyGeos )
624  {
625  if ( !exportWkbToGeos() )
626  {
627  return 0;
628  }
629  }
630 
631  return mGeos;
632 }
633 
634 
636 {
637  QgsConstWkbPtr wkbPtr( asWkb() + 1 ); // ensure that wkb representation exists
638 
639  if ( mGeometry && wkbSize() >= 5 )
640  {
642  wkbPtr >> wkbType;
643  return wkbType;
644  }
645  else
646  {
647  return QGis::WKBUnknown;
648  }
649 }
650 
651 
653 {
654  if ( mDirtyWkb )
655  exportGeosToWkb();
656 
657  switch ( wkbType() )
658  {
659  case QGis::WKBPoint:
660  case QGis::WKBPoint25D:
661  case QGis::WKBMultiPoint:
663  return QGis::Point;
664 
665  case QGis::WKBLineString:
669  return QGis::Line;
670 
671  case QGis::WKBPolygon:
672  case QGis::WKBPolygon25D:
675  return QGis::Polygon;
676 
677  default:
678  return QGis::UnknownGeometry;
679  }
680 }
681 
683 {
684  if ( mDirtyWkb )
685  exportGeosToWkb();
686 
687  return QGis::isMultiType( wkbType() );
688 }
689 
690 void QgsGeometry::fromGeos( GEOSGeometry *geos )
691 {
692  // TODO - make this more heap-friendly
693 
694  if ( mGeos )
695  {
696  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
697  mGeos = 0;
698  }
699 
700  if ( mGeometry )
701  {
702  delete [] mGeometry;
703  mGeometry = 0;
704  }
705 
706  mGeos = geos;
707 
708  mDirtyWkb = true;
709  mDirtyGeos = false;
710 }
711 
712 QgsPoint QgsGeometry::closestVertex( const QgsPoint& point, int& atVertex, int& beforeVertex, int& afterVertex, double& sqrDist )
713 {
714  // TODO: implement with GEOS
715  if ( mDirtyWkb )
716  exportGeosToWkb();
717 
718  if ( !mGeometry )
719  {
720  QgsDebugMsg( "WKB geometry not available!" );
721  return QgsPoint( 0, 0 );
722  }
723 
724  double actdist = std::numeric_limits<double>::max();
725 
726  beforeVertex = -1;
727  afterVertex = -1;
728 
729  QgsWkbPtr wkbPtr( mGeometry + 1 );
731  wkbPtr >> wkbType;
732 
733  QgsPoint p;
734  bool hasZValue = false;
735  int vertexnr = -1;
736  switch ( wkbType )
737  {
738  case QGis::WKBPoint25D:
739  hasZValue = true;
740  //intentional fall-through
741  case QGis::WKBPoint:
742  {
743  double x, y;
744  wkbPtr >> x >> y;
745  p.set( x, y );
746  actdist = point.sqrDist( x, y );
747  vertexnr = 0;
748  break;
749  }
750 
752  hasZValue = true;
753  //intentional fall-through
754  case QGis::WKBLineString:
755  {
756  int nPoints;
757  wkbPtr >> nPoints;
758  for ( int index = 0; index < nPoints; ++index )
759  {
760  double x, y;
761  wkbPtr >> x >> y;
762  if ( hasZValue )
763  wkbPtr += sizeof( double );
764 
765  double dist = point.sqrDist( x, y );
766  if ( dist < actdist )
767  {
768  p.set( x, y );
769  actdist = dist;
770  vertexnr = index;
771 
772  beforeVertex = index - 1;
773  afterVertex = index == nPoints - 1 ? -1 : index + 1;
774  }
775  }
776  break;
777  }
778 
779  case QGis::WKBPolygon25D:
780  hasZValue = true;
781  //intentional fall-through
782  case QGis::WKBPolygon:
783  {
784  int nRings;
785  wkbPtr >> nRings;
786  for ( int index = 0, pointIndex = 0; index < nRings; ++index )
787  {
788  int nPoints;
789  wkbPtr >> nPoints;
790  for ( int index2 = 0; index2 < nPoints; ++index2 )
791  {
792  double x, y;
793  wkbPtr >> x >> y;
794  if ( hasZValue )
795  wkbPtr += sizeof( double );
796 
797  double dist = point.sqrDist( x, y );
798  if ( dist < actdist )
799  {
800  p.set( x, y );
801  actdist = dist;
802  vertexnr = pointIndex;
803 
804  // assign the rubberband indices
805  if ( index2 == 0 )
806  {
807  beforeVertex = pointIndex + ( nPoints - 2 );
808  afterVertex = pointIndex + 1;
809  }
810  else if ( index2 == nPoints - 1 )
811  {
812  beforeVertex = pointIndex - 1;
813  afterVertex = pointIndex - ( nPoints - 2 );
814  }
815  else
816  {
817  beforeVertex = pointIndex - 1;
818  afterVertex = pointIndex + 1;
819  }
820  }
821  ++pointIndex;
822  }
823  }
824  break;
825  }
826 
828  hasZValue = true;
829  //intentional fall-through
830  case QGis::WKBMultiPoint:
831  {
832  int nPoints;
833  wkbPtr >> nPoints;
834  for ( int index = 0; index < nPoints; ++index )
835  {
836  wkbPtr += 1 + sizeof( int ); // skip endian and point type
837 
838  double x, y;
839  wkbPtr >> x >> y;
840  if ( hasZValue )
841  wkbPtr += sizeof( double );
842 
843  double dist = point.sqrDist( x, y );
844  if ( dist < actdist )
845  {
846  p.set( x, y );
847  actdist = dist;
848  vertexnr = index;
849  }
850  }
851  break;
852  }
853 
855  hasZValue = true;
856  //intentional fall-through
858  {
859  int nLines;
860  wkbPtr >> nLines;
861  for ( int index = 0, pointIndex = 0; index < nLines; ++index )
862  {
863  wkbPtr += 1 + sizeof( int );
864 
865  int nPoints;
866  wkbPtr >> nPoints;
867  for ( int index2 = 0; index2 < nPoints; ++index2 )
868  {
869  double x, y;
870  wkbPtr >> x >> y;
871  if ( hasZValue )
872  wkbPtr += sizeof( double );
873 
874  double dist = point.sqrDist( x, y );
875  if ( dist < actdist )
876  {
877  p.set( x, y );
878  actdist = dist;
879  vertexnr = pointIndex;
880 
881  if ( index2 == 0 )//assign the rubber band indices
882  beforeVertex = -1;
883  else
884  beforeVertex = vertexnr - 1;
885 
886  if ( index2 == nPoints - 1 )
887  afterVertex = -1;
888  else
889  afterVertex = vertexnr + 1;
890 
891  }
892  ++pointIndex;
893  }
894  }
895  break;
896  }
897 
899  hasZValue = true;
900  //intentional fall-through
902  {
903  int nPolys;
904  wkbPtr >> nPolys;
905  for ( int index = 0, pointIndex = 0; index < nPolys; ++index )
906  {
907  wkbPtr += 1 + sizeof( int ); //skip endian and polygon type
908  int nRings;
909  wkbPtr >> nRings;
910  for ( int index2 = 0; index2 < nRings; ++index2 )
911  {
912  int nPoints;
913  wkbPtr >> nPoints;
914  for ( int index3 = 0; index3 < nPoints; ++index3 )
915  {
916  double x, y;
917  wkbPtr >> x >> y;
918  if ( hasZValue )
919  wkbPtr += sizeof( double );
920 
921  double dist = point.sqrDist( x, y );
922  if ( dist < actdist )
923  {
924  p.set( x, y );
925  actdist = dist;
926  vertexnr = pointIndex;
927 
928  //assign the rubber band indices
929  if ( index3 == 0 )
930  {
931  beforeVertex = pointIndex + ( nPoints - 2 );
932  afterVertex = pointIndex + 1;
933  }
934  else if ( index3 == nPoints - 1 )
935  {
936  beforeVertex = pointIndex - 1;
937  afterVertex = pointIndex - ( nPoints - 2 );
938  }
939  else
940  {
941  beforeVertex = pointIndex - 1;
942  afterVertex = pointIndex + 1;
943  }
944  }
945  ++pointIndex;
946  }
947  }
948  }
949  break;
950  }
951 
952  default:
953  break;
954  }
955 
956  sqrDist = actdist;
957  atVertex = vertexnr;
958  return p;
959 }
960 
961 void QgsGeometry::adjacentVertices( int atVertex, int& beforeVertex, int& afterVertex )
962 {
963  // TODO: implement with GEOS
964  if ( mDirtyWkb )
965  exportGeosToWkb();
966 
967  beforeVertex = -1;
968  afterVertex = -1;
969 
970  if ( !mGeometry )
971  {
972  QgsDebugMsg( "WKB geometry not available!" );
973  return;
974  }
975 
976  if ( atVertex < 0 )
977  return;
978 
980  bool hasZValue = false;
981  QgsWkbPtr wkbPtr( mGeometry + 1 );
982  wkbPtr >> wkbType;
983 
984  switch ( wkbType )
985  {
986  case QGis::WKBPoint:
987  {
988  // NOOP - Points do not have adjacent verticies
989  break;
990  }
991 
993  case QGis::WKBLineString:
994  {
995  int nPoints;
996  wkbPtr >> nPoints;
997 
998  if ( atVertex >= nPoints )
999  return;
1000 
1001  const int index = atVertex;
1002 
1003  // assign the rubber band indices
1004 
1005  beforeVertex = index - 1;
1006 
1007  if ( index == nPoints - 1 )
1008  afterVertex = -1;
1009  else
1010  afterVertex = index + 1;
1011 
1012  break;
1013  }
1014 
1015  case QGis::WKBPolygon25D:
1016  hasZValue = true;
1017  //intentional fall-through
1018  case QGis::WKBPolygon:
1019  {
1020  int nRings;
1021  wkbPtr >> nRings;
1022 
1023  for ( int index0 = 0, pointIndex = 0; index0 < nRings; ++index0 )
1024  {
1025  int nPoints;
1026  wkbPtr >> nPoints;
1027  for ( int index1 = 0; index1 < nPoints; ++index1 )
1028  {
1029  wkbPtr += ( hasZValue ? 3 : 2 ) * sizeof( double );
1030 
1031  if ( pointIndex == atVertex )
1032  {
1033  if ( index1 == 0 )
1034  {
1035  beforeVertex = pointIndex + ( nPoints - 2 );
1036  afterVertex = pointIndex + 1;
1037  }
1038  else if ( index1 == nPoints - 1 )
1039  {
1040  beforeVertex = pointIndex - 1;
1041  afterVertex = pointIndex - ( nPoints - 2 );
1042  }
1043  else
1044  {
1045  beforeVertex = pointIndex - 1;
1046  afterVertex = pointIndex + 1;
1047  }
1048  }
1049 
1050  ++pointIndex;
1051  }
1052  }
1053  break;
1054  }
1055 
1057  case QGis::WKBMultiPoint:
1058  {
1059  // NOOP - Points do not have adjacent verticies
1060  break;
1061  }
1062 
1064  hasZValue = true;
1065  //intentional fall-through
1067  {
1068  int nLines;
1069  wkbPtr >> nLines;
1070  for ( int index0 = 0, pointIndex = 0; index0 < nLines; ++index0 )
1071  {
1072  wkbPtr += 1 + sizeof( int );
1073 
1074  int nPoints;
1075  wkbPtr >> nPoints;
1076 
1077  for ( int index1 = 0; index1 < nPoints; ++index1 )
1078  {
1079  wkbPtr += ( hasZValue ? 3 : 2 ) * sizeof( double );
1080 
1081  if ( pointIndex == atVertex )
1082  {
1083  // Found the vertex of the linestring we were looking for.
1084  if ( index1 == 0 )
1085  beforeVertex = -1;
1086  else
1087  beforeVertex = pointIndex - 1;
1088 
1089  if ( index1 == nPoints - 1 )
1090  afterVertex = -1;
1091  else
1092  afterVertex = pointIndex + 1;
1093  }
1094 
1095  ++pointIndex;
1096  }
1097  }
1098 
1099  break;
1100  }
1101 
1103  hasZValue = true;
1104  //intentional fall-through
1105  case QGis::WKBMultiPolygon:
1106  {
1107  int nPolys;
1108  wkbPtr >> nPolys;
1109  for ( int index0 = 0, pointIndex = 0; index0 < nPolys; ++index0 )
1110  {
1111  wkbPtr += 1 + sizeof( int ); //skip endian and polygon type
1112  int nRings;
1113  wkbPtr >> nRings;
1114 
1115  for ( int index1 = 0; index1 < nRings; ++index1 )
1116  {
1117  int nPoints;
1118  wkbPtr >> nPoints;
1119  for ( int index2 = 0; index2 < nPoints; ++index2 )
1120  {
1121  wkbPtr += ( hasZValue ? 3 : 2 ) * sizeof( double );
1122 
1123  if ( pointIndex == atVertex )
1124  {
1125  // Found the vertex of the linear-ring of the polygon we were looking for.
1126  // assign the rubber band indices
1127 
1128  if ( index2 == 0 )
1129  {
1130  beforeVertex = pointIndex + ( nPoints - 2 );
1131  afterVertex = pointIndex + 1;
1132  }
1133  else if ( index2 == nPoints - 1 )
1134  {
1135  beforeVertex = pointIndex - 1;
1136  afterVertex = pointIndex - ( nPoints - 2 );
1137  }
1138  else
1139  {
1140  beforeVertex = pointIndex - 1;
1141  afterVertex = pointIndex + 1;
1142  }
1143  }
1144  ++pointIndex;
1145  }
1146  }
1147  }
1148 
1149  break;
1150  }
1151 
1152  default:
1153  break;
1154  } // switch (wkbType)
1155 }
1156 
1157 bool QgsGeometry::insertVertex( double x, double y,
1158  int beforeVertex,
1159  const GEOSCoordSequence *old_sequence,
1160  GEOSCoordSequence **new_sequence )
1161 {
1162  // Bounds checking
1163  if ( beforeVertex < 0 )
1164  {
1165  *new_sequence = 0;
1166  return false;
1167  }
1168 
1169  unsigned int numPoints;
1170  GEOSCoordSeq_getSize_r( geosinit.ctxt, old_sequence, &numPoints );
1171 
1172  *new_sequence = GEOSCoordSeq_create_r( geosinit.ctxt, numPoints + 1, 2 );
1173  if ( !*new_sequence )
1174  return false;
1175 
1176  bool inserted = false;
1177  for ( unsigned int i = 0, j = 0; i < numPoints; i++, j++ )
1178  {
1179  // Do we insert the new vertex here?
1180  if ( beforeVertex == static_cast<int>( i ) )
1181  {
1182  GEOSCoordSeq_setX_r( geosinit.ctxt, *new_sequence, j, x );
1183  GEOSCoordSeq_setY_r( geosinit.ctxt, *new_sequence, j, y );
1184  j++;
1185  inserted = true;
1186  }
1187 
1188  double aX, aY;
1189  GEOSCoordSeq_getX_r( geosinit.ctxt, old_sequence, i, &aX );
1190  GEOSCoordSeq_getY_r( geosinit.ctxt, old_sequence, i, &aY );
1191 
1192  GEOSCoordSeq_setX_r( geosinit.ctxt, *new_sequence, j, aX );
1193  GEOSCoordSeq_setY_r( geosinit.ctxt, *new_sequence, j, aY );
1194  }
1195 
1196  if ( !inserted )
1197  {
1198  // The beforeVertex is greater than the actual number of vertices
1199  // in the geometry - append it.
1200  GEOSCoordSeq_setX_r( geosinit.ctxt, *new_sequence, numPoints, x );
1201  GEOSCoordSeq_setY_r( geosinit.ctxt, *new_sequence, numPoints, y );
1202  }
1203 
1204  // TODO: Check that the sequence is still simple, e.g. with GEOS_GEOM::Geometry->isSimple()
1205 
1206  return inserted;
1207 }
1208 
1209 bool QgsGeometry::moveVertex( QgsWkbPtr &wkbPtr, const double &x, const double &y, int atVertex, bool hasZValue, int &pointIndex, bool isRing )
1210 {
1211  int nPoints;
1212  wkbPtr >> nPoints;
1213 
1214  const int ps = ( hasZValue ? 3 : 2 ) * sizeof( double );
1215 
1216  // Not this linestring/ring?
1217  if ( atVertex >= pointIndex + nPoints )
1218  {
1219  wkbPtr += ps * nPoints;
1220  pointIndex += nPoints;
1221  return false;
1222  }
1223 
1224  if ( isRing && atVertex == pointIndex + nPoints - 1 )
1225  atVertex = pointIndex;
1226 
1227  // Goto point in this linestring/ring
1228  wkbPtr += ps * ( atVertex - pointIndex );
1229  wkbPtr << x << y;
1230  if ( hasZValue )
1231  wkbPtr << 0.0;
1232 
1233  if ( isRing && atVertex == pointIndex )
1234  {
1235  wkbPtr += ps * ( nPoints - 2 );
1236  wkbPtr << x << y;
1237  if ( hasZValue )
1238  wkbPtr << 0.0;
1239  }
1240 
1241  return true;
1242 }
1243 
1244 bool QgsGeometry::moveVertex( double x, double y, int atVertex )
1245 {
1246  if ( atVertex < 0 )
1247  return false;
1248 
1249  if ( mDirtyWkb )
1250  exportGeosToWkb();
1251 
1252  if ( !mGeometry )
1253  {
1254  QgsDebugMsg( "WKB geometry not available!" );
1255  return false;
1256  }
1257 
1259  bool hasZValue = false;
1260  QgsWkbPtr wkbPtr( mGeometry + 1 );
1261  wkbPtr >> wkbType;
1262 
1263  switch ( wkbType )
1264  {
1265  case QGis::WKBPoint25D:
1266  hasZValue = true;
1267  //intentional fall-through
1268  case QGis::WKBPoint:
1269  {
1270  if ( atVertex != 0 )
1271  return false;
1272 
1273  wkbPtr << x << y;
1274  mDirtyGeos = true;
1275  return true;
1276  }
1277 
1279  hasZValue = true;
1280  //intentional fall-through
1281  case QGis::WKBLineString:
1282  {
1283  int pointIndex = 0;
1284  if ( moveVertex( wkbPtr, x, y, atVertex, hasZValue, pointIndex, false ) )
1285  {
1286  mDirtyGeos = true;
1287  return true;
1288  }
1289 
1290  return false;
1291  }
1292 
1294  hasZValue = true;
1295  //intentional fall-through
1296  case QGis::WKBMultiPoint:
1297  {
1298  int nPoints;
1299  wkbPtr >> nPoints;
1300 
1301  if ( atVertex < nPoints )
1302  {
1303  wkbPtr += atVertex * ( 1 + sizeof( int ) + ( hasZValue ? 3 : 2 ) * sizeof( double ) ) + 1 + sizeof( int );
1304  wkbPtr << x << y;
1305  if ( hasZValue )
1306  wkbPtr << 0.0;
1307  mDirtyGeos = true;
1308  return true;
1309  }
1310  else
1311  {
1312  return false;
1313  }
1314  }
1315 
1317  hasZValue = true;
1318  //intentional fall-through
1320  {
1321  int nLines;
1322  wkbPtr >> nLines;
1323 
1324  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
1325  {
1326  wkbPtr += 1 + sizeof( int );
1327  if ( moveVertex( wkbPtr, x, y, atVertex, hasZValue, pointIndex, false ) )
1328  {
1329  mDirtyGeos = true;
1330  return true;
1331  }
1332  }
1333 
1334  return false;
1335  }
1336 
1337  case QGis::WKBPolygon25D:
1338  hasZValue = true;
1339  //intentional fall-through
1340  case QGis::WKBPolygon:
1341  {
1342  int nLines;
1343  wkbPtr >> nLines;
1344 
1345  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
1346  {
1347  if ( moveVertex( wkbPtr, x, y, atVertex, hasZValue, pointIndex, true ) )
1348  {
1349  mDirtyGeos = true;
1350  return true;
1351  }
1352  }
1353  return false;
1354  }
1355 
1357  hasZValue = true;
1358  //intentional fall-through
1359  case QGis::WKBMultiPolygon:
1360  {
1361  int nPolygons;
1362  wkbPtr >> nPolygons;
1363  for ( int polynr = 0, pointIndex = 0; polynr < nPolygons; ++polynr )
1364  {
1365  wkbPtr += 1 + sizeof( int ); // skip endian and polygon type
1366 
1367  int nRings;
1368  wkbPtr >> nRings;
1369 
1370  for ( int ringnr = 0; ringnr < nRings; ++ringnr )
1371  {
1372  if ( moveVertex( wkbPtr, x, y, atVertex, hasZValue, pointIndex, true ) )
1373  {
1374  mDirtyGeos = true;
1375  return true;
1376  }
1377  }
1378  }
1379  return false;
1380  }
1381 
1382  default:
1383  return false;
1384  }
1385 }
1386 
1387 // copy vertices from srcPtr to dstPtr and skip/delete one vertex
1388 // @param srcPtr ring/part starting with number of points (adjusted in each call)
1389 // @param dstPtr ring/part to copy to (adjusted in each call)
1390 // @param atVertex index of vertex to skip
1391 // @param hasZValue points have 3 elements
1392 // @param pointIndex reference to index of first ring/part vertex in overall object (adjusted in each call)
1393 // @param isRing srcPtr points to a ring
1394 // @param lastItem last ring/part, atVertex after this one must be wrong
1395 // @return
1396 // 0 no delete was done
1397 // 1 "normal" delete was done
1398 // 2 last element of the ring/part was deleted
1399 int QgsGeometry::deleteVertex( QgsConstWkbPtr &srcPtr, QgsWkbPtr &dstPtr, int atVertex, bool hasZValue, int &pointIndex, bool isRing, bool lastItem )
1400 {
1401  QgsDebugMsg( QString( "atVertex:%1 hasZValue:%2 pointIndex:%3 isRing:%4" ).arg( atVertex ).arg( hasZValue ).arg( pointIndex ).arg( isRing ) );
1402  const int ps = ( hasZValue ? 3 : 2 ) * sizeof( double );
1403  int nPoints;
1404  srcPtr >> nPoints;
1405 
1406  // copy complete ring/part if vertex is in a following one
1407  if ( atVertex < pointIndex || atVertex >= pointIndex + nPoints )
1408  {
1409  // atVertex does not exist
1410  if ( lastItem && atVertex >= pointIndex + nPoints )
1411  return 0;
1412 
1413  dstPtr << nPoints;
1414 
1415  int len = nPoints * ps;
1416  memcpy( dstPtr, srcPtr, len );
1417  dstPtr += len;
1418  srcPtr += len;
1419  pointIndex += nPoints;
1420  return 0;
1421  }
1422 
1423  // delete the first vertex of a ring instead of the last
1424  if ( isRing && atVertex == pointIndex + nPoints - 1 )
1425  atVertex = pointIndex;
1426 
1427  if ( nPoints == ( isRing ? 2 : 1 ) )
1428  {
1429  // last point of the part/ring is deleted
1430  // skip the whole part/ring
1431  srcPtr += nPoints * ps;
1432  pointIndex += nPoints;
1433  return 2;
1434  }
1435 
1436  dstPtr << nPoints - 1;
1437 
1438  // copy ring before vertex
1439  int len = ( atVertex - pointIndex ) * ps;
1440  if ( len > 0 )
1441  {
1442  memcpy( dstPtr, srcPtr, len );
1443  dstPtr += len;
1444  srcPtr += len;
1445  }
1446 
1447  // skip deleted vertex
1448  srcPtr += ps;
1449 
1450  // copy reset of ring
1451  len = ( pointIndex + nPoints - atVertex - 1 ) * ps;
1452 
1453  // save position of vertex, if we delete the first vertex of a ring
1454  const unsigned char *first = 0;
1455  if ( isRing && atVertex == pointIndex )
1456  {
1457  len -= ps;
1458  first = srcPtr;
1459  }
1460 
1461  if ( len > 0 )
1462  {
1463  memcpy( dstPtr, srcPtr, len );
1464  dstPtr += len;
1465  srcPtr += len;
1466  }
1467 
1468  // copy new first vertex instead of the old last, if we deleted the original first vertex
1469  if ( first )
1470  {
1471  memcpy( dstPtr, first, ps );
1472  dstPtr += ps;
1473  srcPtr += ps;
1474  }
1475 
1476  pointIndex += nPoints;
1477 
1478  return 1;
1479 }
1480 
1481 bool QgsGeometry::deleteVertex( int atVertex )
1482 {
1483  QgsDebugMsg( QString( "atVertex:%1" ).arg( atVertex ) );
1484  if ( atVertex < 0 )
1485  return false;
1486 
1487  if ( mDirtyWkb )
1488  exportGeosToWkb();
1489 
1490  if ( !mGeometry )
1491  {
1492  QgsDebugMsg( "WKB geometry not available!" );
1493  return false;
1494  }
1495 
1496  QgsConstWkbPtr srcPtr( mGeometry );
1497  char endianness;
1499  srcPtr >> endianness >> wkbType;
1500 
1501  bool hasZValue = QGis::wkbDimensions( wkbType ) == 3;
1502 
1503  int ps = ( hasZValue ? 3 : 2 ) * sizeof( double );
1504  if ( QGis::flatType( wkbType ) == QGis::WKBMultiPoint )
1505  ps += 1 + sizeof( int );
1506 
1507  unsigned char *dstBuffer = new unsigned char[mGeometrySize - ps];
1508  QgsWkbPtr dstPtr( dstBuffer );
1509  dstPtr << endianness << wkbType;
1510 
1511  bool deleted = false;
1512  switch ( wkbType )
1513  {
1514  case QGis::WKBPoint25D:
1515  case QGis::WKBPoint:
1516  break; //cannot remove the only point vertex
1517 
1519  case QGis::WKBLineString:
1520  {
1521  int pointIndex = 0;
1522  int res = deleteVertex( srcPtr, dstPtr, atVertex, hasZValue, pointIndex, false, true );
1523  if ( res == 2 )
1524  {
1525  // Linestring with 0 points
1526  dstPtr << 0;
1527  }
1528 
1529  deleted = res != 0;
1530  break;
1531  }
1532 
1533  case QGis::WKBPolygon25D:
1534  case QGis::WKBPolygon:
1535  {
1536  int nRings;
1537  srcPtr >> nRings;
1538  QgsWkbPtr ptrN( dstPtr );
1539  dstPtr << nRings;
1540 
1541  for ( int ringnr = 0, pointIndex = 0; ringnr < nRings; ++ringnr )
1542  {
1543  int res = deleteVertex( srcPtr, dstPtr, atVertex, hasZValue, pointIndex, true, ringnr == nRings - 1 );
1544  if ( res == 2 )
1545  ptrN << nRings - 1;
1546 
1547  deleted |= res != 0;
1548  }
1549 
1550  break;
1551  }
1552 
1554  case QGis::WKBMultiPoint:
1555  {
1556  int nPoints;
1557  srcPtr >> nPoints;
1558 
1559  if ( atVertex < nPoints )
1560  {
1561  dstPtr << nPoints - 1;
1562 
1563  int len = ps * atVertex;
1564  if ( len > 0 )
1565  {
1566  memcpy( dstPtr, srcPtr, len );
1567  srcPtr += len;
1568  dstPtr += len;
1569  }
1570 
1571  srcPtr += ps;
1572 
1573  len = ps * ( nPoints - atVertex - 1 );
1574  if ( len > 0 )
1575  {
1576  memcpy( dstPtr, srcPtr, len );
1577  srcPtr += len;
1578  dstPtr += len;
1579  }
1580 
1581  deleted = true;
1582  }
1583 
1584  break;
1585  }
1586 
1589  {
1590  int nLines;
1591  srcPtr >> nLines;
1592  QgsWkbPtr ptrN( dstPtr );
1593  dstPtr << nLines;
1594 
1595  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
1596  {
1597  QgsWkbPtr saveDstPtr( dstPtr );
1598  srcPtr >> endianness >> wkbType;
1599  dstPtr << endianness << wkbType;
1600 
1601  int res = deleteVertex( srcPtr, dstPtr, atVertex, hasZValue, pointIndex, false, linenr == nLines - 1 );
1602  if ( res == 2 )
1603  {
1604  // line string was completely removed
1605  ptrN << nLines - 1;
1606  dstPtr = saveDstPtr;
1607  }
1608 
1609  deleted |= res != 0;
1610  }
1611 
1612  break;
1613  }
1614 
1616  case QGis::WKBMultiPolygon:
1617  {
1618  int nPolys;
1619  srcPtr >> nPolys;
1620  QgsWkbPtr ptrNPolys( dstPtr );
1621  dstPtr << nPolys;
1622 
1623  for ( int polynr = 0, pointIndex = 0; polynr < nPolys; ++polynr )
1624  {
1625  int nRings;
1626  srcPtr >> endianness >> wkbType >> nRings;
1627  QgsWkbPtr saveDstPolyPtr( dstPtr );
1628  dstPtr << endianness << wkbType;
1629  QgsWkbPtr ptrNRings( dstPtr );
1630  dstPtr << nRings;
1631 
1632  for ( int ringnr = 0; ringnr < nRings; ++ringnr )
1633  {
1634  int res = deleteVertex( srcPtr, dstPtr, atVertex, hasZValue, pointIndex, true, polynr == nPolys - 1 && ringnr == nRings - 1 );
1635  if ( res == 2 )
1636  {
1637  // ring was completely removed
1638  if ( nRings == 1 )
1639  {
1640  // last ring => remove polygon
1641  ptrNPolys << nPolys - 1;
1642  dstPtr = saveDstPolyPtr;
1643  }
1644  else
1645  {
1646  ptrNRings << nRings - 1;
1647  }
1648  }
1649 
1650  deleted |= res != 0;
1651  }
1652  }
1653  break;
1654  }
1655 
1656  case QGis::WKBNoGeometry:
1657  case QGis::WKBUnknown:
1658  break;
1659  }
1660 
1661  if ( deleted )
1662  {
1663  delete [] mGeometry;
1664  mGeometry = dstBuffer;
1665  mGeometrySize -= ps;
1666  mDirtyGeos = true;
1667  return true;
1668  }
1669  else
1670  {
1671  delete [] dstBuffer;
1672  return false;
1673  }
1674 }
1675 
1676 bool QgsGeometry::insertVertex( QgsConstWkbPtr &srcPtr, QgsWkbPtr &dstPtr, int beforeVertex, const double &x, const double &y, bool hasZValue, int &pointIndex, bool isRing )
1677 {
1678  int nPoints;
1679  srcPtr >> nPoints;
1680 
1681  bool insertHere = beforeVertex >= pointIndex && beforeVertex < pointIndex + nPoints;
1682 
1683  int len;
1684  if ( insertHere )
1685  {
1686  dstPtr << nPoints + 1;
1687  len = ( hasZValue ? 3 : 2 ) * ( beforeVertex - pointIndex ) * sizeof( double );
1688  if ( len > 0 )
1689  {
1690  memcpy( dstPtr, srcPtr, len );
1691  srcPtr += len;
1692  dstPtr += len;
1693  }
1694 
1695  dstPtr << x << y;
1696  if ( hasZValue )
1697  dstPtr << 0.0;
1698 
1699  len = ( hasZValue ? 3 : 2 ) * ( pointIndex + nPoints - beforeVertex ) * sizeof( double );
1700  if ( isRing && beforeVertex == pointIndex )
1701  len -= ( hasZValue ? 3 : 2 ) * sizeof( double );
1702  }
1703  else
1704  {
1705  dstPtr << nPoints;
1706  len = ( hasZValue ? 3 : 2 ) * nPoints * sizeof( double );
1707  }
1708 
1709  memcpy( dstPtr, srcPtr, len );
1710  srcPtr += len;
1711  dstPtr += len;
1712 
1713  if ( isRing && beforeVertex == pointIndex )
1714  {
1715  dstPtr << x << y;
1716  if ( hasZValue )
1717  dstPtr << 0.0;
1718  }
1719 
1720  pointIndex += nPoints;
1721  return insertHere;
1722 }
1723 
1724 bool QgsGeometry::insertVertex( double x, double y, int beforeVertex )
1725 {
1726  // TODO: implement with GEOS
1727  if ( mDirtyWkb )
1728  exportGeosToWkb();
1729 
1730  if ( !mGeometry )
1731  {
1732  QgsDebugMsg( "WKB geometry not available!" );
1733  return false;
1734  }
1735 
1736  if ( beforeVertex < 0 )
1737  return false;
1738 
1739  QgsConstWkbPtr srcPtr( mGeometry );
1740  char endianness;
1742  srcPtr >> endianness >> wkbType;
1743 
1744  bool hasZValue = QGis::wkbDimensions( wkbType ) == 3;
1745 
1746  int ps = ( hasZValue ? 3 : 2 ) * sizeof( double );
1747  if ( QGis::flatType( wkbType ) == QGis::WKBMultiPoint )
1748  ps += 1 + sizeof( int );
1749 
1750  unsigned char *dstBuffer = new unsigned char[mGeometrySize + ps];
1751  QgsWkbPtr dstPtr( dstBuffer );
1752  dstPtr << endianness << wkbType;
1753 
1754  bool inserted = false;
1755  switch ( wkbType )
1756  {
1757  case QGis::WKBPoint25D:
1758  case QGis::WKBPoint: //cannot insert a vertex before another one on point types
1759  break;
1760 
1762  case QGis::WKBLineString:
1763  {
1764  int pointIndex = 0;
1765  inserted = insertVertex( srcPtr, dstPtr, beforeVertex, x, y, hasZValue, pointIndex, false );
1766  break;
1767  }
1768 
1769  case QGis::WKBPolygon25D:
1770  case QGis::WKBPolygon:
1771  {
1772  int nRings;
1773  srcPtr >> nRings;
1774  dstPtr << nRings;
1775 
1776  for ( int ringnr = 0, pointIndex = 0; ringnr < nRings; ++ringnr )
1777  inserted |= insertVertex( srcPtr, dstPtr, beforeVertex, x, y, hasZValue, pointIndex, true );
1778 
1779  break;
1780  }
1781 
1783  case QGis::WKBMultiPoint:
1784  {
1785  int nPoints;
1786  srcPtr >> nPoints;
1787 
1788  if ( beforeVertex <= nPoints )
1789  {
1790  dstPtr << nPoints + 1;
1791 
1792  int len = ps * beforeVertex;
1793  if ( len > 0 )
1794  {
1795  memcpy( dstPtr, srcPtr, len );
1796  srcPtr += len;
1797  dstPtr += len;
1798  }
1799 
1800  dstPtr << endianness << ( hasZValue ? QGis::WKBPoint25D : QGis::WKBPoint ) << x << y;
1801  if ( hasZValue )
1802  dstPtr << 0.0;
1803 
1804  len = ps * ( nPoints - beforeVertex );
1805  if ( len > 0 )
1806  memcpy( dstPtr, srcPtr, len );
1807 
1808  inserted = true;
1809  }
1810 
1811  break;
1812  }
1813 
1816  {
1817  int nLines;
1818  srcPtr >> nLines;
1819  dstPtr << nLines;
1820 
1821  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
1822  {
1823  srcPtr >> endianness >> wkbType;
1824  dstPtr << endianness << wkbType;
1825  inserted |= insertVertex( srcPtr, dstPtr, beforeVertex, x, y, hasZValue, pointIndex, false );
1826  }
1827  break;
1828  }
1829 
1831  case QGis::WKBMultiPolygon:
1832  {
1833  int nPolys;
1834  srcPtr >> nPolys;
1835  dstPtr << nPolys;
1836 
1837  for ( int polynr = 0, pointIndex = 0; polynr < nPolys; ++polynr )
1838  {
1839  int nRings;
1840  srcPtr >> endianness >> wkbType >> nRings;
1841  dstPtr << endianness << wkbType << nRings;
1842 
1843  for ( int ringnr = 0; ringnr < nRings; ++ringnr )
1844  inserted |= insertVertex( srcPtr, dstPtr, beforeVertex, x, y, hasZValue, pointIndex, true );
1845  }
1846  break;
1847  }
1848 
1849  case QGis::WKBNoGeometry:
1850  case QGis::WKBUnknown:
1851  break;
1852  }
1853 
1854  if ( inserted )
1855  {
1856  delete [] mGeometry;
1857  mGeometry = dstBuffer;
1858  mGeometrySize += ps;
1859  mDirtyGeos = true;
1860  return true;
1861  }
1862  else
1863  {
1864  delete [] dstBuffer;
1865  return false;
1866  }
1867 }
1868 
1870 {
1871  if ( atVertex < 0 )
1872  return QgsPoint( 0, 0 );
1873 
1874  if ( mDirtyWkb )
1875  exportGeosToWkb();
1876 
1877  if ( !mGeometry )
1878  {
1879  QgsDebugMsg( "WKB geometry not available!" );
1880  return QgsPoint( 0, 0 );
1881  }
1882 
1883  QgsConstWkbPtr wkbPtr( mGeometry + 1 );
1885  wkbPtr >> wkbType;
1886 
1887  bool hasZValue = false;
1888  switch ( wkbType )
1889  {
1890  case QGis::WKBPoint25D:
1891  case QGis::WKBPoint:
1892  {
1893  if ( atVertex != 0 )
1894  return QgsPoint( 0, 0 );
1895 
1896  double x, y;
1897  wkbPtr >> x >> y;
1898 
1899  return QgsPoint( x, y );
1900  }
1901 
1903  hasZValue = true;
1904  //intentional fall-through
1905  case QGis::WKBLineString:
1906  {
1907  // get number of points in the line
1908  int nPoints;
1909  wkbPtr >> nPoints;
1910 
1911  if ( atVertex >= nPoints )
1912  return QgsPoint( 0, 0 );
1913 
1914  // copy the vertex coordinates
1915  wkbPtr += atVertex * ( hasZValue ? 3 : 2 ) * sizeof( double );
1916 
1917  double x, y;
1918  wkbPtr >> x >> y;
1919 
1920  return QgsPoint( x, y );
1921  }
1922 
1923  case QGis::WKBPolygon25D:
1924  hasZValue = true;
1925  //intentional fall-through
1926  case QGis::WKBPolygon:
1927  {
1928  int nRings;
1929  wkbPtr >> nRings;
1930 
1931  for ( int ringnr = 0, pointIndex = 0; ringnr < nRings; ++ringnr )
1932  {
1933  int nPoints;
1934  wkbPtr >> nPoints;
1935 
1936  if ( atVertex >= pointIndex + nPoints )
1937  {
1938  wkbPtr += nPoints * ( hasZValue ? 3 : 2 ) * sizeof( double );
1939  pointIndex += nPoints;
1940  continue;
1941  }
1942 
1943  wkbPtr += ( atVertex - pointIndex ) * ( hasZValue ? 3 : 2 ) * sizeof( double );
1944 
1945  double x, y;
1946  wkbPtr >> x >> y;
1947  return QgsPoint( x, y );
1948  }
1949 
1950  return QgsPoint( 0, 0 );
1951  }
1952 
1954  hasZValue = true;
1955  //intentional fall-through
1956  case QGis::WKBMultiPoint:
1957  {
1958  // get number of points in the line
1959  int nPoints;
1960  wkbPtr >> nPoints;
1961 
1962  if ( atVertex >= nPoints )
1963  return QgsPoint( 0, 0 );
1964 
1965  wkbPtr += atVertex * ( 1 + sizeof( int ) + ( hasZValue ? 3 : 2 ) * sizeof( double ) ) + 1 + sizeof( int );
1966 
1967  double x, y;
1968  wkbPtr >> x >> y;
1969  return QgsPoint( x, y );
1970  }
1971 
1973  hasZValue = true;
1974  //intentional fall-through
1976  {
1977  int nLines;
1978  wkbPtr >> nLines;
1979 
1980  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
1981  {
1982  wkbPtr += 1 + sizeof( int );
1983  int nPoints;
1984  wkbPtr >> nPoints;
1985 
1986  if ( atVertex >= pointIndex + nPoints )
1987  {
1988  wkbPtr += nPoints * ( hasZValue ? 3 : 2 ) * sizeof( double );
1989  pointIndex += nPoints;
1990  continue;
1991  }
1992 
1993  wkbPtr += ( atVertex - pointIndex ) * ( hasZValue ? 3 : 2 ) * sizeof( double );
1994 
1995  double x, y;
1996  wkbPtr >> x >> y;
1997 
1998  return QgsPoint( x, y );
1999  }
2000 
2001  return QgsPoint( 0, 0 );
2002  }
2003 
2005  hasZValue = true;
2006  //intentional fall-through
2007  case QGis::WKBMultiPolygon:
2008  {
2009  int nPolygons;
2010  wkbPtr >> nPolygons;
2011 
2012  for ( int polynr = 0, pointIndex = 0; polynr < nPolygons; ++polynr )
2013  {
2014  wkbPtr += 1 + sizeof( int );
2015 
2016  int nRings;
2017  wkbPtr >> nRings;
2018  for ( int ringnr = 0; ringnr < nRings; ++ringnr )
2019  {
2020  int nPoints;
2021  wkbPtr >> nPoints;
2022 
2023  if ( atVertex >= pointIndex + nPoints )
2024  {
2025  wkbPtr += nPoints * ( hasZValue ? 3 : 2 ) * sizeof( double );
2026  pointIndex += nPoints;
2027  continue;
2028  }
2029 
2030  wkbPtr += ( atVertex - pointIndex ) * ( hasZValue ? 3 : 2 ) * sizeof( double );
2031 
2032  double x, y;
2033  wkbPtr >> x >> y;
2034 
2035  return QgsPoint( x, y );
2036  }
2037  }
2038  return QgsPoint( 0, 0 );
2039  }
2040 
2041  default:
2042  QgsDebugMsg( "error: mGeometry type not recognized" );
2043  return QgsPoint( 0, 0 );
2044  }
2045 }
2046 
2047 double QgsGeometry::sqrDistToVertexAt( QgsPoint& point, int atVertex )
2048 {
2049  QgsPoint pnt = vertexAt( atVertex );
2050  if ( pnt != QgsPoint( 0, 0 ) )
2051  {
2052  QgsDebugMsg( "Exiting with distance to " + pnt.toString() );
2053  return point.sqrDist( pnt );
2054  }
2055  else
2056  {
2057  QgsDebugMsg( "Exiting with std::numeric_limits<double>::max()." );
2058  // probably safest to bail out with a very large number
2060  }
2061 }
2062 
2063 double QgsGeometry::closestVertexWithContext( const QgsPoint& point, int& atVertex )
2064 {
2065  double sqrDist = std::numeric_limits<double>::max();
2066 
2067  try
2068  {
2069  // Initialise some stuff
2070  int closestVertexIndex = 0;
2071 
2072  // set up the GEOS geometry
2073  if ( mDirtyGeos )
2074  exportWkbToGeos();
2075 
2076  if ( !mGeos )
2077  return -1;
2078 
2079  const GEOSGeometry *g = GEOSGetExteriorRing_r( geosinit.ctxt, mGeos );
2080  if ( !g )
2081  return -1;
2082 
2083  const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, g );
2084 
2085  unsigned int n;
2086  GEOSCoordSeq_getSize_r( geosinit.ctxt, sequence, &n );
2087 
2088  for ( unsigned int i = 0; i < n; i++ )
2089  {
2090  double x, y;
2091  GEOSCoordSeq_getX_r( geosinit.ctxt, sequence, i, &x );
2092  GEOSCoordSeq_getY_r( geosinit.ctxt, sequence, i, &y );
2093 
2094  double testDist = point.sqrDist( x, y );
2095  if ( testDist < sqrDist )
2096  {
2097  closestVertexIndex = i;
2098  sqrDist = testDist;
2099  }
2100  }
2101 
2102  atVertex = closestVertexIndex;
2103  }
2104  catch ( GEOSException &e )
2105  {
2106  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2107  return -1;
2108  }
2109 
2110  return sqrDist;
2111 }
2112 
2114  const QgsPoint& point,
2115  QgsPoint& minDistPoint,
2116  int& afterVertex,
2117  double *leftOf,
2118  double epsilon )
2119 {
2120  QgsDebugMsgLevel( "Entering.", 3 );
2121 
2122  // TODO: implement with GEOS
2123  if ( mDirtyWkb ) //convert latest geos to mGeometry
2124  exportGeosToWkb();
2125 
2126  if ( !mGeometry )
2127  {
2128  QgsDebugMsg( "WKB geometry not available!" );
2129  return -1;
2130  }
2131 
2132  QgsWkbPtr wkbPtr( mGeometry + 1 );
2134  wkbPtr >> wkbType;
2135 
2136  // Initialise some stuff
2137  double sqrDist = std::numeric_limits<double>::max();
2138 
2139  QgsPoint distPoint;
2140  int closestSegmentIndex = 0;
2141  bool hasZValue = false;
2142  switch ( wkbType )
2143  {
2144  case QGis::WKBPoint25D:
2145  case QGis::WKBPoint:
2147  case QGis::WKBMultiPoint:
2148  {
2149  // Points have no lines
2150  return -1;
2151  }
2152 
2154  hasZValue = true;
2155  //intentional fall-through
2156  case QGis::WKBLineString:
2157  {
2158  int nPoints;
2159  wkbPtr >> nPoints;
2160 
2161  double prevx = 0.0, prevy = 0.0;
2162  for ( int index = 0; index < nPoints; ++index )
2163  {
2164  double thisx, thisy;
2165  wkbPtr >> thisx >> thisy;
2166  if ( hasZValue )
2167  wkbPtr += sizeof( double );
2168 
2169  if ( index > 0 )
2170  {
2171  double testdist = point.sqrDistToSegment( prevx, prevy, thisx, thisy, distPoint, epsilon );
2172  if ( testdist < sqrDist )
2173  {
2174  closestSegmentIndex = index;
2175  sqrDist = testdist;
2176  minDistPoint = distPoint;
2177  if ( leftOf )
2178  {
2179  *leftOf = QgsGeometry::leftOf( point.x(), point.y(), prevx, prevy, thisx, thisy );
2180  }
2181  }
2182  }
2183 
2184  prevx = thisx;
2185  prevy = thisy;
2186  }
2187  afterVertex = closestSegmentIndex;
2188  break;
2189  }
2190 
2192  hasZValue = true;
2193  //intentional fall-through
2195  {
2196  int nLines;
2197  wkbPtr >> nLines;
2198  for ( int linenr = 0, pointIndex = 0; linenr < nLines; ++linenr )
2199  {
2200  wkbPtr += 1 + sizeof( int );
2201  int nPoints;
2202  wkbPtr >> nPoints;
2203 
2204  double prevx = 0.0, prevy = 0.0;
2205  for ( int pointnr = 0; pointnr < nPoints; ++pointnr )
2206  {
2207  double thisx, thisy;
2208  wkbPtr >> thisx >> thisy;
2209  if ( hasZValue )
2210  wkbPtr += sizeof( double );
2211 
2212  if ( pointnr > 0 )
2213  {
2214  double testdist = point.sqrDistToSegment( prevx, prevy, thisx, thisy, distPoint, epsilon );
2215  if ( testdist < sqrDist )
2216  {
2217  closestSegmentIndex = pointIndex;
2218  sqrDist = testdist;
2219  minDistPoint = distPoint;
2220  if ( leftOf )
2221  {
2222  *leftOf = QgsGeometry::leftOf( point.x(), point.y(), prevx, prevy, thisx, thisy );
2223  }
2224  }
2225  }
2226 
2227  prevx = thisx;
2228  prevy = thisy;
2229  ++pointIndex;
2230  }
2231  }
2232  afterVertex = closestSegmentIndex;
2233  break;
2234  }
2235 
2236  case QGis::WKBPolygon25D:
2237  hasZValue = true;
2238  //intentional fall-through
2239  case QGis::WKBPolygon:
2240  {
2241  int nRings;
2242  wkbPtr >> nRings;
2243 
2244  for ( int ringnr = 0, pointIndex = 0; ringnr < nRings; ++ringnr )//loop over rings
2245  {
2246  int nPoints;
2247  wkbPtr >> nPoints;
2248 
2249  double prevx = 0.0, prevy = 0.0;
2250  for ( int pointnr = 0; pointnr < nPoints; ++pointnr )//loop over points in a ring
2251  {
2252  double thisx, thisy;
2253  wkbPtr >> thisx >> thisy;
2254  if ( hasZValue )
2255  wkbPtr += sizeof( double );
2256 
2257  if ( pointnr > 0 )
2258  {
2259  double testdist = point.sqrDistToSegment( prevx, prevy, thisx, thisy, distPoint, epsilon );
2260  if ( testdist < sqrDist )
2261  {
2262  closestSegmentIndex = pointIndex;
2263  sqrDist = testdist;
2264  minDistPoint = distPoint;
2265  if ( leftOf )
2266  {
2267  *leftOf = QgsGeometry::leftOf( point.x(), point.y(), prevx, prevy, thisx, thisy );
2268  }
2269  }
2270  }
2271 
2272  prevx = thisx;
2273  prevy = thisy;
2274  ++pointIndex;
2275  }
2276  }
2277  afterVertex = closestSegmentIndex;
2278  break;
2279  }
2280 
2282  hasZValue = true;
2283  //intentional fall-through
2284  case QGis::WKBMultiPolygon:
2285  {
2286  int nPolygons;
2287  wkbPtr >> nPolygons;
2288  for ( int polynr = 0, pointIndex = 0; polynr < nPolygons; ++polynr )
2289  {
2290  wkbPtr += 1 + sizeof( int );
2291  int nRings;
2292  wkbPtr >> nRings;
2293  for ( int ringnr = 0; ringnr < nRings; ++ringnr )
2294  {
2295  int nPoints;
2296  wkbPtr >> nPoints;
2297 
2298  double prevx = 0.0, prevy = 0.0;
2299  for ( int pointnr = 0; pointnr < nPoints; ++pointnr )
2300  {
2301  double thisx, thisy;
2302  wkbPtr >> thisx >> thisy;
2303  if ( hasZValue )
2304  wkbPtr += sizeof( double );
2305 
2306  if ( pointnr > 0 )
2307  {
2308  double testdist = point.sqrDistToSegment( prevx, prevy, thisx, thisy, distPoint, epsilon );
2309  if ( testdist < sqrDist )
2310  {
2311  closestSegmentIndex = pointIndex;
2312  sqrDist = testdist;
2313  minDistPoint = distPoint;
2314  if ( leftOf )
2315  {
2316  *leftOf = QgsGeometry::leftOf( point.x(), point.y(), prevx, prevy, thisx, thisy );
2317  }
2318  }
2319  }
2320 
2321  prevx = thisx;
2322  prevy = thisy;
2323  ++pointIndex;
2324  }
2325  }
2326  }
2327  afterVertex = closestSegmentIndex;
2328  break;
2329  }
2330 
2331  case QGis::WKBUnknown:
2332  default:
2333  return -1;
2334  break;
2335  } // switch (wkbType)
2336 
2337  QgsDebugMsgLevel( QString( "Exiting with nearest point %1, dist %2." )
2338  .arg( point.toString() ).arg( sqrDist ), 3 );
2339 
2340  return sqrDist;
2341 }
2342 
2343 int QgsGeometry::addRing( const QList<QgsPoint>& ring )
2344 {
2345  //bail out if this geometry is not polygon/multipolygon
2346  if ( type() != QGis::Polygon )
2347  return 1;
2348 
2349  //test for invalid geometries
2350  if ( ring.size() < 4 )
2351  return 3;
2352 
2353  //ring must be closed
2354  if ( ring.first() != ring.last() )
2355  return 2;
2356 
2357  //create geos geometry from wkb if not already there
2358  if ( mDirtyGeos )
2359  {
2360  exportWkbToGeos();
2361  }
2362 
2363  if ( !mGeos )
2364  {
2365  return 6;
2366  }
2367 
2368  int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
2369 
2370  //Fill GEOS Polygons of the feature into list
2371  QVector<const GEOSGeometry*> polygonList;
2372 
2373  if ( wkbType() == QGis::WKBPolygon )
2374  {
2375  if ( type != GEOS_POLYGON )
2376  return 1;
2377 
2378  polygonList << mGeos;
2379  }
2380  else if ( wkbType() == QGis::WKBMultiPolygon )
2381  {
2382  if ( type != GEOS_MULTIPOLYGON )
2383  return 1;
2384 
2385  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); ++i )
2386  polygonList << GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
2387  }
2388 
2389  //create new ring
2390  GEOSGeometry *newRing = 0;
2391  GEOSGeometry *newRingPolygon = 0;
2392 
2393  try
2394  {
2395  newRing = createGeosLinearRing( ring.toVector() );
2396  if ( !GEOSisValid_r( geosinit.ctxt, newRing ) )
2397  {
2398  throwGEOSException( "ring is invalid" );
2399  }
2400 
2401  newRingPolygon = createGeosPolygon( newRing );
2402  if ( !GEOSisValid_r( geosinit.ctxt, newRingPolygon ) )
2403  {
2404  throwGEOSException( "ring is invalid" );
2405  }
2406  }
2407  catch ( GEOSException &e )
2408  {
2409  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2410 
2411  if ( newRingPolygon )
2412  GEOSGeom_destroy_r( geosinit.ctxt, newRingPolygon );
2413  else if ( newRing )
2414  GEOSGeom_destroy_r( geosinit.ctxt, newRing );
2415 
2416  return 3;
2417  }
2418 
2419  QVector<GEOSGeometry*> rings;
2420 
2421  int i;
2422  for ( i = 0; i < polygonList.size(); i++ )
2423  {
2424  for ( int j = 0; j < rings.size(); j++ )
2425  GEOSGeom_destroy_r( geosinit.ctxt, rings[j] );
2426  rings.clear();
2427 
2428  GEOSGeometry *shellRing = 0;
2429  GEOSGeometry *shell = 0;
2430  try
2431  {
2432  shellRing = GEOSGeom_clone_r( geosinit.ctxt, GEOSGetExteriorRing_r( geosinit.ctxt, polygonList[i] ) );
2433  shell = createGeosPolygon( shellRing );
2434 
2435  if ( !GEOSWithin_r( geosinit.ctxt, newRingPolygon, shell ) )
2436  {
2437  GEOSGeom_destroy_r( geosinit.ctxt, shell );
2438  continue;
2439  }
2440  }
2441  catch ( GEOSException &e )
2442  {
2443  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2444 
2445  if ( shell )
2446  GEOSGeom_destroy_r( geosinit.ctxt, shell );
2447  else if ( shellRing )
2448  GEOSGeom_destroy_r( geosinit.ctxt, shellRing );
2449 
2450  GEOSGeom_destroy_r( geosinit.ctxt, newRingPolygon );
2451 
2452  return 4;
2453  }
2454 
2455  // add outer ring
2456  rings << GEOSGeom_clone_r( geosinit.ctxt, shellRing );
2457 
2458  GEOSGeom_destroy_r( geosinit.ctxt, shell );
2459 
2460  // check inner rings
2461  int n = GEOSGetNumInteriorRings_r( geosinit.ctxt, polygonList[i] );
2462 
2463  int j;
2464  for ( j = 0; j < n; j++ )
2465  {
2466  GEOSGeometry *holeRing = 0;
2467  GEOSGeometry *hole = 0;
2468  try
2469  {
2470  holeRing = GEOSGeom_clone_r( geosinit.ctxt, GEOSGetInteriorRingN_r( geosinit.ctxt, polygonList[i], j ) );
2471  hole = createGeosPolygon( holeRing );
2472 
2473  if ( !GEOSDisjoint_r( geosinit.ctxt, hole, newRingPolygon ) )
2474  {
2475  GEOSGeom_destroy_r( geosinit.ctxt, hole );
2476  break;
2477  }
2478  }
2479  catch ( GEOSException &e )
2480  {
2481  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2482 
2483  if ( hole )
2484  GEOSGeom_destroy_r( geosinit.ctxt, hole );
2485  else if ( holeRing )
2486  GEOSGeom_destroy_r( geosinit.ctxt, holeRing );
2487 
2488  break;
2489  }
2490 
2491  rings << GEOSGeom_clone_r( geosinit.ctxt, holeRing );
2492  GEOSGeom_destroy_r( geosinit.ctxt, hole );
2493  }
2494 
2495  if ( j == n )
2496  // this is it...
2497  break;
2498  }
2499 
2500  if ( i == polygonList.size() )
2501  {
2502  // clear rings
2503  for ( int j = 0; j < rings.size(); j++ )
2504  GEOSGeom_destroy_r( geosinit.ctxt, rings[j] );
2505  rings.clear();
2506 
2507  GEOSGeom_destroy_r( geosinit.ctxt, newRingPolygon );
2508 
2509  // no containing polygon found
2510  return 5;
2511  }
2512 
2513  rings << GEOSGeom_clone_r( geosinit.ctxt, newRing );
2514  GEOSGeom_destroy_r( geosinit.ctxt, newRingPolygon );
2515 
2516  GEOSGeometry *newPolygon = createGeosPolygon( rings );
2517 
2518  if ( wkbType() == QGis::WKBPolygon )
2519  {
2520  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
2521  mGeos = newPolygon;
2522  }
2523  else if ( wkbType() == QGis::WKBMultiPolygon )
2524  {
2525  QVector<GEOSGeometry*> newPolygons;
2526 
2527  for ( int j = 0; j < polygonList.size(); j++ )
2528  {
2529  newPolygons << ( i == j ? newPolygon : GEOSGeom_clone_r( geosinit.ctxt, polygonList[j] ) );
2530  }
2531 
2532  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
2533  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, newPolygons );
2534  }
2535 
2536  mDirtyWkb = true;
2537  mDirtyGeos = false;
2538  return 0;
2539 }
2540 
2541 int QgsGeometry::addPart( const QList<QgsPoint> &points, QGis::GeometryType geomType )
2542 {
2543  if ( geomType == QGis::UnknownGeometry )
2544  {
2545  geomType = type();
2546  }
2547 
2548  switch ( geomType )
2549  {
2550  case QGis::Point:
2551  // only one part at a time
2552  if ( points.size() != 1 )
2553  {
2554  QgsDebugMsg( "expected 1 point: " + QString::number( points.size() ) );
2555  return 2;
2556  }
2557  break;
2558 
2559  case QGis::Line:
2560  // line needs to have at least two points
2561  if ( points.size() < 2 )
2562  {
2563  QgsDebugMsg( "line must at least have two points: " + QString::number( points.size() ) );
2564  return 2;
2565  }
2566  break;
2567 
2568  case QGis::Polygon:
2569  // polygon needs to have at least three distinct points and must be closed
2570  if ( points.size() < 4 )
2571  {
2572  QgsDebugMsg( "polygon must at least have three distinct points and must be closed: " + QString::number( points.size() ) );
2573  return 2;
2574  }
2575 
2576  // Polygon must be closed
2577  if ( points.first() != points.last() )
2578  {
2579  QgsDebugMsg( "polygon not closed" );
2580  return 2;
2581  }
2582  break;
2583 
2584  default:
2585  QgsDebugMsg( "unsupported geometry type: " + QString::number( geomType ) );
2586  return 2;
2587  }
2588 
2589  GEOSGeometry *newPart = 0;
2590 
2591  switch ( geomType )
2592  {
2593  case QGis::Point:
2594  newPart = createGeosPoint( points[0] );
2595  break;
2596 
2597  case QGis::Line:
2598  newPart = createGeosLineString( points.toVector() );
2599  break;
2600 
2601  case QGis::Polygon:
2602  {
2603  //create new polygon from ring
2604  GEOSGeometry *newRing = 0;
2605 
2606  try
2607  {
2608  newRing = createGeosLinearRing( points.toVector() );
2609  if ( !GEOSisValid_r( geosinit.ctxt, newRing ) )
2610  throw GEOSException( "ring invalid" );
2611 
2612  newPart = createGeosPolygon( newRing );
2613  }
2614  catch ( GEOSException &e )
2615  {
2616  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2617 
2618  if ( newRing )
2619  GEOSGeom_destroy_r( geosinit.ctxt, newRing );
2620 
2621  return 2;
2622  }
2623  }
2624  break;
2625 
2626  default:
2627  QgsDebugMsg( "unsupported type: " + QString::number( type() ) );
2628  return 2;
2629  }
2630 
2631  if ( type() == QGis::UnknownGeometry )
2632  {
2633  fromGeos( newPart );
2634  return 0;
2635  }
2636  return addPart( newPart );
2637 }
2638 
2640 {
2641  if ( !newPart )
2642  return 4;
2643 
2644  const GEOSGeometry * geosPart = newPart->asGeos();
2645  return addPart( GEOSGeom_clone_r( geosinit.ctxt, geosPart ) );
2646 }
2647 
2648 int QgsGeometry::addPart( GEOSGeometry *newPart )
2649 {
2650  QGis::GeometryType geomType = type();
2651 
2652  if ( !isMultipart() && !convertToMultiType() )
2653  {
2654  QgsDebugMsg( "could not convert to multipart" );
2655  return 1;
2656  }
2657 
2658  //create geos geometry from wkb if not already there
2659  if ( mDirtyGeos )
2660  {
2661  exportWkbToGeos();
2662  }
2663 
2664  if ( !mGeos )
2665  {
2666  QgsDebugMsg( "GEOS geometry not available!" );
2667  return 4;
2668  }
2669 
2670  int geosType = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
2671 
2672  Q_ASSERT( newPart );
2673 
2674  try
2675  {
2676  if ( !GEOSisValid_r( geosinit.ctxt, newPart ) )
2677  throw GEOSException( "new part geometry invalid" );
2678  }
2679  catch ( GEOSException &e )
2680  {
2681  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
2682 
2683  if ( newPart )
2684  GEOSGeom_destroy_r( geosinit.ctxt, newPart );
2685 
2686  QgsDebugMsg( "part invalid: " + e.what() );
2687  return 2;
2688  }
2689 
2690  QVector<GEOSGeometry*> parts;
2691 
2692  //create new multipolygon
2693  int n = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
2694  int i;
2695  for ( i = 0; i < n; ++i )
2696  {
2697  const GEOSGeometry *partN = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
2698 
2699  if ( geomType == QGis::Polygon && GEOSOverlaps_r( geosinit.ctxt, partN, newPart ) )
2700  //bail out if new polygon overlaps with existing ones
2701  break;
2702 
2703  parts << GEOSGeom_clone_r( geosinit.ctxt, partN );
2704  }
2705 
2706  if ( i < n )
2707  {
2708  // bailed out
2709  for ( int i = 0; i < parts.size(); i++ )
2710  GEOSGeom_destroy_r( geosinit.ctxt, parts[i] );
2711 
2712  QgsDebugMsg( "new polygon part overlaps" );
2713  return 3;
2714  }
2715 
2716  int nPartGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, newPart );
2717  for ( int i = 0; i < nPartGeoms; ++i )
2718  {
2719  parts << GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, newPart, i ) );
2720  }
2721  GEOSGeom_destroy_r( geosinit.ctxt, newPart );
2722 
2723  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
2724 
2725  mGeos = createGeosCollection( geosType, parts );
2726 
2727  mDirtyWkb = true;
2728  mDirtyGeos = false;
2729 
2730  return 0;
2731 }
2732 
2733 int QgsGeometry::transform( const QTransform& t )
2734 {
2735  if ( mDirtyWkb )
2736  exportGeosToWkb();
2737 
2738  if ( !mGeometry )
2739  {
2740  QgsDebugMsg( "WKB geometry not available!" );
2741  return 1;
2742  }
2743 
2744  bool hasZValue = false;
2745  QgsWkbPtr wkbPtr( mGeometry + 1 );
2747  wkbPtr >> wkbType;
2748 
2749  switch ( wkbType )
2750  {
2751  case QGis::WKBPoint25D:
2752  case QGis::WKBPoint:
2753  {
2754  transformVertex( wkbPtr, t, hasZValue );
2755  }
2756  break;
2757 
2759  hasZValue = true;
2760  //intentional fall-through
2761  case QGis::WKBLineString:
2762  {
2763  int nPoints;
2764  wkbPtr >> nPoints;
2765  for ( int index = 0; index < nPoints; ++index )
2766  transformVertex( wkbPtr, t, hasZValue );
2767 
2768  break;
2769  }
2770 
2771  case QGis::WKBPolygon25D:
2772  hasZValue = true;
2773  //intentional fall-through
2774  case QGis::WKBPolygon:
2775  {
2776  int nRings;
2777  wkbPtr >> nRings;
2778  for ( int index = 0; index < nRings; ++index )
2779  {
2780  int nPoints;
2781  wkbPtr >> nPoints;
2782  for ( int index2 = 0; index2 < nPoints; ++index2 )
2783  transformVertex( wkbPtr, t, hasZValue );
2784 
2785  }
2786  break;
2787  }
2788 
2790  hasZValue = true;
2791  //intentional fall-through
2792  case QGis::WKBMultiPoint:
2793  {
2794  int nPoints;
2795  wkbPtr >> nPoints;
2796  for ( int index = 0; index < nPoints; ++index )
2797  {
2798  wkbPtr += 1 + sizeof( int );
2799  transformVertex( wkbPtr, t, hasZValue );
2800  }
2801  break;
2802  }
2803 
2805  hasZValue = true;
2806  //intentional fall-through
2808  {
2809  int nLines;
2810  wkbPtr >> nLines;
2811  for ( int index = 0; index < nLines; ++index )
2812  {
2813  wkbPtr += 1 + sizeof( int );
2814  int nPoints;
2815  wkbPtr >> nPoints;
2816  for ( int index2 = 0; index2 < nPoints; ++index2 )
2817  transformVertex( wkbPtr, t, hasZValue );
2818 
2819  }
2820  break;
2821  }
2822 
2824  hasZValue = true;
2825  //intentional fall-through
2826  case QGis::WKBMultiPolygon:
2827  {
2828  int nPolys;
2829  wkbPtr >> nPolys;
2830  for ( int index = 0; index < nPolys; ++index )
2831  {
2832  wkbPtr += 1 + sizeof( int ); //skip endian and polygon type
2833  int nRings;
2834  wkbPtr >> nRings;
2835  for ( int index2 = 0; index2 < nRings; ++index2 )
2836  {
2837  int nPoints;
2838  wkbPtr >> nPoints;
2839  for ( int index3 = 0; index3 < nPoints; ++index3 )
2840  transformVertex( wkbPtr, t, hasZValue );
2841 
2842  }
2843  }
2844  }
2845 
2846  default:
2847  break;
2848  }
2849  mDirtyGeos = true;
2850  return 0;
2851 }
2852 
2853 int QgsGeometry::translate( double dx, double dy )
2854 {
2855  return transform( QTransform::fromTranslate( dx, dy ) );
2856 }
2857 
2858 int QgsGeometry::rotate( double rotation, const QgsPoint& center )
2859 {
2860  QTransform t = QTransform::fromTranslate( center.x(), center.y() );
2861  t.rotate( -rotation );
2862  t.translate( -center.x(), -center.y() );
2863  return transform( t );
2864 }
2865 
2867 {
2868  if ( mDirtyWkb )
2869  exportGeosToWkb();
2870 
2871  if ( !mGeometry )
2872  {
2873  QgsDebugMsg( "WKB geometry not available!" );
2874  return 1;
2875  }
2876 
2877  bool hasZValue = false;
2878  QgsWkbPtr wkbPtr( mGeometry + 1 );
2880  wkbPtr >> wkbType;
2881 
2882  switch ( wkbType )
2883  {
2884  case QGis::WKBPoint25D:
2885  case QGis::WKBPoint:
2886  {
2887  transformVertex( wkbPtr, ct, hasZValue );
2888  }
2889  break;
2890 
2892  hasZValue = true;
2893  //intentional fall-through
2894  case QGis::WKBLineString:
2895  {
2896  int nPoints;
2897  wkbPtr >> nPoints;
2898  for ( int index = 0; index < nPoints; ++index )
2899  transformVertex( wkbPtr, ct, hasZValue );
2900 
2901  break;
2902  }
2903 
2904  case QGis::WKBPolygon25D:
2905  hasZValue = true;
2906  //intentional fall-through
2907  case QGis::WKBPolygon:
2908  {
2909  int nRings;
2910  wkbPtr >> nRings;
2911  for ( int index = 0; index < nRings; ++index )
2912  {
2913  int nPoints;
2914  wkbPtr >> nPoints;
2915  for ( int index2 = 0; index2 < nPoints; ++index2 )
2916  transformVertex( wkbPtr, ct, hasZValue );
2917 
2918  }
2919  break;
2920  }
2921 
2923  hasZValue = true;
2924  //intentional fall-through
2925  case QGis::WKBMultiPoint:
2926  {
2927  int nPoints;
2928  wkbPtr >> nPoints;
2929  for ( int index = 0; index < nPoints; ++index )
2930  {
2931  wkbPtr += 1 + sizeof( int );
2932  transformVertex( wkbPtr, ct, hasZValue );
2933  }
2934  break;
2935  }
2936 
2938  hasZValue = true;
2939  //intentional fall-through
2941  {
2942  int nLines;
2943  wkbPtr >> nLines;
2944  for ( int index = 0; index < nLines; ++index )
2945  {
2946  wkbPtr += 1 + sizeof( int );
2947  int nPoints;
2948  wkbPtr >> nPoints;
2949  for ( int index2 = 0; index2 < nPoints; ++index2 )
2950  transformVertex( wkbPtr, ct, hasZValue );
2951 
2952  }
2953  break;
2954  }
2955 
2957  hasZValue = true;
2958  //intentional fall-through
2959  case QGis::WKBMultiPolygon:
2960  {
2961  int nPolys;
2962  wkbPtr >> nPolys;
2963  for ( int index = 0; index < nPolys; ++index )
2964  {
2965  wkbPtr += 1 + sizeof( int ); //skip endian and polygon type
2966  int nRings;
2967  wkbPtr >> nRings;
2968  for ( int index2 = 0; index2 < nRings; ++index2 )
2969  {
2970  int nPoints;
2971  wkbPtr >> nPoints;
2972  for ( int index3 = 0; index3 < nPoints; ++index3 )
2973  transformVertex( wkbPtr, ct, hasZValue );
2974 
2975  }
2976  }
2977  }
2978 
2979  default:
2980  break;
2981  }
2982  mDirtyGeos = true;
2983  return 0;
2984 }
2985 
2986 int QgsGeometry::splitGeometry( const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint> &topologyTestPoints )
2987 {
2988  int returnCode = 0;
2989 
2990  //return if this type is point/multipoint
2991  if ( type() == QGis::Point )
2992  {
2993  return 1; //cannot split points
2994  }
2995 
2996  //make sure, mGeos and mWkb are there and up-to-date
2997  if ( mDirtyWkb )
2998  exportGeosToWkb();
2999 
3000  if ( mDirtyGeos )
3001  exportWkbToGeos();
3002 
3003  if ( !mGeos )
3004  return 1;
3005 
3006  if ( !GEOSisValid_r( geosinit.ctxt, mGeos ) )
3007  return 7;
3008 
3009  //make sure splitLine is valid
3010  if (( type() == QGis::Line && splitLine.size() < 1 ) ||
3011  ( type() == QGis::Polygon && splitLine.size() < 2 ) )
3012  return 1;
3013 
3014  newGeometries.clear();
3015 
3016  try
3017  {
3018  GEOSGeometry* splitLineGeos;
3019  if ( splitLine.size() > 1 )
3020  {
3021  splitLineGeos = createGeosLineString( splitLine.toVector() );
3022  }
3023  else if ( splitLine.size() == 1 )
3024  {
3025  splitLineGeos = createGeosPoint( splitLine.at( 0 ) );
3026  }
3027  else
3028  {
3029  return 1;
3030  }
3031  if ( !GEOSisValid_r( geosinit.ctxt, splitLineGeos ) || !GEOSisSimple_r( geosinit.ctxt, splitLineGeos ) )
3032  {
3033  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
3034  return 1;
3035  }
3036 
3037  if ( topological )
3038  {
3039  //find out candidate points for topological corrections
3040  if ( topologicalTestPointsSplit( splitLineGeos, topologyTestPoints ) != 0 )
3041  return 1;
3042  }
3043 
3044  //call split function depending on geometry type
3045  if ( type() == QGis::Line )
3046  {
3047  returnCode = splitLinearGeometry( splitLineGeos, newGeometries );
3048  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
3049  }
3050  else if ( type() == QGis::Polygon )
3051  {
3052  returnCode = splitPolygonGeometry( splitLineGeos, newGeometries );
3053  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
3054  }
3055  else
3056  {
3057  return 1;
3058  }
3059  }
3060  CATCH_GEOS( 2 )
3061 
3062  return returnCode;
3063 }
3064 
3066 int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine )
3067 {
3068  if ( reshapeWithLine.size() < 2 )
3069  return 1;
3070 
3071  if ( type() == QGis::Point )
3072  return 1; //cannot reshape points
3073 
3074  GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() );
3075 
3076  //make sure this geos geometry is up-to-date
3077  if ( mDirtyGeos )
3078  exportWkbToGeos();
3079 
3080  if ( !mGeos )
3081  return 1;
3082 
3083  //single or multi?
3084  int numGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
3085  if ( numGeoms == -1 )
3086  return 1;
3087 
3088  bool isMultiGeom = false;
3089  int geosTypeId = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
3090  if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
3091  isMultiGeom = true;
3092 
3093  bool isLine = ( type() == QGis::Line );
3094 
3095  //polygon or multipolygon?
3096  if ( !isMultiGeom )
3097  {
3098  GEOSGeometry* reshapedGeometry;
3099  if ( isLine )
3100  reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos );
3101  else
3102  reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos );
3103 
3104  GEOSGeom_destroy_r( geosinit.ctxt, reshapeLineGeos );
3105  if ( reshapedGeometry )
3106  {
3107  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
3108  mGeos = reshapedGeometry;
3109  mDirtyWkb = true;
3110  return 0;
3111  }
3112  else
3113  {
3114  return 1;
3115  }
3116  }
3117  else
3118  {
3119  //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
3120  bool reshapeTookPlace = false;
3121 
3122  GEOSGeometry* currentReshapeGeometry = 0;
3123  GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms];
3124 
3125  for ( int i = 0; i < numGeoms; ++i )
3126  {
3127  if ( isLine )
3128  currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ), reshapeLineGeos );
3129  else
3130  currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ), reshapeLineGeos );
3131 
3132  if ( currentReshapeGeometry )
3133  {
3134  newGeoms[i] = currentReshapeGeometry;
3135  reshapeTookPlace = true;
3136  }
3137  else
3138  {
3139  newGeoms[i] = GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) );
3140  }
3141  }
3142  GEOSGeom_destroy_r( geosinit.ctxt, reshapeLineGeos );
3143 
3144  GEOSGeometry* newMultiGeom = 0;
3145  if ( isLine )
3146  {
3147  newMultiGeom = GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, newGeoms, numGeoms );
3148  }
3149  else //multipolygon
3150  {
3151  newMultiGeom = GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTIPOLYGON, newGeoms, numGeoms );
3152  }
3153 
3154  delete[] newGeoms;
3155  if ( !newMultiGeom )
3156  return 3;
3157 
3158  if ( reshapeTookPlace )
3159  {
3160  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
3161  mGeos = newMultiGeom;
3162  mDirtyWkb = true;
3163  return 0;
3164  }
3165  else
3166  {
3167  GEOSGeom_destroy_r( geosinit.ctxt, newMultiGeom );
3168  return 1;
3169  }
3170  }
3171 }
3172 
3174 {
3175  //make sure geos geometry is up to date
3176  if ( !other )
3177  return 1;
3178 
3179  if ( mDirtyGeos )
3180  exportWkbToGeos();
3181 
3182  if ( !mGeos )
3183  return 1;
3184 
3185  if ( !GEOSisValid_r( geosinit.ctxt, mGeos ) )
3186  return 2;
3187 
3188  if ( !GEOSisSimple_r( geosinit.ctxt, mGeos ) )
3189  return 3;
3190 
3191  //convert other geometry to geos
3192  if ( other->mDirtyGeos )
3193  other->exportWkbToGeos();
3194 
3195  if ( !other->mGeos )
3196  return 4;
3197 
3198  //make geometry::difference
3199  try
3200  {
3201  if ( GEOSIntersects_r( geosinit.ctxt, mGeos, other->mGeos ) )
3202  {
3203  //check if multitype before and after
3204  bool multiType = isMultipart();
3205 
3206  mGeos = GEOSDifference_r( geosinit.ctxt, mGeos, other->mGeos );
3207  mDirtyWkb = true;
3208 
3209  if ( multiType && !isMultipart() )
3210  {
3212  exportWkbToGeos();
3213  }
3214  }
3215  else
3216  {
3217  return 0; //nothing to do
3218  }
3219  }
3220  CATCH_GEOS( 5 )
3221 
3222  if ( !mGeos )
3223  {
3224  mDirtyGeos = true;
3225  return 6;
3226  }
3227 
3228  return 0;
3229 }
3230 
3232 {
3233  double xmin = std::numeric_limits<double>::max();
3234  double ymin = std::numeric_limits<double>::max();
3235  double xmax = -std::numeric_limits<double>::max();
3236  double ymax = -std::numeric_limits<double>::max();
3237 
3238  // TODO: implement with GEOS
3239  if ( mDirtyWkb )
3240  exportGeosToWkb();
3241 
3242  if ( !mGeometry )
3243  {
3244  QgsDebugMsg( "WKB geometry not available!" );
3245  // Return minimal QgsRectangle
3246  QgsRectangle invalidRect;
3247  invalidRect.setMinimal();
3248  return invalidRect;
3249  }
3250 
3251  bool hasZValue = false;
3252  QgsWkbPtr wkbPtr( mGeometry + 1 );
3254  wkbPtr >> wkbType;
3255 
3256  // consider endian when fetching feature type
3257  switch ( wkbType )
3258  {
3259  case QGis::WKBPoint25D:
3260  case QGis::WKBPoint:
3261  {
3262  double x, y;
3263  wkbPtr >> x >> y;
3264 
3265  if ( x < xmin )
3266  xmin = x;
3267 
3268  if ( x > xmax )
3269  xmax = x;
3270 
3271  if ( y < ymin )
3272  ymin = y;
3273 
3274  if ( y > ymax )
3275  ymax = y;
3276 
3277  }
3278  break;
3279 
3281  hasZValue = true;
3282  //intentional fall-through
3283  case QGis::WKBMultiPoint:
3284  {
3285  int nPoints;
3286  wkbPtr >> nPoints;
3287  for ( int idx = 0; idx < nPoints; idx++ )
3288  {
3289  wkbPtr += 1 + sizeof( int );
3290 
3291  double x, y;
3292  wkbPtr >> x >> y;
3293  if ( hasZValue )
3294  wkbPtr += sizeof( double );
3295 
3296  if ( x < xmin )
3297  xmin = x;
3298 
3299  if ( x > xmax )
3300  xmax = x;
3301 
3302  if ( y < ymin )
3303  ymin = y;
3304 
3305  if ( y > ymax )
3306  ymax = y;
3307 
3308  }
3309  break;
3310  }
3312  hasZValue = true;
3313  //intentional fall-through
3314  case QGis::WKBLineString:
3315  {
3316  // get number of points in the line
3317  int nPoints;
3318  wkbPtr >> nPoints;
3319  for ( int idx = 0; idx < nPoints; idx++ )
3320  {
3321  double x, y;
3322  wkbPtr >> x >> y;
3323  if ( hasZValue )
3324  wkbPtr += sizeof( double );
3325 
3326  if ( x < xmin )
3327  xmin = x;
3328 
3329  if ( x > xmax )
3330  xmax = x;
3331 
3332  if ( y < ymin )
3333  ymin = y;
3334 
3335  if ( y > ymax )
3336  ymax = y;
3337 
3338  }
3339  break;
3340  }
3342  hasZValue = true;
3343  //intentional fall-through
3345  {
3346  int nLines;
3347  wkbPtr >> nLines;
3348  for ( int jdx = 0; jdx < nLines; jdx++ )
3349  {
3350  // each of these is a wbklinestring so must handle as such
3351  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3352  int nPoints;
3353  wkbPtr >> nPoints;
3354  for ( int idx = 0; idx < nPoints; idx++ )
3355  {
3356  double x, y;
3357  wkbPtr >> x >> y;
3358  if ( hasZValue )
3359  wkbPtr += sizeof( double );
3360 
3361  if ( x < xmin )
3362  xmin = x;
3363 
3364  if ( x > xmax )
3365  xmax = x;
3366 
3367  if ( y < ymin )
3368  ymin = y;
3369 
3370  if ( y > ymax )
3371  ymax = y;
3372 
3373  }
3374  }
3375  break;
3376  }
3377  case QGis::WKBPolygon25D:
3378  hasZValue = true;
3379  //intentional fall-through
3380  case QGis::WKBPolygon:
3381  {
3382  // get number of rings in the polygon
3383  int nRings;
3384  wkbPtr >> nRings;
3385  for ( int idx = 0; idx < nRings; idx++ )
3386  {
3387  // get number of points in the ring
3388  int nPoints;
3389  wkbPtr >> nPoints;
3390  for ( int jdx = 0; jdx < nPoints; jdx++ )
3391  {
3392  // add points to a point array for drawing the polygon
3393  double x, y;
3394  wkbPtr >> x >> y;
3395  if ( hasZValue )
3396  wkbPtr += sizeof( double );
3397 
3398  if ( x < xmin )
3399  xmin = x;
3400 
3401  if ( x > xmax )
3402  xmax = x;
3403 
3404  if ( y < ymin )
3405  ymin = y;
3406 
3407  if ( y > ymax )
3408  ymax = y;
3409 
3410  }
3411  }
3412  break;
3413  }
3415  hasZValue = true;
3416  //intentional fall-through
3417  case QGis::WKBMultiPolygon:
3418  {
3419  // get the number of polygons
3420  int nPolygons;
3421  wkbPtr >> nPolygons;
3422  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3423  {
3424  //skip the endian and mGeometry type info and
3425  // get number of rings in the polygon
3426  wkbPtr += 1 + sizeof( int );
3427  int nRings;
3428  wkbPtr >> nRings;
3429  for ( int idx = 0; idx < nRings; idx++ )
3430  {
3431  // get number of points in the ring
3432  int nPoints;
3433  wkbPtr >> nPoints;
3434  for ( int jdx = 0; jdx < nPoints; jdx++ )
3435  {
3436  // add points to a point array for drawing the polygon
3437  double x, y;
3438  wkbPtr >> x >> y;
3439  if ( hasZValue )
3440  wkbPtr += sizeof( double );
3441 
3442  if ( x < xmin )
3443  xmin = x;
3444 
3445  if ( x > xmax )
3446  xmax = x;
3447 
3448  if ( y < ymin )
3449  ymin = y;
3450 
3451  if ( y > ymax )
3452  ymax = y;
3453 
3454  }
3455  }
3456  }
3457  break;
3458  }
3459 
3460  default:
3461  QgsDebugMsg( QString( "Unknown WkbType %1 ENCOUNTERED" ).arg( wkbType ) );
3462  return QgsRectangle( 0, 0, 0, 0 );
3463  break;
3464 
3465  }
3466  return QgsRectangle( xmin, ymin, xmax, ymax );
3467 }
3468 
3470 {
3471  QgsGeometry* g = fromRect( r );
3472  bool res = intersects( g );
3473  delete g;
3474  return res;
3475 }
3476 
3477 bool QgsGeometry::intersects( const QgsGeometry* geometry ) const
3478 {
3479  if ( !geometry )
3480  return false;
3481 
3482  try // geos might throw exception on error
3483  {
3484  // ensure that both geometries have geos geometry
3485  exportWkbToGeos();
3486  geometry->exportWkbToGeos();
3487 
3488  if ( !mGeos || !geometry->mGeos )
3489  {
3490  QgsDebugMsg( "GEOS geometry not available!" );
3491  return false;
3492  }
3493 
3494  return GEOSIntersects_r( geosinit.ctxt, mGeos, geometry->mGeos );
3495  }
3496  CATCH_GEOS( false )
3497 }
3498 
3499 bool QgsGeometry::contains( const QgsPoint* p ) const
3500 {
3501  exportWkbToGeos();
3502 
3503  if ( !p )
3504  {
3505  QgsDebugMsg( "pointer p is 0" );
3506  return false;
3507  }
3508 
3509  if ( !mGeos )
3510  {
3511  QgsDebugMsg( "GEOS geometry not available!" );
3512  return false;
3513  }
3514 
3515  GEOSGeometry *geosPoint = 0;
3516 
3517  bool returnval = false;
3518 
3519  try
3520  {
3521  geosPoint = createGeosPoint( *p );
3522  returnval = GEOSContains_r( geosinit.ctxt, mGeos, geosPoint );
3523  }
3524  catch ( GEOSException &e )
3525  {
3526  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
3527  returnval = false;
3528  }
3529 
3530  if ( geosPoint )
3531  GEOSGeom_destroy_r( geosinit.ctxt, geosPoint );
3532 
3533  return returnval;
3534 }
3535 
3536 bool QgsGeometry::geosRelOp(
3537  char( *op )( GEOSContextHandle_t handle, const GEOSGeometry*, const GEOSGeometry * ),
3538  const QgsGeometry *a,
3539  const QgsGeometry *b )
3540 {
3541  if ( !a || !b )
3542  return false;
3543 
3544  try // geos might throw exception on error
3545  {
3546  // ensure that both geometries have geos geometry
3547  a->exportWkbToGeos();
3548  b->exportWkbToGeos();
3549 
3550  if ( !a->mGeos || !b->mGeos )
3551  {
3552  QgsDebugMsg( "GEOS geometry not available!" );
3553  return false;
3554  }
3555  return op( geosinit.ctxt, a->mGeos, b->mGeos );
3556  }
3557  CATCH_GEOS( false )
3558 }
3559 
3560 bool QgsGeometry::contains( const QgsGeometry* geometry ) const
3561 {
3562  return geosRelOp( GEOSContains_r, this, geometry );
3563 }
3564 
3565 bool QgsGeometry::disjoint( const QgsGeometry* geometry ) const
3566 {
3567  return geosRelOp( GEOSDisjoint_r, this, geometry );
3568 }
3569 
3570 bool QgsGeometry::equals( const QgsGeometry* geometry ) const
3571 {
3572  return geosRelOp( GEOSEquals_r, this, geometry );
3573 }
3574 
3575 bool QgsGeometry::touches( const QgsGeometry* geometry ) const
3576 {
3577  return geosRelOp( GEOSTouches_r, this, geometry );
3578 }
3579 
3580 bool QgsGeometry::overlaps( const QgsGeometry* geometry ) const
3581 {
3582  return geosRelOp( GEOSOverlaps_r, this, geometry );
3583 }
3584 
3585 bool QgsGeometry::within( const QgsGeometry* geometry ) const
3586 {
3587  return geosRelOp( GEOSWithin_r, this, geometry );
3588 }
3589 
3590 bool QgsGeometry::crosses( const QgsGeometry* geometry ) const
3591 {
3592  return geosRelOp( GEOSCrosses_r, this, geometry );
3593 }
3594 
3595 QString QgsGeometry::exportToWkt( const int &precision ) const
3596 {
3597  QgsDebugMsg( "entered." );
3598 
3599  // TODO: implement with GEOS
3600  if ( mDirtyWkb )
3601  {
3602  exportGeosToWkb();
3603  }
3604 
3605  if ( !mGeometry || wkbSize() < 5 )
3606  {
3607  QgsDebugMsg( "WKB geometry not available or too short!" );
3608  return QString::null;
3609  }
3610 
3611  bool hasZValue = false;
3612  QgsWkbPtr wkbPtr( mGeometry + 1 );
3614  wkbPtr >> wkbType;
3615 
3616  QString wkt;
3617 
3618  switch ( wkbType )
3619  {
3620  case QGis::WKBPoint25D:
3621  case QGis::WKBPoint:
3622  {
3623  double x, y;
3624  wkbPtr >> x >> y;
3625  wkt += "POINT(" + qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision ) + ")";
3626  return wkt;
3627  }
3628 
3630  hasZValue = true;
3631  //intentional fall-through
3632  case QGis::WKBLineString:
3633  {
3634  int nPoints;
3635  wkbPtr >> nPoints;
3636 
3637  wkt += "LINESTRING(";
3638  // get number of points in the line
3639  for ( int idx = 0; idx < nPoints; ++idx )
3640  {
3641  double x, y;
3642  wkbPtr >> x >> y;
3643  if ( hasZValue )
3644  wkbPtr += sizeof( double );
3645 
3646  if ( idx != 0 )
3647  wkt += ", ";
3648 
3649  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3650  }
3651  wkt += ")";
3652  return wkt;
3653  }
3654 
3655  case QGis::WKBPolygon25D:
3656  hasZValue = true;
3657  //intentional fall-through
3658  case QGis::WKBPolygon:
3659  {
3660  wkt += "POLYGON(";
3661  // get number of rings in the polygon
3662  int nRings;
3663  wkbPtr >> nRings;
3664  if ( nRings == 0 ) // sanity check for zero rings in polygon
3665  return QString();
3666 
3667  for ( int idx = 0; idx < nRings; idx++ )
3668  {
3669  if ( idx != 0 )
3670  wkt += ",";
3671 
3672  wkt += "(";
3673  // get number of points in the ring
3674  int nPoints;
3675  wkbPtr >> nPoints;
3676 
3677  for ( int jdx = 0; jdx < nPoints; jdx++ )
3678  {
3679  if ( jdx != 0 )
3680  wkt += ",";
3681 
3682  double x, y;
3683  wkbPtr >> x >> y;
3684  if ( hasZValue )
3685  wkbPtr += sizeof( double );
3686 
3687  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3688  }
3689  wkt += ")";
3690  }
3691  wkt += ")";
3692  return wkt;
3693  }
3694 
3696  hasZValue = true;
3697  //intentional fall-through
3698  case QGis::WKBMultiPoint:
3699  {
3700  int nPoints;
3701  wkbPtr >> nPoints;
3702 
3703  wkt += "MULTIPOINT(";
3704  for ( int idx = 0; idx < nPoints; ++idx )
3705  {
3706  wkbPtr += 1 + sizeof( int );
3707  if ( idx != 0 )
3708  wkt += ", ";
3709 
3710  double x, y;
3711  wkbPtr >> x >> y;
3712  if ( hasZValue )
3713  wkbPtr += sizeof( double );
3714 
3715  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3716  }
3717  wkt += ")";
3718  return wkt;
3719  }
3720 
3722  hasZValue = true;
3723  //intentional fall-through
3725  {
3726  int nLines;
3727  wkbPtr >> nLines;
3728 
3729  wkt += "MULTILINESTRING(";
3730  for ( int jdx = 0; jdx < nLines; jdx++ )
3731  {
3732  if ( jdx != 0 )
3733  wkt += ", ";
3734 
3735  wkt += "(";
3736  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3737  int nPoints;
3738  wkbPtr >> nPoints;
3739  for ( int idx = 0; idx < nPoints; idx++ )
3740  {
3741  if ( idx != 0 )
3742  wkt += ", ";
3743 
3744  double x, y;
3745  wkbPtr >> x >> y;
3746  if ( hasZValue )
3747  wkbPtr += sizeof( double );
3748 
3749  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3750  }
3751  wkt += ")";
3752  }
3753  wkt += ")";
3754  return wkt;
3755  }
3756 
3758  hasZValue = true;
3759  //intentional fall-through
3760  case QGis::WKBMultiPolygon:
3761  {
3762  int nPolygons;
3763  wkbPtr >> nPolygons;
3764 
3765  wkt += "MULTIPOLYGON(";
3766  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3767  {
3768  if ( kdx != 0 )
3769  wkt += ",";
3770 
3771  wkt += "(";
3772  wkbPtr += 1 + sizeof( int );
3773  int nRings;
3774  wkbPtr >> nRings;
3775  for ( int idx = 0; idx < nRings; idx++ )
3776  {
3777  if ( idx != 0 )
3778  wkt += ",";
3779 
3780  wkt += "(";
3781  int nPoints;
3782  wkbPtr >> nPoints;
3783  for ( int jdx = 0; jdx < nPoints; jdx++ )
3784  {
3785  if ( jdx != 0 )
3786  wkt += ",";
3787 
3788  double x, y;
3789  wkbPtr >> x >> y;
3790  if ( hasZValue )
3791  wkbPtr += sizeof( double );
3792 
3793  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3794  }
3795  wkt += ")";
3796  }
3797  wkt += ")";
3798  }
3799  wkt += ")";
3800  return wkt;
3801  }
3802 
3803  default:
3804  QgsDebugMsg( "error: mGeometry type not recognized" );
3805  return QString::null;
3806  }
3807 }
3808 
3809 QString QgsGeometry::exportToGeoJSON( const int &precision ) const
3810 {
3811  QgsDebugMsg( "entered." );
3812 
3813  // TODO: implement with GEOS
3814  if ( mDirtyWkb )
3815  exportGeosToWkb();
3816 
3817  if ( !mGeometry )
3818  {
3819  QgsDebugMsg( "WKB geometry not available!" );
3820  return QString::null;
3821  }
3822 
3823  QgsWkbPtr wkbPtr( mGeometry + 1 );
3825  wkbPtr >> wkbType;
3826  bool hasZValue = false;
3827 
3828  QString wkt;
3829 
3830  switch ( wkbType )
3831  {
3832  case QGis::WKBPoint25D:
3833  case QGis::WKBPoint:
3834  {
3835 
3836  double x, y;
3837  wkbPtr >> x >> y;
3838 
3839  wkt += "{ \"type\": \"Point\", \"coordinates\": ["
3840  + qgsDoubleToString( x, precision ) + ", "
3841  + qgsDoubleToString( y, precision )
3842  + "] }";
3843  return wkt;
3844  }
3845 
3847  hasZValue = true;
3848  //intentional fall-through
3849  case QGis::WKBLineString:
3850  {
3851 
3852  wkt += "{ \"type\": \"LineString\", \"coordinates\": [ ";
3853  // get number of points in the line
3854  int nPoints;
3855  wkbPtr >> nPoints;
3856  for ( int idx = 0; idx < nPoints; ++idx )
3857  {
3858  if ( idx != 0 )
3859  wkt += ", ";
3860 
3861  double x, y;
3862  wkbPtr >> x >> y;
3863  if ( hasZValue )
3864  wkbPtr += sizeof( double );
3865 
3866  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3867  }
3868  wkt += " ] }";
3869  return wkt;
3870  }
3871 
3872  case QGis::WKBPolygon25D:
3873  hasZValue = true;
3874  //intentional fall-through
3875  case QGis::WKBPolygon:
3876  {
3877 
3878  wkt += "{ \"type\": \"Polygon\", \"coordinates\": [ ";
3879  // get number of rings in the polygon
3880  int nRings;
3881  wkbPtr >> nRings;
3882  if ( nRings == 0 ) // sanity check for zero rings in polygon
3883  return QString();
3884 
3885  for ( int idx = 0; idx < nRings; idx++ )
3886  {
3887  if ( idx != 0 )
3888  wkt += ", ";
3889 
3890  wkt += "[ ";
3891  // get number of points in the ring
3892  int nPoints;
3893  wkbPtr >> nPoints;
3894 
3895  for ( int jdx = 0; jdx < nPoints; jdx++ )
3896  {
3897  if ( jdx != 0 )
3898  wkt += ", ";
3899 
3900  double x, y;
3901  wkbPtr >> x >> y;
3902  if ( hasZValue )
3903  wkbPtr += sizeof( double );
3904 
3905  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3906  }
3907  wkt += " ]";
3908  }
3909  wkt += " ] }";
3910  return wkt;
3911  }
3912 
3914  hasZValue = true;
3915  //intentional fall-through
3916  case QGis::WKBMultiPoint:
3917  {
3918  wkt += "{ \"type\": \"MultiPoint\", \"coordinates\": [ ";
3919  int nPoints;
3920  wkbPtr >> nPoints;
3921  for ( int idx = 0; idx < nPoints; ++idx )
3922  {
3923  wkbPtr += 1 + sizeof( int );
3924  if ( idx != 0 )
3925  wkt += ", ";
3926 
3927  double x, y;
3928  wkbPtr >> x >> y;
3929  if ( hasZValue )
3930  wkbPtr += sizeof( double );
3931 
3932  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3933  }
3934  wkt += " ] }";
3935  return wkt;
3936  }
3937 
3939  hasZValue = true;
3940  //intentional fall-through
3942  {
3943  wkt += "{ \"type\": \"MultiLineString\", \"coordinates\": [ ";
3944 
3945  int nLines;
3946  wkbPtr >> nLines;
3947  for ( int jdx = 0; jdx < nLines; jdx++ )
3948  {
3949  if ( jdx != 0 )
3950  wkt += ", ";
3951 
3952  wkt += "[ ";
3953  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3954 
3955  int nPoints;
3956  wkbPtr >> nPoints;
3957  for ( int idx = 0; idx < nPoints; idx++ )
3958  {
3959  if ( idx != 0 )
3960  wkt += ", ";
3961 
3962  double x, y;
3963  wkbPtr >> x >> y;
3964  if ( hasZValue )
3965  wkbPtr += sizeof( double );
3966 
3967  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3968  }
3969  wkt += " ]";
3970  }
3971  wkt += " ] }";
3972  return wkt;
3973  }
3974 
3976  hasZValue = true;
3977  //intentional fall-through
3978  case QGis::WKBMultiPolygon:
3979  {
3980 
3981  wkt += "{ \"type\": \"MultiPolygon\", \"coordinates\": [ ";
3982 
3983  int nPolygons;
3984  wkbPtr >> nPolygons;
3985  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3986  {
3987  if ( kdx != 0 )
3988  wkt += ", ";
3989 
3990  wkt += "[ ";
3991 
3992  wkbPtr += 1 + sizeof( int );
3993 
3994  int nRings;
3995  wkbPtr >> nRings;
3996  for ( int idx = 0; idx < nRings; idx++ )
3997  {
3998  if ( idx != 0 )
3999  wkt += ", ";
4000 
4001  wkt += "[ ";
4002 
4003  int nPoints;
4004  wkbPtr >> nPoints;
4005  for ( int jdx = 0; jdx < nPoints; jdx++ )
4006  {
4007  if ( jdx != 0 )
4008  wkt += ", ";
4009 
4010  double x, y;
4011  wkbPtr >> x >> y;
4012  if ( hasZValue )
4013  wkbPtr += sizeof( double );
4014 
4015  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
4016  }
4017  wkt += " ]";
4018  }
4019  wkt += " ]";
4020  }
4021  wkt += " ] }";
4022  return wkt;
4023  }
4024 
4025  default:
4026  QgsDebugMsg( "error: mGeometry type not recognized" );
4027  return QString::null;
4028  }
4029 }
4030 
4031 bool QgsGeometry::exportWkbToGeos() const
4032 {
4033  QgsDebugMsgLevel( "entered.", 3 );
4034 
4035  if ( !mDirtyGeos )
4036  {
4037  // No need to convert again
4038  return true;
4039  }
4040 
4041  if ( mGeos )
4042  {
4043  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
4044  mGeos = 0;
4045  }
4046 
4047  // this probably shouldn't return true
4048  if ( !mGeometry )
4049  {
4050  // no WKB => no GEOS
4051  mDirtyGeos = false;
4052  return true;
4053  }
4054 
4055  bool hasZValue = false;
4056 
4057  QgsWkbPtr wkbPtr( mGeometry + 1 );
4059  wkbPtr >> wkbType;
4060 
4061  try
4062  {
4063  switch ( wkbType )
4064  {
4065  case QGis::WKBPoint25D:
4066  case QGis::WKBPoint:
4067  {
4068  double x, y;
4069  wkbPtr >> x >> y;
4070  mGeos = createGeosPoint( QgsPoint( x, y ) );
4071  mDirtyGeos = false;
4072  break;
4073  }
4074 
4076  hasZValue = true;
4077  //intentional fall-through
4078  case QGis::WKBMultiPoint:
4079  {
4080  QVector<GEOSGeometry *> points;
4081 
4082  int nPoints;
4083  wkbPtr >> nPoints;
4084  for ( int idx = 0; idx < nPoints; idx++ )
4085  {
4086  double x, y;
4087  wkbPtr += 1 + sizeof( int );
4088  wkbPtr >> x >> y;
4089  if ( hasZValue )
4090  wkbPtr += sizeof( double );
4091 
4092  points << createGeosPoint( QgsPoint( x, y ) );
4093  }
4094  mGeos = createGeosCollection( GEOS_MULTIPOINT, points );
4095  mDirtyGeos = false;
4096  break;
4097  }
4098 
4100  hasZValue = true;
4101  //intentional fall-through
4102  case QGis::WKBLineString:
4103  {
4104  QgsPolyline sequence;
4105 
4106  int nPoints;
4107  wkbPtr >> nPoints;
4108  for ( int idx = 0; idx < nPoints; idx++ )
4109  {
4110  double x, y;
4111  wkbPtr >> x >> y;
4112  if ( hasZValue )
4113  wkbPtr += sizeof( double );
4114 
4115  sequence << QgsPoint( x, y );
4116  }
4117  mDirtyGeos = false;
4118  mGeos = createGeosLineString( sequence );
4119  break;
4120  }
4121 
4123  hasZValue = true;
4124  //intentional fall-through
4126  {
4127  QVector<GEOSGeometry*> lines;
4128  int nLines;
4129  wkbPtr >> nLines;
4130  for ( int jdx = 0; jdx < nLines; jdx++ )
4131  {
4132  QgsPolyline sequence;
4133 
4134  // each of these is a wbklinestring so must handle as such
4135  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
4136  int nPoints;
4137  wkbPtr >> nPoints;
4138  for ( int idx = 0; idx < nPoints; idx++ )
4139  {
4140  double x, y;
4141  wkbPtr >> x >> y;
4142  if ( hasZValue )
4143  wkbPtr += sizeof( double );
4144 
4145  sequence << QgsPoint( x, y );
4146  }
4147 
4148  // ignore invalid parts, it can come from ST_Simplify operations
4149  if ( sequence.count() > 1 )
4150  lines << createGeosLineString( sequence );
4151  }
4152  mGeos = createGeosCollection( GEOS_MULTILINESTRING, lines );
4153  mDirtyGeos = false;
4154  break;
4155  }
4156 
4157  case QGis::WKBPolygon25D:
4158  hasZValue = true;
4159  //intentional fall-through
4160  case QGis::WKBPolygon:
4161  {
4162  // get number of rings in the polygon
4163  int nRings;
4164  wkbPtr >> nRings;
4165 
4166  QVector<GEOSGeometry*> rings;
4167 
4168  for ( int idx = 0; idx < nRings; idx++ )
4169  {
4170  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4171 
4172  QgsPolyline sequence;
4173 
4174  // get number of points in the ring
4175  int nPoints;
4176  wkbPtr >> nPoints;
4177  for ( int jdx = 0; jdx < nPoints; jdx++ )
4178  {
4179  // add points to a point array for drawing the polygon
4180  double x, y;
4181  wkbPtr >> x >> y;
4182  if ( hasZValue )
4183  wkbPtr += sizeof( double );
4184 
4185  sequence << QgsPoint( x, y );
4186  }
4187 
4188  GEOSGeometry *ring = createGeosLinearRing( sequence );
4189  if ( ring ) rings << ring;
4190  }
4191  mGeos = createGeosPolygon( rings );
4192  mDirtyGeos = false;
4193  break;
4194  }
4195 
4197  hasZValue = true;
4198  //intentional fall-through
4199  case QGis::WKBMultiPolygon:
4200  {
4201  QVector<GEOSGeometry*> polygons;
4202 
4203  // get the number of polygons
4204  int nPolygons;
4205  wkbPtr >> nPolygons;
4206 
4207  for ( int kdx = 0; kdx < nPolygons; kdx++ )
4208  {
4209  //QgsDebugMsg("Polygon nr: "+QString::number(kdx));
4210  QVector<GEOSGeometry*> rings;
4211 
4212  //skip the endian and mGeometry type info and
4213  // get number of rings in the polygon
4214  wkbPtr += 1 + sizeof( int );
4215  int numRings;
4216  wkbPtr >> numRings;
4217  for ( int idx = 0; idx < numRings; idx++ )
4218  {
4219  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4220 
4221  QgsPolyline sequence;
4222 
4223  // get number of points in the ring
4224  int nPoints;
4225  wkbPtr >> nPoints;
4226  for ( int jdx = 0; jdx < nPoints; jdx++ )
4227  {
4228  // add points to a point array for drawing the polygon
4229  double x, y;
4230  wkbPtr >> x >> y;
4231  if ( hasZValue )
4232  wkbPtr += sizeof( double );
4233 
4234  sequence << QgsPoint( x, y );
4235  }
4236 
4237  GEOSGeometry *ring = createGeosLinearRing( sequence );
4238  if ( ring ) rings << ring;
4239  }
4240 
4241  GEOSGeometry *polygon = createGeosPolygon( rings );
4242  if ( polygon ) polygons << polygon;
4243  }
4244  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
4245  mDirtyGeos = false;
4246  break;
4247  }
4248 
4249  default:
4250  return false;
4251  }
4252  }
4253  CATCH_GEOS( false )
4254 
4255  return true;
4256 }
4257 
4258 bool QgsGeometry::exportGeosToWkb() const
4259 {
4260  //QgsDebugMsg("entered.");
4261  if ( !mDirtyWkb )
4262  {
4263  // No need to convert again
4264  return true;
4265  }
4266 
4267  // clear the WKB, ready to replace with the new one
4268  if ( mGeometry )
4269  {
4270  delete [] mGeometry;
4271  mGeometry = 0;
4272  }
4273 
4274  if ( !mGeos )
4275  {
4276  // GEOS is null, therefore WKB is null.
4277  mDirtyWkb = false;
4278  return true;
4279  }
4280 
4281  // set up byteOrder
4282  char byteOrder = QgsApplication::endian();
4283 
4284  switch ( GEOSGeomTypeId_r( geosinit.ctxt, mGeos ) )
4285  {
4286  case GEOS_POINT: // a point
4287  {
4288  mGeometrySize = 1 +
4289  sizeof( int ) +
4290  2 * sizeof( double );
4291  mGeometry = new unsigned char[mGeometrySize];
4292 
4293 
4294  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4295 
4296  double x, y;
4297  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, 0, &x );
4298  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, 0, &y );
4299 
4300  QgsWkbPtr wkbPtr( mGeometry );
4301  wkbPtr << byteOrder << QGis::WKBPoint << x << y;
4302 
4303  mDirtyWkb = false;
4304  return true;
4305  } // case GEOS_GEOM::GEOS_POINT
4306 
4307  case GEOS_LINESTRING: // a linestring
4308  {
4309  //QgsDebugMsg("Got a geos::GEOS_LINESTRING.");
4310 
4311  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4312  unsigned int nPoints;
4313  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &nPoints );
4314 
4315  // allocate some space for the WKB
4316  mGeometrySize = 1 + // sizeof(byte)
4317  sizeof( int ) +
4318  sizeof( int ) +
4319  (( sizeof( double ) +
4320  sizeof( double ) ) * nPoints );
4321 
4322  mGeometry = new unsigned char[mGeometrySize];
4323  QgsWkbPtr wkbPtr( mGeometry );
4324 
4325  wkbPtr << byteOrder << QGis::WKBLineString << nPoints;
4326 
4327  const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4328 
4329  // assign points
4330  for ( unsigned int n = 0; n < nPoints; n++ )
4331  {
4332  double x, y;
4333  GEOSCoordSeq_getX_r( geosinit.ctxt, sequence, n, &x );
4334  GEOSCoordSeq_getY_r( geosinit.ctxt, sequence, n, &y );
4335  wkbPtr << x << y;
4336  }
4337 
4338  mDirtyWkb = false;
4339  return true;
4340 
4341  // TODO: Deal with endian-ness
4342  } // case GEOS_GEOM::GEOS_LINESTRING
4343 
4344  case GEOS_LINEARRING: // a linear ring (linestring with 1st point == last point)
4345  {
4346  // TODO
4347  break;
4348  } // case GEOS_GEOM::GEOS_LINEARRING
4349 
4350  case GEOS_POLYGON: // a polygon
4351  {
4352  int nPointsInRing = 0;
4353 
4354  //first calculate the geometry size
4355  int geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of rings
4356  const GEOSGeometry *theRing = GEOSGetExteriorRing_r( geosinit.ctxt, mGeos );
4357  if ( theRing )
4358  {
4359  geometrySize += sizeof( int );
4360  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4361  }
4362  for ( int i = 0; i < GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ); ++i )
4363  {
4364  geometrySize += sizeof( int ); //number of points in ring
4365  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, mGeos, i );
4366  if ( theRing )
4367  {
4368  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4369  }
4370  }
4371 
4372  mGeometry = new unsigned char[geometrySize];
4373  mGeometrySize = geometrySize;
4374 
4375  //then fill the geometry itself into the wkb
4376  QgsWkbPtr wkbPtr( mGeometry );
4377 
4378  int nRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ) + 1;
4379  wkbPtr << byteOrder << QGis::WKBPolygon << nRings;
4380 
4381  //exterior ring first
4382  theRing = GEOSGetExteriorRing_r( geosinit.ctxt, mGeos );
4383  if ( theRing )
4384  {
4385  nPointsInRing = getNumGeosPoints( theRing );
4386 
4387  wkbPtr << nPointsInRing;
4388 
4389  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4390  unsigned int n;
4391  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &n );
4392 
4393  for ( unsigned int j = 0; j < n; ++j )
4394  {
4395  double x, y;
4396  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4397  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4398  wkbPtr << x << y;
4399  }
4400  }
4401 
4402  //interior rings after
4403  for ( int i = 0; i < GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ); i++ )
4404  {
4405  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, mGeos, i );
4406 
4407  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4408 
4409  unsigned int nPointsInRing;
4410  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &nPointsInRing );
4411  wkbPtr << nPointsInRing;
4412 
4413  for ( unsigned int j = 0; j < nPointsInRing; j++ )
4414  {
4415  double x, y;
4416  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4417  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4418  wkbPtr << x << y;
4419  }
4420  }
4421  mDirtyWkb = false;
4422  return true;
4423  } // case GEOS_GEOM::GEOS_POLYGON
4424  break;
4425 
4426  case GEOS_MULTIPOINT: // a collection of points
4427  {
4428  // determine size of geometry
4429  int geometrySize = 1 + 2 * sizeof( int );
4430  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4431  {
4432  geometrySize += 1 + sizeof( int ) + 2 * sizeof( double );
4433  }
4434 
4435  mGeometry = new unsigned char[geometrySize];
4436  mGeometrySize = geometrySize;
4437 
4438  QgsWkbPtr wkbPtr( mGeometry );
4439  int numPoints = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4440 
4441  wkbPtr << byteOrder << QGis::WKBMultiPoint << numPoints;
4442 
4443  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4444  {
4445  //copy endian and point type
4446  wkbPtr << byteOrder << QGis::WKBPoint;
4447 
4448  const GEOSGeometry *currentPoint = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4449 
4450  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentPoint );
4451 
4452  double x, y;
4453  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, 0, &x );
4454  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, 0, &y );
4455  wkbPtr << x << y;
4456  }
4457  mDirtyWkb = false;
4458  return true;
4459  } // case GEOS_GEOM::GEOS_MULTIPOINT
4460 
4461  case GEOS_MULTILINESTRING: // a collection of linestrings
4462  {
4463  // determine size of geometry
4464  int geometrySize = 1 + 2 * sizeof( int );
4465  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4466  {
4467  geometrySize += 1 + 2 * sizeof( int );
4468  geometrySize += getNumGeosPoints( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) ) * 2 * sizeof( double );
4469  }
4470 
4471  mGeometry = new unsigned char[geometrySize];
4472  mGeometrySize = geometrySize;
4473 
4474  QgsWkbPtr wkbPtr( mGeometry );
4475 
4476  int numLines = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4477  wkbPtr << byteOrder << QGis::WKBMultiLineString << numLines;
4478 
4479  //loop over lines
4480  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4481  {
4482  //endian and type WKBLineString
4483  wkbPtr << byteOrder << QGis::WKBLineString;
4484 
4485  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) );
4486 
4487  //line size
4488  unsigned int lineSize;
4489  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &lineSize );
4490  wkbPtr << lineSize;
4491 
4492  //vertex coordinates
4493  for ( unsigned int j = 0; j < lineSize; ++j )
4494  {
4495  double x, y;
4496  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4497  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4498  wkbPtr << x << y;
4499  }
4500  }
4501  mDirtyWkb = false;
4502  return true;
4503  } // case GEOS_GEOM::GEOS_MULTILINESTRING
4504 
4505  case GEOS_MULTIPOLYGON: // a collection of polygons
4506  {
4507  //first determine size of geometry
4508  int geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of polygons
4509  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4510  {
4511  const GEOSGeometry *thePoly = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4512  geometrySize += 1 + 2 * sizeof( int ); //endian, type, number of rings
4513  //exterior ring
4514  geometrySize += sizeof( int ); //number of points in exterior ring
4515  const GEOSGeometry *exRing = GEOSGetExteriorRing_r( geosinit.ctxt, thePoly );
4516  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( exRing );
4517 
4518  const GEOSGeometry *intRing = 0;
4519  for ( int j = 0; j < GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ); j++ )
4520  {
4521  geometrySize += sizeof( int ); //number of points in ring
4522  intRing = GEOSGetInteriorRingN_r( geosinit.ctxt, thePoly, j );
4523  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( intRing );
4524  }
4525  }
4526 
4527  mGeometry = new unsigned char[geometrySize];
4528  mGeometrySize = geometrySize;
4529 
4530  QgsWkbPtr wkbPtr( mGeometry );
4531  int numPolygons = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4532  wkbPtr << byteOrder << QGis::WKBMultiPolygon << numPolygons;
4533 
4534  //loop over polygons
4535  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4536  {
4537  const GEOSGeometry *thePoly = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4538  int numRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ) + 1;
4539 
4540  //exterior ring
4541  const GEOSGeometry *theRing = GEOSGetExteriorRing_r( geosinit.ctxt, thePoly );
4542  int nPointsInRing = getNumGeosPoints( theRing );
4543 
4544  wkbPtr << byteOrder << QGis::WKBPolygon << numRings << nPointsInRing;
4545 
4546  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4547 
4548  for ( int k = 0; k < nPointsInRing; ++k )
4549  {
4550  double x, y;
4551  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, k, &x );
4552  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, k, &y );
4553  wkbPtr << x << y;
4554  }
4555 
4556  //interior rings
4557  for ( int j = 0; j < GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ); j++ )
4558  {
4559  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, thePoly, j );
4560 
4561  int nPointsInRing = getNumGeosPoints( theRing );
4562  wkbPtr << nPointsInRing;
4563 
4564  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4565 
4566  for ( int k = 0; k < nPointsInRing; ++k )
4567  {
4568  double x, y;
4569  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, k, &x );
4570  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, k, &y );
4571  wkbPtr << x << y;
4572  }
4573  }
4574  }
4575  mDirtyWkb = false;
4576  return true;
4577  } // case GEOS_GEOM::GEOS_MULTIPOLYGON
4578 
4579  case GEOS_GEOMETRYCOLLECTION: // a collection of heterogeneus geometries
4580  {
4581  // TODO
4582  QgsDebugMsg( "geometry collection - not supported" );
4583  break;
4584  } // case GEOS_GEOM::GEOS_GEOMETRYCOLLECTION
4585 
4586  } // switch (mGeos->getGeometryTypeId())
4587 
4588  return false;
4589 }
4590 
4592 {
4593  switch ( destType )
4594  {
4595  case QGis::Point:
4596  return convertToPoint( destMultipart );
4597 
4598  case QGis::Line:
4599  return convertToLine( destMultipart );
4600 
4601  case QGis::Polygon:
4602  return convertToPolygon( destMultipart );
4603 
4604  default:
4605  return 0;
4606  }
4607 }
4608 
4610 {
4611  // TODO: implement with GEOS
4612  if ( mDirtyWkb )
4613  {
4614  exportGeosToWkb();
4615  }
4616 
4617  if ( !mGeometry )
4618  {
4619  return false;
4620  }
4621 
4622  QGis::WkbType geomType = wkbType();
4623 
4624  if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D ||
4625  geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D ||
4626  geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown )
4627  {
4628  return false; //no need to convert
4629  }
4630 
4631  size_t newGeomSize = mGeometrySize + 1 + 2 * sizeof( int ); //endian: 1, multitype: sizeof(int), number of geometries: sizeof(int)
4632  unsigned char* newGeometry = new unsigned char[newGeomSize];
4633 
4634  //copy endian
4635  char byteOrder = QgsApplication::endian();
4636 
4637  QgsWkbPtr wkbPtr( newGeometry );
4638  wkbPtr << byteOrder;
4639 
4640  //copy wkbtype
4641  //todo
4642  QGis::WkbType newMultiType;
4643  switch ( geomType )
4644  {
4645  case QGis::WKBPoint:
4646  newMultiType = QGis::WKBMultiPoint;
4647  break;
4648  case QGis::WKBPoint25D:
4649  newMultiType = QGis::WKBMultiPoint25D;
4650  break;
4651  case QGis::WKBLineString:
4652  newMultiType = QGis::WKBMultiLineString;
4653  break;
4655  newMultiType = QGis::WKBMultiLineString25D;
4656  break;
4657  case QGis::WKBPolygon:
4658  newMultiType = QGis::WKBMultiPolygon;
4659  break;
4660  case QGis::WKBPolygon25D:
4661  newMultiType = QGis::WKBMultiPolygon25D;
4662  break;
4663  default:
4664  delete [] newGeometry;
4665  return false;
4666  }
4667 
4668  wkbPtr << newMultiType << 1;
4669 
4670  //copy the existing single geometry
4671  memcpy( wkbPtr, mGeometry, mGeometrySize );
4672 
4673  delete [] mGeometry;
4674  mGeometry = newGeometry;
4675  mGeometrySize = newGeomSize;
4676  mDirtyGeos = true;
4677  return true;
4678 }
4679 
4680 void QgsGeometry::transformVertex( QgsWkbPtr &wkbPtr, const QTransform& trans, bool hasZValue )
4681 {
4682  double x, y, rotated_x, rotated_y;
4683 
4684  QgsWkbPtr tmp = wkbPtr;
4685  tmp >> x >> y;
4686  trans.map( x, y, &rotated_x, &rotated_y );
4687  wkbPtr << rotated_x << rotated_y;
4688 
4689  if ( hasZValue )
4690  wkbPtr += sizeof( double );
4691 }
4692 
4693 void QgsGeometry::transformVertex( QgsWkbPtr &wkbPtr, const QgsCoordinateTransform& ct, bool hasZValue )
4694 {
4695  double x, y, z = 0.0;
4696 
4697  QgsWkbPtr tmp = wkbPtr;
4698  tmp >> x >> y;
4699  ct.transformInPlace( x, y, z );
4700  wkbPtr << x << y;
4701 
4702  if ( hasZValue )
4703  wkbPtr += sizeof( double );
4704 }
4705 
4706 GEOSGeometry* QgsGeometry::linePointDifference( GEOSGeometry* GEOSsplitPoint )
4707 {
4708  int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
4709  QgsMultiPolyline multiLine;
4710 
4711  if ( type == GEOS_MULTILINESTRING )
4712  multiLine = asMultiPolyline();
4713  else if ( type == GEOS_LINESTRING )
4714  multiLine = QgsMultiPolyline() << asPolyline();
4715  else
4716  return 0;
4717 
4718  // GEOSsplitPoint will be deleted in the caller, so make a clone
4719  QgsGeometry* geosPoint = fromGeosGeom( GEOSGeom_clone_r( geosinit.ctxt, GEOSsplitPoint ) );
4720  QgsPoint splitPoint = geosPoint->asPoint();
4721  delete geosPoint;
4722 
4723  QgsMultiPolyline lines;
4724  QgsPolyline line;
4725  QgsPolyline newline;
4726 
4727  //For each part
4728  for ( int i = 0; i < multiLine.size() ; ++i )
4729  {
4730  line = multiLine[i];
4731  newline = QgsPolyline();
4732  newline.append( line[0] );
4733  //For each segment
4734  for ( int j = 1; j < line.size() - 1 ; ++j )
4735  {
4736  newline.append( line[j] );
4737  if ( line[j] == splitPoint )
4738  {
4739  lines.append( newline );
4740  newline = QgsPolyline();
4741  newline.append( line[j] );
4742  }
4743  }
4744  newline.append( line.last() );
4745  lines.append( newline );
4746  }
4747  QgsGeometry* splitLines = fromMultiPolyline( lines );
4748  GEOSGeometry* splitGeom = GEOSGeom_clone_r( geosinit.ctxt, splitLines->asGeos() );
4749  delete splitLines;
4750 
4751  return splitGeom;
4752 }
4753 
4754 int QgsGeometry::splitLinearGeometry( GEOSGeometry *splitLine, QList<QgsGeometry*>& newGeometries )
4755 {
4756  if ( !splitLine )
4757  return 2;
4758 
4759  if ( mDirtyGeos )
4760  exportWkbToGeos();
4761 
4762  if ( !mGeos )
4763  return 5;
4764 
4765  //first test if linestring intersects geometry. If not, return straight away
4766  if ( !GEOSIntersects_r( geosinit.ctxt, splitLine, mGeos ) )
4767  return 1;
4768 
4769  //check that split line has no linear intersection
4770  int linearIntersect = GEOSRelatePattern_r( geosinit.ctxt, mGeos, splitLine, "1********" );
4771  if ( linearIntersect > 0 )
4772  return 3;
4773 
4774  int splitGeomType = GEOSGeomTypeId_r( geosinit.ctxt, splitLine );
4775 
4776  GEOSGeometry* splitGeom;
4777  if ( splitGeomType == GEOS_POINT )
4778  {
4779  splitGeom = linePointDifference( splitLine );
4780  }
4781  else
4782  {
4783  splitGeom = GEOSDifference_r( geosinit.ctxt, mGeos, splitLine );
4784  }
4785  QVector<GEOSGeometry*> lineGeoms;
4786 
4787  int splitType = GEOSGeomTypeId_r( geosinit.ctxt, splitGeom );
4788  if ( splitType == GEOS_MULTILINESTRING )
4789  {
4790  int nGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, splitGeom );
4791  for ( int i = 0; i < nGeoms; ++i )
4792  lineGeoms << GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, splitGeom, i ) );
4793 
4794  }
4795  else
4796  {
4797  lineGeoms << GEOSGeom_clone_r( geosinit.ctxt, splitGeom );
4798  }
4799 
4800  mergeGeometriesMultiTypeSplit( lineGeoms );
4801 
4802  if ( lineGeoms.size() > 0 )
4803  {
4804  fromGeos( lineGeoms[0] );
4805  }
4806 
4807  for ( int i = 1; i < lineGeoms.size(); ++i )
4808  {
4809  newGeometries << fromGeosGeom( lineGeoms[i] );
4810  }
4811 
4812  GEOSGeom_destroy_r( geosinit.ctxt, splitGeom );
4813  return 0;
4814 }
4815 
4816 int QgsGeometry::splitPolygonGeometry( GEOSGeometry* splitLine, QList<QgsGeometry*>& newGeometries )
4817 {
4818  if ( !splitLine )
4819  return 2;
4820 
4821  if ( mDirtyGeos )
4822  exportWkbToGeos();
4823 
4824  if ( !mGeos )
4825  return 5;
4826 
4827  //first test if linestring intersects geometry. If not, return straight away
4828  if ( !GEOSIntersects_r( geosinit.ctxt, splitLine, mGeos ) )
4829  return 1;
4830 
4831  //first union all the polygon rings together (to get them noded, see JTS developer guide)
4832  GEOSGeometry *nodedGeometry = nodeGeometries( splitLine, mGeos );
4833  if ( !nodedGeometry )
4834  return 2; //an error occured during noding
4835 
4836  GEOSGeometry *polygons = GEOSPolygonize_r( geosinit.ctxt, &nodedGeometry, 1 );
4837  if ( !polygons || numberOfGeometries( polygons ) == 0 )
4838  {
4839  if ( polygons )
4840  GEOSGeom_destroy_r( geosinit.ctxt, polygons );
4841 
4842  GEOSGeom_destroy_r( geosinit.ctxt, nodedGeometry );
4843 
4844  return 4;
4845  }
4846 
4847  GEOSGeom_destroy_r( geosinit.ctxt, nodedGeometry );
4848 
4849  //test every polygon if contained in original geometry
4850  //include in result if yes
4851  QVector<GEOSGeometry*> testedGeometries;
4852  GEOSGeometry *intersectGeometry = 0;
4853 
4854  //ratio intersect geometry / geometry. This should be close to 1
4855  //if the polygon belongs to the input geometry
4856 
4857  for ( int i = 0; i < numberOfGeometries( polygons ); i++ )
4858  {
4859  const GEOSGeometry *polygon = GEOSGetGeometryN_r( geosinit.ctxt, polygons, i );
4860  intersectGeometry = GEOSIntersection_r( geosinit.ctxt, mGeos, polygon );
4861  if ( !intersectGeometry )
4862  {
4863  QgsDebugMsg( "intersectGeometry is NULL" );
4864  continue;
4865  }
4866 
4867  double intersectionArea;
4868  GEOSArea_r( geosinit.ctxt, intersectGeometry, &intersectionArea );
4869 
4870  double polygonArea;
4871  GEOSArea_r( geosinit.ctxt, polygon, &polygonArea );
4872 
4873  const double areaRatio = intersectionArea / polygonArea;
4874  if ( areaRatio > 0.99 && areaRatio < 1.01 )
4875  testedGeometries << GEOSGeom_clone_r( geosinit.ctxt, polygon );
4876 
4877  GEOSGeom_destroy_r( geosinit.ctxt, intersectGeometry );
4878  }
4879 
4880  bool splitDone = true;
4881  int nGeometriesThis = numberOfGeometries( mGeos ); //original number of geometries
4882  if ( testedGeometries.size() == nGeometriesThis )
4883  {
4884  splitDone = false;
4885  }
4886 
4887  mergeGeometriesMultiTypeSplit( testedGeometries );
4888 
4889  //no split done, preserve original geometry
4890  if ( !splitDone )
4891  {
4892  for ( int i = 0; i < testedGeometries.size(); ++i )
4893  {
4894  GEOSGeom_destroy_r( geosinit.ctxt, testedGeometries[i] );
4895  }
4896  return 1;
4897  }
4898  else if ( testedGeometries.size() > 0 ) //split successfull
4899  {
4900  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
4901  mGeos = testedGeometries[0];
4902  mDirtyWkb = true;
4903  }
4904 
4905  int i;
4906  for ( i = 1; i < testedGeometries.size() && GEOSisValid_r( geosinit.ctxt, testedGeometries[i] ); ++i )
4907  ;
4908 
4909  if ( i < testedGeometries.size() )
4910  {
4911  for ( i = 0; i < testedGeometries.size(); ++i )
4912  GEOSGeom_destroy_r( geosinit.ctxt, testedGeometries[i] );
4913 
4914  return 3;
4915  }
4916 
4917  for ( i = 1; i < testedGeometries.size(); ++i )
4918  newGeometries << fromGeosGeom( testedGeometries[i] );
4919 
4920  GEOSGeom_destroy_r( geosinit.ctxt, polygons );
4921  return 0;
4922 }
4923 
4924 GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos )
4925 {
4926  //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
4927  int nIntersections = 0;
4928  int lastIntersectingRing = -2;
4929  const GEOSGeometry* lastIntersectingGeom = 0;
4930 
4931  int nRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, polygon );
4932  if ( nRings < 0 )
4933  return 0;
4934 
4935  //does outer ring intersect?
4936  const GEOSGeometry* outerRing = GEOSGetExteriorRing_r( geosinit.ctxt, polygon );
4937  if ( GEOSIntersects_r( geosinit.ctxt, outerRing, reshapeLineGeos ) == 1 )
4938  {
4939  ++nIntersections;
4940  lastIntersectingRing = -1;
4941  lastIntersectingGeom = outerRing;
4942  }
4943 
4944  //do inner rings intersect?
4945  const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
4946 
4947  try
4948  {
4949  for ( int i = 0; i < nRings; ++i )
4950  {
4951  innerRings[i] = GEOSGetInteriorRingN_r( geosinit.ctxt, polygon, i );
4952  if ( GEOSIntersects_r( geosinit.ctxt, innerRings[i], reshapeLineGeos ) == 1 )
4953  {
4954  ++nIntersections;
4955  lastIntersectingRing = i;
4956  lastIntersectingGeom = innerRings[i];
4957  }
4958  }
4959  }
4960  catch ( GEOSException &e )
4961  {
4962  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
4963  nIntersections = 0;
4964  }
4965 
4966  if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
4967  {
4968  delete [] innerRings;
4969  return 0;
4970  }
4971 
4972  //we have one intersecting ring, let's try to reshape it
4973  GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos );
4974  if ( !reshapeResult )
4975  {
4976  delete [] innerRings;
4977  return 0;
4978  }
4979 
4980  //if reshaping took place, we need to reassemble the polygon and its rings
4981  GEOSGeometry* newRing = 0;
4982  const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, reshapeResult );
4983  GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone_r( geosinit.ctxt, reshapeSequence );
4984 
4985  GEOSGeom_destroy_r( geosinit.ctxt, reshapeResult );
4986 
4987  try
4988  {
4989  newRing = GEOSGeom_createLinearRing_r( geosinit.ctxt, newCoordSequence );
4990  }
4991  catch ( GEOSException &e )
4992  {
4993  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
4994  }
4995 
4996  if ( !newRing )
4997  {
4998  delete [] innerRings;
4999  return 0;
5000  }
5001 
5002  GEOSGeometry* newOuterRing = 0;
5003  if ( lastIntersectingRing == -1 )
5004  newOuterRing = newRing;
5005  else
5006  newOuterRing = GEOSGeom_clone_r( geosinit.ctxt, outerRing );
5007 
5008  //check if all the rings are still inside the outer boundary
5009  QList<GEOSGeometry*> ringList;
5010  if ( nRings > 0 )
5011  {
5012  GEOSGeometry* outerRingPoly = GEOSGeom_createPolygon_r( geosinit.ctxt, GEOSGeom_clone_r( geosinit.ctxt, newOuterRing ), 0, 0 );
5013  if ( outerRingPoly )
5014  {
5015  GEOSGeometry* currentRing = 0;
5016  for ( int i = 0; i < nRings; ++i )
5017  {
5018  if ( lastIntersectingRing == i )
5019  currentRing = newRing;
5020  else
5021  currentRing = GEOSGeom_clone_r( geosinit.ctxt, innerRings[i] );
5022 
5023  //possibly a ring is no longer contained in the result polygon after reshape
5024  if ( GEOSContains_r( geosinit.ctxt, outerRingPoly, currentRing ) == 1 )
5025  ringList.push_back( currentRing );
5026  else
5027  GEOSGeom_destroy_r( geosinit.ctxt, currentRing );
5028  }
5029  }
5030  GEOSGeom_destroy_r( geosinit.ctxt, outerRingPoly );
5031  }
5032 
5033  GEOSGeometry** newInnerRings = new GEOSGeometry*[ringList.size()];
5034  for ( int i = 0; i < ringList.size(); ++i )
5035  newInnerRings[i] = ringList.at( i );
5036 
5037  delete [] innerRings;
5038 
5039  GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon_r( geosinit.ctxt, newOuterRing, newInnerRings, ringList.size() );
5040  delete[] newInnerRings;
5041 
5042  return reshapedPolygon;
5043 }
5044 
5045 GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos )
5046 {
5047  if ( !line || !reshapeLineGeos )
5048  return 0;
5049 
5050  bool atLeastTwoIntersections = false;
5051 
5052  try
5053  {
5054  //make sure there are at least two intersection between line and reshape geometry
5055  GEOSGeometry* intersectGeom = GEOSIntersection_r( geosinit.ctxt, line, reshapeLineGeos );
5056  if ( intersectGeom )
5057  {
5058  atLeastTwoIntersections = ( GEOSGeomTypeId_r( geosinit.ctxt, intersectGeom ) == GEOS_MULTIPOINT && GEOSGetNumGeometries_r( geosinit.ctxt, intersectGeom ) > 1 );
5059  GEOSGeom_destroy_r( geosinit.ctxt, intersectGeom );
5060  }
5061  }
5062  catch ( GEOSException &e )
5063  {
5064  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
5065  atLeastTwoIntersections = false;
5066  }
5067 
5068  if ( !atLeastTwoIntersections )
5069  return 0;
5070 
5071  //begin and end point of original line
5072  const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, line );
5073  if ( !lineCoordSeq )
5074  return 0;
5075 
5076  unsigned int lineCoordSeqSize;
5077  if ( GEOSCoordSeq_getSize_r( geosinit.ctxt, lineCoordSeq, &lineCoordSeqSize ) == 0 )
5078  return 0;
5079 
5080  if ( lineCoordSeqSize < 2 )
5081  return 0;
5082 
5083  //first and last vertex of line
5084  double x1, y1, x2, y2;
5085  GEOSCoordSeq_getX_r( geosinit.ctxt, lineCoordSeq, 0, &x1 );
5086  GEOSCoordSeq_getY_r( geosinit.ctxt, lineCoordSeq, 0, &y1 );
5087  GEOSCoordSeq_getX_r( geosinit.ctxt, lineCoordSeq, lineCoordSeqSize - 1, &x2 );
5088  GEOSCoordSeq_getY_r( geosinit.ctxt, lineCoordSeq, lineCoordSeqSize - 1, &y2 );
5089  GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) );
5090  GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) );
5091 
5092  bool isRing = false;
5093  if ( GEOSGeomTypeId_r( geosinit.ctxt, line ) == GEOS_LINEARRING || GEOSEquals_r( geosinit.ctxt, beginLineVertex, endLineVertex ) == 1 )
5094  isRing = true;
5095 
5096 //node line and reshape line
5097  GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line );
5098  if ( !nodedGeometry )
5099  {
5100  GEOSGeom_destroy_r( geosinit.ctxt, beginLineVertex );
5101  GEOSGeom_destroy_r( geosinit.ctxt, endLineVertex );
5102  return 0;
5103  }
5104 
5105  //and merge them together
5106  GEOSGeometry *mergedLines = GEOSLineMerge_r( geosinit.ctxt, nodedGeometry );
5107  GEOSGeom_destroy_r( geosinit.ctxt, nodedGeometry );
5108  if ( !mergedLines )
5109  {
5110  GEOSGeom_destroy_r( geosinit.ctxt, beginLineVertex );
5111  GEOSGeom_destroy_r( geosinit.ctxt, endLineVertex );
5112  return 0;
5113  }
5114 
5115  int numMergedLines = GEOSGetNumGeometries_r( geosinit.ctxt, mergedLines );
5116  if ( numMergedLines < 2 ) //some special cases. Normally it is >2
5117  {
5118  GEOSGeom_destroy_r( geosinit.ctxt, beginLineVertex );
5119  GEOSGeom_destroy_r( geosinit.ctxt, endLineVertex );
5120  if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
5121  return GEOSGeom_clone_r( geosinit.ctxt, reshapeLineGeos );
5122  else
5123  return 0;
5124  }
5125 
5126  QList<GEOSGeometry*> resultLineParts; //collection with the line segments that will be contained in result
5127  QList<GEOSGeometry*> probableParts; //parts where we can decide on inclusion only after going through all the candidates
5128 
5129  for ( int i = 0; i < numMergedLines; ++i )
5130  {
5131  const GEOSGeometry* currentGeom;
5132 
5133  currentGeom = GEOSGetGeometryN_r( geosinit.ctxt, mergedLines, i );
5134  const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentGeom );
5135  unsigned int currentCoordSeqSize;
5136  GEOSCoordSeq_getSize_r( geosinit.ctxt, currentCoordSeq, &currentCoordSeqSize );
5137  if ( currentCoordSeqSize < 2 )
5138  continue;
5139 
5140  //get the two endpoints of the current line merge result
5141  double xBegin, xEnd, yBegin, yEnd;
5142  GEOSCoordSeq_getX_r( geosinit.ctxt, currentCoordSeq, 0, &xBegin );
5143  GEOSCoordSeq_getY_r( geosinit.ctxt, currentCoordSeq, 0, &yBegin );
5144  GEOSCoordSeq_getX_r( geosinit.ctxt, currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
5145  GEOSCoordSeq_getY_r( geosinit.ctxt, currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
5146  GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) );
5147  GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) );
5148 
5149  //check how many endpoints of the line merge result are on the (original) line
5150  int nEndpointsOnOriginalLine = 0;
5151  if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 )
5152  nEndpointsOnOriginalLine += 1;
5153 
5154  if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 )
5155  nEndpointsOnOriginalLine += 1;
5156 
5157  //check how many endpoints equal the endpoints of the original line
5158  int nEndpointsSameAsOriginalLine = 0;
5159  if ( GEOSEquals_r( geosinit.ctxt, beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals_r( geosinit.ctxt, beginCurrentGeomVertex, endLineVertex ) == 1 )
5160  nEndpointsSameAsOriginalLine += 1;
5161 
5162  if ( GEOSEquals_r( geosinit.ctxt, endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals_r( geosinit.ctxt, endCurrentGeomVertex, endLineVertex ) == 1 )
5163  nEndpointsSameAsOriginalLine += 1;
5164 
5165  //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
5166  bool currentGeomOverlapsOriginalGeom = false;
5167  bool currentGeomOverlapsReshapeLine = false;
5168  if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 )
5169  currentGeomOverlapsOriginalGeom = true;
5170 
5171  if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
5172  currentGeomOverlapsReshapeLine = true;
5173 
5174  //logic to decide if this part belongs to the result
5175  if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
5176  {
5177  resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
5178  }
5179  //for closed rings, we take one segment from the candidate list
5180  else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
5181  {
5182  probableParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
5183  }
5184  else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
5185  {
5186  resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
5187  }
5188  else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
5189  {
5190  resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
5191  }
5192  else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
5193  {
5194  resultLineParts.push_back( GEOSGeom_clone_r( geosinit.ctxt, currentGeom ) );
5195  }
5196 
5197  GEOSGeom_destroy_r( geosinit.ctxt, beginCurrentGeomVertex );
5198  GEOSGeom_destroy_r( geosinit.ctxt, endCurrentGeomVertex );
5199  }
5200 
5201  //add the longest segment from the probable list for rings (only used for polygon rings)
5202  if ( isRing && probableParts.size() > 0 )
5203  {
5204  GEOSGeometry* maxGeom = 0; //the longest geometry in the probabla list
5205  GEOSGeometry* currentGeom = 0;
5206  double maxLength = -DBL_MAX;
5207  double currentLength = 0;
5208  for ( int i = 0; i < probableParts.size(); ++i )
5209  {
5210  currentGeom = probableParts.at( i );
5211  GEOSLength_r( geosinit.ctxt, currentGeom, &currentLength );
5212  if ( currentLength > maxLength )
5213  {
5214  maxLength = currentLength;
5215  GEOSGeom_destroy_r( geosinit.ctxt, maxGeom );
5216  maxGeom = currentGeom;
5217  }
5218  else
5219  {
5220  GEOSGeom_destroy_r( geosinit.ctxt, currentGeom );
5221  }
5222  }
5223  resultLineParts.push_back( maxGeom );
5224  }
5225 
5226  GEOSGeom_destroy_r( geosinit.ctxt, beginLineVertex );
5227  GEOSGeom_destroy_r( geosinit.ctxt, endLineVertex );
5228  GEOSGeom_destroy_r( geosinit.ctxt, mergedLines );
5229 
5230  GEOSGeometry* result = 0;
5231  if ( resultLineParts.size() < 1 )
5232  return 0;
5233 
5234  if ( resultLineParts.size() == 1 ) //the whole result was reshaped
5235  {
5236  result = resultLineParts[0];
5237  }
5238  else //>1
5239  {
5240  GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
5241  for ( int i = 0; i < resultLineParts.size(); ++i )
5242  {
5243  lineArray[i] = resultLineParts[i];
5244  }
5245 
5246  //create multiline from resultLineParts
5247  GEOSGeometry* multiLineGeom = GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, lineArray, resultLineParts.size() );
5248  delete [] lineArray;
5249 
5250  //then do a linemerge with the newly combined partstrings
5251  result = GEOSLineMerge_r( geosinit.ctxt, multiLineGeom );
5252  GEOSGeom_destroy_r( geosinit.ctxt, multiLineGeom );
5253  }
5254 
5255  //now test if the result is a linestring. Otherwise something went wrong
5256  if ( GEOSGeomTypeId_r( geosinit.ctxt, result ) != GEOS_LINESTRING )
5257  {
5258  GEOSGeom_destroy_r( geosinit.ctxt, result );
5259  return 0;
5260  }
5261 
5262  return result;
5263 }
5264 
5265 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const
5266 {
5267  //Find out the intersection points between splitLineGeos and this geometry.
5268  //These points need to be tested for topological correctness by the calling function
5269  //if topological editing is enabled
5270 
5271  testPoints.clear();
5272  GEOSGeometry* intersectionGeom = GEOSIntersection_r( geosinit.ctxt, mGeos, splitLine );
5273  if ( !intersectionGeom )
5274  return 1;
5275 
5276  bool simple = false;
5277  int nIntersectGeoms = 1;
5278  if ( GEOSGeomTypeId_r( geosinit.ctxt, intersectionGeom ) == GEOS_LINESTRING || GEOSGeomTypeId_r( geosinit.ctxt, intersectionGeom ) == GEOS_POINT )
5279  simple = true;
5280 
5281  if ( !simple )
5282  nIntersectGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, intersectionGeom );
5283 
5284  for ( int i = 0; i < nIntersectGeoms; ++i )
5285  {
5286  const GEOSGeometry* currentIntersectGeom;
5287  if ( simple )
5288  currentIntersectGeom = intersectionGeom;
5289  else
5290  currentIntersectGeom = GEOSGetGeometryN_r( geosinit.ctxt, intersectionGeom, i );
5291 
5292  const GEOSCoordSequence* lineSequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentIntersectGeom );
5293  unsigned int sequenceSize = 0;
5294  double x, y;
5295  if ( GEOSCoordSeq_getSize_r( geosinit.ctxt, lineSequence, &sequenceSize ) != 0 )
5296  {
5297  for ( unsigned int i = 0; i < sequenceSize; ++i )
5298  {
5299  if ( GEOSCoordSeq_getX_r( geosinit.ctxt, lineSequence, i, &x ) != 0 )
5300  {
5301  if ( GEOSCoordSeq_getY_r( geosinit.ctxt, lineSequence, i, &y ) != 0 )
5302  {
5303  testPoints.push_back( QgsPoint( x, y ) );
5304  }
5305  }
5306  }
5307  }
5308  }
5309  GEOSGeom_destroy_r( geosinit.ctxt, intersectionGeom );
5310  return 0;
5311 }
5312 
5313 GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
5314 {
5315  if ( !splitLine || !geom )
5316  return 0;
5317 
5318  if ( GEOSGeomTypeId_r( geosinit.ctxt, geom ) == GEOS_POLYGON || GEOSGeomTypeId_r( geosinit.ctxt, geom ) == GEOS_MULTIPOLYGON )
5319  {
5320  GEOSGeometry *geometryBoundary = GEOSBoundary_r( geosinit.ctxt, geom );
5321  GEOSGeometry *unionGeometry = GEOSUnion_r( geosinit.ctxt, splitLine, geometryBoundary );
5322  GEOSGeom_destroy_r( geosinit.ctxt, geometryBoundary );
5323  return unionGeometry;
5324  }
5325  else
5326  {
5327  return GEOSUnion_r( geosinit.ctxt, splitLine, geom );
5328  }
5329 }
5330 
5331 int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 )
5332 {
5333  if ( !line1 || !line2 )
5334  {
5335  return -1;
5336  }
5337 
5338  double bufferDistance = pow( 10.0L, geomDigits( line2 ) - 11 );
5339 
5340  GEOSGeometry* bufferGeom = GEOSBuffer_r( geosinit.ctxt, line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS );
5341  if ( !bufferGeom )
5342  return -2;
5343 
5344  GEOSGeometry* intersectionGeom = GEOSIntersection_r( geosinit.ctxt, bufferGeom, line1 );
5345 
5346  //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
5347  double intersectGeomLength;
5348  double line1Length;
5349 
5350  GEOSLength_r( geosinit.ctxt, intersectionGeom, &intersectGeomLength );
5351  GEOSLength_r( geosinit.ctxt, line1, &line1Length );
5352 
5353  GEOSGeom_destroy_r( geosinit.ctxt, bufferGeom );
5354  GEOSGeom_destroy_r( geosinit.ctxt, intersectionGeom );
5355 
5356  double intersectRatio = line1Length / intersectGeomLength;
5357  if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
5358  return 1;
5359 
5360  return 0;
5361 }
5362 
5363 int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line )
5364 {
5365  if ( !point || !line )
5366  return -1;
5367 
5368  double bufferDistance = pow( 10.0L, geomDigits( line ) - 11 );
5369 
5370  GEOSGeometry* lineBuffer = GEOSBuffer_r( geosinit.ctxt, line, bufferDistance, 8 );
5371  if ( !lineBuffer )
5372  return -2;
5373 
5374  bool contained = false;
5375  if ( GEOSContains_r( geosinit.ctxt, lineBuffer, point ) == 1 )
5376  contained = true;
5377 
5378  GEOSGeom_destroy_r( geosinit.ctxt, lineBuffer );
5379  return contained;
5380 }
5381 
5382 int QgsGeometry::geomDigits( const GEOSGeometry* geom )
5383 {
5384  GEOSGeometry* bbox = GEOSEnvelope_r( geosinit.ctxt, geom );
5385  if ( !bbox )
5386  return -1;
5387 
5388  const GEOSGeometry* bBoxRing = GEOSGetExteriorRing_r( geosinit.ctxt, bbox );
5389  if ( !bBoxRing )
5390  return -1;
5391 
5392  const GEOSCoordSequence* bBoxCoordSeq = GEOSGeom_getCoordSeq_r( geosinit.ctxt, bBoxRing );
5393 
5394  if ( !bBoxCoordSeq )
5395  return -1;
5396 
5397  unsigned int nCoords = 0;
5398  if ( !GEOSCoordSeq_getSize_r( geosinit.ctxt, bBoxCoordSeq, &nCoords ) )
5399  return -1;
5400 
5401  int maxDigits = -1;
5402  for ( unsigned int i = 0; i < nCoords - 1; ++i )
5403  {
5404  double t;
5405  GEOSCoordSeq_getX_r( geosinit.ctxt, bBoxCoordSeq, i, &t );
5406 
5407  int digits;
5408  digits = ceil( log10( fabs( t ) ) );
5409  if ( digits > maxDigits )
5410  maxDigits = digits;
5411 
5412  GEOSCoordSeq_getY_r( geosinit.ctxt, bBoxCoordSeq, i, &t );
5413  digits = ceil( log10( fabs( t ) ) );
5414  if ( digits > maxDigits )
5415  maxDigits = digits;
5416  }
5417 
5418  return maxDigits;
5419 }
5420 
5421 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const
5422 {
5423  if ( !g )
5424  return 0;
5425 
5426  int geometryType = GEOSGeomTypeId_r( geosinit.ctxt, g );
5427  if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
5428  || geometryType == GEOS_POLYGON )
5429  return 1;
5430 
5431  //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
5432  return GEOSGetNumGeometries_r( geosinit.ctxt, g );
5433 }
5434 
5435 int QgsGeometry::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry*>& splitResult )
5436 {
5437  if ( mDirtyGeos )
5438  exportWkbToGeos();
5439 
5440  if ( !mGeos )
5441  return 1;
5442 
5443  //convert mGeos to geometry collection
5444  int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
5445  if ( type != GEOS_GEOMETRYCOLLECTION &&
5446  type != GEOS_MULTILINESTRING &&
5447  type != GEOS_MULTIPOLYGON &&
5448  type != GEOS_MULTIPOINT )
5449  return 0;
5450 
5451  QVector<GEOSGeometry*> copyList = splitResult;
5452  splitResult.clear();
5453 
5454  //collect all the geometries that belong to the initial multifeature
5455  QVector<GEOSGeometry*> unionGeom;
5456 
5457  for ( int i = 0; i < copyList.size(); ++i )
5458  {
5459  //is this geometry a part of the original multitype?
5460  bool isPart = false;
5461  for ( int j = 0; j < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); j++ )
5462  {
5463  if ( GEOSEquals_r( geosinit.ctxt, copyList[i], GEOSGetGeometryN_r( geosinit.ctxt, mGeos, j ) ) )
5464  {
5465  isPart = true;
5466  break;
5467  }
5468  }
5469 
5470  if ( isPart )
5471  {
5472  unionGeom << copyList[i];
5473  }
5474  else
5475  {
5476  QVector<GEOSGeometry*> geomVector;
5477  geomVector << copyList[i];
5478 
5479  if ( type == GEOS_MULTILINESTRING )
5480  splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector );
5481  else if ( type == GEOS_MULTIPOLYGON )
5482  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
5483  else
5484  GEOSGeom_destroy_r( geosinit.ctxt, copyList[i] );
5485  }
5486  }
5487 
5488  //make multifeature out of unionGeom
5489  if ( unionGeom.size() > 0 )
5490  {
5491  if ( type == GEOS_MULTILINESTRING )
5492  splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom );
5493  else if ( type == GEOS_MULTIPOLYGON )
5494  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom );
5495  }
5496  else
5497  {
5498  unionGeom.clear();
5499  }
5500 
5501  return 0;
5502 }
5503 
5504 QgsPoint QgsGeometry::asPoint( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5505 {
5506  wkbPtr += 1 + sizeof( int );
5507 
5508  double x, y;
5509  wkbPtr >> x >> y;
5510  if ( hasZValue )
5511  wkbPtr += sizeof( double );
5512 
5513  return QgsPoint( x, y );
5514 }
5515 
5516 QgsPolyline QgsGeometry::asPolyline( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5517 {
5518  wkbPtr += 1 + sizeof( int );
5519 
5520  unsigned int nPoints;
5521  wkbPtr >> nPoints;
5522 
5523  QgsPolyline line( nPoints );
5524 
5525  // Extract the points from the WKB format into the x and y vectors.
5526  for ( uint i = 0; i < nPoints; ++i )
5527  {
5528  double x, y;
5529  wkbPtr >> x >> y;
5530  if ( hasZValue )
5531  wkbPtr += sizeof( double );
5532 
5533  line[i] = QgsPoint( x, y );
5534  }
5535 
5536  return line;
5537 }
5538 
5539 QgsPolygon QgsGeometry::asPolygon( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5540 {
5541  wkbPtr += 1 + sizeof( int );
5542 
5543  // get number of rings in the polygon
5544  unsigned int numRings;
5545  wkbPtr >> numRings;
5546 
5547  if ( numRings == 0 ) // sanity check for zero rings in polygon
5548  return QgsPolygon();
5549 
5550  QgsPolygon rings( numRings );
5551 
5552  for ( uint idx = 0; idx < numRings; idx++ )
5553  {
5554  int nPoints;
5555  wkbPtr >> nPoints;
5556 
5557  QgsPolyline ring( nPoints );
5558 
5559  for ( int jdx = 0; jdx < nPoints; jdx++ )
5560  {
5561  double x, y;
5562  wkbPtr >> x >> y;
5563  if ( hasZValue )
5564  wkbPtr += sizeof( double );
5565 
5566  ring[jdx] = QgsPoint( x, y );
5567  }
5568 
5569  rings[idx] = ring;
5570  }
5571 
5572  return rings;
5573 }
5574 
5576 {
5577  QGis::WkbType type = wkbType();
5578  if ( type != QGis::WKBPoint && type != QGis::WKBPoint25D )
5579  return QgsPoint( 0, 0 );
5580 
5581  QgsConstWkbPtr wkbPtr( mGeometry );
5582  return asPoint( wkbPtr, type == QGis::WKBPoint25D );
5583 }
5584 
5586 {
5587  QGis::WkbType type = wkbType();
5588  if ( type != QGis::WKBLineString && type != QGis::WKBLineString25D )
5589  return QgsPolyline();
5590 
5591  QgsConstWkbPtr wkbPtr( mGeometry );
5592  return asPolyline( wkbPtr, type == QGis::WKBLineString25D );
5593 }
5594 
5596 {
5597  QGis::WkbType type = wkbType();
5598  if ( type != QGis::WKBPolygon && type != QGis::WKBPolygon25D )
5599  return QgsPolygon();
5600 
5601  QgsConstWkbPtr wkbPtr( mGeometry );
5602  return asPolygon( wkbPtr, type == QGis::WKBPolygon25D );
5603 }
5604 
5606 {
5607  QGis::WkbType type = wkbType();
5608  if ( type != QGis::WKBMultiPoint && type != QGis::WKBMultiPoint25D )
5609  return QgsMultiPoint();
5610 
5611  bool hasZValue = ( type == QGis::WKBMultiPoint25D );
5612 
5613  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5614 
5615  int nPoints;
5616  wkbPtr >> nPoints;
5617 
5618  QgsMultiPoint points( nPoints );
5619  for ( int i = 0; i < nPoints; i++ )
5620  {
5621  points[i] = asPoint( wkbPtr, hasZValue );
5622  }
5623 
5624  return points;
5625 }
5626 
5628 {
5629  QGis::WkbType type = wkbType();
5630  if ( type != QGis::WKBMultiLineString && type != QGis::WKBMultiLineString25D )
5631  return QgsMultiPolyline();
5632 
5633  bool hasZValue = ( type == QGis::WKBMultiLineString25D );
5634 
5635  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5636 
5637  int numLineStrings;
5638  wkbPtr >> numLineStrings;
5639 
5640  QgsMultiPolyline lines( numLineStrings );
5641 
5642  for ( int i = 0; i < numLineStrings; i++ )
5643  lines[i] = asPolyline( wkbPtr, hasZValue );
5644 
5645  return lines;
5646 }
5647 
5649 {
5650  QGis::WkbType type = wkbType();
5651  if ( type != QGis::WKBMultiPolygon && type != QGis::WKBMultiPolygon25D )
5652  return QgsMultiPolygon();
5653 
5654  bool hasZValue = ( type == QGis::WKBMultiPolygon25D );
5655 
5656  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5657 
5658  int numPolygons;
5659  wkbPtr >> numPolygons;
5660 
5661  QgsMultiPolygon polygons( numPolygons );
5662 
5663  for ( int i = 0; i < numPolygons; i++ )
5664  polygons[i] = asPolygon( wkbPtr, hasZValue );
5665 
5666  return polygons;
5667 }
5668 
5670 {
5671  if ( mDirtyGeos )
5672  exportWkbToGeos();
5673 
5674  if ( !mGeos )
5675  return -1.0;
5676 
5677  double area;
5678 
5679  try
5680  {
5681  if ( GEOSArea_r( geosinit.ctxt, mGeos, &area ) == 0 )
5682  return -1.0;
5683  }
5684  CATCH_GEOS( -1.0 )
5685 
5686  return area;
5687 }
5688 
5690 {
5691  if ( mDirtyGeos )
5692  exportWkbToGeos();
5693 
5694  if ( !mGeos )
5695  return -1.0;
5696 
5697  double length;
5698 
5699  try
5700  {
5701  if ( GEOSLength_r( geosinit.ctxt, mGeos, &length ) == 0 )
5702  return -1.0;
5703  }
5704  CATCH_GEOS( -1.0 )
5705 
5706  return length;
5707 }
5709 {
5710  if ( mDirtyGeos )
5711  exportWkbToGeos();
5712 
5713  if ( geom.mDirtyGeos )
5714  geom.exportWkbToGeos();
5715 
5716  if ( !mGeos || !geom.mGeos )
5717  return -1.0;
5718 
5719  double dist = -1.0;
5720 
5721  try
5722  {
5723  GEOSDistance_r( geosinit.ctxt, mGeos, geom.mGeos, &dist );
5724  }
5725  CATCH_GEOS( -1.0 )
5726 
5727  return dist;
5728 }
5729 
5730 QgsGeometry* QgsGeometry::buffer( double distance, int segments )
5731 {
5732  if ( mDirtyGeos )
5733  exportWkbToGeos();
5734 
5735  if ( !mGeos )
5736  return 0;
5737 
5738  try
5739  {
5740  return fromGeosGeom( GEOSBuffer_r( geosinit.ctxt, mGeos, distance, segments ) );
5741  }
5742  CATCH_GEOS( 0 )
5743 }
5744 
5745 QgsGeometry*QgsGeometry::buffer( double distance, int segments, int endCapStyle, int joinStyle, double mitreLimit )
5746 {
5747 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
5748  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)))
5749  if ( mDirtyGeos )
5750  exportWkbToGeos();
5751 
5752  if ( !mGeos )
5753  return 0;
5754 
5755  try
5756  {
5757  return fromGeosGeom( GEOSBufferWithStyle_r( geosinit.ctxt, mGeos, distance, segments, endCapStyle, joinStyle, mitreLimit ) );
5758  }
5759  CATCH_GEOS( 0 )
5760 #else
5761  return 0;
5762 #endif
5763 }
5764 
5765 QgsGeometry* QgsGeometry::offsetCurve( double distance, int segments, int joinStyle, double mitreLimit )
5766 {
5767 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
5768  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)))
5769  if ( mDirtyGeos )
5770  exportWkbToGeos();
5771 
5772  if ( !mGeos || this->type() != QGis::Line )
5773  return 0;
5774 
5775  try
5776  {
5777  return fromGeosGeom( GEOSOffsetCurve_r( geosinit.ctxt, mGeos, distance, segments, joinStyle, mitreLimit ) );
5778  }
5779  CATCH_GEOS( 0 )
5780 #else
5781  return 0;
5782 #endif
5783 }
5784 
5786 {
5787  if ( mDirtyGeos )
5788  exportWkbToGeos();
5789 
5790  if ( !mGeos )
5791  return 0;
5792 
5793  try
5794  {
5795  return fromGeosGeom( GEOSTopologyPreserveSimplify_r( geosinit.ctxt, mGeos, tolerance ) );
5796  }
5797  CATCH_GEOS( 0 )
5798 }
5799 
5801 {
5802  if ( mDirtyGeos )
5803  exportWkbToGeos();
5804 
5805  if ( !mGeos )
5806  return 0;
5807 
5808  try
5809  {
5810  return fromGeosGeom( GEOSGetCentroid_r( geosinit.ctxt, mGeos ) );
5811  }
5812  CATCH_GEOS( 0 )
5813 }
5814 
5816 {
5817  if ( mDirtyGeos )
5818  exportWkbToGeos();
5819 
5820  if ( !mGeos )
5821  return 0;
5822 
5823  try
5824  {
5825  return fromGeosGeom( GEOSPointOnSurface_r( geosinit.ctxt, mGeos ) );
5826  }
5827  CATCH_GEOS( 0 )
5828 }
5829 
5831 {
5832  if ( mDirtyGeos )
5833  exportWkbToGeos();
5834 
5835  if ( !mGeos )
5836  return 0;
5837 
5838  try
5839  {
5840  return fromGeosGeom( GEOSConvexHull_r( geosinit.ctxt, mGeos ) );
5841  }
5842  CATCH_GEOS( 0 )
5843 }
5844 
5846 {
5847 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
5848  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=2)))
5849  if ( mDirtyGeos )
5850  exportWkbToGeos();
5851 
5852  if ( !mGeos )
5853  return 0;
5854 
5855  try
5856  {
5857  return fromGeosGeom( GEOSInterpolate_r( geosinit.ctxt, mGeos, distance ) );
5858  }
5859  CATCH_GEOS( 0 )
5860 #else
5861  QgsMessageLog::logMessage( QObject::tr( "GEOS prior to 3.2 doesn't support GEOSInterpolate" ), QObject::tr( "GEOS" ) );
5862  return NULL;
5863 #endif
5864 }
5865 
5867 {
5868  if ( !geometry )
5869  return NULL;
5870 
5871  if ( mDirtyGeos )
5872  exportWkbToGeos();
5873 
5874  if ( geometry->mDirtyGeos )
5875  geometry->exportWkbToGeos();
5876 
5877  if ( !mGeos || !geometry->mGeos )
5878  return 0;
5879 
5880  try
5881  {
5882  return fromGeosGeom( GEOSIntersection_r( geosinit.ctxt, mGeos, geometry->mGeos ) );
5883  }
5884  CATCH_GEOS( 0 )
5885 }
5886 
5888 {
5889  if ( !geometry )
5890  return NULL;
5891 
5892  if ( mDirtyGeos )
5893  exportWkbToGeos();
5894 
5895  if ( geometry->mDirtyGeos )
5896  geometry->exportWkbToGeos();
5897 
5898  if ( !mGeos || !geometry->mGeos )
5899  return 0;
5900 
5901  try
5902  {
5903  GEOSGeometry* unionGeom = GEOSUnion_r( geosinit.ctxt, mGeos, geometry->mGeos );
5904  if ( !unionGeom )
5905  return 0;
5906 
5907  if ( type() == QGis::Line )
5908  {
5909  GEOSGeometry* mergedGeom = GEOSLineMerge_r( geosinit.ctxt, unionGeom );
5910  if ( mergedGeom )
5911  {
5912  GEOSGeom_destroy_r( geosinit.ctxt, unionGeom );
5913  unionGeom = mergedGeom;
5914  }
5915  }
5916  return fromGeosGeom( unionGeom );
5917  }
5918  CATCH_GEOS( new QgsGeometry( *this ) ) //return this geometry if union not possible
5919 }
5920 
5922 {
5923  if ( !geometry )
5924  return NULL;
5925 
5926  if ( mDirtyGeos )
5927  exportWkbToGeos();
5928 
5929  if ( geometry->mDirtyGeos )
5930  geometry->exportWkbToGeos();
5931 
5932  if ( !mGeos || !geometry->mGeos )
5933  return 0;
5934 
5935  try
5936  {
5937  return fromGeosGeom( GEOSDifference_r( geosinit.ctxt, mGeos, geometry->mGeos ) );
5938  }
5939  CATCH_GEOS( 0 )
5940 }
5941 
5943 {
5944  if ( !geometry )
5945  return NULL;
5946 
5947  if ( mDirtyGeos )
5948  exportWkbToGeos();
5949 
5950  if ( geometry->mDirtyGeos )
5951  geometry->exportWkbToGeos();
5952 
5953  if ( !mGeos || !geometry->mGeos )
5954  return 0;
5955 
5956  try
5957  {
5958  return fromGeosGeom( GEOSSymDifference_r( geosinit.ctxt, mGeos, geometry->mGeos ) );
5959  }
5960  CATCH_GEOS( 0 )
5961 }
5962 
5963 QList<QgsGeometry*> QgsGeometry::asGeometryCollection() const
5964 {
5965  if ( mDirtyGeos )
5966  exportWkbToGeos();
5967 
5968  if ( !mGeos )
5969  return QList<QgsGeometry*>();
5970 
5971  int type = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
5972  QgsDebugMsg( "geom type: " + QString::number( type ) );
5973 
5974  QList<QgsGeometry*> geomCollection;
5975 
5976  if ( type != GEOS_MULTIPOINT &&
5977  type != GEOS_MULTILINESTRING &&
5978  type != GEOS_MULTIPOLYGON &&
5979  type != GEOS_GEOMETRYCOLLECTION )
5980  {
5981  // we have a single-part geometry - put there a copy of this one
5982  geomCollection.append( new QgsGeometry( *this ) );
5983  return geomCollection;
5984  }
5985 
5986  int count = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
5987  QgsDebugMsg( "geom count: " + QString::number( count ) );
5988 
5989  for ( int i = 0; i < count; ++i )
5990  {
5991  const GEOSGeometry * geometry = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
5992  geomCollection.append( fromGeosGeom( GEOSGeom_clone_r( geosinit.ctxt, geometry ) ) );
5993  }
5994 
5995  return geomCollection;
5996 }
5997 
5998 QPointF QgsGeometry::asQPointF() const
5999 {
6000  QgsPoint point = asPoint();
6001  return point.toQPointF();
6002 }
6003 
6004 QPolygonF QgsGeometry::asQPolygonF() const
6005 {
6006  QPolygonF result;
6007  QgsPolyline polyline;
6008  QGis::WkbType type = wkbType();
6009  if ( type == QGis::WKBLineString || type == QGis::WKBLineString25D )
6010  {
6011  polyline = asPolyline();
6012  }
6013  else if ( type == QGis::WKBPolygon || type == QGis::WKBPolygon25D )
6014  {
6015  QgsPolygon polygon = asPolygon();
6016  if ( polygon.size() < 1 )
6017  return result;
6018  polyline = polygon.at( 0 );
6019  }
6020  else
6021  {
6022  return result;
6023  }
6024 
6025  QgsPolyline::const_iterator lineIt = polyline.constBegin();
6026  for ( ; lineIt != polyline.constEnd(); ++lineIt )
6027  {
6028  result << lineIt->toQPointF();
6029  }
6030  return result;
6031 }
6032 
6033 bool QgsGeometry::deleteRing( int ringNum, int partNum )
6034 {
6035  if ( ringNum <= 0 || partNum < 0 )
6036  return false;
6037 
6038  switch ( wkbType() )
6039  {
6040  case QGis::WKBPolygon25D:
6041  case QGis::WKBPolygon:
6042  {
6043  if ( partNum != 0 )
6044  return false;
6045 
6046  QgsPolygon polygon = asPolygon();
6047  if ( ringNum >= polygon.count() )
6048  return false;
6049 
6050  polygon.remove( ringNum );
6051 
6052  QgsGeometry* g2 = QgsGeometry::fromPolygon( polygon );
6053  *this = *g2;
6054  delete g2;
6055  return true;
6056  }
6057 
6059  case QGis::WKBMultiPolygon:
6060  {
6061  QgsMultiPolygon mpolygon = asMultiPolygon();
6062 
6063  if ( partNum >= mpolygon.count() )
6064  return false;
6065 
6066  if ( ringNum >= mpolygon[partNum].count() )
6067  return false;
6068 
6069  mpolygon[partNum].remove( ringNum );
6070 
6071  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
6072  *this = *g2;
6073  delete g2;
6074  return true;
6075  }
6076 
6077  default:
6078  return false; // only makes sense with polygons and multipolygons
6079  }
6080 }
6081 
6082 bool QgsGeometry::deletePart( int partNum )
6083 {
6084  if ( partNum < 0 )
6085  return false;
6086 
6087  switch ( wkbType() )
6088  {
6090  case QGis::WKBMultiPoint:
6091  {
6092  QgsMultiPoint mpoint = asMultiPoint();
6093 
6094  if ( partNum >= mpoint.size() || mpoint.size() == 1 )
6095  return false;
6096 
6097  mpoint.remove( partNum );
6098 
6099  QgsGeometry* g2 = QgsGeometry::fromMultiPoint( mpoint );
6100  *this = *g2;
6101  delete g2;
6102  break;
6103  }
6104 
6107  {
6109 
6110  if ( partNum >= mline.size() || mline.size() == 1 )
6111  return false;
6112 
6113  mline.remove( partNum );
6114 
6116  *this = *g2;
6117  delete g2;
6118  break;
6119  }
6120 
6122  case QGis::WKBMultiPolygon:
6123  {
6124  QgsMultiPolygon mpolygon = asMultiPolygon();
6125 
6126  if ( partNum >= mpolygon.size() || mpolygon.size() == 1 )
6127  return false;
6128 
6129  mpolygon.remove( partNum );
6130 
6131  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
6132  *this = *g2;
6133  delete g2;
6134  break;
6135  }
6136 
6137  default:
6138  // single part geometries are ignored
6139  return false;
6140  }
6141 
6142  return true;
6143 }
6144 
6147 static GEOSGeometry* _makeUnion( QList<GEOSGeometry*> geoms )
6148 {
6149 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && (((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)) || (GEOS_VERSION_MAJOR>3))
6150  GEOSGeometry* geomCollection = 0;
6151  geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geoms.toVector() );
6152  GEOSGeometry* geomUnion = GEOSUnaryUnion_r( geosinit.ctxt, geomCollection );
6153  GEOSGeom_destroy_r( geosinit.ctxt, geomCollection );
6154  return geomUnion;
6155 #else
6156  GEOSGeometry* geomCollection = geoms.takeFirst();
6157 
6158  while ( !geoms.isEmpty() )
6159  {
6160  GEOSGeometry* g = geoms.takeFirst();
6161  GEOSGeometry* geomCollectionNew = GEOSUnion_r( geosinit.ctxt, geomCollection, g );
6162  GEOSGeom_destroy_r( geosinit.ctxt, geomCollection );
6163  GEOSGeom_destroy_r( geosinit.ctxt, g );
6164  geomCollection = geomCollectionNew;
6165  }
6166 
6167  return geomCollection;
6168 #endif
6169 }
6170 
6171 int QgsGeometry::avoidIntersections( QMap<QgsVectorLayer*, QSet< QgsFeatureId > > ignoreFeatures )
6172 {
6173  int returnValue = 0;
6174 
6175  //check if g has polygon type
6176  if ( type() != QGis::Polygon )
6177  return 1;
6178 
6179  QGis::WkbType geomTypeBeforeModification = wkbType();
6180 
6181  //read avoid intersections list from project properties
6182  bool listReadOk;
6183  QStringList avoidIntersectionsList = QgsProject::instance()->readListEntry( "Digitizing", "/AvoidIntersectionsList", QStringList(), &listReadOk );
6184  if ( !listReadOk )
6185  return true; //no intersections stored in project does not mean error
6186 
6187  QList<GEOSGeometry*> nearGeometries;
6188 
6189  //go through list, convert each layer to vector layer and call QgsVectorLayer::removePolygonIntersections for each
6190  QgsVectorLayer* currentLayer = 0;
6191  QStringList::const_iterator aIt = avoidIntersectionsList.constBegin();
6192  for ( ; aIt != avoidIntersectionsList.constEnd(); ++aIt )
6193  {
6194  currentLayer = dynamic_cast<QgsVectorLayer*>( QgsMapLayerRegistry::instance()->mapLayer( *aIt ) );
6195  if ( currentLayer )
6196  {
6197  QgsFeatureIds ignoreIds;
6198  QMap<QgsVectorLayer*, QSet<qint64> >::const_iterator ignoreIt = ignoreFeatures.find( currentLayer );
6199  if ( ignoreIt != ignoreFeatures.constEnd() )
6200  ignoreIds = ignoreIt.value();
6201 
6204  .setSubsetOfAttributes( QgsAttributeList() ) );
6205  QgsFeature f;
6206  while ( fi.nextFeature( f ) )
6207  {
6208  if ( ignoreIds.contains( f.id() ) )
6209  continue;
6210 
6211  if ( !f.geometry() )
6212  continue;
6213 
6214  nearGeometries << GEOSGeom_clone_r( geosinit.ctxt, f.geometry()->asGeos() );
6215  }
6216  }
6217  }
6218 
6219  if ( nearGeometries.isEmpty() )
6220  return 0;
6221 
6222  GEOSGeometry* nearGeometriesUnion = 0;
6223  GEOSGeometry* geomWithoutIntersections = 0;
6224  try
6225  {
6226  nearGeometriesUnion = _makeUnion( nearGeometries );
6227  geomWithoutIntersections = GEOSDifference_r( geosinit.ctxt, asGeos(), nearGeometriesUnion );
6228 
6229  fromGeos( geomWithoutIntersections );
6230 
6231  GEOSGeom_destroy_r( geosinit.ctxt, nearGeometriesUnion );
6232  }
6233  catch ( GEOSException &e )
6234  {
6235  if ( nearGeometriesUnion )
6236  GEOSGeom_destroy_r( geosinit.ctxt, nearGeometriesUnion );
6237  if ( geomWithoutIntersections )
6238  GEOSGeom_destroy_r( geosinit.ctxt, geomWithoutIntersections );
6239 
6240  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
6241  return 3;
6242  }
6243 
6244  //make sure the geometry still has the same type (e.g. no change from polygon to multipolygon)
6245  if ( wkbType() != geomTypeBeforeModification )
6246  return 2;
6247 
6248  return returnValue;
6249 }
6250 
6251 void QgsGeometry::validateGeometry( QList<Error> &errors )
6252 {
6254 }
6255 
6257 {
6258  try
6259  {
6260  const GEOSGeometry *g = asGeos();
6261 
6262  if ( !g )
6263  return false;
6264 
6265  return GEOSisValid_r( geosinit.ctxt, g );
6266  }
6267  catch ( GEOSException &e )
6268  {
6269  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
6270  return false;
6271  }
6272 }
6273 
6275 {
6276  return geosRelOp( GEOSEquals_r, this, &g );
6277 }
6278 
6280 {
6281  try
6282  {
6283  const GEOSGeometry *g = asGeos();
6284 
6285  if ( !g )
6286  return false;
6287 
6288  return GEOSisEmpty_r( geosinit.ctxt, g );
6289  }
6290  catch ( GEOSException &e )
6291  {
6292  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
6293  return false;
6294  }
6295 }
6296 
6297 double QgsGeometry::leftOf( double x, double y, double& x1, double& y1, double& x2, double& y2 )
6298 {
6299  double f1 = x - x1;
6300  double f2 = y2 - y1;
6301  double f3 = y - y1;
6302  double f4 = x2 - x1;
6303  return f1*f2 - f3*f4;
6304 }
6305 
6306 QgsGeometry* QgsGeometry::convertToPoint( bool destMultipart )
6307 {
6308  switch ( type() )
6309  {
6310  case QGis::Point:
6311  {
6312  bool srcIsMultipart = isMultipart();
6313 
6314  if (( destMultipart && srcIsMultipart ) ||
6315  ( !destMultipart && !srcIsMultipart ) )
6316  {
6317  // return a copy of the same geom
6318  return new QgsGeometry( *this );
6319  }
6320  if ( destMultipart )
6321  {
6322  // layer is multipart => make a multipoint with a single point
6323  return fromMultiPoint( QgsMultiPoint() << asPoint() );
6324  }
6325  else
6326  {
6327  // destination is singlepart => make a single part if possible
6328  QgsMultiPoint multiPoint = asMultiPoint();
6329  if ( multiPoint.count() == 1 )
6330  {
6331  return fromPoint( multiPoint[0] );
6332  }
6333  }
6334  return 0;
6335  }
6336 
6337  case QGis::Line:
6338  {
6339  // only possible if destination is multipart
6340  if ( !destMultipart )
6341  return 0;
6342 
6343  // input geometry is multipart
6344  if ( isMultipart() )
6345  {
6346  QgsMultiPolyline multiLine = asMultiPolyline();
6347  QgsMultiPoint multiPoint;
6348  for ( QgsMultiPolyline::const_iterator multiLineIt = multiLine.constBegin(); multiLineIt != multiLine.constEnd(); ++multiLineIt )
6349  for ( QgsPolyline::const_iterator lineIt = ( *multiLineIt ).constBegin(); lineIt != ( *multiLineIt ).constEnd(); ++lineIt )
6350  multiPoint << *lineIt;
6351  return fromMultiPoint( multiPoint );
6352  }
6353  // input geometry is not multipart: copy directly the line into a multipoint
6354  else
6355  {
6356  QgsPolyline line = asPolyline();
6357  if ( !line.isEmpty() )
6358  return fromMultiPoint( line );
6359  }
6360  return 0;
6361  }
6362 
6363  case QGis::Polygon:
6364  {
6365  // can only transform if destination is multipoint
6366  if ( !destMultipart )
6367  return 0;
6368 
6369  // input geometry is multipart: make a multipoint from multipolygon
6370  if ( isMultipart() )
6371  {
6372  QgsMultiPolygon multiPolygon = asMultiPolygon();
6373  QgsMultiPoint multiPoint;
6374  for ( QgsMultiPolygon::const_iterator polygonIt = multiPolygon.constBegin(); polygonIt != multiPolygon.constEnd(); ++polygonIt )
6375  for ( QgsMultiPolyline::const_iterator multiLineIt = ( *polygonIt ).constBegin(); multiLineIt != ( *polygonIt ).constEnd(); ++multiLineIt )
6376  for ( QgsPolyline::const_iterator lineIt = ( *multiLineIt ).constBegin(); lineIt != ( *multiLineIt ).constEnd(); ++lineIt )
6377  multiPoint << *lineIt;
6378  return fromMultiPoint( multiPoint );
6379  }
6380  // input geometry is not multipart: make a multipoint from polygon
6381  else
6382  {
6383  QgsPolygon polygon = asPolygon();
6384  QgsMultiPoint multiPoint;
6385  for ( QgsMultiPolyline::const_iterator multiLineIt = polygon.constBegin(); multiLineIt != polygon.constEnd(); ++multiLineIt )
6386  for ( QgsPolyline::const_iterator lineIt = ( *multiLineIt ).constBegin(); lineIt != ( *multiLineIt ).constEnd(); ++lineIt )
6387  multiPoint << *lineIt;
6388  return fromMultiPoint( multiPoint );
6389  }
6390  }
6391 
6392  default:
6393  return 0;
6394  }
6395 }
6396 
6397 QgsGeometry* QgsGeometry::convertToLine( bool destMultipart )
6398 {
6399  switch ( type() )
6400  {
6401  case QGis::Point:
6402  {
6403  if ( !isMultipart() )
6404  return 0;
6405 
6406  QgsMultiPoint multiPoint = asMultiPoint();
6407  if ( multiPoint.count() < 2 )
6408  return 0;
6409 
6410  if ( destMultipart )
6411  return fromMultiPolyline( QgsMultiPolyline() << multiPoint );
6412  else
6413  return fromPolyline( multiPoint );
6414  }
6415 
6416  case QGis::Line:
6417  {
6418  bool srcIsMultipart = isMultipart();
6419 
6420  if (( destMultipart && srcIsMultipart ) ||
6421  ( !destMultipart && ! srcIsMultipart ) )
6422  {
6423  // return a copy of the same geom
6424  return new QgsGeometry( *this );
6425  }
6426  if ( destMultipart )
6427  {
6428  // destination is multipart => makes a multipoint with a single line
6429  QgsPolyline line = asPolyline();
6430  if ( !line.isEmpty() )
6431  return fromMultiPolyline( QgsMultiPolyline() << line );
6432  }
6433  else
6434  {
6435  // destination is singlepart => make a single part if possible
6436  QgsMultiPolyline multiLine = asMultiPolyline();
6437  if ( multiLine.count() == 1 )
6438  return fromPolyline( multiLine[0] );
6439  }
6440  return 0;
6441  }
6442 
6443  case QGis::Polygon:
6444  {
6445  // input geometry is multipolygon
6446  if ( isMultipart() )
6447  {
6448  QgsMultiPolygon multiPolygon = asMultiPolygon();
6449  QgsMultiPolyline multiLine;
6450  for ( QgsMultiPolygon::const_iterator polygonIt = multiPolygon.constBegin(); polygonIt != multiPolygon.constEnd(); ++polygonIt )
6451  for ( QgsMultiPolyline::const_iterator multiLineIt = ( *polygonIt ).constBegin(); multiLineIt != ( *polygonIt ).constEnd(); ++multiLineIt )
6452  multiLine << *multiLineIt;
6453 
6454  if ( destMultipart )
6455  {
6456  // destination is multipart
6457  return fromMultiPolyline( multiLine );
6458  }
6459  else if ( multiLine.count() == 1 )
6460  {
6461  // destination is singlepart => make a single part if possible
6462  return fromPolyline( multiLine[0] );
6463  }
6464  }
6465  // input geometry is single polygon
6466  else
6467  {
6468  QgsPolygon polygon = asPolygon();
6469  // if polygon has rings
6470  if ( polygon.count() > 1 )
6471  {
6472  // cannot fit a polygon with rings in a single line layer
6473  // TODO: would it be better to remove rings?
6474  if ( destMultipart )
6475  {
6476  QgsPolygon polygon = asPolygon();
6477  QgsMultiPolyline multiLine;
6478  for ( QgsMultiPolyline::const_iterator multiLineIt = polygon.constBegin(); multiLineIt != polygon.constEnd(); ++multiLineIt )
6479  multiLine << *multiLineIt;
6480  return fromMultiPolyline( multiLine );
6481  }
6482  }
6483  // no rings
6484  else if ( polygon.count() == 1 )
6485  {
6486  if ( destMultipart )
6487  {
6488  return fromMultiPolyline( polygon );
6489  }
6490  else
6491  {
6492  return fromPolyline( polygon[0] );
6493  }
6494  }
6495  }
6496  return 0;
6497  }
6498 
6499  default:
6500  return 0;
6501  }
6502 }
6503 
6504 QgsGeometry* QgsGeometry::convertToPolygon( bool destMultipart )
6505 {
6506  switch ( type() )
6507  {
6508  case QGis::Point:
6509  {
6510  if ( !isMultipart() )
6511  return 0;
6512 
6513  QgsMultiPoint multiPoint = asMultiPoint();
6514  if ( multiPoint.count() < 3 )
6515  return 0;
6516 
6517  if ( multiPoint.last() != multiPoint.first() )
6518  multiPoint << multiPoint.first();
6519 
6520  QgsPolygon polygon = QgsPolygon() << multiPoint;
6521  if ( destMultipart )
6522  return fromMultiPolygon( QgsMultiPolygon() << polygon );
6523  else
6524  return fromPolygon( polygon );
6525  }
6526 
6527  case QGis::Line:
6528  {
6529  // input geometry is multiline
6530  if ( isMultipart() )
6531  {
6532  QgsMultiPolyline multiLine = asMultiPolyline();
6533  QgsMultiPolygon multiPolygon;
6534  for ( QgsMultiPolyline::iterator multiLineIt = multiLine.begin(); multiLineIt != multiLine.end(); ++multiLineIt )
6535  {
6536  // do not create polygon for a 1 segment line
6537  if (( *multiLineIt ).count() < 3 )
6538  return 0;
6539  if (( *multiLineIt ).count() == 3 && ( *multiLineIt ).first() == ( *multiLineIt ).last() )
6540  return 0;
6541 
6542  // add closing node
6543  if (( *multiLineIt ).first() != ( *multiLineIt ).last() )
6544  *multiLineIt << ( *multiLineIt ).first();
6545  multiPolygon << ( QgsPolygon() << *multiLineIt );
6546  }
6547  // check that polygons were inserted
6548  if ( !multiPolygon.isEmpty() )
6549  {
6550  if ( destMultipart )
6551  {
6552  return fromMultiPolygon( multiPolygon );
6553  }
6554  else if ( multiPolygon.count() == 1 )
6555  {
6556  // destination is singlepart => make a single part if possible
6557  return fromPolygon( multiPolygon[0] );
6558  }
6559  }
6560  }
6561  // input geometry is single line
6562  else
6563  {
6564  QgsPolyline line = asPolyline();
6565 
6566  // do not create polygon for a 1 segment line
6567  if ( line.count() < 3 )
6568  return 0;
6569  if ( line.count() == 3 && line.first() == line.last() )
6570  return 0;
6571 
6572  // add closing node
6573  if ( line.first() != line.last() )
6574  line << line.first();
6575 
6576  // destination is multipart
6577  if ( destMultipart )
6578  {
6579  return fromMultiPolygon( QgsMultiPolygon() << ( QgsPolygon() << line ) );
6580  }
6581  else
6582  {
6583  return fromPolygon( QgsPolygon() << line );
6584  }
6585  }
6586  return 0;
6587  }
6588 
6589  case QGis::Polygon:
6590  {
6591  bool srcIsMultipart = isMultipart();
6592 
6593  if (( destMultipart && srcIsMultipart ) ||
6594  ( !destMultipart && ! srcIsMultipart ) )
6595  {
6596  // return a copy of the same geom
6597  return new QgsGeometry( *this );
6598  }
6599  if ( destMultipart )
6600  {
6601  // destination is multipart => makes a multipoint with a single polygon
6602  QgsPolygon polygon = asPolygon();
6603  if ( !polygon.isEmpty() )
6604  return fromMultiPolygon( QgsMultiPolygon() << polygon );
6605  }
6606  else
6607  {
6608  QgsMultiPolygon multiPolygon = asMultiPolygon();
6609  if ( multiPolygon.count() == 1 )
6610  {
6611  // destination is singlepart => make a single part if possible
6612  return fromPolygon( multiPolygon[0] );
6613  }
6614  }
6615  return 0;
6616  }
6617 
6618  default:
6619  return 0;
6620  }
6621 }
6622 
6623 QgsGeometry *QgsGeometry::unaryUnion( const QList<QgsGeometry*> &geometryList )
6624 {
6625  QList<GEOSGeometry*> geoms;
6626  foreach ( QgsGeometry* g, geometryList )
6627  {
6628  geoms.append( GEOSGeom_clone_r( geosinit.ctxt, g->asGeos() ) );
6629  }
6630  GEOSGeometry *geomUnion = _makeUnion( geoms );
6631  QgsGeometry *ret = new QgsGeometry();
6632  ret->fromGeos( geomUnion );
6633  return ret;
6634 }
QgsFeatureId id() const
Get the feature id for this feature.
Definition: qgsfeature.cpp:100
static QgsGeometry * fromQPolygonF(const QPolygonF &polygon)
Construct geometry from a QPolygonF.
Wrapper for iterator of features from vector data provider or vector layer.
GEOSException(QString theMsg)
Definition: qgsgeometry.cpp:55
static unsigned index
A rectangle specified with double values.
Definition: qgsrectangle.h:35
static void validateGeometry(QgsGeometry *g, QList< QgsGeometry::Error > &errors)
Validate geometry and produce a list of geometry errors.
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
void adjacentVertices(int atVertex, int &beforeVertex, int &afterVertex)
Returns the indexes of the vertices before and after the given vertex index.
int makeDifference(QgsGeometry *other)
Changes this geometry such that it does not intersect the other geometry.
double length()
get length of geometry using GEOS
int reshapeGeometry(const QList< QgsPoint > &reshapeWithLine)
Replaces a part of this geometry with another line.
Use exact geometry intersection (slower) instead of bounding boxes.
QPolygonF asQPolygonF() const
Return contents of the geometry as a QPolygonF.
QgsGeometry * combine(QgsGeometry *geometry)
Returns a geometry representing all the points in this geometry and other (a union geometry operation...
size_t wkbSize() const
Returns the size of the WKB in asWkb().
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:188
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeature.h:317
QgsGeometry * convertToType(QGis::GeometryType destType, bool destMultipart=false)
try to convert the geometry to the requested type
QgsMultiPolyline asMultiPolyline() const
return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
QVector< QgsPoint > QgsPolyline
polyline is represented as a vector of points
Definition: qgsgeometry.h:33
QPointF asQPointF() const
Return contents of the geometry as a QPointF if wkbType is WKBPoint, otherwise returns a null QPointF...
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
Definition: qgsfeature.cpp:112
QgsPolygon asPolygon() const
return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list ...
QList< QgsGeometry * > asGeometryCollection() const
return contents of the geometry as a list of geometries
static QgsMapLayerRegistry * instance()
Definition: qgssingleton.h:23
static GEOSGeometry * createGeosCollection(int typeId, QVector< GEOSGeometry * > geoms)
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON)
Searches for the closest segment of geometry to the given point.
bool moveVertex(double x, double y, int atVertex)
Moves the vertex at the given position number and item (first number is index 0) to the given coordin...
QGis::GeometryType type() const
Returns type of the vector.
double closestVertexWithContext(const QgsPoint &point, int &atVertex)
Searches for the closest vertex in this geometry to the given point.
bool crosses(const QgsGeometry *geometry) const
Test for if geometry crosses another (uses GEOS)
int addPart(const QList< QgsPoint > &points, QGis::GeometryType geomType=QGis::UnknownGeometry)
Adds a new island polygon to a multipolygon feature.
GeometryType
Definition: qgis.h:155
QgsGeometry()
Constructor.
QString what()
Definition: qgsgeometry.cpp:80
static GEOSGeometry * createGeosPolygon(const QVector< GEOSGeometry * > &rings)
WkbType
Used for symbology operations.
Definition: qgis.h:53
The feature class encapsulates a single feature including its id, geometry and a list of field/values...
Definition: qgsfeature.h:113
double area()
get area of geometry using GEOS
bool insertVertex(double x, double y, int beforeVertex)
Insert a new vertex before the given vertex index, ring and item (first number is index 0) If the req...
QgsPoint closestVertex(const QgsPoint &point, int &atVertex, int &beforeVertex, int &afterVertex, double &sqrDist)
Returns the vertex closest to the given point, the corresponding vertex index, squared distance snap ...
QgsPoint vertexAt(int atVertex)
Returns coordinates of a vertex.
static endian_t endian()
Returns whether this machine uses big or little endian.
QgsGeometry & operator=(QgsGeometry const &rhs)
assignments will prompt a deep copy of the object
double sqrDist(double x, double y) const
Returns the squared distance between this point and x,y.
Definition: qgspoint.cpp:333
GEOSException(const GEOSException &rhs)
Definition: qgsgeometry.cpp:69
QgsGeometry * difference(QgsGeometry *geometry)
Returns a geometry representing the points making up this geometry that do not make up other...
double x() const
Definition: qgspoint.h:126
bool deleteRing(int ringNum, int partNum=0)
delete a ring in polygon or multipolygon.
const GEOSGeometry * asGeos() const
Returns a geos geometry.
static void throwGEOSException(const char *fmt,...)
Definition: qgsgeometry.cpp:92
QgsGeometry * centroid()
Returns the center of mass of a geometry.
static void logMessage(QString message, QString tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
QgsMultiPolygon asMultiPolygon() const
return contents of the geometry as a multi polygon if wkbType is WKBMultiPolygon, otherwise an empty ...
QgsGeometry * offsetCurve(double distance, int segments, int joinStyle, double mitreLimit)
Returns an offset line at a given distance and side from an input line.
bool deletePart(int partNum)
delete part identified by the part number
bool isGeosValid()
check validity using GEOS
static WkbType flatType(WkbType type)
Definition: qgis.h:99
int splitGeometry(const QList< QgsPoint > &splitLine, QList< QgsGeometry * > &newGeometries, bool topological, QList< QgsPoint > &topologyTestPoints)
Splits this geometry according to a given line.
void transformInPlace(double &x, double &y, double &z, TransformDirection direction=ForwardTransform) const
bool contains(const QgsPoint *p) const
Test for containment of a point (uses GEOS)
QStringList readListEntry(const QString &scope, const QString &key, QStringList def=QStringList(), bool *ok=0) const
key value accessors
QgsGeometry * buffer(double distance, int segments)
Returns a buffer region around this geometry having the given width and with a specified number of se...
double yMinimum() const
Get the y minimum value (bottom side of rectangle)
Definition: qgsrectangle.h:193
static GEOSCoordSequence * createGeosCoordSequence(const QgsPolyline &points)
bool isGeosEmpty()
check if geometry is empty using GEOS
double xMaximum() const
Get the x maximum value (right side of rectangle)
Definition: qgsrectangle.h:178
#define QgsDebugMsgLevel(str, level)
Definition: qgslogger.h:34
QgsGeometry * convexHull()
Returns the smallest convex polygon that contains all the points in the geometry. ...
bool overlaps(const QgsGeometry *geometry) const
Test for if geometry overlaps another (uses GEOS)
QgsGeometry * simplify(double tolerance)
Returns a simplified version of this geometry using a specified tolerance value.
bool deleteVertex(int atVertex)
Deletes the vertex at the given position number and item (first number is index 0) Returns false if a...
double sqrDistToSegment(double x1, double y1, double x2, double y2, QgsPoint &minDistPoint, double epsilon=DEFAULT_SEGMENT_EPSILON) const
Returns the minimum distance between this point and a segment.
Definition: qgspoint.cpp:414
int avoidIntersections(QMap< QgsVectorLayer *, QSet< QgsFeatureId > > ignoreFeatures=(QMap< QgsVectorLayer *, QSet< QgsFeatureId > >()))
Modifies geometry to avoid intersections with the layers specified in project properties.
This class wraps a request for features to a vector layer (or directly its vector data provider)...
QVector< QgsPolygon > QgsMultiPolygon
a collection of QgsPolygons that share a common collection of attributes
Definition: qgsgeometry.h:48
static GEOSGeometry * createGeosLinearRing(const QgsPolyline &polyline)
QList< int > QgsAttributeList
static GEOSGeometry * createGeosLineString(const QgsPolyline &polyline)
QVector< QgsPoint > QgsMultiPoint
a collection of QgsPoints that share a common collection of attributes
Definition: qgsgeometry.h:42
void fromGeos(GEOSGeometry *geos)
Set the geometry, feeding in a geometry in GEOS format.
QString toString() const
String representation of the point (x,y)
Definition: qgspoint.cpp:126
static GEOSGeometry * createGeosPoint(const double x, const double y)
QGis::WkbType wkbType() const
Returns type of wkb (point / linestring / polygon etc.)
GEOSContextHandle_t ctxt
QVector< QgsPolyline > QgsPolygon
polygon: first item of the list is outer ring, inner rings (if any) start from second item ...
Definition: qgsgeometry.h:39
void validateGeometry(QList< Error > &errors)
Validate geometry and produce a list of geometry errors.