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