QGIS API Documentation  2.9.0-Master
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
qgsgeometry.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgeometry.cpp - Geometry (stored as Open Geospatial Consortium WKB)
3  -------------------------------------------------------------------
4 Date : 02 May 2005
5 Copyright : (C) 2005 by Brendan Morley
6 email : morb at ozemail dot com dot au
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include <limits>
17 #include <cstdarg>
18 #include <cstdio>
19 #include <cmath>
20 
21 #include "qgis.h"
22 #include "qgsgeometry.h"
23 #include "qgsapplication.h"
24 #include "qgslogger.h"
25 #include "qgsmessagelog.h"
26 #include "qgspoint.h"
27 #include "qgsrectangle.h"
28 
29 #include "qgsmaplayerregistry.h"
30 #include "qgsvectorlayer.h"
31 #include "qgsproject.h"
32 #include "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 )
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 )
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 
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 )
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 )
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 )
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 
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 
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 
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 }
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 )
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 )
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 )
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 
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 )
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.geometry() )
6321  continue;
6322 
6323  nearGeometries << GEOSGeom_clone_r( geosinit.ctxt, f.geometry()->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 )
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 )
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 )
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 )
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 }
QgsFeatureId id() const
Get the feature id for this feature.
Definition: qgsfeature.cpp:100
static QgsGeometry * fromQPolygonF(const QPolygonF &polygon)
Construct geometry from a QPolygonF.
Wrapper for iterator of features from vector data provider or vector layer.
GEOSException(QString theMsg)
Definition: qgsgeometry.cpp:54
static unsigned index
A rectangle specified with double values.
Definition: qgsrectangle.h:35
static void validateGeometry(QgsGeometry *g, QList< QgsGeometry::Error > &errors)
Validate geometry and produce a list of geometry errors.
void setMinimal()
Set a rectangle so that min corner is at max and max corner is at min.
void adjacentVertices(int atVertex, int &beforeVertex, int &afterVertex)
Returns the indexes of the vertices before and after the given vertex index.
int makeDifference(QgsGeometry *other)
Changes this geometry such that it does not intersect the other geometry.
double length()
get length of geometry using GEOS
int reshapeGeometry(const QList< QgsPoint > &reshapeWithLine)
Replaces a part of this geometry with another line.
Use exact geometry intersection (slower) instead of bounding boxes.
QPolygonF asQPolygonF() const
Return contents of the geometry as a QPolygonF.
QgsGeometry * combine(QgsGeometry *geometry)
Returns a geometry representing all the points in this geometry and other (a union geometry operation...
size_t wkbSize() const
Returns the size of the WKB in asWkb().
double yMaximum() const
Get the y maximum value (top side of rectangle)
Definition: qgsrectangle.h:188
#define QgsDebugMsg(str)
Definition: qgslogger.h:33
QSet< QgsFeatureId > QgsFeatureIds
Definition: qgsfeature.h:317
QgsGeometry * convertToType(QGis::GeometryType destType, bool destMultipart=false)
try to convert the geometry to the requested type
QgsMultiPolyline asMultiPolyline() const
return contents of the geometry as a multi linestring if wkbType is WKBMultiLineString, otherwise an empty list
QgsFeatureIterator getFeatures(const QgsFeatureRequest &request=QgsFeatureRequest())
Query the provider for features specified in request.
QVector< QgsPoint > QgsPolyline
polyline is represented as a vector of points
Definition: qgsgeometry.h:33
QPointF asQPointF() const
Return contents of the geometry as a QPointF if wkbType is WKBPoint, otherwise returns a null QPointF...
QgsGeometry * geometry() const
Get the geometry object associated with this feature.
Definition: qgsfeature.cpp:112
QgsPolygon asPolygon() const
return contents of the geometry as a polygon if wkbType is WKBPolygon, otherwise an empty list ...
QList< QgsGeometry * > asGeometryCollection() const
return contents of the geometry as a list of geometries
static QgsMapLayerRegistry * instance()
Definition: qgssingleton.h:23
static GEOSGeometry * createGeosCollection(int typeId, QVector< GEOSGeometry * > geoms)
double closestSegmentWithContext(const QgsPoint &point, QgsPoint &minDistPoint, int &afterVertex, double *leftOf=0, double epsilon=DEFAULT_SEGMENT_EPSILON)
Searches for the closest segment of geometry to the given point.
bool moveVertex(double x, double y, int atVertex)
Moves the vertex at the given position number and item (first number is index 0) to the given coordin...
QGis::GeometryType type() const
Returns type of the vector.
double closestVertexWithContext(const QgsPoint &point, int &atVertex)
Searches for the closest vertex in this geometry to the given point.
bool crosses(const QgsGeometry *geometry) const
Test for if geometry crosses another (uses GEOS)
int addPart(const QList< QgsPoint > &points, QGis::GeometryType geomType=QGis::UnknownGeometry)
Adds a new island polygon to a multipolygon feature.