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