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_WS_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  QgsDebugMsg( "Entering." );
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  QgsDebugMsg( QString( "Exiting with nearest point %1, dist %2." )
2314  .arg( point.toString() ).arg( sqrDist ) );
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::translate( double dx, double dy )
2710 {
2711  if ( mDirtyWkb )
2712  exportGeosToWkb();
2713 
2714  if ( !mGeometry )
2715  {
2716  QgsDebugMsg( "WKB geometry not available!" );
2717  return 1;
2718  }
2719 
2720  QgsWkbPtr wkbPtr( mGeometry + 1 );
2722  wkbPtr >> wkbType;
2723 
2724  bool hasZValue = false;
2725  switch ( wkbType )
2726  {
2727  case QGis::WKBPoint25D:
2728  case QGis::WKBPoint:
2729  {
2730  translateVertex( wkbPtr, dx, dy, 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  translateVertex( wkbPtr, dx, dy, 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  translateVertex( wkbPtr, dx, dy, hasZValue );
2758 
2759  }
2760  break;
2761  }
2762 
2764  hasZValue = true;
2765  case QGis::WKBMultiPoint:
2766  {
2767  int nPoints;
2768  wkbPtr >> nPoints;
2769 
2770  for ( int index = 0; index < nPoints; ++index )
2771  {
2772  wkbPtr += 1 + sizeof( int );
2773  translateVertex( wkbPtr, dx, dy, hasZValue );
2774  }
2775  break;
2776  }
2777 
2779  hasZValue = true;
2781  {
2782  int nLines;
2783  wkbPtr >> nLines;
2784  for ( int index = 0; index < nLines; ++index )
2785  {
2786  wkbPtr += 1 + sizeof( int );
2787  int nPoints;
2788  wkbPtr >> nPoints;
2789  for ( int index2 = 0; index2 < nPoints; ++index2 )
2790  translateVertex( wkbPtr, dx, dy, hasZValue );
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 
2805  int nRings;
2806  wkbPtr >> nRings;
2807 
2808  for ( int index2 = 0; index2 < nRings; ++index2 )
2809  {
2810  int nPoints;
2811  wkbPtr >> nPoints;
2812  for ( int index3 = 0; index3 < nPoints; ++index3 )
2813  translateVertex( wkbPtr, dx, dy, hasZValue );
2814  }
2815  }
2816  break;
2817  }
2818 
2819  default:
2820  break;
2821  }
2822  mDirtyGeos = true;
2823  return 0;
2824 }
2825 
2827 {
2828  if ( mDirtyWkb )
2829  exportGeosToWkb();
2830 
2831  if ( !mGeometry )
2832  {
2833  QgsDebugMsg( "WKB geometry not available!" );
2834  return 1;
2835  }
2836 
2837  bool hasZValue = false;
2838  QgsWkbPtr wkbPtr( mGeometry + 1 );
2840  wkbPtr >> wkbType;
2841 
2842  switch ( wkbType )
2843  {
2844  case QGis::WKBPoint25D:
2845  case QGis::WKBPoint:
2846  {
2847  transformVertex( wkbPtr, ct, hasZValue );
2848  }
2849  break;
2850 
2852  hasZValue = true;
2853  case QGis::WKBLineString:
2854  {
2855  int nPoints;
2856  wkbPtr >> nPoints;
2857  for ( int index = 0; index < nPoints; ++index )
2858  transformVertex( wkbPtr, ct, hasZValue );
2859 
2860  break;
2861  }
2862 
2863  case QGis::WKBPolygon25D:
2864  hasZValue = true;
2865  case QGis::WKBPolygon:
2866  {
2867  int nRings;
2868  wkbPtr >> nRings;
2869  for ( int index = 0; index < nRings; ++index )
2870  {
2871  int nPoints;
2872  wkbPtr >> nPoints;
2873  for ( int index2 = 0; index2 < nPoints; ++index2 )
2874  transformVertex( wkbPtr, ct, hasZValue );
2875 
2876  }
2877  break;
2878  }
2879 
2881  hasZValue = true;
2882  case QGis::WKBMultiPoint:
2883  {
2884  int nPoints;
2885  wkbPtr >> nPoints;
2886  for ( int index = 0; index < nPoints; ++index )
2887  {
2888  wkbPtr += 1 + sizeof( int );
2889  transformVertex( wkbPtr, ct, hasZValue );
2890  }
2891  break;
2892  }
2893 
2895  hasZValue = true;
2897  {
2898  int nLines;
2899  wkbPtr >> nLines;
2900  for ( int index = 0; index < nLines; ++index )
2901  {
2902  wkbPtr += 1 + sizeof( int );
2903  int nPoints;
2904  wkbPtr >> nPoints;
2905  for ( int index2 = 0; index2 < nPoints; ++index2 )
2906  transformVertex( wkbPtr, ct, hasZValue );
2907 
2908  }
2909  break;
2910  }
2911 
2913  hasZValue = true;
2914  case QGis::WKBMultiPolygon:
2915  {
2916  int nPolys;
2917  wkbPtr >> nPolys;
2918  for ( int index = 0; index < nPolys; ++index )
2919  {
2920  wkbPtr += 1 + sizeof( int ); //skip endian and polygon type
2921  int nRings;
2922  wkbPtr >> nRings;
2923  for ( int index2 = 0; index2 < nRings; ++index2 )
2924  {
2925  int nPoints;
2926  wkbPtr >> nPoints;
2927  for ( int index3 = 0; index3 < nPoints; ++index3 )
2928  transformVertex( wkbPtr, ct, hasZValue );
2929 
2930  }
2931  }
2932  }
2933 
2934  default:
2935  break;
2936  }
2937  mDirtyGeos = true;
2938  return 0;
2939 }
2940 
2941 int QgsGeometry::splitGeometry( const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint> &topologyTestPoints )
2942 {
2943  int returnCode = 0;
2944 
2945  //return if this type is point/multipoint
2946  if ( type() == QGis::Point )
2947  {
2948  return 1; //cannot split points
2949  }
2950 
2951  //make sure, mGeos and mWkb are there and up-to-date
2952  if ( mDirtyWkb )
2953  exportGeosToWkb();
2954 
2955  if ( mDirtyGeos )
2956  exportWkbToGeos();
2957 
2958  if ( !mGeos )
2959  return 1;
2960 
2961  if ( !GEOSisValid_r( geosinit.ctxt, mGeos ) )
2962  return 7;
2963 
2964  //make sure splitLine is valid
2965  if (( type() == QGis::Line && splitLine.size() < 1 ) ||
2966  ( type() == QGis::Polygon && splitLine.size() < 2 ) )
2967  return 1;
2968 
2969  newGeometries.clear();
2970 
2971  try
2972  {
2973  GEOSGeometry* splitLineGeos;
2974  if ( splitLine.size() > 1 )
2975  {
2976  splitLineGeos = createGeosLineString( splitLine.toVector() );
2977  }
2978  else if ( splitLine.size() == 1 )
2979  {
2980  splitLineGeos = createGeosPoint( splitLine.at( 0 ) );
2981  }
2982  else
2983  {
2984  return 1;
2985  }
2986  if ( !GEOSisValid_r( geosinit.ctxt, splitLineGeos ) || !GEOSisSimple_r( geosinit.ctxt, splitLineGeos ) )
2987  {
2988  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
2989  return 1;
2990  }
2991 
2992  if ( topological )
2993  {
2994  //find out candidate points for topological corrections
2995  if ( topologicalTestPointsSplit( splitLineGeos, topologyTestPoints ) != 0 )
2996  return 1;
2997  }
2998 
2999  //call split function depending on geometry type
3000  if ( type() == QGis::Line )
3001  {
3002  returnCode = splitLinearGeometry( splitLineGeos, newGeometries );
3003  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
3004  }
3005  else if ( type() == QGis::Polygon )
3006  {
3007  returnCode = splitPolygonGeometry( splitLineGeos, newGeometries );
3008  GEOSGeom_destroy_r( geosinit.ctxt, splitLineGeos );
3009  }
3010  else
3011  {
3012  return 1;
3013  }
3014  }
3015  CATCH_GEOS( 2 )
3016 
3017  return returnCode;
3018 }
3019 
3021 int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine )
3022 {
3023  if ( reshapeWithLine.size() < 2 )
3024  return 1;
3025 
3026  if ( type() == QGis::Point )
3027  return 1; //cannot reshape points
3028 
3029  GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() );
3030 
3031  //make sure this geos geometry is up-to-date
3032  if ( mDirtyGeos )
3033  exportWkbToGeos();
3034 
3035  if ( !mGeos )
3036  return 1;
3037 
3038  //single or multi?
3039  int numGeoms = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
3040  if ( numGeoms == -1 )
3041  return 1;
3042 
3043  bool isMultiGeom = false;
3044  int geosTypeId = GEOSGeomTypeId_r( geosinit.ctxt, mGeos );
3045  if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON )
3046  isMultiGeom = true;
3047 
3048  bool isLine = ( type() == QGis::Line );
3049 
3050  //polygon or multipolygon?
3051  if ( !isMultiGeom )
3052  {
3053  GEOSGeometry* reshapedGeometry;
3054  if ( isLine )
3055  reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos );
3056  else
3057  reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos );
3058 
3059  GEOSGeom_destroy_r( geosinit.ctxt, reshapeLineGeos );
3060  if ( reshapedGeometry )
3061  {
3062  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
3063  mGeos = reshapedGeometry;
3064  mDirtyWkb = true;
3065  return 0;
3066  }
3067  else
3068  {
3069  return 1;
3070  }
3071  }
3072  else
3073  {
3074  //call reshape for each geometry part and replace mGeos with new geometry if reshape took place
3075  bool reshapeTookPlace = false;
3076 
3077  GEOSGeometry* currentReshapeGeometry = 0;
3078  GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms];
3079 
3080  for ( int i = 0; i < numGeoms; ++i )
3081  {
3082  if ( isLine )
3083  currentReshapeGeometry = reshapeLine( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ), reshapeLineGeos );
3084  else
3085  currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ), reshapeLineGeos );
3086 
3087  if ( currentReshapeGeometry )
3088  {
3089  newGeoms[i] = currentReshapeGeometry;
3090  reshapeTookPlace = true;
3091  }
3092  else
3093  {
3094  newGeoms[i] = GEOSGeom_clone_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) );
3095  }
3096  }
3097  GEOSGeom_destroy_r( geosinit.ctxt, reshapeLineGeos );
3098 
3099  GEOSGeometry* newMultiGeom = 0;
3100  if ( isLine )
3101  {
3102  newMultiGeom = GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTILINESTRING, newGeoms, numGeoms );
3103  }
3104  else //multipolygon
3105  {
3106  newMultiGeom = GEOSGeom_createCollection_r( geosinit.ctxt, GEOS_MULTIPOLYGON, newGeoms, numGeoms );
3107  }
3108 
3109  delete[] newGeoms;
3110  if ( !newMultiGeom )
3111  return 3;
3112 
3113  if ( reshapeTookPlace )
3114  {
3115  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
3116  mGeos = newMultiGeom;
3117  mDirtyWkb = true;
3118  return 0;
3119  }
3120  else
3121  {
3122  GEOSGeom_destroy_r( geosinit.ctxt, newMultiGeom );
3123  return 1;
3124  }
3125  }
3126 }
3127 
3129 {
3130  //make sure geos geometry is up to date
3131  if ( !other )
3132  return 1;
3133 
3134  if ( mDirtyGeos )
3135  exportWkbToGeos();
3136 
3137  if ( !mGeos )
3138  return 1;
3139 
3140  if ( !GEOSisValid_r( geosinit.ctxt, mGeos ) )
3141  return 2;
3142 
3143  if ( !GEOSisSimple_r( geosinit.ctxt, mGeos ) )
3144  return 3;
3145 
3146  //convert other geometry to geos
3147  if ( other->mDirtyGeos )
3148  other->exportWkbToGeos();
3149 
3150  if ( !other->mGeos )
3151  return 4;
3152 
3153  //make geometry::difference
3154  try
3155  {
3156  if ( GEOSIntersects_r( geosinit.ctxt, mGeos, other->mGeos ) )
3157  {
3158  //check if multitype before and after
3159  bool multiType = isMultipart();
3160 
3161  mGeos = GEOSDifference_r( geosinit.ctxt, mGeos, other->mGeos );
3162  mDirtyWkb = true;
3163 
3164  if ( multiType && !isMultipart() )
3165  {
3167  exportWkbToGeos();
3168  }
3169  }
3170  else
3171  {
3172  return 0; //nothing to do
3173  }
3174  }
3175  CATCH_GEOS( 5 )
3176 
3177  if ( !mGeos )
3178  {
3179  mDirtyGeos = true;
3180  return 6;
3181  }
3182 
3183  return 0;
3184 }
3185 
3187 {
3188  double xmin = std::numeric_limits<double>::max();
3189  double ymin = std::numeric_limits<double>::max();
3190  double xmax = -std::numeric_limits<double>::max();
3191  double ymax = -std::numeric_limits<double>::max();
3192 
3193  // TODO: implement with GEOS
3194  if ( mDirtyWkb )
3195  exportGeosToWkb();
3196 
3197  if ( !mGeometry )
3198  {
3199  QgsDebugMsg( "WKB geometry not available!" );
3200  // Return minimal QgsRectangle
3201  QgsRectangle invalidRect;
3202  invalidRect.setMinimal();
3203  return invalidRect;
3204  }
3205 
3206  bool hasZValue = false;
3207  QgsWkbPtr wkbPtr( mGeometry + 1 );
3209  wkbPtr >> wkbType;
3210 
3211  // consider endian when fetching feature type
3212  switch ( wkbType )
3213  {
3214  case QGis::WKBPoint25D:
3215  case QGis::WKBPoint:
3216  {
3217  double x, y;
3218  wkbPtr >> x >> y;
3219 
3220  if ( x < xmin )
3221  xmin = x;
3222 
3223  if ( x > xmax )
3224  xmax = x;
3225 
3226  if ( y < ymin )
3227  ymin = y;
3228 
3229  if ( y > ymax )
3230  ymax = y;
3231 
3232  }
3233  break;
3234 
3236  hasZValue = true;
3237  case QGis::WKBMultiPoint:
3238  {
3239  int nPoints;
3240  wkbPtr >> nPoints;
3241  for ( int idx = 0; idx < nPoints; idx++ )
3242  {
3243  wkbPtr += 1 + sizeof( int );
3244 
3245  double x, y;
3246  wkbPtr >> x >> y;
3247  if ( hasZValue )
3248  wkbPtr += sizeof( double );
3249 
3250  if ( x < xmin )
3251  xmin = x;
3252 
3253  if ( x > xmax )
3254  xmax = x;
3255 
3256  if ( y < ymin )
3257  ymin = y;
3258 
3259  if ( y > ymax )
3260  ymax = y;
3261 
3262  }
3263  break;
3264  }
3266  hasZValue = true;
3267  case QGis::WKBLineString:
3268  {
3269  // get number of points in the line
3270  int nPoints;
3271  wkbPtr >> nPoints;
3272  for ( int idx = 0; idx < nPoints; idx++ )
3273  {
3274  double x, y;
3275  wkbPtr >> x >> y;
3276  if ( hasZValue )
3277  wkbPtr += sizeof( double );
3278 
3279  if ( x < xmin )
3280  xmin = x;
3281 
3282  if ( x > xmax )
3283  xmax = x;
3284 
3285  if ( y < ymin )
3286  ymin = y;
3287 
3288  if ( y > ymax )
3289  ymax = y;
3290 
3291  }
3292  break;
3293  }
3295  hasZValue = true;
3297  {
3298  int nLines;
3299  wkbPtr >> nLines;
3300  for ( int jdx = 0; jdx < nLines; jdx++ )
3301  {
3302  // each of these is a wbklinestring so must handle as such
3303  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3304  int nPoints;
3305  wkbPtr >> nPoints;
3306  for ( int idx = 0; idx < nPoints; idx++ )
3307  {
3308  double x, y;
3309  wkbPtr >> x >> y;
3310  if ( hasZValue )
3311  wkbPtr += sizeof( double );
3312 
3313  if ( x < xmin )
3314  xmin = x;
3315 
3316  if ( x > xmax )
3317  xmax = x;
3318 
3319  if ( y < ymin )
3320  ymin = y;
3321 
3322  if ( y > ymax )
3323  ymax = y;
3324 
3325  }
3326  }
3327  break;
3328  }
3329  case QGis::WKBPolygon25D:
3330  hasZValue = true;
3331  case QGis::WKBPolygon:
3332  {
3333  // get number of rings in the polygon
3334  int nRings;
3335  wkbPtr >> nRings;
3336  for ( int idx = 0; idx < nRings; idx++ )
3337  {
3338  // get number of points in the ring
3339  int nPoints;
3340  wkbPtr >> nPoints;
3341  for ( int jdx = 0; jdx < nPoints; jdx++ )
3342  {
3343  // add points to a point array for drawing the polygon
3344  double x, y;
3345  wkbPtr >> x >> y;
3346  if ( hasZValue )
3347  wkbPtr += sizeof( double );
3348 
3349  if ( x < xmin )
3350  xmin = x;
3351 
3352  if ( x > xmax )
3353  xmax = x;
3354 
3355  if ( y < ymin )
3356  ymin = y;
3357 
3358  if ( y > ymax )
3359  ymax = y;
3360 
3361  }
3362  }
3363  break;
3364  }
3366  hasZValue = true;
3367  case QGis::WKBMultiPolygon:
3368  {
3369  // get the number of polygons
3370  int nPolygons;
3371  wkbPtr >> nPolygons;
3372  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3373  {
3374  //skip the endian and mGeometry type info and
3375  // get number of rings in the polygon
3376  wkbPtr += 1 + sizeof( int );
3377  int nRings;
3378  wkbPtr >> nRings;
3379  for ( int idx = 0; idx < nRings; idx++ )
3380  {
3381  // get number of points in the ring
3382  int nPoints;
3383  wkbPtr >> nPoints;
3384  for ( int jdx = 0; jdx < nPoints; jdx++ )
3385  {
3386  // add points to a point array for drawing the polygon
3387  double x, y;
3388  wkbPtr >> x >> y;
3389  if ( hasZValue )
3390  wkbPtr += sizeof( double );
3391 
3392  if ( x < xmin )
3393  xmin = x;
3394 
3395  if ( x > xmax )
3396  xmax = x;
3397 
3398  if ( y < ymin )
3399  ymin = y;
3400 
3401  if ( y > ymax )
3402  ymax = y;
3403 
3404  }
3405  }
3406  }
3407  break;
3408  }
3409 
3410  default:
3411  QgsDebugMsg( QString( "Unknown WkbType %1 ENCOUNTERED" ).arg( wkbType ) );
3412  return QgsRectangle( 0, 0, 0, 0 );
3413  break;
3414 
3415  }
3416  return QgsRectangle( xmin, ymin, xmax, ymax );
3417 }
3418 
3420 {
3421  QgsGeometry* g = fromRect( r );
3422  bool res = intersects( g );
3423  delete g;
3424  return res;
3425 }
3426 
3427 bool QgsGeometry::intersects( const QgsGeometry* geometry ) const
3428 {
3429  if ( !geometry )
3430  return false;
3431 
3432  try // geos might throw exception on error
3433  {
3434  // ensure that both geometries have geos geometry
3435  exportWkbToGeos();
3436  geometry->exportWkbToGeos();
3437 
3438  if ( !mGeos || !geometry->mGeos )
3439  {
3440  QgsDebugMsg( "GEOS geometry not available!" );
3441  return false;
3442  }
3443 
3444  return GEOSIntersects_r( geosinit.ctxt, mGeos, geometry->mGeos );
3445  }
3446  CATCH_GEOS( false )
3447 }
3448 
3449 bool QgsGeometry::contains( const QgsPoint* p ) const
3450 {
3451  exportWkbToGeos();
3452 
3453  if ( !p )
3454  {
3455  QgsDebugMsg( "pointer p is 0" );
3456  return false;
3457  }
3458 
3459  if ( !mGeos )
3460  {
3461  QgsDebugMsg( "GEOS geometry not available!" );
3462  return false;
3463  }
3464 
3465  GEOSGeometry *geosPoint = 0;
3466 
3467  bool returnval = false;
3468 
3469  try
3470  {
3471  geosPoint = createGeosPoint( *p );
3472  returnval = GEOSContains_r( geosinit.ctxt, mGeos, geosPoint );
3473  }
3474  catch ( GEOSException &e )
3475  {
3476  QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
3477  returnval = false;
3478  }
3479 
3480  if ( geosPoint )
3481  GEOSGeom_destroy_r( geosinit.ctxt, geosPoint );
3482 
3483  return returnval;
3484 }
3485 
3486 bool QgsGeometry::geosRelOp(
3487  char( *op )( GEOSContextHandle_t handle, const GEOSGeometry*, const GEOSGeometry * ),
3488  const QgsGeometry *a,
3489  const QgsGeometry *b )
3490 {
3491  if ( !a || !b )
3492  return false;
3493 
3494  try // geos might throw exception on error
3495  {
3496  // ensure that both geometries have geos geometry
3497  a->exportWkbToGeos();
3498  b->exportWkbToGeos();
3499 
3500  if ( !a->mGeos || !b->mGeos )
3501  {
3502  QgsDebugMsg( "GEOS geometry not available!" );
3503  return false;
3504  }
3505  return op( geosinit.ctxt, a->mGeos, b->mGeos );
3506  }
3507  CATCH_GEOS( false )
3508 }
3509 
3510 bool QgsGeometry::contains( const QgsGeometry* geometry ) const
3511 {
3512  return geosRelOp( GEOSContains_r, this, geometry );
3513 }
3514 
3515 bool QgsGeometry::disjoint( const QgsGeometry* geometry ) const
3516 {
3517  return geosRelOp( GEOSDisjoint_r, this, geometry );
3518 }
3519 
3520 bool QgsGeometry::equals( const QgsGeometry* geometry ) const
3521 {
3522  return geosRelOp( GEOSEquals_r, this, geometry );
3523 }
3524 
3525 bool QgsGeometry::touches( const QgsGeometry* geometry ) const
3526 {
3527  return geosRelOp( GEOSTouches_r, this, geometry );
3528 }
3529 
3530 bool QgsGeometry::overlaps( const QgsGeometry* geometry ) const
3531 {
3532  return geosRelOp( GEOSOverlaps_r, this, geometry );
3533 }
3534 
3535 bool QgsGeometry::within( const QgsGeometry* geometry ) const
3536 {
3537  return geosRelOp( GEOSWithin_r, this, geometry );
3538 }
3539 
3540 bool QgsGeometry::crosses( const QgsGeometry* geometry ) const
3541 {
3542  return geosRelOp( GEOSCrosses_r, this, geometry );
3543 }
3544 
3545 QString QgsGeometry::exportToWkt( const int &precision ) const
3546 {
3547  QgsDebugMsg( "entered." );
3548 
3549  // TODO: implement with GEOS
3550  if ( mDirtyWkb )
3551  {
3552  exportGeosToWkb();
3553  }
3554 
3555  if ( !mGeometry || wkbSize() < 5 )
3556  {
3557  QgsDebugMsg( "WKB geometry not available or too short!" );
3558  return QString::null;
3559  }
3560 
3561  bool hasZValue = false;
3562  QgsWkbPtr wkbPtr( mGeometry + 1 );
3564  wkbPtr >> wkbType;
3565 
3566  QString wkt;
3567 
3568  switch ( wkbType )
3569  {
3570  case QGis::WKBPoint25D:
3571  case QGis::WKBPoint:
3572  {
3573  double x, y;
3574  wkbPtr >> x >> y;
3575  wkt += "POINT(" + qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision ) + ")";
3576  return wkt;
3577  }
3578 
3580  hasZValue = true;
3581  case QGis::WKBLineString:
3582  {
3583  int nPoints;
3584  wkbPtr >> nPoints;
3585 
3586  wkt += "LINESTRING(";
3587  // get number of points in the line
3588  for ( int idx = 0; idx < nPoints; ++idx )
3589  {
3590  double x, y;
3591  wkbPtr >> x >> y;
3592  if ( hasZValue )
3593  wkbPtr += sizeof( double );
3594 
3595  if ( idx != 0 )
3596  wkt += ", ";
3597 
3598  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3599  }
3600  wkt += ")";
3601  return wkt;
3602  }
3603 
3604  case QGis::WKBPolygon25D:
3605  hasZValue = true;
3606  case QGis::WKBPolygon:
3607  {
3608  wkt += "POLYGON(";
3609  // get number of rings in the polygon
3610  int nRings;
3611  wkbPtr >> nRings;
3612  if ( nRings == 0 ) // sanity check for zero rings in polygon
3613  return QString();
3614 
3615  for ( int idx = 0; idx < nRings; idx++ )
3616  {
3617  if ( idx != 0 )
3618  wkt += ",";
3619 
3620  wkt += "(";
3621  // get number of points in the ring
3622  int nPoints;
3623  wkbPtr >> nPoints;
3624 
3625  for ( int jdx = 0; jdx < nPoints; jdx++ )
3626  {
3627  if ( jdx != 0 )
3628  wkt += ",";
3629 
3630  double x, y;
3631  wkbPtr >> x >> y;
3632  if ( hasZValue )
3633  wkbPtr += sizeof( double );
3634 
3635  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3636  }
3637  wkt += ")";
3638  }
3639  wkt += ")";
3640  return wkt;
3641  }
3642 
3644  hasZValue = true;
3645  case QGis::WKBMultiPoint:
3646  {
3647  int nPoints;
3648  wkbPtr >> nPoints;
3649 
3650  wkt += "MULTIPOINT(";
3651  for ( int idx = 0; idx < nPoints; ++idx )
3652  {
3653  wkbPtr += 1 + sizeof( int );
3654  if ( idx != 0 )
3655  wkt += ", ";
3656 
3657  double x, y;
3658  wkbPtr >> x >> y;
3659  if ( hasZValue )
3660  wkbPtr += sizeof( double );
3661 
3662  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3663  }
3664  wkt += ")";
3665  return wkt;
3666  }
3667 
3669  hasZValue = true;
3671  {
3672  int nLines;
3673  wkbPtr >> nLines;
3674 
3675  wkt += "MULTILINESTRING(";
3676  for ( int jdx = 0; jdx < nLines; jdx++ )
3677  {
3678  if ( jdx != 0 )
3679  wkt += ", ";
3680 
3681  wkt += "(";
3682  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3683  int nPoints;
3684  wkbPtr >> nPoints;
3685  for ( int idx = 0; idx < nPoints; idx++ )
3686  {
3687  if ( idx != 0 )
3688  wkt += ", ";
3689 
3690  double x, y;
3691  wkbPtr >> x >> y;
3692  if ( hasZValue )
3693  wkbPtr += sizeof( double );
3694 
3695  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3696  }
3697  wkt += ")";
3698  }
3699  wkt += ")";
3700  return wkt;
3701  }
3702 
3704  hasZValue = true;
3705  case QGis::WKBMultiPolygon:
3706  {
3707  int nPolygons;
3708  wkbPtr >> nPolygons;
3709 
3710  wkt += "MULTIPOLYGON(";
3711  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3712  {
3713  if ( kdx != 0 )
3714  wkt += ",";
3715 
3716  wkt += "(";
3717  wkbPtr += 1 + sizeof( int );
3718  int nRings;
3719  wkbPtr >> nRings;
3720  for ( int idx = 0; idx < nRings; idx++ )
3721  {
3722  if ( idx != 0 )
3723  wkt += ",";
3724 
3725  wkt += "(";
3726  int nPoints;
3727  wkbPtr >> nPoints;
3728  for ( int jdx = 0; jdx < nPoints; jdx++ )
3729  {
3730  if ( jdx != 0 )
3731  wkt += ",";
3732 
3733  double x, y;
3734  wkbPtr >> x >> y;
3735  if ( hasZValue )
3736  wkbPtr += sizeof( double );
3737 
3738  wkt += qgsDoubleToString( x, precision ) + " " + qgsDoubleToString( y, precision );
3739  }
3740  wkt += ")";
3741  }
3742  wkt += ")";
3743  }
3744  wkt += ")";
3745  return wkt;
3746  }
3747 
3748  default:
3749  QgsDebugMsg( "error: mGeometry type not recognized" );
3750  return QString::null;
3751  }
3752 }
3753 
3754 QString QgsGeometry::exportToGeoJSON( const int &precision ) const
3755 {
3756  QgsDebugMsg( "entered." );
3757 
3758  // TODO: implement with GEOS
3759  if ( mDirtyWkb )
3760  exportGeosToWkb();
3761 
3762  if ( !mGeometry )
3763  {
3764  QgsDebugMsg( "WKB geometry not available!" );
3765  return QString::null;
3766  }
3767 
3768  QgsWkbPtr wkbPtr( mGeometry + 1 );
3770  wkbPtr >> wkbType;
3771  bool hasZValue = false;
3772 
3773  QString wkt;
3774 
3775  switch ( wkbType )
3776  {
3777  case QGis::WKBPoint25D:
3778  case QGis::WKBPoint:
3779  {
3780 
3781  double x, y;
3782  wkbPtr >> x >> y;
3783 
3784  wkt += "{ \"type\": \"Point\", \"coordinates\": ["
3785  + qgsDoubleToString( x, precision ) + ", "
3786  + qgsDoubleToString( y, precision )
3787  + "] }";
3788  return wkt;
3789  }
3790 
3792  hasZValue = true;
3793  case QGis::WKBLineString:
3794  {
3795 
3796  wkt += "{ \"type\": \"LineString\", \"coordinates\": [ ";
3797  // get number of points in the line
3798  int nPoints;
3799  wkbPtr >> nPoints;
3800  for ( int idx = 0; idx < nPoints; ++idx )
3801  {
3802  if ( idx != 0 )
3803  wkt += ", ";
3804 
3805  double x, y;
3806  wkbPtr >> x >> y;
3807  if ( hasZValue )
3808  wkbPtr += sizeof( double );
3809 
3810  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3811  }
3812  wkt += " ] }";
3813  return wkt;
3814  }
3815 
3816  case QGis::WKBPolygon25D:
3817  hasZValue = true;
3818  case QGis::WKBPolygon:
3819  {
3820 
3821  wkt += "{ \"type\": \"Polygon\", \"coordinates\": [ ";
3822  // get number of rings in the polygon
3823  int nRings;
3824  wkbPtr >> nRings;
3825  if ( nRings == 0 ) // sanity check for zero rings in polygon
3826  return QString();
3827 
3828  for ( int idx = 0; idx < nRings; idx++ )
3829  {
3830  if ( idx != 0 )
3831  wkt += ", ";
3832 
3833  wkt += "[ ";
3834  // get number of points in the ring
3835  int nPoints;
3836  wkbPtr >> nPoints;
3837 
3838  for ( int jdx = 0; jdx < nPoints; jdx++ )
3839  {
3840  if ( jdx != 0 )
3841  wkt += ", ";
3842 
3843  double x, y;
3844  wkbPtr >> x >> y;
3845  if ( hasZValue )
3846  wkbPtr += sizeof( double );
3847 
3848  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3849  }
3850  wkt += " ]";
3851  }
3852  wkt += " ] }";
3853  return wkt;
3854  }
3855 
3857  hasZValue = true;
3858  case QGis::WKBMultiPoint:
3859  {
3860  wkt += "{ \"type\": \"MultiPoint\", \"coordinates\": [ ";
3861  int nPoints;
3862  wkbPtr >> nPoints;
3863  for ( int idx = 0; idx < nPoints; ++idx )
3864  {
3865  wkbPtr += 1 + sizeof( int );
3866  if ( idx != 0 )
3867  wkt += ", ";
3868 
3869  double x, y;
3870  wkbPtr >> x >> y;
3871  if ( hasZValue )
3872  wkbPtr += sizeof( double );
3873 
3874  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3875  }
3876  wkt += " ] }";
3877  return wkt;
3878  }
3879 
3881  hasZValue = true;
3883  {
3884  wkt += "{ \"type\": \"MultiLineString\", \"coordinates\": [ ";
3885 
3886  int nLines;
3887  wkbPtr >> nLines;
3888  for ( int jdx = 0; jdx < nLines; jdx++ )
3889  {
3890  if ( jdx != 0 )
3891  wkt += ", ";
3892 
3893  wkt += "[ ";
3894  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
3895 
3896  int nPoints;
3897  wkbPtr >> nPoints;
3898  for ( int idx = 0; idx < nPoints; idx++ )
3899  {
3900  if ( idx != 0 )
3901  wkt += ", ";
3902 
3903  double x, y;
3904  wkbPtr >> x >> y;
3905  if ( hasZValue )
3906  wkbPtr += sizeof( double );
3907 
3908  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3909  }
3910  wkt += " ]";
3911  }
3912  wkt += " ] }";
3913  return wkt;
3914  }
3915 
3917  hasZValue = true;
3918  case QGis::WKBMultiPolygon:
3919  {
3920 
3921  wkt += "{ \"type\": \"MultiPolygon\", \"coordinates\": [ ";
3922 
3923  int nPolygons;
3924  wkbPtr >> nPolygons;
3925  for ( int kdx = 0; kdx < nPolygons; kdx++ )
3926  {
3927  if ( kdx != 0 )
3928  wkt += ", ";
3929 
3930  wkt += "[ ";
3931 
3932  wkbPtr += 1 + sizeof( int );
3933 
3934  int nRings;
3935  wkbPtr >> nRings;
3936  for ( int idx = 0; idx < nRings; idx++ )
3937  {
3938  if ( idx != 0 )
3939  wkt += ", ";
3940 
3941  wkt += "[ ";
3942 
3943  int nPoints;
3944  wkbPtr >> nPoints;
3945  for ( int jdx = 0; jdx < nPoints; jdx++ )
3946  {
3947  if ( jdx != 0 )
3948  wkt += ", ";
3949 
3950  double x, y;
3951  wkbPtr >> x >> y;
3952  if ( hasZValue )
3953  wkbPtr += sizeof( double );
3954 
3955  wkt += "[" + qgsDoubleToString( x, precision ) + ", " + qgsDoubleToString( y, precision ) + "]";
3956  }
3957  wkt += " ]";
3958  }
3959  wkt += " ]";
3960  }
3961  wkt += " ] }";
3962  return wkt;
3963  }
3964 
3965  default:
3966  QgsDebugMsg( "error: mGeometry type not recognized" );
3967  return QString::null;
3968  }
3969 }
3970 
3971 bool QgsGeometry::exportWkbToGeos() const
3972 {
3973  QgsDebugMsgLevel( "entered.", 3 );
3974 
3975  if ( !mDirtyGeos )
3976  {
3977  // No need to convert again
3978  return true;
3979  }
3980 
3981  if ( mGeos )
3982  {
3983  GEOSGeom_destroy_r( geosinit.ctxt, mGeos );
3984  mGeos = 0;
3985  }
3986 
3987  // this probably shouldn't return true
3988  if ( !mGeometry )
3989  {
3990  // no WKB => no GEOS
3991  mDirtyGeos = false;
3992  return true;
3993  }
3994 
3995  bool hasZValue = false;
3996 
3997  QgsWkbPtr wkbPtr( mGeometry + 1 );
3999  wkbPtr >> wkbType;
4000 
4001  try
4002  {
4003  switch ( wkbType )
4004  {
4005  case QGis::WKBPoint25D:
4006  case QGis::WKBPoint:
4007  {
4008  double x, y;
4009  wkbPtr >> x >> y;
4010  mGeos = createGeosPoint( QgsPoint( x, y ) );
4011  mDirtyGeos = false;
4012  break;
4013  }
4014 
4016  hasZValue = true;
4017  case QGis::WKBMultiPoint:
4018  {
4019  QVector<GEOSGeometry *> points;
4020 
4021  int nPoints;
4022  wkbPtr >> nPoints;
4023  for ( int idx = 0; idx < nPoints; idx++ )
4024  {
4025  double x, y;
4026  wkbPtr += 1 + sizeof( int );
4027  wkbPtr >> x >> y;
4028  if ( hasZValue )
4029  wkbPtr += sizeof( double );
4030 
4031  points << createGeosPoint( QgsPoint( x, y ) );
4032  }
4033  mGeos = createGeosCollection( GEOS_MULTIPOINT, points );
4034  mDirtyGeos = false;
4035  break;
4036  }
4037 
4039  hasZValue = true;
4040  case QGis::WKBLineString:
4041  {
4042  QgsPolyline sequence;
4043 
4044  int nPoints;
4045  wkbPtr >> nPoints;
4046  for ( int idx = 0; idx < nPoints; idx++ )
4047  {
4048  double x, y;
4049  wkbPtr >> x >> y;
4050  if ( hasZValue )
4051  wkbPtr += sizeof( double );
4052 
4053  sequence << QgsPoint( x, y );
4054  }
4055  mDirtyGeos = false;
4056  mGeos = createGeosLineString( sequence );
4057  break;
4058  }
4059 
4061  hasZValue = true;
4063  {
4064  QVector<GEOSGeometry*> lines;
4065  int nLines;
4066  wkbPtr >> nLines;
4067  for ( int jdx = 0; jdx < nLines; jdx++ )
4068  {
4069  QgsPolyline sequence;
4070 
4071  // each of these is a wbklinestring so must handle as such
4072  wkbPtr += 1 + sizeof( int ); // skip type since we know its 2
4073  int nPoints;
4074  wkbPtr >> nPoints;
4075  for ( int idx = 0; idx < nPoints; idx++ )
4076  {
4077  double x, y;
4078  wkbPtr >> x >> y;
4079  if ( hasZValue )
4080  wkbPtr += sizeof( double );
4081 
4082  sequence << QgsPoint( x, y );
4083  }
4084 
4085  // ignore invalid parts, it can come from ST_Simplify operations
4086  if ( sequence.count() > 1 )
4087  lines << createGeosLineString( sequence );
4088  }
4089  mGeos = createGeosCollection( GEOS_MULTILINESTRING, lines );
4090  mDirtyGeos = false;
4091  break;
4092  }
4093 
4094  case QGis::WKBPolygon25D:
4095  hasZValue = true;
4096  case QGis::WKBPolygon:
4097  {
4098  // get number of rings in the polygon
4099  int nRings;
4100  wkbPtr >> nRings;
4101 
4102  QVector<GEOSGeometry*> rings;
4103 
4104  for ( int idx = 0; idx < nRings; idx++ )
4105  {
4106  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4107 
4108  QgsPolyline sequence;
4109 
4110  // get number of points in the ring
4111  int nPoints;
4112  wkbPtr >> nPoints;
4113  for ( int jdx = 0; jdx < nPoints; jdx++ )
4114  {
4115  // add points to a point array for drawing the polygon
4116  double x, y;
4117  wkbPtr >> x >> y;
4118  if ( hasZValue )
4119  wkbPtr += sizeof( double );
4120 
4121  sequence << QgsPoint( x, y );
4122  }
4123 
4124  GEOSGeometry *ring = createGeosLinearRing( sequence );
4125  if ( ring ) rings << ring;
4126  }
4127  mGeos = createGeosPolygon( rings );
4128  mDirtyGeos = false;
4129  break;
4130  }
4131 
4133  hasZValue = true;
4134  case QGis::WKBMultiPolygon:
4135  {
4136  QVector<GEOSGeometry*> polygons;
4137 
4138  // get the number of polygons
4139  int nPolygons;
4140  wkbPtr >> nPolygons;
4141 
4142  for ( int kdx = 0; kdx < nPolygons; kdx++ )
4143  {
4144  //QgsDebugMsg("Polygon nr: "+QString::number(kdx));
4145  QVector<GEOSGeometry*> rings;
4146 
4147  //skip the endian and mGeometry type info and
4148  // get number of rings in the polygon
4149  wkbPtr += 1 + sizeof( int );
4150  int numRings;
4151  wkbPtr >> numRings;
4152  for ( int idx = 0; idx < numRings; idx++ )
4153  {
4154  //QgsDebugMsg("Ring nr: "+QString::number(idx));
4155 
4156  QgsPolyline sequence;
4157 
4158  // get number of points in the ring
4159  int nPoints;
4160  wkbPtr >> nPoints;
4161  for ( int jdx = 0; jdx < nPoints; jdx++ )
4162  {
4163  // add points to a point array for drawing the polygon
4164  double x, y;
4165  wkbPtr >> x >> y;
4166  if ( hasZValue )
4167  wkbPtr += sizeof( double );
4168 
4169  sequence << QgsPoint( x, y );
4170  }
4171 
4172  GEOSGeometry *ring = createGeosLinearRing( sequence );
4173  if ( ring ) rings << ring;
4174  }
4175 
4176  GEOSGeometry *polygon = createGeosPolygon( rings );
4177  if ( polygon ) polygons << polygon;
4178  }
4179  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
4180  mDirtyGeos = false;
4181  break;
4182  }
4183 
4184  default:
4185  return false;
4186  }
4187  }
4188  CATCH_GEOS( false )
4189 
4190  return true;
4191 }
4192 
4193 bool QgsGeometry::exportGeosToWkb() const
4194 {
4195  //QgsDebugMsg("entered.");
4196  if ( !mDirtyWkb )
4197  {
4198  // No need to convert again
4199  return true;
4200  }
4201 
4202  // clear the WKB, ready to replace with the new one
4203  if ( mGeometry )
4204  {
4205  delete [] mGeometry;
4206  mGeometry = 0;
4207  }
4208 
4209  if ( !mGeos )
4210  {
4211  // GEOS is null, therefore WKB is null.
4212  mDirtyWkb = false;
4213  return true;
4214  }
4215 
4216  // set up byteOrder
4217  char byteOrder = QgsApplication::endian();
4218 
4219  switch ( GEOSGeomTypeId_r( geosinit.ctxt, mGeos ) )
4220  {
4221  case GEOS_POINT: // a point
4222  {
4223  mGeometrySize = 1 +
4224  sizeof( int ) +
4225  2 * sizeof( double );
4226  mGeometry = new unsigned char[mGeometrySize];
4227 
4228 
4229  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4230 
4231  double x, y;
4232  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, 0, &x );
4233  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, 0, &y );
4234 
4235  QgsWkbPtr wkbPtr( mGeometry );
4236  wkbPtr << byteOrder << QGis::WKBPoint << x << y;
4237 
4238  mDirtyWkb = false;
4239  return true;
4240  } // case GEOS_GEOM::GEOS_POINT
4241 
4242  case GEOS_LINESTRING: // a linestring
4243  {
4244  //QgsDebugMsg("Got a geos::GEOS_LINESTRING.");
4245 
4246  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4247  unsigned int nPoints;
4248  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &nPoints );
4249 
4250  // allocate some space for the WKB
4251  mGeometrySize = 1 + // sizeof(byte)
4252  sizeof( int ) +
4253  sizeof( int ) +
4254  (( sizeof( double ) +
4255  sizeof( double ) ) * nPoints );
4256 
4257  mGeometry = new unsigned char[mGeometrySize];
4258  QgsWkbPtr wkbPtr( mGeometry );
4259 
4260  wkbPtr << byteOrder << QGis::WKBLineString << nPoints;
4261 
4262  const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq_r( geosinit.ctxt, mGeos );
4263 
4264  // assign points
4265  for ( unsigned int n = 0; n < nPoints; n++ )
4266  {
4267  double x, y;
4268  GEOSCoordSeq_getX_r( geosinit.ctxt, sequence, n, &x );
4269  GEOSCoordSeq_getY_r( geosinit.ctxt, sequence, n, &y );
4270  wkbPtr << x << y;
4271  }
4272 
4273  mDirtyWkb = false;
4274  return true;
4275 
4276  // TODO: Deal with endian-ness
4277  } // case GEOS_GEOM::GEOS_LINESTRING
4278 
4279  case GEOS_LINEARRING: // a linear ring (linestring with 1st point == last point)
4280  {
4281  // TODO
4282  break;
4283  } // case GEOS_GEOM::GEOS_LINEARRING
4284 
4285  case GEOS_POLYGON: // a polygon
4286  {
4287  int nPointsInRing = 0;
4288 
4289  //first calculate the geometry size
4290  int geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of rings
4291  const GEOSGeometry *theRing = GEOSGetExteriorRing_r( geosinit.ctxt, mGeos );
4292  if ( theRing )
4293  {
4294  geometrySize += sizeof( int );
4295  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4296  }
4297  for ( int i = 0; i < GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ); ++i )
4298  {
4299  geometrySize += sizeof( int ); //number of points in ring
4300  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, mGeos, i );
4301  if ( theRing )
4302  {
4303  geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double );
4304  }
4305  }
4306 
4307  mGeometry = new unsigned char[geometrySize];
4308  mGeometrySize = geometrySize;
4309 
4310  //then fill the geometry itself into the wkb
4311  QgsWkbPtr wkbPtr( mGeometry );
4312 
4313  int nRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ) + 1;
4314  wkbPtr << byteOrder << QGis::WKBPolygon << nRings;
4315 
4316  //exterior ring first
4317  theRing = GEOSGetExteriorRing_r( geosinit.ctxt, mGeos );
4318  if ( theRing )
4319  {
4320  nPointsInRing = getNumGeosPoints( theRing );
4321 
4322  wkbPtr << nPointsInRing;
4323 
4324  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4325  unsigned int n;
4326  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &n );
4327 
4328  for ( unsigned int j = 0; j < n; ++j )
4329  {
4330  double x, y;
4331  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4332  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4333  wkbPtr << x << y;
4334  }
4335  }
4336 
4337  //interior rings after
4338  for ( int i = 0; i < GEOSGetNumInteriorRings_r( geosinit.ctxt, mGeos ); i++ )
4339  {
4340  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, mGeos, i );
4341 
4342  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4343 
4344  unsigned int nPointsInRing;
4345  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &nPointsInRing );
4346  wkbPtr << nPointsInRing;
4347 
4348  for ( unsigned int j = 0; j < nPointsInRing; j++ )
4349  {
4350  double x, y;
4351  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4352  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4353  wkbPtr << x << y;
4354  }
4355  }
4356  mDirtyWkb = false;
4357  return true;
4358  } // case GEOS_GEOM::GEOS_POLYGON
4359  break;
4360 
4361  case GEOS_MULTIPOINT: // a collection of points
4362  {
4363  // determine size of geometry
4364  int geometrySize = 1 + 2 * sizeof( int );
4365  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4366  {
4367  geometrySize += 1 + sizeof( int ) + 2 * sizeof( double );
4368  }
4369 
4370  mGeometry = new unsigned char[geometrySize];
4371  mGeometrySize = geometrySize;
4372 
4373  QgsWkbPtr wkbPtr( mGeometry );
4374  int numPoints = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4375 
4376  wkbPtr << byteOrder << QGis::WKBMultiPoint << numPoints;
4377 
4378  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4379  {
4380  //copy endian and point type
4381  wkbPtr << byteOrder << QGis::WKBPoint;
4382 
4383  const GEOSGeometry *currentPoint = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4384 
4385  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, currentPoint );
4386 
4387  double x, y;
4388  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, 0, &x );
4389  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, 0, &y );
4390  wkbPtr << x << y;
4391  }
4392  mDirtyWkb = false;
4393  return true;
4394  } // case GEOS_GEOM::GEOS_MULTIPOINT
4395 
4396  case GEOS_MULTILINESTRING: // a collection of linestrings
4397  {
4398  // determine size of geometry
4399  int geometrySize = 1 + 2 * sizeof( int );
4400  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4401  {
4402  geometrySize += 1 + 2 * sizeof( int );
4403  geometrySize += getNumGeosPoints( GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) ) * 2 * sizeof( double );
4404  }
4405 
4406  mGeometry = new unsigned char[geometrySize];
4407  mGeometrySize = geometrySize;
4408 
4409  QgsWkbPtr wkbPtr( mGeometry );
4410 
4411  int numLines = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4412  wkbPtr << byteOrder << QGis::WKBMultiLineString << numLines;
4413 
4414  //loop over lines
4415  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4416  {
4417  //endian and type WKBLineString
4418  wkbPtr << byteOrder << QGis::WKBLineString;
4419 
4420  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i ) );
4421 
4422  //line size
4423  unsigned int lineSize;
4424  GEOSCoordSeq_getSize_r( geosinit.ctxt, cs, &lineSize );
4425  wkbPtr << lineSize;
4426 
4427  //vertex coordinates
4428  for ( unsigned int j = 0; j < lineSize; ++j )
4429  {
4430  double x, y;
4431  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, j, &x );
4432  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, j, &y );
4433  wkbPtr << x << y;
4434  }
4435  }
4436  mDirtyWkb = false;
4437  return true;
4438  } // case GEOS_GEOM::GEOS_MULTILINESTRING
4439 
4440  case GEOS_MULTIPOLYGON: // a collection of polygons
4441  {
4442  //first determine size of geometry
4443  int geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of polygons
4444  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4445  {
4446  const GEOSGeometry *thePoly = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4447  geometrySize += 1 + 2 * sizeof( int ); //endian, type, number of rings
4448  //exterior ring
4449  geometrySize += sizeof( int ); //number of points in exterior ring
4450  const GEOSGeometry *exRing = GEOSGetExteriorRing_r( geosinit.ctxt, thePoly );
4451  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( exRing );
4452 
4453  const GEOSGeometry *intRing = 0;
4454  for ( int j = 0; j < GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ); j++ )
4455  {
4456  geometrySize += sizeof( int ); //number of points in ring
4457  intRing = GEOSGetInteriorRingN_r( geosinit.ctxt, thePoly, j );
4458  geometrySize += 2 * sizeof( double ) * getNumGeosPoints( intRing );
4459  }
4460  }
4461 
4462  mGeometry = new unsigned char[geometrySize];
4463  mGeometrySize = geometrySize;
4464 
4465  QgsWkbPtr wkbPtr( mGeometry );
4466  int numPolygons = GEOSGetNumGeometries_r( geosinit.ctxt, mGeos );
4467  wkbPtr << byteOrder << QGis::WKBMultiPolygon << numPolygons;
4468 
4469  //loop over polygons
4470  for ( int i = 0; i < GEOSGetNumGeometries_r( geosinit.ctxt, mGeos ); i++ )
4471  {
4472  const GEOSGeometry *thePoly = GEOSGetGeometryN_r( geosinit.ctxt, mGeos, i );
4473  int numRings = GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ) + 1;
4474 
4475  //exterior ring
4476  const GEOSGeometry *theRing = GEOSGetExteriorRing_r( geosinit.ctxt, thePoly );
4477  int nPointsInRing = getNumGeosPoints( theRing );
4478 
4479  wkbPtr << byteOrder << QGis::WKBPolygon << numRings << nPointsInRing;
4480 
4481  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4482 
4483  for ( int k = 0; k < nPointsInRing; ++k )
4484  {
4485  double x, y;
4486  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, k, &x );
4487  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, k, &y );
4488  wkbPtr << x << y;
4489  }
4490 
4491  //interior rings
4492  for ( int j = 0; j < GEOSGetNumInteriorRings_r( geosinit.ctxt, thePoly ); j++ )
4493  {
4494  theRing = GEOSGetInteriorRingN_r( geosinit.ctxt, thePoly, j );
4495 
4496  int nPointsInRing = getNumGeosPoints( theRing );
4497  wkbPtr << nPointsInRing;
4498 
4499  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq_r( geosinit.ctxt, theRing );
4500 
4501  for ( int k = 0; k < nPointsInRing; ++k )
4502  {
4503  double x, y;
4504  GEOSCoordSeq_getX_r( geosinit.ctxt, cs, k, &x );
4505  GEOSCoordSeq_getY_r( geosinit.ctxt, cs, k, &y );
4506  wkbPtr << x << y;
4507  }
4508  }
4509  }
4510  mDirtyWkb = false;
4511  return true;
4512  } // case GEOS_GEOM::GEOS_MULTIPOLYGON
4513 
4514  case GEOS_GEOMETRYCOLLECTION: // a collection of heterogeneus geometries
4515  {
4516  // TODO
4517  QgsDebugMsg( "geometry collection - not supported" );
4518  break;
4519  } // case GEOS_GEOM::GEOS_GEOMETRYCOLLECTION
4520 
4521  } // switch (mGeos->getGeometryTypeId())
4522 
4523  return false;
4524 }
4525 
4527 {
4528  switch ( destType )
4529  {
4530  case QGis::Point:
4531  return convertToPoint( destMultipart );
4532 
4533  case QGis::Line:
4534  return convertToLine( destMultipart );
4535 
4536  case QGis::Polygon:
4537  return convertToPolygon( destMultipart );
4538 
4539  default:
4540  return 0;
4541  }
4542 }
4543 
4545 {
4546  // TODO: implement with GEOS
4547  if ( mDirtyWkb )
4548  {
4549  exportGeosToWkb();
4550  }
4551 
4552  if ( !mGeometry )
4553  {
4554  return false;
4555  }
4556 
4557  QGis::WkbType geomType = wkbType();
4558 
4559  if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D ||
4560  geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D ||
4561  geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown )
4562  {
4563  return false; //no need to convert
4564  }
4565 
4566  size_t newGeomSize = mGeometrySize + 1 + 2 * sizeof( int ); //endian: 1, multitype: sizeof(int), number of geometries: sizeof(int)
4567  unsigned char* newGeometry = new unsigned char[newGeomSize];
4568 
4569  //copy endian
4570  char byteOrder = QgsApplication::endian();
4571 
4572  QgsWkbPtr wkbPtr( newGeometry );
4573  wkbPtr << byteOrder;
4574 
4575  //copy wkbtype
4576  //todo
4577  QGis::WkbType newMultiType;
4578  switch ( geomType )
4579  {
4580  case QGis::WKBPoint:
4581  newMultiType = QGis::WKBMultiPoint;
4582  break;
4583  case QGis::WKBPoint25D:
4584  newMultiType = QGis::WKBMultiPoint25D;
4585  break;
4586  case QGis::WKBLineString:
4587  newMultiType = QGis::WKBMultiLineString;
4588  break;
4590  newMultiType = QGis::WKBMultiLineString25D;
4591  break;
4592  case QGis::WKBPolygon:
4593  newMultiType = QGis::WKBMultiPolygon;
4594  break;
4595  case QGis::WKBPolygon25D:
4596  newMultiType = QGis::WKBMultiPolygon25D;
4597  break;
4598  default:
4599  delete [] newGeometry;
4600  return false;
4601  }
4602 
4603  wkbPtr << newMultiType << 1;
4604 
4605  //copy the existing single geometry
4606  memcpy( wkbPtr, mGeometry, mGeometrySize );
4607 
4608  delete [] mGeometry;
4609  mGeometry = newGeometry;
4610  mGeometrySize = newGeomSize;
4611  mDirtyGeos = true;
4612  return true;
4613 }
4614 
4615 void QgsGeometry::translateVertex( QgsWkbPtr &wkbPtr, double dx, double dy, bool hasZValue )
4616 {
4617  double x, y, translated_x, translated_y;
4618 
4619  //x-coordinate
4620  memcpy( &x, wkbPtr, sizeof( double ) );
4621  translated_x = x + dx;
4622  memcpy( wkbPtr, &translated_x, sizeof( double ) );
4623  wkbPtr += sizeof( double );
4624 
4625  //y-coordinate
4626  memcpy( &y, wkbPtr, sizeof( double ) );
4627  translated_y = y + dy;
4628  memcpy( wkbPtr, &translated_y, sizeof( double ) );
4629  wkbPtr += sizeof( double );
4630 
4631  if ( hasZValue )
4632  wkbPtr += sizeof( double );
4633 }
4634 
4635 void QgsGeometry::transformVertex( QgsWkbPtr &wkbPtr, const QgsCoordinateTransform& ct, bool hasZValue )
4636 {
4637  double x, y, z;
4638 
4639  memcpy( &x, wkbPtr, sizeof( double ) );
4640  memcpy( &y, wkbPtr + sizeof( double ), sizeof( double ) );
4641  z = 0.0; // Ignore Z for now.
4642 
4643  ct.transformInPlace( x, y, z );
4644 
4645  // new coordinate
4646  wkbPtr << x << y;
4647  if ( hasZValue )
4648  wkbPtr += sizeof( double );
4649 
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 GEOSGeometry * createGeosCollection(int typeId, QVector< GEOSGeometry * > geoms)
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON)
Searches for the closest segment of geometry to the given point.
bool moveVertex(double x, double y, int atVertex)
Moves the vertex at the given position number and item (first number is index 0) to the given coordin...
double closestVertexWithContext(const QgsPoint &point, int &atVertex)
Searches for the closest vertex in this geometry to the given point.
QGis::GeometryType type()
Returns type of the vector.
bool crosses(const QgsGeometry *geometry) const
Test for if geometry crosses another (uses GEOS)
int addPart(const QList< QgsPoint > &points, QGis::GeometryType geomType=QGis::UnknownGeometry)
Adds a new island polygon to a multipolygon feature.
GeometryType
Definition: qgis.h:155
QgsGeometry()
Constructor.
QString what()
Definition: qgsgeometry.cpp: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:313
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 ...
QgsRectangle boundingBox()
Returns the bounding box of this feature.
#define CATCH_GEOS(r)
Definition: qgsgeometry.cpp:45
<