QGIS API Documentation  2.2.0-Valmiera
 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  // TODO: implement with GEOS
4429  if ( mDirtyWkb )
4430  {
4431  exportGeosToWkb();
4432  }
4433 
4434  if ( !mGeometry )
4435  {
4436  return false;
4437  }
4438 
4439  QGis::WkbType geomType = wkbType();
4440 
4441  if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D ||
4442  geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D ||
4443  geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown )
4444  {
4445  return false; //no need to convert
4446  }
4447 
4448  size_t newGeomSize = mGeometrySize + 1 + 2 * sizeof( int ); //endian: 1, multitype: sizeof(int), number of geometries: sizeof(int)
4449  unsigned char* newGeometry = new unsigned char[newGeomSize];
4450 
4451  //copy endian
4452  char byteOrder = QgsApplication::endian();
4453 
4454  QgsWkbPtr wkbPtr( newGeometry );
4455  wkbPtr << byteOrder;
4456 
4457  //copy wkbtype
4458  //todo
4459  QGis::WkbType newMultiType;
4460  switch ( geomType )
4461  {
4462  case QGis::WKBPoint:
4463  newMultiType = QGis::WKBMultiPoint;
4464  break;
4465  case QGis::WKBPoint25D:
4466  newMultiType = QGis::WKBMultiPoint25D;
4467  break;
4468  case QGis::WKBLineString:
4469  newMultiType = QGis::WKBMultiLineString;
4470  break;
4472  newMultiType = QGis::WKBMultiLineString25D;
4473  break;
4474  case QGis::WKBPolygon:
4475  newMultiType = QGis::WKBMultiPolygon;
4476  break;
4477  case QGis::WKBPolygon25D:
4478  newMultiType = QGis::WKBMultiPolygon25D;
4479  break;
4480  default:
4481  delete [] newGeometry;
4482  return false;
4483  }
4484 
4485  wkbPtr << newMultiType << 1;
4486 
4487  //copy the existing single geometry
4488  memcpy( wkbPtr, mGeometry, mGeometrySize );
4489 
4490  delete [] mGeometry;
4491  mGeometry = newGeometry;
4492  mGeometrySize = newGeomSize;
4493  mDirtyGeos = true;
4494  return true;
4495 }
4496 
4497 void QgsGeometry::translateVertex( QgsWkbPtr &wkbPtr, double dx, double dy, bool hasZValue )
4498 {
4499  double x, y, translated_x, translated_y;
4500 
4501  //x-coordinate
4502  memcpy( &x, wkbPtr, sizeof( double ) );
4503  translated_x = x + dx;
4504  memcpy( wkbPtr, &translated_x, sizeof( double ) );
4505  wkbPtr += sizeof( double );
4506 
4507  //y-coordinate
4508  memcpy( &y, wkbPtr, sizeof( double ) );
4509  translated_y = y + dy;
4510  memcpy( wkbPtr, &translated_y, sizeof( double ) );
4511  wkbPtr += sizeof( double );
4512 
4513  if ( hasZValue )
4514  wkbPtr += sizeof( double );
4515 }
4516 
4517 void QgsGeometry::transformVertex( QgsWkbPtr &wkbPtr, const QgsCoordinateTransform& ct, bool hasZValue )
4518 {
4519  double x, y, z;
4520 
4521  memcpy( &x, wkbPtr, sizeof( double ) );
4522  memcpy( &y, wkbPtr + sizeof( double ), sizeof( double ) );
4523  z = 0.0; // Ignore Z for now.
4524 
4525  ct.transformInPlace( x, y, z );
4526 
4527  // new coordinate
4528  wkbPtr << x << y;
4529  if ( hasZValue )
4530  wkbPtr += sizeof( double );
4531 
4532 }
4533 
4534 int QgsGeometry::splitLinearGeometry( GEOSGeometry *splitLine, QList<QgsGeometry*>& newGeometries )
4535 {
4536  if ( !splitLine )
4537  return 2;
4538 
4539  if ( mDirtyGeos )
4540  exportWkbToGeos();
4541 
4542  if ( !mGeos )
4543  return 5;
4544 
4545  //first test if linestring intersects geometry. If not, return straight away
4546  if ( !GEOSIntersects( splitLine, mGeos ) )
4547  return 1;
4548 
4549  //check that split line has no linear intersection
4550  int linearIntersect = GEOSRelatePattern( mGeos, splitLine, "1********" );
4551  if ( linearIntersect > 0 )
4552  return 3;
4553 
4554  GEOSGeometry* splitGeom = GEOSDifference( mGeos, splitLine );
4555  QVector<GEOSGeometry*> lineGeoms;
4556 
4557  int splitType = GEOSGeomTypeId( splitGeom );
4558  if ( splitType == GEOS_MULTILINESTRING )
4559  {
4560  int nGeoms = GEOSGetNumGeometries( splitGeom );
4561  for ( int i = 0; i < nGeoms; ++i )
4562  lineGeoms << GEOSGeom_clone( GEOSGetGeometryN( splitGeom, i ) );
4563 
4564  }
4565  else
4566  {
4567  lineGeoms << GEOSGeom_clone( splitGeom );
4568  }
4569 
4570  mergeGeometriesMultiTypeSplit( lineGeoms );
4571 
4572  if ( lineGeoms.size() > 0 )
4573  {
4574  GEOSGeom_destroy( mGeos );
4575  mGeos = lineGeoms[0];
4576  mDirtyWkb = true;
4577  }
4578 
4579  for ( int i = 1; i < lineGeoms.size(); ++i )
4580  {
4581  newGeometries << fromGeosGeom( lineGeoms[i] );
4582  }
4583 
4584  GEOSGeom_destroy( splitGeom );
4585  return 0;
4586 }
4587 
4588 int QgsGeometry::splitPolygonGeometry( GEOSGeometry* splitLine, QList<QgsGeometry*>& newGeometries )
4589 {
4590  if ( !splitLine )
4591  return 2;
4592 
4593  if ( mDirtyGeos )
4594  exportWkbToGeos();
4595 
4596  if ( !mGeos )
4597  return 5;
4598 
4599  //first test if linestring intersects geometry. If not, return straight away
4600  if ( !GEOSIntersects( splitLine, mGeos ) )
4601  return 1;
4602 
4603  //first union all the polygon rings together (to get them noded, see JTS developer guide)
4604  GEOSGeometry *nodedGeometry = nodeGeometries( splitLine, mGeos );
4605  if ( !nodedGeometry )
4606  return 2; //an error occured during noding
4607 
4608  GEOSGeometry *polygons = GEOSPolygonize( &nodedGeometry, 1 );
4609  if ( !polygons || numberOfGeometries( polygons ) == 0 )
4610  {
4611  if ( polygons )
4612  GEOSGeom_destroy( polygons );
4613 
4614  GEOSGeom_destroy( nodedGeometry );
4615 
4616  return 4;
4617  }
4618 
4619  GEOSGeom_destroy( nodedGeometry );
4620 
4621  //test every polygon if contained in original geometry
4622  //include in result if yes
4623  QVector<GEOSGeometry*> testedGeometries;
4624  GEOSGeometry *intersectGeometry = 0;
4625 
4626  //ratio intersect geometry / geometry. This should be close to 1
4627  //if the polygon belongs to the input geometry
4628 
4629  for ( int i = 0; i < numberOfGeometries( polygons ); i++ )
4630  {
4631  const GEOSGeometry *polygon = GEOSGetGeometryN( polygons, i );
4632  intersectGeometry = GEOSIntersection( mGeos, polygon );
4633  if ( !intersectGeometry )
4634  {
4635  QgsDebugMsg( "intersectGeometry is NULL" );
4636  continue;
4637  }
4638 
4639  double intersectionArea;
4640  GEOSArea( intersectGeometry, &intersectionArea );
4641 
4642  double polygonArea;
4643  GEOSArea( polygon, &polygonArea );
4644 
4645  const double areaRatio = intersectionArea / polygonArea;
4646  if ( areaRatio > 0.99 && areaRatio < 1.01 )
4647  testedGeometries << GEOSGeom_clone( polygon );
4648 
4649  GEOSGeom_destroy( intersectGeometry );
4650  }
4651 
4652  bool splitDone = true;
4653  int nGeometriesThis = numberOfGeometries( mGeos ); //original number of geometries
4654  if ( testedGeometries.size() == nGeometriesThis )
4655  {
4656  splitDone = false;
4657  }
4658 
4659  mergeGeometriesMultiTypeSplit( testedGeometries );
4660 
4661  //no split done, preserve original geometry
4662  if ( !splitDone )
4663  {
4664  for ( int i = 0; i < testedGeometries.size(); ++i )
4665  {
4666  GEOSGeom_destroy( testedGeometries[i] );
4667  }
4668  return 1;
4669  }
4670  else if ( testedGeometries.size() > 0 ) //split successfull
4671  {
4672  GEOSGeom_destroy( mGeos );
4673  mGeos = testedGeometries[0];
4674  mDirtyWkb = true;
4675  }
4676 
4677  int i;
4678  for ( i = 1; i < testedGeometries.size() && GEOSisValid( testedGeometries[i] ); ++i )
4679  ;
4680 
4681  if ( i < testedGeometries.size() )
4682  {
4683  for ( i = 0; i < testedGeometries.size(); ++i )
4684  GEOSGeom_destroy( testedGeometries[i] );
4685 
4686  return 3;
4687  }
4688 
4689  for ( i = 1; i < testedGeometries.size(); ++i )
4690  newGeometries << fromGeosGeom( testedGeometries[i] );
4691 
4692  GEOSGeom_destroy( polygons );
4693  return 0;
4694 }
4695 
4696 GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos )
4697 {
4698  //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line
4699  int nIntersections = 0;
4700  int lastIntersectingRing = -2;
4701  const GEOSGeometry* lastIntersectingGeom = 0;
4702 
4703  int nRings = GEOSGetNumInteriorRings( polygon );
4704  if ( nRings < 0 )
4705  return 0;
4706 
4707  //does outer ring intersect?
4708  const GEOSGeometry* outerRing = GEOSGetExteriorRing( polygon );
4709  if ( GEOSIntersects( outerRing, reshapeLineGeos ) == 1 )
4710  {
4711  ++nIntersections;
4712  lastIntersectingRing = -1;
4713  lastIntersectingGeom = outerRing;
4714  }
4715 
4716  //do inner rings intersect?
4717  const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings];
4718 
4719  try
4720  {
4721  for ( int i = 0; i < nRings; ++i )
4722  {
4723  innerRings[i] = GEOSGetInteriorRingN( polygon, i );
4724  if ( GEOSIntersects( innerRings[i], reshapeLineGeos ) == 1 )
4725  {
4726  ++nIntersections;
4727  lastIntersectingRing = i;
4728  lastIntersectingGeom = innerRings[i];
4729  }
4730  }
4731  }
4732  catch ( GEOSException &e )
4733  {
4734  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
4735  nIntersections = 0;
4736  }
4737 
4738  if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring
4739  {
4740  delete [] innerRings;
4741  return 0;
4742  }
4743 
4744  //we have one intersecting ring, let's try to reshape it
4745  GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos );
4746  if ( !reshapeResult )
4747  {
4748  delete [] innerRings;
4749  return 0;
4750  }
4751 
4752  //if reshaping took place, we need to reassemble the polygon and its rings
4753  GEOSGeometry* newRing = 0;
4754  const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq( reshapeResult );
4755  GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone( reshapeSequence );
4756 
4757  GEOSGeom_destroy( reshapeResult );
4758 
4759  newRing = GEOSGeom_createLinearRing( newCoordSequence );
4760  if ( !newRing )
4761  {
4762  delete [] innerRings;
4763  return 0;
4764  }
4765 
4766  GEOSGeometry* newOuterRing = 0;
4767  if ( lastIntersectingRing == -1 )
4768  newOuterRing = newRing;
4769  else
4770  newOuterRing = GEOSGeom_clone( outerRing );
4771 
4772  //check if all the rings are still inside the outer boundary
4773  QList<GEOSGeometry*> ringList;
4774  if ( nRings > 0 )
4775  {
4776  GEOSGeometry* outerRingPoly = GEOSGeom_createPolygon( GEOSGeom_clone( newOuterRing ), 0, 0 );
4777  if ( outerRingPoly )
4778  {
4779  GEOSGeometry* currentRing = 0;
4780  for ( int i = 0; i < nRings; ++i )
4781  {
4782  if ( lastIntersectingRing == i )
4783  currentRing = newRing;
4784  else
4785  currentRing = GEOSGeom_clone( innerRings[i] );
4786 
4787  //possibly a ring is no longer contained in the result polygon after reshape
4788  if ( GEOSContains( outerRingPoly, currentRing ) == 1 )
4789  ringList.push_back( currentRing );
4790  else
4791  GEOSGeom_destroy( currentRing );
4792  }
4793  }
4794  GEOSGeom_destroy( outerRingPoly );
4795  }
4796 
4797  GEOSGeometry** newInnerRings = new GEOSGeometry*[ringList.size()];
4798  for ( int i = 0; i < ringList.size(); ++i )
4799  newInnerRings[i] = ringList.at( i );
4800 
4801  delete [] innerRings;
4802 
4803  GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon( newOuterRing, newInnerRings, ringList.size() );
4804  delete[] newInnerRings;
4805 
4806  return reshapedPolygon;
4807 }
4808 
4809 GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos )
4810 {
4811  if ( !line || !reshapeLineGeos )
4812  return 0;
4813 
4814  bool atLeastTwoIntersections = false;
4815 
4816  try
4817  {
4818  //make sure there are at least two intersection between line and reshape geometry
4819  GEOSGeometry* intersectGeom = GEOSIntersection( line, reshapeLineGeos );
4820  if ( intersectGeom )
4821  {
4822  atLeastTwoIntersections = ( GEOSGeomTypeId( intersectGeom ) == GEOS_MULTIPOINT && GEOSGetNumGeometries( intersectGeom ) > 1 );
4823  GEOSGeom_destroy( intersectGeom );
4824  }
4825  }
4826  catch ( GEOSException &e )
4827  {
4828  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
4829  atLeastTwoIntersections = false;
4830  }
4831 
4832  if ( !atLeastTwoIntersections )
4833  return 0;
4834 
4835  //begin and end point of original line
4836  const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq( line );
4837  if ( !lineCoordSeq )
4838  return 0;
4839 
4840  unsigned int lineCoordSeqSize;
4841  if ( GEOSCoordSeq_getSize( lineCoordSeq, &lineCoordSeqSize ) == 0 )
4842  return 0;
4843 
4844  if ( lineCoordSeqSize < 2 )
4845  return 0;
4846 
4847  //first and last vertex of line
4848  double x1, y1, x2, y2;
4849  GEOSCoordSeq_getX( lineCoordSeq, 0, &x1 );
4850  GEOSCoordSeq_getY( lineCoordSeq, 0, &y1 );
4851  GEOSCoordSeq_getX( lineCoordSeq, lineCoordSeqSize - 1, &x2 );
4852  GEOSCoordSeq_getY( lineCoordSeq, lineCoordSeqSize - 1, &y2 );
4853  GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) );
4854  GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) );
4855 
4856  bool isRing = false;
4857  if ( GEOSGeomTypeId( line ) == GEOS_LINEARRING || GEOSEquals( beginLineVertex, endLineVertex ) == 1 )
4858  isRing = true;
4859 
4860 //node line and reshape line
4861  GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line );
4862  if ( !nodedGeometry )
4863  {
4864  GEOSGeom_destroy( beginLineVertex );
4865  GEOSGeom_destroy( endLineVertex );
4866  return 0;
4867  }
4868 
4869  //and merge them together
4870  GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry );
4871  GEOSGeom_destroy( nodedGeometry );
4872  if ( !mergedLines )
4873  {
4874  GEOSGeom_destroy( beginLineVertex );
4875  GEOSGeom_destroy( endLineVertex );
4876  return 0;
4877  }
4878 
4879  int numMergedLines = GEOSGetNumGeometries( mergedLines );
4880  if ( numMergedLines < 2 ) //some special cases. Normally it is >2
4881  {
4882  GEOSGeom_destroy( beginLineVertex );
4883  GEOSGeom_destroy( endLineVertex );
4884  if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline
4885  return GEOSGeom_clone( reshapeLineGeos );
4886  else
4887  return 0;
4888  }
4889 
4890  QList<GEOSGeometry*> resultLineParts; //collection with the line segments that will be contained in result
4891  QList<GEOSGeometry*> probableParts; //parts where we can decide on inclusion only after going through all the candidates
4892 
4893  for ( int i = 0; i < numMergedLines; ++i )
4894  {
4895  const GEOSGeometry* currentGeom;
4896 
4897  currentGeom = GEOSGetGeometryN( mergedLines, i );
4898  const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq( currentGeom );
4899  unsigned int currentCoordSeqSize;
4900  GEOSCoordSeq_getSize( currentCoordSeq, &currentCoordSeqSize );
4901  if ( currentCoordSeqSize < 2 )
4902  continue;
4903 
4904  //get the two endpoints of the current line merge result
4905  double xBegin, xEnd, yBegin, yEnd;
4906  GEOSCoordSeq_getX( currentCoordSeq, 0, &xBegin );
4907  GEOSCoordSeq_getY( currentCoordSeq, 0, &yBegin );
4908  GEOSCoordSeq_getX( currentCoordSeq, currentCoordSeqSize - 1, &xEnd );
4909  GEOSCoordSeq_getY( currentCoordSeq, currentCoordSeqSize - 1, &yEnd );
4910  GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) );
4911  GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) );
4912 
4913  //check how many endpoints of the line merge result are on the (original) line
4914  int nEndpointsOnOriginalLine = 0;
4915  if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 )
4916  nEndpointsOnOriginalLine += 1;
4917 
4918  if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 )
4919  nEndpointsOnOriginalLine += 1;
4920 
4921  //check how many endpoints equal the endpoints of the original line
4922  int nEndpointsSameAsOriginalLine = 0;
4923  if ( GEOSEquals( beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( beginCurrentGeomVertex, endLineVertex ) == 1 )
4924  nEndpointsSameAsOriginalLine += 1;
4925 
4926  if ( GEOSEquals( endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( endCurrentGeomVertex, endLineVertex ) == 1 )
4927  nEndpointsSameAsOriginalLine += 1;
4928 
4929  //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings)
4930  bool currentGeomOverlapsOriginalGeom = false;
4931  bool currentGeomOverlapsReshapeLine = false;
4932  if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 )
4933  currentGeomOverlapsOriginalGeom = true;
4934 
4935  if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 )
4936  currentGeomOverlapsReshapeLine = true;
4937 
4938  //logic to decide if this part belongs to the result
4939  if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
4940  {
4941  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
4942  }
4943  //for closed rings, we take one segment from the candidate list
4944  else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom )
4945  {
4946  probableParts.push_back( GEOSGeom_clone( currentGeom ) );
4947  }
4948  else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
4949  {
4950  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
4951  }
4952  else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom )
4953  {
4954  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
4955  }
4956  else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine )
4957  {
4958  resultLineParts.push_back( GEOSGeom_clone( currentGeom ) );
4959  }
4960 
4961  GEOSGeom_destroy( beginCurrentGeomVertex );
4962  GEOSGeom_destroy( endCurrentGeomVertex );
4963  }
4964 
4965  //add the longest segment from the probable list for rings (only used for polygon rings)
4966  if ( isRing && probableParts.size() > 0 )
4967  {
4968  GEOSGeometry* maxGeom = 0; //the longest geometry in the probabla list
4969  GEOSGeometry* currentGeom = 0;
4970  double maxLength = -DBL_MAX;
4971  double currentLength = 0;
4972  for ( int i = 0; i < probableParts.size(); ++i )
4973  {
4974  currentGeom = probableParts.at( i );
4975  GEOSLength( currentGeom, &currentLength );
4976  if ( currentLength > maxLength )
4977  {
4978  maxLength = currentLength;
4979  GEOSGeom_destroy( maxGeom );
4980  maxGeom = currentGeom;
4981  }
4982  else
4983  {
4984  GEOSGeom_destroy( currentGeom );
4985  }
4986  }
4987  resultLineParts.push_back( maxGeom );
4988  }
4989 
4990  GEOSGeom_destroy( beginLineVertex );
4991  GEOSGeom_destroy( endLineVertex );
4992  GEOSGeom_destroy( mergedLines );
4993 
4994  GEOSGeometry* result = 0;
4995  if ( resultLineParts.size() < 1 )
4996  return 0;
4997 
4998  if ( resultLineParts.size() == 1 ) //the whole result was reshaped
4999  {
5000  result = resultLineParts[0];
5001  }
5002  else //>1
5003  {
5004  GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()];
5005  for ( int i = 0; i < resultLineParts.size(); ++i )
5006  {
5007  lineArray[i] = resultLineParts[i];
5008  }
5009 
5010  //create multiline from resultLineParts
5011  GEOSGeometry* multiLineGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, lineArray, resultLineParts.size() );
5012  delete [] lineArray;
5013 
5014  //then do a linemerge with the newly combined partstrings
5015  result = GEOSLineMerge( multiLineGeom );
5016  GEOSGeom_destroy( multiLineGeom );
5017  }
5018 
5019  //now test if the result is a linestring. Otherwise something went wrong
5020  if ( GEOSGeomTypeId( result ) != GEOS_LINESTRING )
5021  {
5022  GEOSGeom_destroy( result );
5023  return 0;
5024  }
5025 
5026  return result;
5027 }
5028 
5029 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const
5030 {
5031  //Find out the intersection points between splitLineGeos and this geometry.
5032  //These points need to be tested for topological correctness by the calling function
5033  //if topological editing is enabled
5034 
5035  testPoints.clear();
5036  GEOSGeometry* intersectionGeom = GEOSIntersection( mGeos, splitLine );
5037  if ( !intersectionGeom )
5038  return 1;
5039 
5040  bool simple = false;
5041  int nIntersectGeoms = 1;
5042  if ( GEOSGeomTypeId( intersectionGeom ) == GEOS_LINESTRING || GEOSGeomTypeId( intersectionGeom ) == GEOS_POINT )
5043  simple = true;
5044 
5045  if ( !simple )
5046  nIntersectGeoms = GEOSGetNumGeometries( intersectionGeom );
5047 
5048  for ( int i = 0; i < nIntersectGeoms; ++i )
5049  {
5050  const GEOSGeometry* currentIntersectGeom;
5051  if ( simple )
5052  currentIntersectGeom = intersectionGeom;
5053  else
5054  currentIntersectGeom = GEOSGetGeometryN( intersectionGeom, i );
5055 
5056  const GEOSCoordSequence* lineSequence = GEOSGeom_getCoordSeq( currentIntersectGeom );
5057  unsigned int sequenceSize = 0;
5058  double x, y;
5059  if ( GEOSCoordSeq_getSize( lineSequence, &sequenceSize ) != 0 )
5060  {
5061  for ( unsigned int i = 0; i < sequenceSize; ++i )
5062  {
5063  if ( GEOSCoordSeq_getX( lineSequence, i, &x ) != 0 )
5064  {
5065  if ( GEOSCoordSeq_getY( lineSequence, i, &y ) != 0 )
5066  {
5067  testPoints.push_back( QgsPoint( x, y ) );
5068  }
5069  }
5070  }
5071  }
5072  }
5073  GEOSGeom_destroy( intersectionGeom );
5074  return 0;
5075 }
5076 
5077 GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom )
5078 {
5079  if ( !splitLine || !geom )
5080  return 0;
5081 
5082  GEOSGeometry *geometryBoundary = 0;
5083  if ( GEOSGeomTypeId( geom ) == GEOS_POLYGON || GEOSGeomTypeId( geom ) == GEOS_MULTIPOLYGON )
5084  geometryBoundary = GEOSBoundary( geom );
5085  else
5086  geometryBoundary = GEOSGeom_clone( geom );
5087 
5088  GEOSGeometry *splitLineClone = GEOSGeom_clone( splitLine );
5089  GEOSGeometry *unionGeometry = GEOSUnion( splitLineClone, geometryBoundary );
5090  GEOSGeom_destroy( splitLineClone );
5091 
5092  GEOSGeom_destroy( geometryBoundary );
5093  return unionGeometry;
5094 }
5095 
5096 int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 )
5097 {
5098  if ( !line1 || !line2 )
5099  {
5100  return -1;
5101  }
5102 
5103  double bufferDistance = 0.00001;
5104  if ( geomInDegrees( line2 ) ) //use more accurate tolerance for degrees
5105  bufferDistance = 0.00000001;
5106 
5107  GEOSGeometry* bufferGeom = GEOSBuffer( line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS );
5108  if ( !bufferGeom )
5109  return -2;
5110 
5111  GEOSGeometry* intersectionGeom = GEOSIntersection( bufferGeom, line1 );
5112 
5113  //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2)
5114  double intersectGeomLength;
5115  double line1Length;
5116 
5117  GEOSLength( intersectionGeom, &intersectGeomLength );
5118  GEOSLength( line1, &line1Length );
5119 
5120  GEOSGeom_destroy( bufferGeom );
5121  GEOSGeom_destroy( intersectionGeom );
5122 
5123  double intersectRatio = line1Length / intersectGeomLength;
5124  if ( intersectRatio > 0.9 && intersectRatio < 1.1 )
5125  return 1;
5126 
5127  return 0;
5128 }
5129 
5130 int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line )
5131 {
5132  if ( !point || !line )
5133  return -1;
5134 
5135  double bufferDistance = 0.000001;
5136  if ( geomInDegrees( line ) )
5137  bufferDistance = 0.00000001;
5138 
5139  GEOSGeometry* lineBuffer = GEOSBuffer( line, bufferDistance, 8 );
5140  if ( !lineBuffer )
5141  return -2;
5142 
5143  bool contained = false;
5144  if ( GEOSContains( lineBuffer, point ) == 1 )
5145  contained = true;
5146 
5147  GEOSGeom_destroy( lineBuffer );
5148  return contained;
5149 }
5150 
5151 bool QgsGeometry::geomInDegrees( const GEOSGeometry* geom )
5152 {
5153  GEOSGeometry* bbox = GEOSEnvelope( geom );
5154  if ( !bbox )
5155  return false;
5156 
5157  const GEOSGeometry* bBoxRing = GEOSGetExteriorRing( bbox );
5158  if ( !bBoxRing )
5159  return false;
5160 
5161  const GEOSCoordSequence* bBoxCoordSeq = GEOSGeom_getCoordSeq( bBoxRing );
5162 
5163  if ( !bBoxCoordSeq )
5164  return false;
5165 
5166  unsigned int nCoords = 0;
5167  if ( !GEOSCoordSeq_getSize( bBoxCoordSeq, &nCoords ) )
5168  return false;
5169 
5170  double x, y;
5171  for ( unsigned int i = 0; i < ( nCoords - 1 ); ++i )
5172  {
5173  GEOSCoordSeq_getX( bBoxCoordSeq, i, &x );
5174  if ( x > 180 || x < -180 )
5175  return false;
5176 
5177  GEOSCoordSeq_getY( bBoxCoordSeq, i, &y );
5178  if ( y > 90 || y < -90 )
5179  return false;
5180 
5181  }
5182 
5183  return true;
5184 }
5185 
5186 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const
5187 {
5188  if ( !g )
5189  return 0;
5190 
5191  int geometryType = GEOSGeomTypeId( g );
5192  if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING
5193  || geometryType == GEOS_POLYGON )
5194  return 1;
5195 
5196  //calling GEOSGetNumGeometries is save for multi types and collections also in geos2
5197  return GEOSGetNumGeometries( g );
5198 }
5199 
5200 int QgsGeometry::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry*>& splitResult )
5201 {
5202  if ( mDirtyGeos )
5203  exportWkbToGeos();
5204 
5205  if ( !mGeos )
5206  return 1;
5207 
5208  //convert mGeos to geometry collection
5209  int type = GEOSGeomTypeId( mGeos );
5210  if ( type != GEOS_GEOMETRYCOLLECTION &&
5211  type != GEOS_MULTILINESTRING &&
5212  type != GEOS_MULTIPOLYGON &&
5213  type != GEOS_MULTIPOINT )
5214  return 0;
5215 
5216  QVector<GEOSGeometry*> copyList = splitResult;
5217  splitResult.clear();
5218 
5219  //collect all the geometries that belong to the initial multifeature
5220  QVector<GEOSGeometry*> unionGeom;
5221 
5222  for ( int i = 0; i < copyList.size(); ++i )
5223  {
5224  //is this geometry a part of the original multitype?
5225  bool isPart = false;
5226  for ( int j = 0; j < GEOSGetNumGeometries( mGeos ); j++ )
5227  {
5228  if ( GEOSEquals( copyList[i], GEOSGetGeometryN( mGeos, j ) ) )
5229  {
5230  isPart = true;
5231  break;
5232  }
5233  }
5234 
5235  if ( isPart )
5236  {
5237  unionGeom << copyList[i];
5238  }
5239  else
5240  {
5241  QVector<GEOSGeometry*> geomVector;
5242  geomVector << copyList[i];
5243 
5244  if ( type == GEOS_MULTILINESTRING )
5245  splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector );
5246  else if ( type == GEOS_MULTIPOLYGON )
5247  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector );
5248  else
5249  GEOSGeom_destroy( copyList[i] );
5250  }
5251  }
5252 
5253  //make multifeature out of unionGeom
5254  if ( unionGeom.size() > 0 )
5255  {
5256  if ( type == GEOS_MULTILINESTRING )
5257  splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom );
5258  else if ( type == GEOS_MULTIPOLYGON )
5259  splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom );
5260  }
5261  else
5262  {
5263  unionGeom.clear();
5264  }
5265 
5266  return 0;
5267 }
5268 
5269 QgsPoint QgsGeometry::asPoint( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5270 {
5271  wkbPtr += 1 + sizeof( int );
5272 
5273  double x, y;
5274  wkbPtr >> x >> y;
5275  if ( hasZValue )
5276  wkbPtr += sizeof( double );
5277 
5278  return QgsPoint( x, y );
5279 }
5280 
5281 QgsPolyline QgsGeometry::asPolyline( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5282 {
5283  wkbPtr += 1 + sizeof( int );
5284 
5285  unsigned int nPoints;
5286  wkbPtr >> nPoints;
5287 
5288  QgsPolyline line( nPoints );
5289 
5290  // Extract the points from the WKB format into the x and y vectors.
5291  for ( uint i = 0; i < nPoints; ++i )
5292  {
5293  double x, y;
5294  wkbPtr >> x >> y;
5295  if ( hasZValue )
5296  wkbPtr += sizeof( double );
5297 
5298  line[i] = QgsPoint( x, y );
5299  }
5300 
5301  return line;
5302 }
5303 
5304 QgsPolygon QgsGeometry::asPolygon( QgsConstWkbPtr &wkbPtr, bool hasZValue ) const
5305 {
5306  wkbPtr += 1 + sizeof( int );
5307 
5308  // get number of rings in the polygon
5309  unsigned int numRings;
5310  wkbPtr >> numRings;
5311 
5312  if ( numRings == 0 ) // sanity check for zero rings in polygon
5313  return QgsPolygon();
5314 
5315  QgsPolygon rings( numRings );
5316 
5317  for ( uint idx = 0; idx < numRings; idx++ )
5318  {
5319  int nPoints;
5320  wkbPtr >> nPoints;
5321 
5322  QgsPolyline ring( nPoints );
5323 
5324  for ( int jdx = 0; jdx < nPoints; jdx++ )
5325  {
5326  double x, y;
5327  wkbPtr >> x >> y;
5328  if ( hasZValue )
5329  wkbPtr += sizeof( double );
5330 
5331  ring[jdx] = QgsPoint( x, y );
5332  }
5333 
5334  rings[idx] = ring;
5335  }
5336 
5337  return rings;
5338 }
5339 
5341 {
5342  QGis::WkbType type = wkbType();
5343  if ( type != QGis::WKBPoint && type != QGis::WKBPoint25D )
5344  return QgsPoint( 0, 0 );
5345 
5346  QgsConstWkbPtr wkbPtr( mGeometry );
5347  return asPoint( wkbPtr, type == QGis::WKBPoint25D );
5348 }
5349 
5351 {
5352  QGis::WkbType type = wkbType();
5353  if ( type != QGis::WKBLineString && type != QGis::WKBLineString25D )
5354  return QgsPolyline();
5355 
5356  QgsConstWkbPtr wkbPtr( mGeometry );
5357  return asPolyline( wkbPtr, type == QGis::WKBLineString25D );
5358 }
5359 
5361 {
5362  QGis::WkbType type = wkbType();
5363  if ( type != QGis::WKBPolygon && type != QGis::WKBPolygon25D )
5364  return QgsPolygon();
5365 
5366  QgsConstWkbPtr wkbPtr( mGeometry );
5367  return asPolygon( wkbPtr, type == QGis::WKBPolygon25D );
5368 }
5369 
5371 {
5372  QGis::WkbType type = wkbType();
5373  if ( type != QGis::WKBMultiPoint && type != QGis::WKBMultiPoint25D )
5374  return QgsMultiPoint();
5375 
5376  bool hasZValue = ( type == QGis::WKBMultiPoint25D );
5377 
5378  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5379 
5380  int nPoints;
5381  wkbPtr >> nPoints;
5382 
5383  QgsMultiPoint points( nPoints );
5384  for ( int i = 0; i < nPoints; i++ )
5385  {
5386  points[i] = asPoint( wkbPtr, hasZValue );
5387  }
5388 
5389  return points;
5390 }
5391 
5393 {
5394  QGis::WkbType type = wkbType();
5395  if ( type != QGis::WKBMultiLineString && type != QGis::WKBMultiLineString25D )
5396  return QgsMultiPolyline();
5397 
5398  bool hasZValue = ( type == QGis::WKBMultiLineString25D );
5399 
5400  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5401 
5402  int numLineStrings;
5403  wkbPtr >> numLineStrings;
5404 
5405  QgsMultiPolyline lines( numLineStrings );
5406 
5407  for ( int i = 0; i < numLineStrings; i++ )
5408  lines[i] = asPolyline( wkbPtr, hasZValue );
5409 
5410  return lines;
5411 }
5412 
5414 {
5415  QGis::WkbType type = wkbType();
5416  if ( type != QGis::WKBMultiPolygon && type != QGis::WKBMultiPolygon25D )
5417  return QgsMultiPolygon();
5418 
5419  bool hasZValue = ( type == QGis::WKBMultiPolygon25D );
5420 
5421  QgsConstWkbPtr wkbPtr( mGeometry + 1 + sizeof( int ) );
5422 
5423  int numPolygons;
5424  wkbPtr >> numPolygons;
5425 
5426  QgsMultiPolygon polygons( numPolygons );
5427 
5428  for ( int i = 0; i < numPolygons; i++ )
5429  polygons[i] = asPolygon( wkbPtr, hasZValue );
5430 
5431  return polygons;
5432 }
5433 
5435 {
5436  if ( mDirtyGeos )
5437  exportWkbToGeos();
5438 
5439  if ( !mGeos )
5440  return -1.0;
5441 
5442  double area;
5443 
5444  try
5445  {
5446  if ( GEOSArea( mGeos, &area ) == 0 )
5447  return -1.0;
5448  }
5449  CATCH_GEOS( -1.0 )
5450 
5451  return area;
5452 }
5453 
5455 {
5456  if ( mDirtyGeos )
5457  exportWkbToGeos();
5458 
5459  if ( !mGeos )
5460  return -1.0;
5461 
5462  double length;
5463 
5464  try
5465  {
5466  if ( GEOSLength( mGeos, &length ) == 0 )
5467  return -1.0;
5468  }
5469  CATCH_GEOS( -1.0 )
5470 
5471  return length;
5472 }
5474 {
5475  if ( mDirtyGeos )
5476  exportWkbToGeos();
5477 
5478  if ( geom.mDirtyGeos )
5479  geom.exportWkbToGeos();
5480 
5481  if ( !mGeos || !geom.mGeos )
5482  return -1.0;
5483 
5484  double dist = -1.0;
5485 
5486  try
5487  {
5488  GEOSDistance( mGeos, geom.mGeos, &dist );
5489  }
5490  CATCH_GEOS( -1.0 )
5491 
5492  return dist;
5493 }
5494 
5495 QgsGeometry* QgsGeometry::buffer( double distance, int segments )
5496 {
5497  if ( mDirtyGeos )
5498  exportWkbToGeos();
5499 
5500  if ( !mGeos )
5501  return 0;
5502 
5503  try
5504  {
5505  return fromGeosGeom( GEOSBuffer( mGeos, distance, segments ) );
5506  }
5507  CATCH_GEOS( 0 )
5508 }
5509 
5511 {
5512  if ( mDirtyGeos )
5513  exportWkbToGeos();
5514 
5515  if ( !mGeos )
5516  return 0;
5517 
5518  try
5519  {
5520  return fromGeosGeom( GEOSTopologyPreserveSimplify( mGeos, tolerance ) );
5521  }
5522  CATCH_GEOS( 0 )
5523 }
5524 
5526 {
5527  if ( mDirtyGeos )
5528  exportWkbToGeos();
5529 
5530  if ( !mGeos )
5531  return 0;
5532 
5533  try
5534  {
5535  return fromGeosGeom( GEOSGetCentroid( mGeos ) );
5536  }
5537  CATCH_GEOS( 0 )
5538 }
5539 
5541 {
5542  if ( mDirtyGeos )
5543  exportWkbToGeos();
5544 
5545  if ( !mGeos )
5546  return 0;
5547 
5548  try
5549  {
5550  return fromGeosGeom( GEOSConvexHull( mGeos ) );
5551  }
5552  CATCH_GEOS( 0 )
5553 }
5554 
5556 {
5557 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \
5558  ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=2)))
5559  if ( mDirtyGeos )
5560  exportWkbToGeos();
5561 
5562  if ( !mGeos )
5563  return 0;
5564 
5565  try
5566  {
5567  return fromGeosGeom( GEOSInterpolate( mGeos, distance ) );
5568  }
5569  CATCH_GEOS( 0 )
5570 #else
5571  QgsMessageLog::logMessage( QObject::tr( "GEOS prior to 3.2 doesn't support GEOSInterpolate" ), QObject::tr( "GEOS" ) );
5572 #endif
5573 }
5574 
5576 {
5577  if ( !geometry )
5578  return NULL;
5579 
5580  if ( mDirtyGeos )
5581  exportWkbToGeos();
5582 
5583  if ( geometry->mDirtyGeos )
5584  geometry->exportWkbToGeos();
5585 
5586  if ( !mGeos || !geometry->mGeos )
5587  return 0;
5588 
5589  try
5590  {
5591  return fromGeosGeom( GEOSIntersection( mGeos, geometry->mGeos ) );
5592  }
5593  CATCH_GEOS( 0 )
5594 }
5595 
5597 {
5598  if ( !geometry )
5599  return NULL;
5600 
5601  if ( mDirtyGeos )
5602  exportWkbToGeos();
5603 
5604  if ( geometry->mDirtyGeos )
5605  geometry->exportWkbToGeos();
5606 
5607  if ( !mGeos || !geometry->mGeos )
5608  return 0;
5609 
5610  try
5611  {
5612  GEOSGeometry* unionGeom = GEOSUnion( mGeos, geometry->mGeos );
5613  if ( !unionGeom )
5614  return 0;
5615 
5616  if ( type() == QGis::Line )
5617  {
5618  GEOSGeometry* mergedGeom = GEOSLineMerge( unionGeom );
5619  if ( mergedGeom )
5620  {
5621  GEOSGeom_destroy( unionGeom );
5622  unionGeom = mergedGeom;
5623  }
5624  }
5625  return fromGeosGeom( unionGeom );
5626  }
5627  CATCH_GEOS( new QgsGeometry( *this ) ) //return this geometry if union not possible
5628 }
5629 
5631 {
5632  if ( !geometry )
5633  return NULL;
5634 
5635  if ( mDirtyGeos )
5636  exportWkbToGeos();
5637 
5638  if ( geometry->mDirtyGeos )
5639  geometry->exportWkbToGeos();
5640 
5641  if ( !mGeos || !geometry->mGeos )
5642  return 0;
5643 
5644  try
5645  {
5646  return fromGeosGeom( GEOSDifference( mGeos, geometry->mGeos ) );
5647  }
5648  CATCH_GEOS( 0 )
5649 }
5650 
5652 {
5653  if ( !geometry )
5654  return NULL;
5655 
5656  if ( mDirtyGeos )
5657  exportWkbToGeos();
5658 
5659  if ( geometry->mDirtyGeos )
5660  geometry->exportWkbToGeos();
5661 
5662  if ( !mGeos || !geometry->mGeos )
5663  return 0;
5664 
5665  try
5666  {
5667  return fromGeosGeom( GEOSSymDifference( mGeos, geometry->mGeos ) );
5668  }
5669  CATCH_GEOS( 0 )
5670 }
5671 
5672 QList<QgsGeometry*> QgsGeometry::asGeometryCollection() const
5673 {
5674  if ( mDirtyGeos )
5675  exportWkbToGeos();
5676 
5677  if ( !mGeos )
5678  return QList<QgsGeometry*>();
5679 
5680  int type = GEOSGeomTypeId( mGeos );
5681  QgsDebugMsg( "geom type: " + QString::number( type ) );
5682 
5683  QList<QgsGeometry*> geomCollection;
5684 
5685  if ( type != GEOS_MULTIPOINT &&
5686  type != GEOS_MULTILINESTRING &&
5687  type != GEOS_MULTIPOLYGON &&
5688  type != GEOS_GEOMETRYCOLLECTION )
5689  {
5690  // we have a single-part geometry - put there a copy of this one
5691  geomCollection.append( new QgsGeometry( *this ) );
5692  return geomCollection;
5693  }
5694 
5695  int count = GEOSGetNumGeometries( mGeos );
5696  QgsDebugMsg( "geom count: " + QString::number( count ) );
5697 
5698  for ( int i = 0; i < count; ++i )
5699  {
5700  const GEOSGeometry * geometry = GEOSGetGeometryN( mGeos, i );
5701  geomCollection.append( fromGeosGeom( GEOSGeom_clone( geometry ) ) );
5702  }
5703 
5704  return geomCollection;
5705 }
5706 
5707 bool QgsGeometry::deleteRing( int ringNum, int partNum )
5708 {
5709  if ( ringNum <= 0 || partNum < 0 )
5710  return false;
5711 
5712  switch ( wkbType() )
5713  {
5714  case QGis::WKBPolygon25D:
5715  case QGis::WKBPolygon:
5716  {
5717  if ( partNum != 0 )
5718  return false;
5719 
5720  QgsPolygon polygon = asPolygon();
5721  if ( ringNum >= polygon.count() )
5722  return false;
5723 
5724  polygon.remove( ringNum );
5725 
5726  QgsGeometry* g2 = QgsGeometry::fromPolygon( polygon );
5727  *this = *g2;
5728  delete g2;
5729  return true;
5730  }
5731 
5733  case QGis::WKBMultiPolygon:
5734  {
5735  QgsMultiPolygon mpolygon = asMultiPolygon();
5736 
5737  if ( partNum >= mpolygon.count() )
5738  return false;
5739 
5740  if ( ringNum >= mpolygon[partNum].count() )
5741  return false;
5742 
5743  mpolygon[partNum].remove( ringNum );
5744 
5745  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
5746  *this = *g2;
5747  delete g2;
5748  return true;
5749  }
5750 
5751  default:
5752  return false; // only makes sense with polygons and multipolygons
5753  }
5754 }
5755 
5756 bool QgsGeometry::deletePart( int partNum )
5757 {
5758  if ( partNum < 0 )
5759  return false;
5760 
5761  switch ( wkbType() )
5762  {
5764  case QGis::WKBMultiPoint:
5765  {
5766  QgsMultiPoint mpoint = asMultiPoint();
5767 
5768  if ( partNum >= mpoint.size() || mpoint.size() == 1 )
5769  return false;
5770 
5771  mpoint.remove( partNum );
5772 
5773  QgsGeometry* g2 = QgsGeometry::fromMultiPoint( mpoint );
5774  *this = *g2;
5775  delete g2;
5776  break;
5777  }
5778 
5781  {
5783 
5784  if ( partNum >= mline.size() || mline.size() == 1 )
5785  return false;
5786 
5787  mline.remove( partNum );
5788 
5790  *this = *g2;
5791  delete g2;
5792  break;
5793  }
5794 
5796  case QGis::WKBMultiPolygon:
5797  {
5798  QgsMultiPolygon mpolygon = asMultiPolygon();
5799 
5800  if ( partNum >= mpolygon.size() || mpolygon.size() == 1 )
5801  return false;
5802 
5803  mpolygon.remove( partNum );
5804 
5805  QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon );
5806  *this = *g2;
5807  delete g2;
5808  break;
5809  }
5810 
5811  default:
5812  // single part geometries are ignored
5813  return false;
5814  }
5815 
5816  return true;
5817 }
5818 
5821 static GEOSGeometry* _makeUnion( QList<GEOSGeometry*> geoms )
5822 {
5823 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && (((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=3)) || (GEOS_VERSION_MAJOR>3))
5824  GEOSGeometry* geomCollection = 0;
5825  geomCollection = createGeosCollection( GEOS_GEOMETRYCOLLECTION, geoms.toVector() );
5826  GEOSGeometry* geomUnion = GEOSUnaryUnion( geomCollection );
5827  GEOSGeom_destroy( geomCollection );
5828  return geomUnion;
5829 #else
5830  GEOSGeometry* geomCollection = geoms.takeFirst();
5831 
5832  while ( !geoms.isEmpty() )
5833  {
5834  GEOSGeometry* g = geoms.takeFirst();
5835  GEOSGeometry* geomCollectionNew = GEOSUnion( geomCollection, g );
5836  GEOSGeom_destroy( geomCollection );
5837  GEOSGeom_destroy( g );
5838  geomCollection = geomCollectionNew;
5839  }
5840 
5841  return geomCollection;
5842 #endif
5843 }
5844 
5845 int QgsGeometry::avoidIntersections( QMap<QgsVectorLayer*, QSet< QgsFeatureId > > ignoreFeatures )
5846 {
5847  int returnValue = 0;
5848 
5849  //check if g has polygon type
5850  if ( type() != QGis::Polygon )
5851  return 1;
5852 
5853  QGis::WkbType geomTypeBeforeModification = wkbType();
5854 
5855  //read avoid intersections list from project properties
5856  bool listReadOk;
5857  QStringList avoidIntersectionsList = QgsProject::instance()->readListEntry( "Digitizing", "/AvoidIntersectionsList", QStringList(), &listReadOk );
5858  if ( !listReadOk )
5859  return true; //no intersections stored in project does not mean error
5860 
5861  QList<GEOSGeometry*> nearGeometries;
5862 
5863  //go through list, convert each layer to vector layer and call QgsVectorLayer::removePolygonIntersections for each
5864  QgsVectorLayer* currentLayer = 0;
5865  QStringList::const_iterator aIt = avoidIntersectionsList.constBegin();
5866  for ( ; aIt != avoidIntersectionsList.constEnd(); ++aIt )
5867  {
5868  currentLayer = dynamic_cast<QgsVectorLayer*>( QgsMapLayerRegistry::instance()->mapLayer( *aIt ) );
5869  if ( currentLayer )
5870  {
5871  QgsFeatureIds ignoreIds;
5872  QMap<QgsVectorLayer*, QSet<qint64> >::const_iterator ignoreIt = ignoreFeatures.find( currentLayer );
5873  if ( ignoreIt != ignoreFeatures.constEnd() )
5874  ignoreIds = ignoreIt.value();
5875 
5878  .setSubsetOfAttributes( QgsAttributeList() ) );
5879  QgsFeature f;
5880  while ( fi.nextFeature( f ) )
5881  {
5882  if ( ignoreIds.contains( f.id() ) )
5883  continue;
5884 
5885  if ( !f.geometry() )
5886  continue;
5887 
5888  nearGeometries << GEOSGeom_clone( f.geometry()->asGeos() );
5889  }
5890  }
5891  }
5892 
5893  if ( nearGeometries.isEmpty() )
5894  return 0;
5895 
5896  GEOSGeometry* nearGeometriesUnion = 0;
5897  GEOSGeometry* geomWithoutIntersections = 0;
5898  try
5899  {
5900  nearGeometriesUnion = _makeUnion( nearGeometries );
5901  geomWithoutIntersections = GEOSDifference( asGeos(), nearGeometriesUnion );
5902 
5903  fromGeos( geomWithoutIntersections );
5904 
5905  GEOSGeom_destroy( nearGeometriesUnion );
5906  }
5907  catch ( GEOSException &e )
5908  {
5909  if ( nearGeometriesUnion )
5910  GEOSGeom_destroy( nearGeometriesUnion );
5911  if ( geomWithoutIntersections )
5912  GEOSGeom_destroy( geomWithoutIntersections );
5913 
5914  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
5915  return 3;
5916  }
5917 
5918  //make sure the geometry still has the same type (e.g. no change from polygon to multipolygon)
5919  if ( wkbType() != geomTypeBeforeModification )
5920  return 2;
5921 
5922  return returnValue;
5923 }
5924 
5925 void QgsGeometry::validateGeometry( QList<Error> &errors )
5926 {
5928 }
5929 
5931 {
5932  try
5933  {
5934  const GEOSGeometry *g = asGeos();
5935 
5936  if ( !g )
5937  return false;
5938 
5939  return GEOSisValid( g );
5940  }
5941  catch ( GEOSException &e )
5942  {
5943  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
5944  return false;
5945  }
5946 }
5947 
5949 {
5950  return geosRelOp( GEOSEquals, this, &g );
5951 }
5952 
5954 {
5955  try
5956  {
5957  const GEOSGeometry *g = asGeos();
5958 
5959  if ( !g )
5960  return false;
5961 
5962  return GEOSisEmpty( g );
5963  }
5964  catch ( GEOSException &e )
5965  {
5966  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
5967  return false;
5968  }
5969 }
5970 
5971 double QgsGeometry::leftOf( double x, double y, double& x1, double& y1, double& x2, double& y2 )
5972 {
5973  double f1 = x - x1;
5974  double f2 = y2 - y1;
5975  double f3 = y - y1;
5976  double f4 = x2 - x1;
5977  return f1*f2 - f3*f4;
5978 }
QgsFeatureId id() const
Get the feature id for this feature.
Definition: qgsfeature.cpp:101
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:497
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:322
bool isMultipart()
Returns true if wkb of the geometry is of WKBMulti* 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:297
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
Definition: qgsfeature.cpp:113
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:491
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
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)
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:500
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:488
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:494