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