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