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