Quantum GIS API Documentation
1.8
|
00001 /*************************************************************************** 00002 qgsgeometry.cpp - Geometry (stored as Open Geospatial Consortium WKB) 00003 ------------------------------------------------------------------- 00004 Date : 02 May 2005 00005 Copyright : (C) 2005 by Brendan Morley 00006 email : morb at ozemail dot com dot au 00007 *************************************************************************** 00008 * * 00009 * This program is free software; you can redistribute it and/or modify * 00010 * it under the terms of the GNU General Public License as published by * 00011 * the Free Software Foundation; either version 2 of the License, or * 00012 * (at your option) any later version. * 00013 * * 00014 ***************************************************************************/ 00015 00016 #include <limits> 00017 #include <cstdarg> 00018 #include <cstdio> 00019 #include <cmath> 00020 00021 #include "qgis.h" 00022 #include "qgsgeometry.h" 00023 #include "qgsapplication.h" 00024 #include "qgslogger.h" 00025 #include "qgsmessagelog.h" 00026 #include "qgspoint.h" 00027 #include "qgsrectangle.h" 00028 00029 #include "qgsmaplayerregistry.h" 00030 #include "qgsvectorlayer.h" 00031 #include "qgsproject.h" 00032 #include "qgsmessagelog.h" 00033 #include "qgsgeometryvalidator.h" 00034 00035 #define DEFAULT_QUADRANT_SEGMENTS 8 00036 00037 #define CATCH_GEOS(r) \ 00038 catch (GEOSException &e) \ 00039 { \ 00040 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr("GEOS") ); \ 00041 return r; \ 00042 } 00043 00044 class GEOSException 00045 { 00046 public: 00047 GEOSException( QString theMsg ) 00048 { 00049 if ( theMsg == "Unknown exception thrown" && lastMsg.isNull() ) 00050 { 00051 msg = theMsg; 00052 } 00053 else 00054 { 00055 msg = theMsg; 00056 lastMsg = msg; 00057 } 00058 } 00059 00060 // copy constructor 00061 GEOSException( const GEOSException &rhs ) 00062 { 00063 *this = rhs; 00064 } 00065 00066 ~GEOSException() 00067 { 00068 if ( lastMsg == msg ) 00069 lastMsg = QString::null; 00070 } 00071 00072 QString what() 00073 { 00074 return msg; 00075 } 00076 00077 private: 00078 QString msg; 00079 static QString lastMsg; 00080 }; 00081 00082 QString GEOSException::lastMsg; 00083 00084 static void throwGEOSException( const char *fmt, ... ) 00085 { 00086 va_list ap; 00087 char buffer[1024]; 00088 00089 va_start( ap, fmt ); 00090 vsnprintf( buffer, sizeof buffer, fmt, ap ); 00091 va_end( ap ); 00092 00093 QgsDebugMsg( QString( "GEOS exception: %1" ).arg( buffer ) ); 00094 00095 throw GEOSException( QString::fromUtf8( buffer ) ); 00096 } 00097 00098 static void printGEOSNotice( const char *fmt, ... ) 00099 { 00100 #if defined(QGISDEBUG) 00101 va_list ap; 00102 char buffer[1024]; 00103 00104 va_start( ap, fmt ); 00105 vsnprintf( buffer, sizeof buffer, fmt, ap ); 00106 va_end( ap ); 00107 00108 QgsDebugMsg( QString( "GEOS notice: %1" ).arg( QString::fromUtf8( buffer ) ) ); 00109 #else 00110 Q_UNUSED( fmt ); 00111 #endif 00112 } 00113 00114 class GEOSInit 00115 { 00116 public: 00117 GEOSInit() 00118 { 00119 initGEOS( printGEOSNotice, throwGEOSException ); 00120 } 00121 00122 ~GEOSInit() 00123 { 00124 finishGEOS(); 00125 } 00126 }; 00127 00128 static GEOSInit geosinit; 00129 00130 00131 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR<3) 00132 #define GEOSGeom_getCoordSeq(g) GEOSGeom_getCoordSeq( (GEOSGeometry *) g ) 00133 #define GEOSGetExteriorRing(g) GEOSGetExteriorRing( (GEOSGeometry *)g ) 00134 #define GEOSGetNumInteriorRings(g) GEOSGetNumInteriorRings( (GEOSGeometry *)g ) 00135 #define GEOSGetInteriorRingN(g,i) GEOSGetInteriorRingN( (GEOSGeometry *)g, i ) 00136 #define GEOSDisjoint(g0,g1) GEOSDisjoint( (GEOSGeometry *)g0, (GEOSGeometry*)g1 ) 00137 #define GEOSIntersection(g0,g1) GEOSIntersection( (GEOSGeometry*) g0, (GEOSGeometry*)g1 ) 00138 #define GEOSBuffer(g, d, s) GEOSBuffer( (GEOSGeometry*) g, d, s ) 00139 #define GEOSArea(g, a) GEOSArea( (GEOSGeometry*) g, a ) 00140 #define GEOSSimplify(g, t) GEOSSimplify( (GEOSGeometry*) g, t ) 00141 #define GEOSGetCentroid(g) GEOSGetCentroid( (GEOSGeometry*) g ) 00142 00143 #define GEOSCoordSeq_getSize(cs,n) GEOSCoordSeq_getSize( (GEOSCoordSequence *) cs, n ) 00144 #define GEOSCoordSeq_getX(cs,i,x) GEOSCoordSeq_getX( (GEOSCoordSequence *)cs, i, x ) 00145 #define GEOSCoordSeq_getY(cs,i,y) GEOSCoordSeq_getY( (GEOSCoordSequence *)cs, i, y ) 00146 00147 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms ); 00148 00149 static GEOSGeometry *cloneGeosGeom( const GEOSGeometry *geom ) 00150 { 00151 // for GEOS < 3.0 we have own cloning function 00152 // because when cloning multipart geometries they're copied into more general geometry collection instance 00153 int type = GEOSGeomTypeId(( GEOSGeometry * ) geom ); 00154 00155 if ( type == GEOS_MULTIPOINT || type == GEOS_MULTILINESTRING || type == GEOS_MULTIPOLYGON ) 00156 { 00157 QVector<GEOSGeometry *> geoms; 00158 00159 try 00160 { 00161 00162 for ( int i = 0; i < GEOSGetNumGeometries(( GEOSGeometry * )geom ); ++i ) 00163 geoms << GEOSGeom_clone(( GEOSGeometry * ) GEOSGetGeometryN(( GEOSGeometry * ) geom, i ) ); 00164 00165 return createGeosCollection( type, geoms ); 00166 } 00167 catch ( GEOSException &e ) 00168 { 00169 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00170 for ( int i = 0; i < geoms.count(); i++ ) 00171 GEOSGeom_destroy( geoms[i] ); 00172 00173 return 0; 00174 } 00175 } 00176 else 00177 { 00178 return GEOSGeom_clone(( GEOSGeometry * ) geom ); 00179 } 00180 } 00181 00182 #define GEOSGeom_clone(g) cloneGeosGeom(g) 00183 #endif 00184 00185 QgsGeometry::QgsGeometry() 00186 : mGeometry( 0 ) 00187 , mGeometrySize( 0 ) 00188 , mGeos( 0 ) 00189 , mDirtyWkb( false ) 00190 , mDirtyGeos( false ) 00191 { 00192 } 00193 00194 QgsGeometry::QgsGeometry( QgsGeometry const & rhs ) 00195 : mGeometry( 0 ) 00196 , mGeometrySize( rhs.mGeometrySize ) 00197 , mDirtyWkb( rhs.mDirtyWkb ) 00198 , mDirtyGeos( rhs.mDirtyGeos ) 00199 { 00200 if ( mGeometrySize && rhs.mGeometry ) 00201 { 00202 mGeometry = new unsigned char[mGeometrySize]; 00203 memcpy( mGeometry, rhs.mGeometry, mGeometrySize ); 00204 } 00205 00206 // deep-copy the GEOS Geometry if appropriate 00207 if ( rhs.mGeos ) 00208 { 00209 mGeos = GEOSGeom_clone( rhs.mGeos ); 00210 } 00211 else 00212 { 00213 mGeos = 0; 00214 } 00215 } 00216 00218 QgsGeometry::~QgsGeometry() 00219 { 00220 if ( mGeometry ) 00221 { 00222 delete [] mGeometry; 00223 } 00224 00225 if ( mGeos ) 00226 { 00227 GEOSGeom_destroy( mGeos ); 00228 } 00229 } 00230 00231 static unsigned int getNumGeosPoints( const GEOSGeometry *geom ) 00232 { 00233 unsigned int n; 00234 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( geom ); 00235 GEOSCoordSeq_getSize( cs, &n ); 00236 return n; 00237 } 00238 00239 static GEOSGeometry *createGeosPoint( const QgsPoint &point ) 00240 { 00241 GEOSCoordSequence *coord = GEOSCoordSeq_create( 1, 2 ); 00242 GEOSCoordSeq_setX( coord, 0, point.x() ); 00243 GEOSCoordSeq_setY( coord, 0, point.y() ); 00244 return GEOSGeom_createPoint( coord ); 00245 } 00246 00247 static GEOSCoordSequence *createGeosCoordSequence( const QgsPolyline& points ) 00248 { 00249 GEOSCoordSequence *coord = 0; 00250 00251 try 00252 { 00253 coord = GEOSCoordSeq_create( points.count(), 2 ); 00254 int i; 00255 for ( i = 0; i < points.count(); i++ ) 00256 { 00257 GEOSCoordSeq_setX( coord, i, points[i].x() ); 00258 GEOSCoordSeq_setY( coord, i, points[i].y() ); 00259 } 00260 return coord; 00261 } 00262 catch ( GEOSException &e ) 00263 { 00264 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00265 /*if ( coord ) 00266 GEOSCoordSeq_destroy( coord );*/ 00267 throw; 00268 } 00269 } 00270 00271 static GEOSGeometry *createGeosCollection( int typeId, QVector<GEOSGeometry*> geoms ) 00272 { 00273 GEOSGeometry **geomarr = new GEOSGeometry*[ geoms.size()]; 00274 if ( !geomarr ) 00275 return 0; 00276 00277 for ( int i = 0; i < geoms.size(); i++ ) 00278 geomarr[i] = geoms[i]; 00279 00280 GEOSGeometry *geom = 0; 00281 00282 try 00283 { 00284 geom = GEOSGeom_createCollection( typeId, geomarr, geoms.size() ); 00285 } 00286 catch ( GEOSException &e ) 00287 { 00288 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00289 } 00290 00291 delete [] geomarr; 00292 00293 return geom; 00294 } 00295 00296 static GEOSGeometry *createGeosLineString( const QgsPolyline& polyline ) 00297 { 00298 GEOSCoordSequence *coord = 0; 00299 00300 try 00301 { 00302 coord = createGeosCoordSequence( polyline ); 00303 return GEOSGeom_createLineString( coord ); 00304 } 00305 catch ( GEOSException &e ) 00306 { 00307 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00308 //MH: for strange reasons, geos3 crashes when removing the coordinate sequence 00309 //if ( coord ) 00310 //GEOSCoordSeq_destroy( coord ); 00311 return 0; 00312 } 00313 } 00314 00315 static GEOSGeometry *createGeosLinearRing( const QgsPolyline& polyline ) 00316 { 00317 GEOSCoordSequence *coord = 0; 00318 00319 if ( polyline.count() == 0 ) 00320 return 0; 00321 00322 try 00323 { 00324 if ( polyline[0] != polyline[polyline.size()-1] ) 00325 { 00326 // Ring not closed 00327 QgsPolyline closed( polyline ); 00328 closed << closed[0]; 00329 coord = createGeosCoordSequence( closed ); 00330 } 00331 else 00332 { 00333 coord = createGeosCoordSequence( polyline ); 00334 } 00335 00336 return GEOSGeom_createLinearRing( coord ); 00337 } 00338 catch ( GEOSException &e ) 00339 { 00340 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00341 /* as MH has noticed ^, this crashes geos 00342 if ( coord ) 00343 GEOSCoordSeq_destroy( coord );*/ 00344 return 0; 00345 } 00346 } 00347 00348 static GEOSGeometry *createGeosPolygon( const QVector<GEOSGeometry*> &rings ) 00349 { 00350 GEOSGeometry *shell = rings[0]; 00351 GEOSGeometry **holes = NULL; 00352 00353 if ( rings.size() > 1 ) 00354 { 00355 holes = new GEOSGeometry*[ rings.size()-1 ]; 00356 if ( !holes ) 00357 return 0; 00358 00359 for ( int i = 0; i < rings.size() - 1; i++ ) 00360 holes[i] = rings[i+1]; 00361 } 00362 00363 GEOSGeometry *geom = GEOSGeom_createPolygon( shell, holes, rings.size() - 1 ); 00364 00365 if ( holes ) 00366 delete [] holes; 00367 00368 return geom; 00369 } 00370 00371 static GEOSGeometry *createGeosPolygon( GEOSGeometry *shell ) 00372 { 00373 return createGeosPolygon( QVector<GEOSGeometry*>() << shell ); 00374 } 00375 00376 static GEOSGeometry *createGeosPolygon( const QgsPolygon& polygon ) 00377 { 00378 if ( polygon.count() == 0 ) 00379 return 0; 00380 00381 QVector<GEOSGeometry *> geoms; 00382 00383 try 00384 { 00385 for ( int i = 0; i < polygon.count(); i++ ) 00386 geoms << createGeosLinearRing( polygon[i] ); 00387 00388 return createGeosPolygon( geoms ); 00389 } 00390 catch ( GEOSException &e ) 00391 { 00392 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00393 for ( int i = 0; i < geoms.count(); i++ ) 00394 GEOSGeom_destroy( geoms[i] ); 00395 return 0; 00396 } 00397 } 00398 00399 static QgsGeometry *fromGeosGeom( GEOSGeometry *geom ) 00400 { 00401 if ( !geom ) 00402 return 0; 00403 00404 QgsGeometry* g = new QgsGeometry; 00405 g->fromGeos( geom ); 00406 return g; 00407 } 00408 00409 QgsGeometry* QgsGeometry::fromWkt( QString wkt ) 00410 { 00411 try 00412 { 00413 #if defined(GEOS_VERSION_MAJOR) && (GEOS_VERSION_MAJOR>=3) 00414 GEOSWKTReader *reader = GEOSWKTReader_create(); 00415 QgsGeometry *g = fromGeosGeom( GEOSWKTReader_read( reader, wkt.toLocal8Bit().data() ) ); 00416 GEOSWKTReader_destroy( reader ); 00417 return g; 00418 #else 00419 return fromGeosGeom( GEOSGeomFromWKT( wkt.toLocal8Bit().data() ) ); 00420 #endif 00421 } 00422 catch ( GEOSException &e ) 00423 { 00424 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00425 return 0; 00426 } 00427 } 00428 00429 QgsGeometry* QgsGeometry::fromPoint( const QgsPoint& point ) 00430 { 00431 return fromGeosGeom( createGeosPoint( point ) ); 00432 } 00433 00434 QgsGeometry* QgsGeometry::fromPolyline( const QgsPolyline& polyline ) 00435 { 00436 return fromGeosGeom( createGeosLineString( polyline ) ); 00437 } 00438 00439 QgsGeometry* QgsGeometry::fromPolygon( const QgsPolygon& polygon ) 00440 { 00441 return fromGeosGeom( createGeosPolygon( polygon ) ); 00442 } 00443 00444 QgsGeometry* QgsGeometry::fromMultiPoint( const QgsMultiPoint& multipoint ) 00445 { 00446 QVector<GEOSGeometry *> geoms; 00447 00448 try 00449 { 00450 for ( int i = 0; i < multipoint.size(); ++i ) 00451 { 00452 geoms << createGeosPoint( multipoint[i] ); 00453 } 00454 00455 return fromGeosGeom( createGeosCollection( GEOS_MULTIPOINT, geoms ) ); 00456 } 00457 catch ( GEOSException &e ) 00458 { 00459 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00460 00461 for ( int i = 0; i < geoms.size(); ++i ) 00462 GEOSGeom_destroy( geoms[i] ); 00463 00464 return 0; 00465 } 00466 } 00467 00468 QgsGeometry* QgsGeometry::fromMultiPolyline( const QgsMultiPolyline& multiline ) 00469 { 00470 QVector<GEOSGeometry *> geoms; 00471 00472 try 00473 { 00474 for ( int i = 0; i < multiline.count(); i++ ) 00475 geoms << createGeosLineString( multiline[i] ); 00476 00477 return fromGeosGeom( createGeosCollection( GEOS_MULTILINESTRING, geoms ) ); 00478 } 00479 catch ( GEOSException &e ) 00480 { 00481 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00482 00483 for ( int i = 0; i < geoms.count(); i++ ) 00484 GEOSGeom_destroy( geoms[i] ); 00485 00486 return 0; 00487 } 00488 } 00489 00490 QgsGeometry* QgsGeometry::fromMultiPolygon( const QgsMultiPolygon& multipoly ) 00491 { 00492 if ( multipoly.count() == 0 ) 00493 return 0; 00494 00495 QVector<GEOSGeometry *> geoms; 00496 00497 try 00498 { 00499 for ( int i = 0; i < multipoly.count(); i++ ) 00500 geoms << createGeosPolygon( multipoly[i] ); 00501 00502 return fromGeosGeom( createGeosCollection( GEOS_MULTIPOLYGON, geoms ) ); 00503 } 00504 catch ( GEOSException &e ) 00505 { 00506 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 00507 00508 for ( int i = 0; i < geoms.count(); i++ ) 00509 GEOSGeom_destroy( geoms[i] ); 00510 00511 return 0; 00512 } 00513 } 00514 00515 QgsGeometry* QgsGeometry::fromRect( const QgsRectangle& rect ) 00516 { 00517 QgsPolyline ring; 00518 ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) ); 00519 ring.append( QgsPoint( rect.xMaximum(), rect.yMinimum() ) ); 00520 ring.append( QgsPoint( rect.xMaximum(), rect.yMaximum() ) ); 00521 ring.append( QgsPoint( rect.xMinimum(), rect.yMaximum() ) ); 00522 ring.append( QgsPoint( rect.xMinimum(), rect.yMinimum() ) ); 00523 00524 QgsPolygon polygon; 00525 polygon.append( ring ); 00526 00527 return fromPolygon( polygon ); 00528 } 00529 00530 QgsGeometry & QgsGeometry::operator=( QgsGeometry const & rhs ) 00531 { 00532 if ( &rhs == this ) 00533 { 00534 return *this; 00535 } 00536 00537 // remove old geometry if it exists 00538 if ( mGeometry ) 00539 { 00540 delete [] mGeometry; 00541 mGeometry = 0; 00542 } 00543 00544 mGeometrySize = rhs.mGeometrySize; 00545 00546 // deep-copy the GEOS Geometry if appropriate 00547 GEOSGeom_destroy( mGeos ); 00548 mGeos = rhs.mGeos ? GEOSGeom_clone( rhs.mGeos ) : 0; 00549 00550 mDirtyGeos = rhs.mDirtyGeos; 00551 mDirtyWkb = rhs.mDirtyWkb; 00552 00553 if ( mGeometrySize && rhs.mGeometry ) 00554 { 00555 mGeometry = new unsigned char[mGeometrySize]; 00556 memcpy( mGeometry, rhs.mGeometry, mGeometrySize ); 00557 } 00558 00559 return *this; 00560 } // QgsGeometry::operator=( QgsGeometry const & rhs ) 00561 00562 00563 void QgsGeometry::fromWkb( unsigned char * wkb, size_t length ) 00564 { 00565 // delete any existing WKB geometry before assigning new one 00566 if ( mGeometry ) 00567 { 00568 delete [] mGeometry; 00569 mGeometry = 0; 00570 } 00571 if ( mGeos ) 00572 { 00573 GEOSGeom_destroy( mGeos ); 00574 mGeos = 0; 00575 } 00576 00577 mGeometry = wkb; 00578 mGeometrySize = length; 00579 00580 mDirtyWkb = false; 00581 mDirtyGeos = true; 00582 } 00583 00584 unsigned char * QgsGeometry::asWkb() 00585 { 00586 if ( mDirtyWkb ) 00587 { 00588 // convert from GEOS 00589 exportGeosToWkb(); 00590 } 00591 00592 return mGeometry; 00593 } 00594 00595 size_t QgsGeometry::wkbSize() 00596 { 00597 if ( mDirtyWkb ) 00598 { 00599 exportGeosToWkb(); 00600 } 00601 00602 return mGeometrySize; 00603 } 00604 00605 GEOSGeometry* QgsGeometry::asGeos() 00606 { 00607 if ( mDirtyGeos ) 00608 { 00609 if ( !exportWkbToGeos() ) 00610 { 00611 return 0; 00612 } 00613 } 00614 return mGeos; 00615 } 00616 00617 00618 QGis::WkbType QgsGeometry::wkbType() 00619 { 00620 unsigned char *geom = asWkb(); // ensure that wkb representation exists 00621 if ( geom ) 00622 { 00623 unsigned int wkbType; 00624 memcpy( &wkbType, ( geom + 1 ), sizeof( wkbType ) ); 00625 return ( QGis::WkbType ) wkbType; 00626 } 00627 else 00628 { 00629 return QGis::WKBUnknown; 00630 } 00631 } 00632 00633 00634 QGis::GeometryType QgsGeometry::type() 00635 { 00636 if ( mDirtyWkb ) 00637 { 00638 // convert from GEOS 00639 exportGeosToWkb(); 00640 } 00641 00642 QGis::WkbType type = wkbType(); 00643 if ( type == QGis::WKBPoint || type == QGis::WKBPoint25D || 00644 type == QGis::WKBMultiPoint || type == QGis::WKBMultiPoint25D ) 00645 return QGis::Point; 00646 if ( type == QGis::WKBLineString || type == QGis::WKBLineString25D || 00647 type == QGis::WKBMultiLineString || type == QGis::WKBMultiLineString25D ) 00648 return QGis::Line; 00649 if ( type == QGis::WKBPolygon || type == QGis::WKBPolygon25D || 00650 type == QGis::WKBMultiPolygon || type == QGis::WKBMultiPolygon25D ) 00651 return QGis::Polygon; 00652 00653 return QGis::UnknownGeometry; 00654 } 00655 00656 bool QgsGeometry::isMultipart() 00657 { 00658 if ( mDirtyWkb ) 00659 { 00660 // convert from GEOS 00661 exportGeosToWkb(); 00662 } 00663 00664 QGis::WkbType type = wkbType(); 00665 if ( type == QGis::WKBMultiPoint || 00666 type == QGis::WKBMultiPoint25D || 00667 type == QGis::WKBMultiLineString || 00668 type == QGis::WKBMultiLineString25D || 00669 type == QGis::WKBMultiPolygon || 00670 type == QGis::WKBMultiPolygon25D ) 00671 return true; 00672 00673 return false; 00674 } 00675 00676 00677 void QgsGeometry::fromGeos( GEOSGeometry* geos ) 00678 { 00679 // TODO - make this more heap-friendly 00680 00681 if ( mGeos ) 00682 { 00683 GEOSGeom_destroy( mGeos ); 00684 mGeos = 0; 00685 } 00686 if ( mGeometry ) 00687 { 00688 delete [] mGeometry; 00689 mGeometry = 0; 00690 } 00691 00692 mGeos = geos; 00693 00694 mDirtyWkb = true; 00695 mDirtyGeos = false; 00696 } 00697 00698 QgsPoint QgsGeometry::closestVertex( const QgsPoint& point, int& atVertex, int& beforeVertex, int& afterVertex, double& sqrDist ) 00699 { 00700 // TODO: implement with GEOS 00701 if ( mDirtyWkb ) 00702 { 00703 exportGeosToWkb(); 00704 } 00705 00706 if ( !mGeometry ) 00707 { 00708 QgsDebugMsg( "WKB geometry not available!" ); 00709 return QgsPoint( 0, 0 ); 00710 } 00711 00712 int vertexnr = -1; 00713 int vertexcounter = 0; 00714 QGis::WkbType wkbType; 00715 double actdist = std::numeric_limits<double>::max(); 00716 double x = 0; 00717 double y = 0; 00718 double *tempx, *tempy; 00719 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) ); 00720 beforeVertex = -1; 00721 afterVertex = -1; 00722 bool hasZValue = false; 00723 00724 switch ( wkbType ) 00725 { 00726 case QGis::WKBPoint25D: 00727 case QGis::WKBPoint: 00728 { 00729 x = *(( double * )( mGeometry + 5 ) ); 00730 y = *(( double * )( mGeometry + 5 + sizeof( double ) ) ); 00731 actdist = point.sqrDist( x, y ); 00732 vertexnr = 0; 00733 break; 00734 } 00735 case QGis::WKBLineString25D: 00736 hasZValue = true; 00737 00738 // fall-through 00739 00740 case QGis::WKBLineString: 00741 { 00742 unsigned char* ptr = mGeometry + 5; 00743 int* npoints = ( int* )ptr; 00744 ptr += sizeof( int ); 00745 for ( int index = 0; index < *npoints; ++index ) 00746 { 00747 tempx = ( double* )ptr; 00748 ptr += sizeof( double ); 00749 tempy = ( double* )ptr; 00750 if ( point.sqrDist( *tempx, *tempy ) < actdist ) 00751 { 00752 x = *tempx; 00753 y = *tempy; 00754 actdist = point.sqrDist( *tempx, *tempy ); 00755 vertexnr = index; 00756 if ( index == 0 )//assign the rubber band indices 00757 { 00758 beforeVertex = -1; 00759 } 00760 else 00761 { 00762 beforeVertex = index - 1; 00763 } 00764 if ( index == ( *npoints - 1 ) ) 00765 { 00766 afterVertex = -1; 00767 } 00768 else 00769 { 00770 afterVertex = index + 1; 00771 } 00772 } 00773 ptr += sizeof( double ); 00774 if ( hasZValue ) //skip z-coordinate for 25D geometries 00775 { 00776 ptr += sizeof( double ); 00777 } 00778 } 00779 break; 00780 } 00781 case QGis::WKBPolygon25D: 00782 hasZValue = true; 00783 case QGis::WKBPolygon: 00784 { 00785 int* nrings = ( int* )( mGeometry + 5 ); 00786 int* npoints; 00787 unsigned char* ptr = mGeometry + 9; 00788 for ( int index = 0; index < *nrings; ++index ) 00789 { 00790 npoints = ( int* )ptr; 00791 ptr += sizeof( int ); 00792 for ( int index2 = 0; index2 < *npoints; ++index2 ) 00793 { 00794 tempx = ( double* )ptr; 00795 ptr += sizeof( double ); 00796 tempy = ( double* )ptr; 00797 if ( point.sqrDist( *tempx, *tempy ) < actdist ) 00798 { 00799 x = *tempx; 00800 y = *tempy; 00801 actdist = point.sqrDist( *tempx, *tempy ); 00802 vertexnr = vertexcounter; 00803 //assign the rubber band indices 00804 if ( index2 == 0 ) 00805 { 00806 beforeVertex = vertexcounter + ( *npoints - 2 ); 00807 afterVertex = vertexcounter + 1; 00808 } 00809 else if ( index2 == ( *npoints - 1 ) ) 00810 { 00811 beforeVertex = vertexcounter - 1; 00812 afterVertex = vertexcounter - ( *npoints - 2 ); 00813 } 00814 else 00815 { 00816 beforeVertex = vertexcounter - 1; 00817 afterVertex = vertexcounter + 1; 00818 } 00819 } 00820 ptr += sizeof( double ); 00821 if ( hasZValue ) //skip z-coordinate for 25D geometries 00822 { 00823 ptr += sizeof( double ); 00824 } 00825 ++vertexcounter; 00826 } 00827 } 00828 break; 00829 } 00830 case QGis::WKBMultiPoint25D: 00831 hasZValue = true; 00832 case QGis::WKBMultiPoint: 00833 { 00834 unsigned char* ptr = mGeometry + 5; 00835 int* npoints = ( int* )ptr; 00836 ptr += sizeof( int ); 00837 for ( int index = 0; index < *npoints; ++index ) 00838 { 00839 ptr += ( 1 + sizeof( int ) ); //skip endian and point type 00840 tempx = ( double* )ptr; 00841 tempy = ( double* )( ptr + sizeof( double ) ); 00842 if ( point.sqrDist( *tempx, *tempy ) < actdist ) 00843 { 00844 x = *tempx; 00845 y = *tempy; 00846 actdist = point.sqrDist( *tempx, *tempy ); 00847 vertexnr = index; 00848 } 00849 ptr += ( 2 * sizeof( double ) ); 00850 if ( hasZValue ) //skip z-coordinate for 25D geometries 00851 { 00852 ptr += sizeof( double ); 00853 } 00854 } 00855 break; 00856 } 00857 case QGis::WKBMultiLineString25D: 00858 hasZValue = true; 00859 case QGis::WKBMultiLineString: 00860 { 00861 unsigned char* ptr = mGeometry + 5; 00862 int* nlines = ( int* )ptr; 00863 int* npoints = 0; 00864 ptr += sizeof( int ); 00865 for ( int index = 0; index < *nlines; ++index ) 00866 { 00867 ptr += ( sizeof( int ) + 1 ); 00868 npoints = ( int* )ptr; 00869 ptr += sizeof( int ); 00870 for ( int index2 = 0; index2 < *npoints; ++index2 ) 00871 { 00872 tempx = ( double* )ptr; 00873 ptr += sizeof( double ); 00874 tempy = ( double* )ptr; 00875 ptr += sizeof( double ); 00876 if ( point.sqrDist( *tempx, *tempy ) < actdist ) 00877 { 00878 x = *tempx; 00879 y = *tempy; 00880 actdist = point.sqrDist( *tempx, *tempy ); 00881 vertexnr = vertexcounter; 00882 00883 if ( index2 == 0 )//assign the rubber band indices 00884 { 00885 beforeVertex = -1; 00886 } 00887 else 00888 { 00889 beforeVertex = vertexnr - 1; 00890 } 00891 if ( index2 == ( *npoints ) - 1 ) 00892 { 00893 afterVertex = -1; 00894 } 00895 else 00896 { 00897 afterVertex = vertexnr + 1; 00898 } 00899 } 00900 if ( hasZValue ) //skip z-coordinate for 25D geometries 00901 { 00902 ptr += sizeof( double ); 00903 } 00904 ++vertexcounter; 00905 } 00906 } 00907 break; 00908 } 00909 case QGis::WKBMultiPolygon25D: 00910 hasZValue = true; 00911 case QGis::WKBMultiPolygon: 00912 { 00913 unsigned char* ptr = mGeometry + 5; 00914 int* npolys = ( int* )ptr; 00915 int* nrings; 00916 int* npoints; 00917 ptr += sizeof( int ); 00918 for ( int index = 0; index < *npolys; ++index ) 00919 { 00920 ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type 00921 nrings = ( int* )ptr; 00922 ptr += sizeof( int ); 00923 for ( int index2 = 0; index2 < *nrings; ++index2 ) 00924 { 00925 npoints = ( int* )ptr; 00926 ptr += sizeof( int ); 00927 for ( int index3 = 0; index3 < *npoints; ++index3 ) 00928 { 00929 tempx = ( double* )ptr; 00930 ptr += sizeof( double ); 00931 tempy = ( double* )ptr; 00932 if ( point.sqrDist( *tempx, *tempy ) < actdist ) 00933 { 00934 x = *tempx; 00935 y = *tempy; 00936 actdist = point.sqrDist( *tempx, *tempy ); 00937 vertexnr = vertexcounter; 00938 00939 //assign the rubber band indices 00940 if ( index3 == 0 ) 00941 { 00942 beforeVertex = vertexcounter + ( *npoints - 2 ); 00943 afterVertex = vertexcounter + 1; 00944 } 00945 else if ( index3 == ( *npoints - 1 ) ) 00946 { 00947 beforeVertex = vertexcounter - 1; 00948 afterVertex = vertexcounter - ( *npoints - 2 ); 00949 } 00950 else 00951 { 00952 beforeVertex = vertexcounter - 1; 00953 afterVertex = vertexcounter + 1; 00954 } 00955 } 00956 ptr += sizeof( double ); 00957 if ( hasZValue ) //skip z-coordinate for 25D geometries 00958 { 00959 ptr += sizeof( double ); 00960 } 00961 ++vertexcounter; 00962 } 00963 } 00964 } 00965 break; 00966 } 00967 default: 00968 break; 00969 } 00970 sqrDist = actdist; 00971 atVertex = vertexnr; 00972 return QgsPoint( x, y ); 00973 } 00974 00975 00976 void QgsGeometry::adjacentVertices( int atVertex, int& beforeVertex, int& afterVertex ) 00977 { 00978 // TODO: implement with GEOS 00979 if ( mDirtyWkb ) 00980 { 00981 exportGeosToWkb(); 00982 } 00983 00984 beforeVertex = -1; 00985 afterVertex = -1; 00986 00987 if ( !mGeometry ) 00988 { 00989 QgsDebugMsg( "WKB geometry not available!" ); 00990 return; 00991 } 00992 00993 int vertexcounter = 0; 00994 00995 QGis::WkbType wkbType; 00996 bool hasZValue = false; 00997 00998 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) ); 00999 01000 switch ( wkbType ) 01001 { 01002 case QGis::WKBPoint: 01003 { 01004 // NOOP - Points do not have adjacent verticies 01005 break; 01006 } 01007 case QGis::WKBLineString25D: 01008 case QGis::WKBLineString: 01009 { 01010 unsigned char* ptr = mGeometry + 5; 01011 int* npoints = ( int* ) ptr; 01012 01013 const int index = atVertex; 01014 01015 // assign the rubber band indices 01016 01017 if ( index == 0 ) 01018 { 01019 beforeVertex = -1; 01020 } 01021 else 01022 { 01023 beforeVertex = index - 1; 01024 } 01025 01026 if ( index == ( *npoints - 1 ) ) 01027 { 01028 afterVertex = -1; 01029 } 01030 else 01031 { 01032 afterVertex = index + 1; 01033 } 01034 01035 break; 01036 } 01037 case QGis::WKBPolygon25D: 01038 hasZValue = true; 01039 case QGis::WKBPolygon: 01040 { 01041 int* nrings = ( int* )( mGeometry + 5 ); 01042 int* npoints; 01043 unsigned char* ptr = mGeometry + 9; 01044 01045 // Walk through the POLYGON WKB 01046 01047 for ( int index0 = 0; index0 < *nrings; ++index0 ) 01048 { 01049 npoints = ( int* )ptr; 01050 ptr += sizeof( int ); 01051 01052 for ( int index1 = 0; index1 < *npoints; ++index1 ) 01053 { 01054 ptr += sizeof( double ); 01055 ptr += sizeof( double ); 01056 if ( hasZValue ) //skip z-coordinate for 25D geometries 01057 { 01058 ptr += sizeof( double ); 01059 } 01060 if ( vertexcounter == atVertex ) 01061 { 01062 if ( index1 == 0 ) 01063 { 01064 beforeVertex = vertexcounter + ( *npoints - 2 ); 01065 afterVertex = vertexcounter + 1; 01066 } 01067 else if ( index1 == ( *npoints - 1 ) ) 01068 { 01069 beforeVertex = vertexcounter - 1; 01070 afterVertex = vertexcounter - ( *npoints - 2 ); 01071 } 01072 else 01073 { 01074 beforeVertex = vertexcounter - 1; 01075 afterVertex = vertexcounter + 1; 01076 } 01077 } 01078 01079 ++vertexcounter; 01080 } 01081 } 01082 break; 01083 } 01084 case QGis::WKBMultiPoint25D: 01085 case QGis::WKBMultiPoint: 01086 { 01087 // NOOP - Points do not have adjacent verticies 01088 break; 01089 } 01090 01091 case QGis::WKBMultiLineString25D: 01092 hasZValue = true; 01093 case QGis::WKBMultiLineString: 01094 { 01095 unsigned char* ptr = mGeometry + 5; 01096 int* nlines = ( int* )ptr; 01097 int* npoints = 0; 01098 ptr += sizeof( int ); 01099 01100 for ( int index0 = 0; index0 < *nlines; ++index0 ) 01101 { 01102 ptr += ( sizeof( int ) + 1 ); 01103 npoints = ( int* )ptr; 01104 ptr += sizeof( int ); 01105 01106 for ( int index1 = 0; index1 < *npoints; ++index1 ) 01107 { 01108 ptr += sizeof( double ); 01109 ptr += sizeof( double ); 01110 if ( hasZValue ) //skip z-coordinate for 25D geometries 01111 { 01112 ptr += sizeof( double ); 01113 } 01114 01115 if ( vertexcounter == atVertex ) 01116 { 01117 // Found the vertex of the linestring we were looking for. 01118 if ( index1 == 0 ) 01119 { 01120 beforeVertex = -1; 01121 } 01122 else 01123 { 01124 beforeVertex = vertexcounter - 1; 01125 } 01126 if ( index1 == ( *npoints ) - 1 ) 01127 { 01128 afterVertex = -1; 01129 } 01130 else 01131 { 01132 afterVertex = vertexcounter + 1; 01133 } 01134 } 01135 ++vertexcounter; 01136 } 01137 } 01138 break; 01139 } 01140 case QGis::WKBMultiPolygon25D: 01141 hasZValue = true; 01142 case QGis::WKBMultiPolygon: 01143 { 01144 unsigned char* ptr = mGeometry + 5; 01145 int* npolys = ( int* )ptr; 01146 int* nrings; 01147 int* npoints; 01148 ptr += sizeof( int ); 01149 01150 for ( int index0 = 0; index0 < *npolys; ++index0 ) 01151 { 01152 ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type 01153 nrings = ( int* )ptr; 01154 ptr += sizeof( int ); 01155 01156 for ( int index1 = 0; index1 < *nrings; ++index1 ) 01157 { 01158 npoints = ( int* )ptr; 01159 ptr += sizeof( int ); 01160 01161 for ( int index2 = 0; index2 < *npoints; ++index2 ) 01162 { 01163 ptr += sizeof( double ); 01164 ptr += sizeof( double ); 01165 if ( hasZValue ) //skip z-coordinate for 25D geometries 01166 { 01167 ptr += sizeof( double ); 01168 } 01169 if ( vertexcounter == atVertex ) 01170 { 01171 // Found the vertex of the linear-ring of the polygon we were looking for. 01172 // assign the rubber band indices 01173 01174 if ( index2 == 0 ) 01175 { 01176 beforeVertex = vertexcounter + ( *npoints - 2 ); 01177 afterVertex = vertexcounter + 1; 01178 } 01179 else if ( index2 == ( *npoints - 1 ) ) 01180 { 01181 beforeVertex = vertexcounter - 1; 01182 afterVertex = vertexcounter - ( *npoints - 2 ); 01183 } 01184 else 01185 { 01186 beforeVertex = vertexcounter - 1; 01187 afterVertex = vertexcounter + 1; 01188 } 01189 } 01190 ++vertexcounter; 01191 } 01192 } 01193 } 01194 01195 break; 01196 } 01197 01198 default: 01199 break; 01200 } // switch (wkbType) 01201 } 01202 01203 01204 01205 bool QgsGeometry::insertVertex( double x, double y, 01206 int beforeVertex, 01207 const GEOSCoordSequence* old_sequence, 01208 GEOSCoordSequence** new_sequence ) 01209 01210 { 01211 // Bounds checking 01212 if ( beforeVertex < 0 ) 01213 { 01214 *new_sequence = 0; 01215 return false; 01216 } 01217 01218 unsigned int numPoints; 01219 GEOSCoordSeq_getSize( old_sequence, &numPoints ); 01220 01221 *new_sequence = GEOSCoordSeq_create( numPoints + 1, 2 ); 01222 if ( !*new_sequence ) 01223 return false; 01224 01225 bool inserted = false; 01226 for ( unsigned int i = 0, j = 0; i < numPoints; i++, j++ ) 01227 { 01228 // Do we insert the new vertex here? 01229 if ( beforeVertex == static_cast<int>( i ) ) 01230 { 01231 GEOSCoordSeq_setX( *new_sequence, j, x ); 01232 GEOSCoordSeq_setY( *new_sequence, j, y ); 01233 j++; 01234 inserted = true; 01235 } 01236 01237 double aX, aY; 01238 GEOSCoordSeq_getX( old_sequence, i, &aX ); 01239 GEOSCoordSeq_getY( old_sequence, i, &aY ); 01240 01241 GEOSCoordSeq_setX( *new_sequence, j, aX ); 01242 GEOSCoordSeq_setY( *new_sequence, j, aY ); 01243 } 01244 01245 if ( !inserted ) 01246 { 01247 // The beforeVertex is greater than the actual number of vertices 01248 // in the geometry - append it. 01249 GEOSCoordSeq_setX( *new_sequence, numPoints, x ); 01250 GEOSCoordSeq_setY( *new_sequence, numPoints, y ); 01251 } 01252 // TODO: Check that the sequence is still simple, e.g. with GEOS_GEOM::Geometry->isSimple() 01253 01254 return inserted; 01255 } 01256 01257 bool QgsGeometry::moveVertex( double x, double y, int atVertex ) 01258 { 01259 int vertexnr = atVertex; 01260 01261 // TODO: implement with GEOS 01262 if ( mDirtyWkb ) 01263 { 01264 exportGeosToWkb(); 01265 } 01266 01267 if ( !mGeometry ) 01268 { 01269 QgsDebugMsg( "WKB geometry not available!" ); 01270 return false; 01271 } 01272 01273 QGis::WkbType wkbType; 01274 bool hasZValue = false; 01275 unsigned char* ptr = mGeometry + 1; 01276 memcpy( &wkbType, ptr, sizeof( wkbType ) ); 01277 ptr += sizeof( wkbType ); 01278 01279 switch ( wkbType ) 01280 { 01281 case QGis::WKBPoint25D: 01282 case QGis::WKBPoint: 01283 { 01284 if ( vertexnr == 0 ) 01285 { 01286 memcpy( ptr, &x, sizeof( double ) ); 01287 ptr += sizeof( double ); 01288 memcpy( ptr, &y, sizeof( double ) ); 01289 mDirtyGeos = true; 01290 return true; 01291 } 01292 else 01293 { 01294 return false; 01295 } 01296 } 01297 case QGis::WKBMultiPoint25D: 01298 hasZValue = true; 01299 case QGis::WKBMultiPoint: 01300 { 01301 int* nrPoints = ( int* )ptr; 01302 if ( vertexnr > *nrPoints || vertexnr < 0 ) 01303 { 01304 return false; 01305 } 01306 ptr += sizeof( int ); 01307 if ( hasZValue ) 01308 { 01309 ptr += ( 3 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr; 01310 } 01311 else 01312 { 01313 ptr += ( 2 * sizeof( double ) + 1 + sizeof( int ) ) * vertexnr; 01314 } 01315 ptr += ( 1 + sizeof( int ) ); 01316 memcpy( ptr, &x, sizeof( double ) ); 01317 ptr += sizeof( double ); 01318 memcpy( ptr, &y, sizeof( double ) ); 01319 mDirtyGeos = true; 01320 return true; 01321 } 01322 case QGis::WKBLineString25D: 01323 hasZValue = true; 01324 case QGis::WKBLineString: 01325 { 01326 int* nrPoints = ( int* )ptr; 01327 if ( vertexnr > *nrPoints || vertexnr < 0 ) 01328 { 01329 return false; 01330 } 01331 ptr += sizeof( int ); 01332 ptr += 2 * sizeof( double ) * vertexnr; 01333 if ( hasZValue ) 01334 { 01335 ptr += sizeof( double ) * vertexnr; 01336 } 01337 memcpy( ptr, &x, sizeof( double ) ); 01338 ptr += sizeof( double ); 01339 memcpy( ptr, &y, sizeof( double ) ); 01340 mDirtyGeos = true; 01341 return true; 01342 } 01343 case QGis::WKBMultiLineString25D: 01344 hasZValue = true; 01345 case QGis::WKBMultiLineString: 01346 { 01347 int* nrLines = ( int* )ptr; 01348 ptr += sizeof( int ); 01349 int* nrPoints = 0; //numer of points in a line 01350 int pointindex = 0; 01351 for ( int linenr = 0; linenr < *nrLines; ++linenr ) 01352 { 01353 ptr += sizeof( int ) + 1; 01354 nrPoints = ( int* )ptr; 01355 ptr += sizeof( int ); 01356 if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nrPoints ) ) 01357 { 01358 ptr += ( vertexnr - pointindex ) * 2 * sizeof( double ); 01359 if ( hasZValue ) 01360 { 01361 ptr += ( vertexnr - pointindex ) * sizeof( double ); 01362 } 01363 memcpy( ptr, &x, sizeof( double ) ); 01364 memcpy( ptr + sizeof( double ), &y, sizeof( double ) ); 01365 mDirtyGeos = true; 01366 return true; 01367 } 01368 pointindex += ( *nrPoints ); 01369 ptr += 2 * sizeof( double ) * ( *nrPoints ); 01370 if ( hasZValue ) 01371 { 01372 ptr += sizeof( double ) * ( *nrPoints ); 01373 } 01374 } 01375 return false; 01376 } 01377 case QGis::WKBPolygon25D: 01378 hasZValue = true; 01379 case QGis::WKBPolygon: 01380 { 01381 int* nrRings = ( int* )ptr; 01382 ptr += sizeof( int ); 01383 int* nrPoints = 0; //numer of points in a ring 01384 int pointindex = 0; 01385 01386 for ( int ringnr = 0; ringnr < *nrRings; ++ringnr ) 01387 { 01388 nrPoints = ( int* )ptr; 01389 ptr += sizeof( int ); 01390 if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )//move two points 01391 { 01392 memcpy( ptr, &x, sizeof( double ) ); 01393 memcpy( ptr + sizeof( double ), &y, sizeof( double ) ); 01394 if ( hasZValue ) 01395 { 01396 memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) ); 01397 } 01398 else 01399 { 01400 memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) ); 01401 } 01402 if ( hasZValue ) 01403 { 01404 memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) ); 01405 } 01406 else 01407 { 01408 memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) ); 01409 } 01410 mDirtyGeos = true; 01411 return true; 01412 } 01413 else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )//move only one point 01414 { 01415 ptr += 2 * sizeof( double ) * ( vertexnr - pointindex ); 01416 if ( hasZValue ) 01417 { 01418 ptr += sizeof( double ) * ( vertexnr - pointindex ); 01419 } 01420 memcpy( ptr, &x, sizeof( double ) ); 01421 ptr += sizeof( double ); 01422 memcpy( ptr, &y, sizeof( double ) ); 01423 mDirtyGeos = true; 01424 return true; 01425 } 01426 ptr += 2 * sizeof( double ) * ( *nrPoints ); 01427 if ( hasZValue ) 01428 { 01429 ptr += sizeof( double ) * ( *nrPoints ); 01430 } 01431 pointindex += *nrPoints; 01432 } 01433 return false; 01434 } 01435 case QGis::WKBMultiPolygon25D: 01436 hasZValue = true; 01437 case QGis::WKBMultiPolygon: 01438 { 01439 int* nrPolygons = ( int* )ptr; 01440 ptr += sizeof( int ); 01441 int* nrRings = 0; //number of rings in a polygon 01442 int* nrPoints = 0; //number of points in a ring 01443 int pointindex = 0; 01444 01445 for ( int polynr = 0; polynr < *nrPolygons; ++polynr ) 01446 { 01447 ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type 01448 nrRings = ( int* )ptr; 01449 ptr += sizeof( int ); 01450 for ( int ringnr = 0; ringnr < *nrRings; ++ringnr ) 01451 { 01452 nrPoints = ( int* )ptr; 01453 ptr += sizeof( int ); 01454 if ( vertexnr == pointindex || vertexnr == pointindex + ( *nrPoints - 1 ) )//move two points 01455 { 01456 memcpy( ptr, &x, sizeof( double ) ); 01457 memcpy( ptr + sizeof( double ), &y, sizeof( double ) ); 01458 if ( hasZValue ) 01459 { 01460 memcpy( ptr + 3*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) ); 01461 } 01462 else 01463 { 01464 memcpy( ptr + 2*sizeof( double )*( *nrPoints - 1 ), &x, sizeof( double ) ); 01465 } 01466 if ( hasZValue ) 01467 { 01468 memcpy( ptr + sizeof( double ) + 3*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) ); 01469 } 01470 else 01471 { 01472 memcpy( ptr + sizeof( double ) + 2*sizeof( double )*( *nrPoints - 1 ), &y, sizeof( double ) ); 01473 } 01474 mDirtyGeos = true; 01475 return true; 01476 } 01477 else if ( vertexnr > pointindex && vertexnr < pointindex + ( *nrPoints - 1 ) )//move only one point 01478 { 01479 ptr += 2 * sizeof( double ) * ( vertexnr - pointindex ); 01480 if ( hasZValue ) 01481 { 01482 ptr += sizeof( double ) * ( vertexnr - pointindex ); 01483 } 01484 memcpy( ptr, &x, sizeof( double ) ); 01485 ptr += sizeof( double ); 01486 memcpy( ptr, &y, sizeof( double ) ); 01487 mDirtyGeos = true; 01488 return true; 01489 } 01490 ptr += 2 * sizeof( double ) * ( *nrPoints ); 01491 if ( hasZValue ) 01492 { 01493 ptr += sizeof( double ) * ( *nrPoints ); 01494 } 01495 pointindex += *nrPoints; 01496 } 01497 } 01498 return false; 01499 } 01500 01501 default: 01502 return false; 01503 } 01504 } 01505 01506 bool QgsGeometry::deleteVertex( int atVertex ) 01507 { 01508 int vertexnr = atVertex; 01509 bool success = false; 01510 01511 // TODO: implement with GEOS 01512 if ( mDirtyWkb ) 01513 { 01514 exportGeosToWkb(); 01515 } 01516 01517 if ( !mGeometry ) 01518 { 01519 QgsDebugMsg( "WKB geometry not available!" ); 01520 return false; 01521 } 01522 01523 //create a new geometry buffer for the modified geometry 01524 unsigned char* newbuffer; 01525 int pointindex = 0; 01526 QGis::WkbType wkbType; 01527 bool hasZValue = false; 01528 unsigned char* ptr = mGeometry + 1; 01529 memcpy( &wkbType, ptr, sizeof( wkbType ) ); 01530 01531 switch ( wkbType ) 01532 { 01533 case QGis::WKBPoint25D: 01534 case QGis::WKBLineString25D: 01535 case QGis::WKBPolygon25D: 01536 case QGis::WKBMultiPoint25D: 01537 case QGis::WKBMultiLineString25D: 01538 case QGis::WKBMultiPolygon25D: 01539 newbuffer = new unsigned char[mGeometrySize-3*sizeof( double )]; 01540 break; 01541 default: 01542 newbuffer = new unsigned char[mGeometrySize-2*sizeof( double )]; 01543 } 01544 01545 memcpy( newbuffer, mGeometry, 1 + sizeof( int ) ); //endian and type are the same 01546 01547 ptr += sizeof( wkbType ); 01548 unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int ); 01549 01550 switch ( wkbType ) 01551 { 01552 case QGis::WKBPoint25D: 01553 case QGis::WKBPoint: 01554 { 01555 break; //cannot remove the only point vertex 01556 } 01557 case QGis::WKBMultiPoint25D: 01558 hasZValue = true; 01559 case QGis::WKBMultiPoint: 01560 { 01561 //todo 01562 break; 01563 } 01564 case QGis::WKBLineString25D: 01565 hasZValue = true; 01566 case QGis::WKBLineString: 01567 { 01568 int* nPoints = ( int* )ptr; 01569 if (( *nPoints ) < 3 || vertexnr > ( *nPoints ) - 1 || vertexnr < 0 ) //line needs at least 2 vertices 01570 { 01571 delete newbuffer; 01572 return false; 01573 } 01574 int newNPoints = ( *nPoints ) - 1; //new number of points 01575 memcpy( newBufferPtr, &newNPoints, sizeof( int ) ); 01576 ptr += sizeof( int ); 01577 newBufferPtr += sizeof( int ); 01578 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01579 { 01580 if ( vertexnr != pointindex ) 01581 { 01582 memcpy( newBufferPtr, ptr, sizeof( double ) ); 01583 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) ); 01584 newBufferPtr += 2 * sizeof( double ); 01585 if ( hasZValue ) 01586 { 01587 newBufferPtr += sizeof( double ); 01588 } 01589 } 01590 else 01591 { 01592 success = true; 01593 } 01594 ptr += 2 * sizeof( double ); 01595 if ( hasZValue ) 01596 { 01597 ptr += sizeof( double ); 01598 } 01599 ++pointindex; 01600 } 01601 break; 01602 } 01603 case QGis::WKBMultiLineString25D: 01604 hasZValue = true; 01605 case QGis::WKBMultiLineString: 01606 { 01607 int* nLines = ( int* )ptr; 01608 memcpy( newBufferPtr, nLines, sizeof( int ) ); 01609 newBufferPtr += sizeof( int ); 01610 ptr += sizeof( int ); 01611 int* nPoints = 0; //number of points in a line 01612 int pointindex = 0; 01613 for ( int linenr = 0; linenr < *nLines; ++linenr ) 01614 { 01615 memcpy( newBufferPtr, ptr, sizeof( int ) + 1 ); 01616 ptr += ( sizeof( int ) + 1 ); 01617 newBufferPtr += ( sizeof( int ) + 1 ); 01618 nPoints = ( int* )ptr; 01619 ptr += sizeof( int ); 01620 int newNPoint; 01621 01622 //find out if the vertex is in this line 01623 if ( vertexnr >= pointindex && vertexnr < pointindex + ( *nPoints ) ) 01624 { 01625 if ( *nPoints < 3 ) //line needs at least 2 vertices 01626 { 01627 delete newbuffer; 01628 return false; 01629 } 01630 newNPoint = ( *nPoints ) - 1; 01631 } 01632 else 01633 { 01634 newNPoint = *nPoints; 01635 } 01636 memcpy( newBufferPtr, &newNPoint, sizeof( int ) ); 01637 newBufferPtr += sizeof( int ); 01638 01639 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01640 { 01641 if ( vertexnr != pointindex ) 01642 { 01643 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 01644 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 01645 newBufferPtr += 2 * sizeof( double ); 01646 if ( hasZValue ) 01647 { 01648 newBufferPtr += sizeof( double ); 01649 } 01650 } 01651 else 01652 { 01653 success = true; 01654 } 01655 ptr += 2 * sizeof( double ); 01656 if ( hasZValue ) 01657 { 01658 ptr += sizeof( double ); 01659 } 01660 ++pointindex; 01661 } 01662 } 01663 break; 01664 } 01665 case QGis::WKBPolygon25D: 01666 hasZValue = true; 01667 case QGis::WKBPolygon: 01668 { 01669 int* nRings = ( int* )ptr; 01670 memcpy( newBufferPtr, nRings, sizeof( int ) ); 01671 ptr += sizeof( int ); 01672 newBufferPtr += sizeof( int ); 01673 int* nPoints = 0; //number of points in a ring 01674 int pointindex = 0; 01675 01676 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 01677 { 01678 nPoints = ( int* )ptr; 01679 ptr += sizeof( int ); 01680 int newNPoints; 01681 if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )//vertex to delete is in this ring 01682 { 01683 if ( *nPoints < 5 ) //a ring has at least 3 points 01684 { 01685 delete newbuffer; 01686 return false; 01687 } 01688 newNPoints = *nPoints - 1; 01689 } 01690 else 01691 { 01692 newNPoints = *nPoints; 01693 } 01694 memcpy( newBufferPtr, &newNPoints, sizeof( int ) ); 01695 newBufferPtr += sizeof( int ); 01696 01697 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01698 { 01699 if ( vertexnr != pointindex ) 01700 { 01701 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 01702 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 01703 newBufferPtr += 2 * sizeof( double ); 01704 if ( hasZValue ) 01705 { 01706 newBufferPtr += sizeof( double ); 01707 } 01708 } 01709 else 01710 { 01711 if ( pointnr == 0 )//adjust the last point of the ring 01712 { 01713 memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) ); 01714 memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) ); 01715 } 01716 if ( pointnr == ( *nPoints ) - 1 ) 01717 { 01718 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) ); 01719 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) ); 01720 } 01721 success = true; 01722 } 01723 ptr += 2 * sizeof( double ); 01724 if ( hasZValue ) 01725 { 01726 ptr += sizeof( double ); 01727 } 01728 ++pointindex; 01729 } 01730 } 01731 break; 01732 } 01733 case QGis::WKBMultiPolygon25D: 01734 hasZValue = true; 01735 case QGis::WKBMultiPolygon: 01736 { 01737 int* nPolys = ( int* )ptr; 01738 memcpy( newBufferPtr, nPolys, sizeof( int ) ); 01739 newBufferPtr += sizeof( int ); 01740 ptr += sizeof( int ); 01741 int* nRings = 0; 01742 int* nPoints = 0; 01743 int pointindex = 0; 01744 01745 for ( int polynr = 0; polynr < *nPolys; ++polynr ) 01746 { 01747 memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) ); 01748 ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type 01749 newBufferPtr += ( 1 + sizeof( int ) ); 01750 nRings = ( int* )ptr; 01751 memcpy( newBufferPtr, nRings, sizeof( int ) ); 01752 newBufferPtr += sizeof( int ); 01753 ptr += sizeof( int ); 01754 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 01755 { 01756 nPoints = ( int* )ptr; 01757 ptr += sizeof( int ); 01758 int newNPoints; 01759 if ( vertexnr >= pointindex && vertexnr < pointindex + *nPoints )//vertex to delete is in this ring 01760 { 01761 if ( *nPoints < 5 ) //a ring has at least 3 points 01762 { 01763 delete newbuffer; 01764 return false; 01765 } 01766 newNPoints = *nPoints - 1; 01767 } 01768 else 01769 { 01770 newNPoints = *nPoints; 01771 } 01772 memcpy( newBufferPtr, &newNPoints, sizeof( int ) ); 01773 newBufferPtr += sizeof( int ); 01774 01775 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01776 { 01777 if ( vertexnr != pointindex ) 01778 { 01779 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 01780 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 01781 newBufferPtr += 2 * sizeof( double ); 01782 if ( hasZValue ) 01783 { 01784 newBufferPtr += sizeof( double ); 01785 } 01786 } 01787 else 01788 { 01789 if ( pointnr == 0 )//adjust the last point of the ring 01790 { 01791 memcpy( ptr + ( *nPoints - 1 )*2*sizeof( double ), ptr + 2*sizeof( double ), sizeof( double ) ); 01792 memcpy( ptr + sizeof( double ) + ( *nPoints - 1 )*2*sizeof( double ), ptr + 3*sizeof( double ), sizeof( double ) ); 01793 } 01794 if ( pointnr == ( *nPoints ) - 1 ) 01795 { 01796 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ), ptr - 2*sizeof( double ), sizeof( double ) ); 01797 memcpy( newBufferPtr - ( *nPoints - 1 )*2*sizeof( double ) + sizeof( double ), ptr - sizeof( double ), sizeof( double ) ); 01798 } 01799 success = true; 01800 } 01801 ptr += 2 * sizeof( double ); 01802 if ( hasZValue ) 01803 { 01804 ptr += sizeof( double ); 01805 } 01806 ++pointindex; 01807 } 01808 } 01809 } 01810 break; 01811 } 01812 case QGis::WKBNoGeometry: 01813 case QGis::WKBUnknown: 01814 break; 01815 } 01816 if ( success ) 01817 { 01818 delete [] mGeometry; 01819 mGeometry = newbuffer; 01820 mGeometrySize -= ( 2 * sizeof( double ) ); 01821 if ( hasZValue ) 01822 { 01823 mGeometrySize -= sizeof( double ); 01824 } 01825 mDirtyGeos = true; 01826 return true; 01827 } 01828 else 01829 { 01830 delete[] newbuffer; 01831 return false; 01832 } 01833 } 01834 01835 bool QgsGeometry::insertVertex( double x, double y, int beforeVertex ) 01836 { 01837 int vertexnr = beforeVertex; 01838 bool success = false; 01839 01840 // TODO: implement with GEOS 01841 if ( mDirtyWkb ) 01842 { 01843 exportGeosToWkb(); 01844 } 01845 01846 if ( !mGeometry ) 01847 { 01848 QgsDebugMsg( "WKB geometry not available!" ); 01849 return false; 01850 } 01851 01852 //create a new geometry buffer for the modified geometry 01853 unsigned char* newbuffer; 01854 01855 int pointindex = 0; 01856 QGis::WkbType wkbType; 01857 bool hasZValue = false; 01858 01859 unsigned char* ptr = mGeometry + 1; 01860 memcpy( &wkbType, ptr, sizeof( wkbType ) ); 01861 01862 switch ( wkbType ) 01863 { 01864 case QGis::WKBPoint25D: 01865 case QGis::WKBLineString25D: 01866 case QGis::WKBPolygon25D: 01867 case QGis::WKBMultiPoint25D: 01868 case QGis::WKBMultiLineString25D: 01869 case QGis::WKBMultiPolygon25D: 01870 newbuffer = new unsigned char[mGeometrySize+3*sizeof( double )]; 01871 break; 01872 default: 01873 newbuffer = new unsigned char[mGeometrySize+2*sizeof( double )]; 01874 } 01875 memcpy( newbuffer, mGeometry, 1 + sizeof( int ) ); //endian and type are the same 01876 01877 ptr += sizeof( wkbType ); 01878 unsigned char* newBufferPtr = newbuffer + 1 + sizeof( int ); 01879 01880 switch ( wkbType ) 01881 { 01882 case QGis::WKBPoint25D: 01883 case QGis::WKBPoint://cannot insert a vertex before another one on point types 01884 { 01885 delete newbuffer; 01886 return false; 01887 } 01888 case QGis::WKBMultiPoint25D: 01889 hasZValue = true; 01890 case QGis::WKBMultiPoint: 01891 { 01892 //todo 01893 break; 01894 } 01895 case QGis::WKBLineString25D: 01896 hasZValue = true; 01897 case QGis::WKBLineString: 01898 { 01899 int* nPoints = ( int* )ptr; 01900 if ( vertexnr > *nPoints || vertexnr < 0 ) 01901 { 01902 break; 01903 } 01904 int newNPoints = ( *nPoints ) + 1; 01905 memcpy( newBufferPtr, &newNPoints, sizeof( int ) ); 01906 newBufferPtr += sizeof( int ); 01907 ptr += sizeof( int ); 01908 01909 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01910 { 01911 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 01912 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//x 01913 ptr += 2 * sizeof( double ); 01914 if ( hasZValue ) 01915 { 01916 ptr += sizeof( double ); 01917 } 01918 newBufferPtr += 2 * sizeof( double ); 01919 if ( hasZValue ) 01920 { 01921 newBufferPtr += sizeof( double ); 01922 } 01923 ++pointindex; 01924 if ( pointindex == vertexnr ) 01925 { 01926 memcpy( newBufferPtr, &x, sizeof( double ) ); 01927 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) ); 01928 newBufferPtr += 2 * sizeof( double ); 01929 if ( hasZValue ) 01930 { 01931 newBufferPtr += sizeof( double ); 01932 } 01933 success = true; 01934 } 01935 } 01936 break; 01937 } 01938 case QGis::WKBMultiLineString25D: 01939 hasZValue = true; 01940 case QGis::WKBMultiLineString: 01941 { 01942 int* nLines = ( int* )ptr; 01943 int* nPoints = 0; //number of points in a line 01944 ptr += sizeof( int ); 01945 memcpy( newBufferPtr, nLines, sizeof( int ) ); 01946 newBufferPtr += sizeof( int ); 01947 int pointindex = 0; 01948 01949 for ( int linenr = 0; linenr < *nLines; ++linenr ) 01950 { 01951 memcpy( newBufferPtr, ptr, sizeof( int ) + 1 ); 01952 ptr += ( sizeof( int ) + 1 ); 01953 newBufferPtr += ( sizeof( int ) + 1 ); 01954 nPoints = ( int* )ptr; 01955 int newNPoints; 01956 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring 01957 { 01958 newNPoints = ( *nPoints ) + 1; 01959 } 01960 else 01961 { 01962 newNPoints = *nPoints; 01963 } 01964 memcpy( newBufferPtr, &newNPoints, sizeof( double ) ); 01965 newBufferPtr += sizeof( int ); 01966 ptr += sizeof( int ); 01967 01968 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 01969 { 01970 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 01971 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 01972 ptr += 2 * sizeof( double ); 01973 newBufferPtr += 2 * sizeof( double ); 01974 if ( hasZValue ) 01975 { 01976 ptr += sizeof( double ); 01977 newBufferPtr += sizeof( double ); 01978 } 01979 ++pointindex; 01980 if ( pointindex == vertexnr ) 01981 { 01982 memcpy( newBufferPtr, &x, sizeof( double ) ); 01983 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) ); 01984 newBufferPtr += 2 * sizeof( double ); 01985 if ( hasZValue ) 01986 { 01987 newBufferPtr += sizeof( double ); 01988 } 01989 success = true; 01990 } 01991 } 01992 } 01993 break; 01994 } 01995 case QGis::WKBPolygon25D: 01996 hasZValue = true; 01997 case QGis::WKBPolygon: 01998 { 01999 int* nRings = ( int* )ptr; 02000 int* nPoints = 0; //number of points in a ring 02001 ptr += sizeof( int ); 02002 memcpy( newBufferPtr, nRings, sizeof( int ) ); 02003 newBufferPtr += sizeof( int ); 02004 int pointindex = 0; 02005 02006 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 02007 { 02008 nPoints = ( int* )ptr; 02009 int newNPoints; 02010 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring 02011 { 02012 newNPoints = ( *nPoints ) + 1; 02013 } 02014 else 02015 { 02016 newNPoints = *nPoints; 02017 } 02018 memcpy( newBufferPtr, &newNPoints, sizeof( double ) ); 02019 newBufferPtr += sizeof( int ); 02020 ptr += sizeof( int ); 02021 02022 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02023 { 02024 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 02025 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 02026 ptr += 2 * sizeof( double ); 02027 newBufferPtr += 2 * sizeof( double ); 02028 if ( hasZValue ) 02029 { 02030 ptr += sizeof( double ); 02031 newBufferPtr += sizeof( double ); 02032 } 02033 ++pointindex; 02034 if ( pointindex == vertexnr ) 02035 { 02036 memcpy( newBufferPtr, &x, sizeof( double ) ); 02037 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) ); 02038 newBufferPtr += 2 * sizeof( double ); 02039 if ( hasZValue ) 02040 { 02041 newBufferPtr += sizeof( double ); 02042 } 02043 success = true; 02044 } 02045 } 02046 } 02047 break; 02048 } 02049 case QGis::WKBMultiPolygon25D: 02050 hasZValue = true; 02051 case QGis::WKBMultiPolygon: 02052 { 02053 int* nPolys = ( int* )ptr; 02054 int* nRings = 0; //number of rings in a polygon 02055 int* nPoints = 0; //number of points in a ring 02056 memcpy( newBufferPtr, nPolys, sizeof( int ) ); 02057 ptr += sizeof( int ); 02058 newBufferPtr += sizeof( int ); 02059 int pointindex = 0; 02060 02061 for ( int polynr = 0; polynr < *nPolys; ++polynr ) 02062 { 02063 memcpy( newBufferPtr, ptr, ( 1 + sizeof( int ) ) ); 02064 ptr += ( 1 + sizeof( int ) );//skip endian and polygon type 02065 newBufferPtr += ( 1 + sizeof( int ) ); 02066 nRings = ( int* )ptr; 02067 ptr += sizeof( int ); 02068 memcpy( newBufferPtr, nRings, sizeof( int ) ); 02069 newBufferPtr += sizeof( int ); 02070 02071 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 02072 { 02073 nPoints = ( int* )ptr; 02074 int newNPoints; 02075 if ( vertexnr >= pointindex && vertexnr < ( pointindex + ( *nPoints ) ) )//point is in this ring 02076 { 02077 newNPoints = ( *nPoints ) + 1; 02078 } 02079 else 02080 { 02081 newNPoints = *nPoints; 02082 } 02083 memcpy( newBufferPtr, &newNPoints, sizeof( double ) ); 02084 newBufferPtr += sizeof( int ); 02085 ptr += sizeof( int ); 02086 02087 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02088 { 02089 memcpy( newBufferPtr, ptr, sizeof( double ) );//x 02090 memcpy( newBufferPtr + sizeof( double ), ptr + sizeof( double ), sizeof( double ) );//y 02091 ptr += 2 * sizeof( double ); 02092 newBufferPtr += 2 * sizeof( double ); 02093 if ( hasZValue ) 02094 { 02095 ptr += sizeof( double ); 02096 newBufferPtr += sizeof( double ); 02097 } 02098 ++pointindex; 02099 if ( pointindex == vertexnr ) 02100 { 02101 memcpy( newBufferPtr, &x, sizeof( double ) ); 02102 memcpy( newBufferPtr + sizeof( double ), &y, sizeof( double ) ); 02103 newBufferPtr += 2 * sizeof( double ); 02104 if ( hasZValue ) 02105 { 02106 newBufferPtr += sizeof( double ); 02107 } 02108 success = true; 02109 } 02110 } 02111 } 02112 02113 } 02114 break; 02115 } 02116 case QGis::WKBNoGeometry: 02117 case QGis::WKBUnknown: 02118 break; 02119 } 02120 02121 if ( success ) 02122 { 02123 delete [] mGeometry; 02124 mGeometry = newbuffer; 02125 mGeometrySize += 2 * sizeof( double ); 02126 if ( hasZValue ) 02127 { 02128 mGeometrySize += sizeof( double ); 02129 } 02130 mDirtyGeos = true; 02131 return true; 02132 } 02133 else 02134 { 02135 delete newbuffer; 02136 return false; 02137 } 02138 } 02139 02140 QgsPoint QgsGeometry::vertexAt( int atVertex ) 02141 { 02142 double x, y; 02143 02144 if ( mDirtyWkb ) 02145 { 02146 exportGeosToWkb(); 02147 } 02148 02149 if ( !mGeometry ) 02150 { 02151 QgsDebugMsg( "WKB geometry not available!" ); 02152 return QgsPoint( 0, 0 ); 02153 } 02154 02155 QGis::WkbType wkbType; 02156 bool hasZValue = false; 02157 unsigned char* ptr; 02158 02159 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) ); 02160 switch ( wkbType ) 02161 { 02162 case QGis::WKBPoint25D: 02163 case QGis::WKBPoint: 02164 { 02165 if ( atVertex == 0 ) 02166 { 02167 ptr = mGeometry + 1 + sizeof( int ); 02168 memcpy( &x, ptr, sizeof( double ) ); 02169 ptr += sizeof( double ); 02170 memcpy( &y, ptr, sizeof( double ) ); 02171 return QgsPoint( x, y ); 02172 } 02173 else 02174 { 02175 return QgsPoint( 0, 0 ); 02176 } 02177 } 02178 case QGis::WKBLineString25D: 02179 hasZValue = true; 02180 case QGis::WKBLineString: 02181 { 02182 int *nPoints; 02183 // get number of points in the line 02184 ptr = mGeometry + 1 + sizeof( int ); // now at mGeometry.numPoints 02185 nPoints = ( int * ) ptr; 02186 02187 // return error if underflow 02188 if ( 0 > atVertex || *nPoints <= atVertex ) 02189 { 02190 return QgsPoint( 0, 0 ); 02191 } 02192 02193 // copy the vertex coordinates 02194 if ( hasZValue ) 02195 { 02196 ptr = mGeometry + 9 + ( atVertex * 3 * sizeof( double ) ); 02197 } 02198 else 02199 { 02200 ptr = mGeometry + 9 + ( atVertex * 2 * sizeof( double ) ); 02201 } 02202 memcpy( &x, ptr, sizeof( double ) ); 02203 ptr += sizeof( double ); 02204 memcpy( &y, ptr, sizeof( double ) ); 02205 return QgsPoint( x, y ); 02206 } 02207 case QGis::WKBPolygon25D: 02208 hasZValue = true; 02209 case QGis::WKBPolygon: 02210 { 02211 int *nRings; 02212 int *nPoints = 0; 02213 ptr = mGeometry + 1 + sizeof( int ); 02214 nRings = ( int* )ptr; 02215 ptr += sizeof( int ); 02216 int pointindex = 0; 02217 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 02218 { 02219 nPoints = ( int* )ptr; 02220 ptr += sizeof( int ); 02221 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02222 { 02223 if ( pointindex == atVertex ) 02224 { 02225 memcpy( &x, ptr, sizeof( double ) ); 02226 ptr += sizeof( double ); 02227 memcpy( &y, ptr, sizeof( double ) ); 02228 return QgsPoint( x, y ); 02229 } 02230 ptr += 2 * sizeof( double ); 02231 if ( hasZValue ) 02232 { 02233 ptr += sizeof( double ); 02234 } 02235 ++pointindex; 02236 } 02237 } 02238 return QgsPoint( 0, 0 ); 02239 } 02240 case QGis::WKBMultiPoint25D: 02241 hasZValue = true; 02242 case QGis::WKBMultiPoint: 02243 { 02244 ptr = mGeometry + 1 + sizeof( int ); 02245 int* nPoints = ( int* )ptr; 02246 if ( atVertex < 0 || atVertex >= *nPoints ) 02247 { 02248 return QgsPoint( 0, 0 ); 02249 } 02250 if ( hasZValue ) 02251 { 02252 ptr += atVertex * ( 3 * sizeof( double ) + 1 + sizeof( int ) ); 02253 } 02254 else 02255 { 02256 ptr += atVertex * ( 2 * sizeof( double ) + 1 + sizeof( int ) ); 02257 } 02258 ptr += 1 + sizeof( int ); 02259 memcpy( &x, ptr, sizeof( double ) ); 02260 ptr += sizeof( double ); 02261 memcpy( &y, ptr, sizeof( double ) ); 02262 return QgsPoint( x, y ); 02263 } 02264 case QGis::WKBMultiLineString25D: 02265 hasZValue = true; 02266 case QGis::WKBMultiLineString: 02267 { 02268 ptr = mGeometry + 1 + sizeof( int ); 02269 int* nLines = ( int* )ptr; 02270 int* nPoints = 0; //number of points in a line 02271 int pointindex = 0; //global point counter 02272 ptr += sizeof( int ); 02273 for ( int linenr = 0; linenr < *nLines; ++linenr ) 02274 { 02275 ptr += sizeof( int ) + 1; 02276 nPoints = ( int* )ptr; 02277 ptr += sizeof( int ); 02278 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02279 { 02280 if ( pointindex == atVertex ) 02281 { 02282 memcpy( &x, ptr, sizeof( double ) ); 02283 ptr += sizeof( double ); 02284 memcpy( &y, ptr, sizeof( double ) ); 02285 return QgsPoint( x, y ); 02286 } 02287 ptr += 2 * sizeof( double ); 02288 if ( hasZValue ) 02289 { 02290 ptr += sizeof( double ); 02291 } 02292 ++pointindex; 02293 } 02294 } 02295 return QgsPoint( 0, 0 ); 02296 } 02297 case QGis::WKBMultiPolygon25D: 02298 hasZValue = true; 02299 case QGis::WKBMultiPolygon: 02300 { 02301 ptr = mGeometry + 1 + sizeof( int ); 02302 int* nRings = 0;//number of rings in a polygon 02303 int* nPoints = 0;//number of points in a ring 02304 int pointindex = 0; //global point counter 02305 int* nPolygons = ( int* )ptr; 02306 ptr += sizeof( int ); 02307 for ( int polynr = 0; polynr < *nPolygons; ++polynr ) 02308 { 02309 ptr += ( 1 + sizeof( int ) ); //skip endian and polygon type 02310 nRings = ( int* )ptr; 02311 ptr += sizeof( int ); 02312 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 02313 { 02314 nPoints = ( int* )ptr; 02315 ptr += sizeof( int ); 02316 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02317 { 02318 if ( pointindex == atVertex ) 02319 { 02320 memcpy( &x, ptr, sizeof( double ) ); 02321 ptr += sizeof( double ); 02322 memcpy( &y, ptr, sizeof( double ) ); 02323 return QgsPoint( x, y ); 02324 } 02325 ++pointindex; 02326 ptr += 2 * sizeof( double ); 02327 if ( hasZValue ) 02328 { 02329 ptr += sizeof( double ); 02330 } 02331 } 02332 } 02333 } 02334 return QgsPoint( 0, 0 ); 02335 } 02336 default: 02337 QgsDebugMsg( "error: mGeometry type not recognized" ); 02338 return QgsPoint( 0, 0 ); 02339 } 02340 } 02341 02342 02343 double QgsGeometry::sqrDistToVertexAt( QgsPoint& point, int atVertex ) 02344 { 02345 QgsPoint pnt = vertexAt( atVertex ); 02346 if ( pnt != QgsPoint( 0, 0 ) ) 02347 { 02348 QgsDebugMsg( "Exiting with distance to " + pnt.toString() ); 02349 return point.sqrDist( pnt ); 02350 } 02351 else 02352 { 02353 QgsDebugMsg( "Exiting with std::numeric_limits<double>::max()." ); 02354 // probably safest to bail out with a very large number 02355 return std::numeric_limits<double>::max(); 02356 } 02357 } 02358 02359 02360 double QgsGeometry::closestVertexWithContext( const QgsPoint& point, int& atVertex ) 02361 { 02362 double sqrDist = std::numeric_limits<double>::max(); 02363 02364 try 02365 { 02366 // Initialise some stuff 02367 int closestVertexIndex = 0; 02368 02369 // set up the GEOS geometry 02370 if ( !exportWkbToGeos() ) 02371 return -1; 02372 02373 const GEOSGeometry *g = GEOSGetExteriorRing( mGeos ); 02374 if ( !g ) 02375 return -1; 02376 02377 const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( g ); 02378 02379 unsigned int n; 02380 GEOSCoordSeq_getSize( sequence, &n ); 02381 02382 for ( unsigned int i = 0; i < n; i++ ) 02383 { 02384 double x, y; 02385 GEOSCoordSeq_getX( sequence, i, &x ); 02386 GEOSCoordSeq_getY( sequence, i, &y ); 02387 02388 double testDist = point.sqrDist( x, y ); 02389 if ( testDist < sqrDist ) 02390 { 02391 closestVertexIndex = i; 02392 sqrDist = testDist; 02393 } 02394 } 02395 02396 atVertex = closestVertexIndex; 02397 } 02398 catch ( GEOSException &e ) 02399 { 02400 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02401 return -1; 02402 } 02403 02404 return sqrDist; 02405 } 02406 02407 02408 double QgsGeometry::closestSegmentWithContext( 02409 const QgsPoint& point, 02410 QgsPoint& minDistPoint, 02411 int& afterVertex, 02412 double* leftOf, 02413 double epsilon ) 02414 { 02415 QgsDebugMsg( "Entering." ); 02416 QgsPoint distPoint; 02417 02418 QGis::WkbType wkbType; 02419 bool hasZValue = false; 02420 double *thisx = NULL; 02421 double *thisy = NULL; 02422 double *prevx = NULL; 02423 double *prevy = NULL; 02424 double testdist; 02425 int closestSegmentIndex = 0; 02426 02427 // Initialise some stuff 02428 double sqrDist = std::numeric_limits<double>::max(); 02429 02430 // TODO: implement with GEOS 02431 if ( mDirtyWkb ) //convert latest geos to mGeometry 02432 { 02433 exportGeosToWkb(); 02434 } 02435 02436 if ( !mGeometry ) 02437 { 02438 QgsDebugMsg( "WKB geometry not available!" ); 02439 return -1; 02440 } 02441 02442 memcpy( &wkbType, ( mGeometry + 1 ), sizeof( int ) ); 02443 02444 switch ( wkbType ) 02445 { 02446 case QGis::WKBPoint25D: 02447 case QGis::WKBPoint: 02448 case QGis::WKBMultiPoint25D: 02449 case QGis::WKBMultiPoint: 02450 { 02451 // Points have no lines 02452 return -1; 02453 } 02454 case QGis::WKBLineString25D: 02455 hasZValue = true; 02456 case QGis::WKBLineString: 02457 { 02458 unsigned char* ptr = mGeometry + 1 + sizeof( int ); 02459 int* npoints = ( int* ) ptr; 02460 ptr += sizeof( int ); 02461 for ( int index = 0; index < *npoints; ++index ) 02462 { 02463 if ( index > 0 ) 02464 { 02465 prevx = thisx; 02466 prevy = thisy; 02467 } 02468 thisx = ( double* ) ptr; 02469 ptr += sizeof( double ); 02470 thisy = ( double* ) ptr; 02471 02472 if ( index > 0 ) 02473 { 02474 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint, epsilon ) ) < sqrDist ) 02475 { 02476 closestSegmentIndex = index; 02477 sqrDist = testdist; 02478 minDistPoint = distPoint; 02479 if ( leftOf ) 02480 { 02481 *leftOf = QgsGeometry::leftOf( point.x(), point.y(), *prevx, *prevy, *thisx, *thisy ); 02482 } 02483 } 02484 } 02485 ptr += sizeof( double ); 02486 if ( hasZValue ) 02487 { 02488 ptr += sizeof( double ); 02489 } 02490 } 02491 afterVertex = closestSegmentIndex; 02492 break; 02493 } 02494 case QGis::WKBMultiLineString25D: 02495 hasZValue = true; 02496 case QGis::WKBMultiLineString: 02497 { 02498 unsigned char* ptr = mGeometry + 1 + sizeof( int ); 02499 int* nLines = ( int* )ptr; 02500 ptr += sizeof( int ); 02501 int* nPoints = 0; //number of points in a line 02502 int pointindex = 0;//global pointindex 02503 for ( int linenr = 0; linenr < *nLines; ++linenr ) 02504 { 02505 ptr += sizeof( int ) + 1; 02506 nPoints = ( int* )ptr; 02507 ptr += sizeof( int ); 02508 prevx = 0; 02509 prevy = 0; 02510 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02511 { 02512 thisx = ( double* ) ptr; 02513 ptr += sizeof( double ); 02514 thisy = ( double* ) ptr; 02515 ptr += sizeof( double ); 02516 if ( hasZValue ) 02517 { 02518 ptr += sizeof( double ); 02519 } 02520 if ( prevx && prevy ) 02521 { 02522 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint, epsilon ) ) < sqrDist ) 02523 { 02524 closestSegmentIndex = pointindex; 02525 sqrDist = testdist; 02526 minDistPoint = distPoint; 02527 if ( leftOf ) 02528 { 02529 *leftOf = QgsGeometry::leftOf( point.x(), point.y(), *prevx, *prevy, *thisx, *thisy ); 02530 } 02531 } 02532 } 02533 prevx = thisx; 02534 prevy = thisy; 02535 ++pointindex; 02536 } 02537 } 02538 afterVertex = closestSegmentIndex; 02539 break; 02540 } 02541 case QGis::WKBPolygon25D: 02542 hasZValue = true; 02543 case QGis::WKBPolygon: 02544 { 02545 int index = 0; 02546 unsigned char* ptr = mGeometry + 1 + sizeof( int ); 02547 int* nrings = ( int* )ptr; 02548 int* npoints = 0; //number of points in a ring 02549 ptr += sizeof( int ); 02550 for ( int ringnr = 0; ringnr < *nrings; ++ringnr )//loop over rings 02551 { 02552 npoints = ( int* )ptr; 02553 ptr += sizeof( int ); 02554 prevx = 0; 02555 prevy = 0; 02556 for ( int pointnr = 0; pointnr < *npoints; ++pointnr )//loop over points in a ring 02557 { 02558 thisx = ( double* )ptr; 02559 ptr += sizeof( double ); 02560 thisy = ( double* )ptr; 02561 ptr += sizeof( double ); 02562 if ( hasZValue ) 02563 { 02564 ptr += sizeof( double ); 02565 } 02566 if ( prevx && prevy ) 02567 { 02568 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint, epsilon ) ) < sqrDist ) 02569 { 02570 closestSegmentIndex = index; 02571 sqrDist = testdist; 02572 minDistPoint = distPoint; 02573 if ( leftOf ) 02574 { 02575 *leftOf = QgsGeometry::leftOf( point.x(), point.y(), *prevx, *prevy, *thisx, *thisy ); 02576 } 02577 } 02578 } 02579 prevx = thisx; 02580 prevy = thisy; 02581 ++index; 02582 } 02583 } 02584 afterVertex = closestSegmentIndex; 02585 break; 02586 } 02587 case QGis::WKBMultiPolygon25D: 02588 hasZValue = true; 02589 case QGis::WKBMultiPolygon: 02590 { 02591 unsigned char* ptr = mGeometry + 1 + sizeof( int ); 02592 int* nRings = 0; 02593 int* nPoints = 0; 02594 int pointindex = 0; 02595 int* nPolygons = ( int* )ptr; 02596 ptr += sizeof( int ); 02597 for ( int polynr = 0; polynr < *nPolygons; ++polynr ) 02598 { 02599 ptr += ( 1 + sizeof( int ) ); 02600 nRings = ( int* )ptr; 02601 ptr += sizeof( int ); 02602 for ( int ringnr = 0; ringnr < *nRings; ++ringnr ) 02603 { 02604 nPoints = ( int* )ptr; 02605 ptr += sizeof( int ); 02606 prevx = 0; 02607 prevy = 0; 02608 for ( int pointnr = 0; pointnr < *nPoints; ++pointnr ) 02609 { 02610 thisx = ( double* )ptr; 02611 ptr += sizeof( double ); 02612 thisy = ( double* )ptr; 02613 ptr += sizeof( double ); 02614 if ( hasZValue ) 02615 { 02616 ptr += sizeof( double ); 02617 } 02618 if ( prevx && prevy ) 02619 { 02620 if (( testdist = point.sqrDistToSegment( *prevx, *prevy, *thisx, *thisy, distPoint, epsilon ) ) < sqrDist ) 02621 { 02622 closestSegmentIndex = pointindex; 02623 sqrDist = testdist; 02624 minDistPoint = distPoint; 02625 if ( leftOf ) 02626 { 02627 *leftOf = QgsGeometry::leftOf( point.x(), point.y(), *prevx, *prevy, *thisx, *thisy ); 02628 } 02629 } 02630 } 02631 prevx = thisx; 02632 prevy = thisy; 02633 ++pointindex; 02634 } 02635 } 02636 } 02637 afterVertex = closestSegmentIndex; 02638 break; 02639 } 02640 case QGis::WKBUnknown: 02641 default: 02642 return -1; 02643 break; 02644 } // switch (wkbType) 02645 02646 02647 QgsDebugMsg( QString( "Exiting with nearest point %1, dist %2." ) 02648 .arg( point.toString() ).arg( sqrDist ) ); 02649 02650 return sqrDist; 02651 } 02652 02653 int QgsGeometry::addRing( const QList<QgsPoint>& ring ) 02654 { 02655 //bail out if this geometry is not polygon/multipolygon 02656 if ( type() != QGis::Polygon ) 02657 return 1; 02658 02659 //test for invalid geometries 02660 if ( ring.size() < 4 ) 02661 return 3; 02662 02663 //ring must be closed 02664 if ( ring.first() != ring.last() ) 02665 return 2; 02666 02667 //create geos geometry from wkb if not already there 02668 if ( !mGeos || mDirtyGeos ) 02669 if ( !exportWkbToGeos() ) 02670 return 6; 02671 02672 int type = GEOSGeomTypeId( mGeos ); 02673 02674 //Fill GEOS Polygons of the feature into list 02675 QVector<const GEOSGeometry*> polygonList; 02676 02677 if ( wkbType() == QGis::WKBPolygon ) 02678 { 02679 if ( type != GEOS_POLYGON ) 02680 return 1; 02681 02682 polygonList << mGeos; 02683 } 02684 else if ( wkbType() == QGis::WKBMultiPolygon ) 02685 { 02686 if ( type != GEOS_MULTIPOLYGON ) 02687 return 1; 02688 02689 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); ++i ) 02690 polygonList << GEOSGetGeometryN( mGeos, i ); 02691 } 02692 02693 //create new ring 02694 GEOSGeometry *newRing = 0; 02695 GEOSGeometry *newRingPolygon = 0; 02696 02697 try 02698 { 02699 newRing = createGeosLinearRing( ring.toVector() ); 02700 if ( !GEOSisValid( newRing ) ) 02701 { 02702 throwGEOSException( "ring is invalid" ); 02703 } 02704 02705 newRingPolygon = createGeosPolygon( newRing ); 02706 if ( !GEOSisValid( newRingPolygon ) ) 02707 { 02708 throwGEOSException( "ring is invalid" ); 02709 } 02710 } 02711 catch ( GEOSException &e ) 02712 { 02713 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02714 02715 if ( newRingPolygon ) 02716 GEOSGeom_destroy( newRingPolygon ); 02717 else if ( newRing ) 02718 GEOSGeom_destroy( newRing ); 02719 02720 return 3; 02721 } 02722 02723 QVector<GEOSGeometry*> rings; 02724 02725 int i; 02726 for ( i = 0; i < polygonList.size(); i++ ) 02727 { 02728 for ( int j = 0; j < rings.size(); j++ ) 02729 GEOSGeom_destroy( rings[j] ); 02730 rings.clear(); 02731 02732 GEOSGeometry *shellRing = 0; 02733 GEOSGeometry *shell = 0; 02734 try 02735 { 02736 shellRing = GEOSGeom_clone( GEOSGetExteriorRing( polygonList[i] ) ); 02737 shell = createGeosPolygon( shellRing ); 02738 02739 if ( !GEOSWithin( newRingPolygon, shell ) ) 02740 { 02741 GEOSGeom_destroy( shell ); 02742 continue; 02743 } 02744 } 02745 catch ( GEOSException &e ) 02746 { 02747 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02748 02749 if ( shell ) 02750 GEOSGeom_destroy( shell ); 02751 else if ( shellRing ) 02752 GEOSGeom_destroy( shellRing ); 02753 02754 GEOSGeom_destroy( newRingPolygon ); 02755 02756 return 4; 02757 } 02758 02759 // add outer ring 02760 rings << GEOSGeom_clone( shellRing ); 02761 02762 GEOSGeom_destroy( shell ); 02763 02764 // check inner rings 02765 int n = GEOSGetNumInteriorRings( polygonList[i] ); 02766 02767 int j; 02768 for ( j = 0; j < n; j++ ) 02769 { 02770 GEOSGeometry *holeRing = 0; 02771 GEOSGeometry *hole = 0; 02772 try 02773 { 02774 holeRing = GEOSGeom_clone( GEOSGetInteriorRingN( polygonList[i], j ) ); 02775 hole = createGeosPolygon( holeRing ); 02776 02777 if ( !GEOSDisjoint( hole, newRingPolygon ) ) 02778 { 02779 GEOSGeom_destroy( hole ); 02780 break; 02781 } 02782 } 02783 catch ( GEOSException &e ) 02784 { 02785 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02786 02787 if ( hole ) 02788 GEOSGeom_destroy( hole ); 02789 else if ( holeRing ) 02790 GEOSGeom_destroy( holeRing ); 02791 02792 break; 02793 } 02794 02795 rings << GEOSGeom_clone( holeRing ); 02796 GEOSGeom_destroy( hole ); 02797 } 02798 02799 if ( j == n ) 02800 // this is it... 02801 break; 02802 } 02803 02804 if ( i == polygonList.size() ) 02805 { 02806 // clear rings 02807 for ( int j = 0; j < rings.size(); j++ ) 02808 GEOSGeom_destroy( rings[j] ); 02809 rings.clear(); 02810 02811 GEOSGeom_destroy( newRingPolygon ); 02812 02813 // no containing polygon found 02814 return 5; 02815 } 02816 02817 rings << GEOSGeom_clone( newRing ); 02818 GEOSGeom_destroy( newRingPolygon ); 02819 02820 GEOSGeometry *newPolygon = createGeosPolygon( rings ); 02821 02822 if ( wkbType() == QGis::WKBPolygon ) 02823 { 02824 GEOSGeom_destroy( mGeos ); 02825 mGeos = newPolygon; 02826 } 02827 else if ( wkbType() == QGis::WKBMultiPolygon ) 02828 { 02829 QVector<GEOSGeometry*> newPolygons; 02830 02831 for ( int j = 0; j < polygonList.size(); j++ ) 02832 { 02833 newPolygons << ( i == j ? newPolygon : GEOSGeom_clone( polygonList[j] ) ); 02834 } 02835 02836 GEOSGeom_destroy( mGeos ); 02837 mGeos = createGeosCollection( GEOS_MULTIPOLYGON, newPolygons ); 02838 } 02839 02840 mDirtyWkb = true; 02841 mDirtyGeos = false; 02842 return 0; 02843 } 02844 02845 int QgsGeometry::addPart( const QList<QgsPoint> &points ) 02846 { 02847 QGis::GeometryType geomType = type(); 02848 02849 switch ( geomType ) 02850 { 02851 case QGis::Point: 02852 // only one part at a time 02853 if ( points.size() != 1 ) 02854 { 02855 QgsDebugMsg( "expected 1 point: " + QString::number( points.size() ) ); 02856 return 2; 02857 } 02858 break; 02859 02860 case QGis::Line: 02861 // Line needs to have at least two points and must be closed 02862 if ( points.size() < 3 ) 02863 { 02864 QgsDebugMsg( "line must at least have two points: " + QString::number( points.size() ) ); 02865 return 2; 02866 } 02867 break; 02868 02869 case QGis::Polygon: 02870 // Ring needs to have at least three points and must be closed 02871 if ( points.size() < 4 ) 02872 { 02873 QgsDebugMsg( "polygon must at least have three points: " + QString::number( points.size() ) ); 02874 return 2; 02875 } 02876 02877 // ring must be closed 02878 if ( points.first() != points.last() ) 02879 { 02880 QgsDebugMsg( "polygon not closed" ); 02881 return 2; 02882 } 02883 break; 02884 02885 default: 02886 QgsDebugMsg( "unsupported geometry type: " + QString::number( geomType ) ); 02887 return 2; 02888 } 02889 02890 if ( !isMultipart() && !convertToMultiType() ) 02891 { 02892 QgsDebugMsg( "could not convert to multipart" ); 02893 return 1; 02894 } 02895 02896 //create geos geometry from wkb if not already there 02897 if ( !mGeos || mDirtyGeos ) 02898 { 02899 if ( !exportWkbToGeos() ) 02900 { 02901 QgsDebugMsg( "could not export to GEOS geometry" ); 02902 return 4; 02903 } 02904 } 02905 02906 int geosType = GEOSGeomTypeId( mGeos ); 02907 GEOSGeometry *newPart = 0; 02908 02909 switch ( geomType ) 02910 { 02911 case QGis::Point: 02912 newPart = createGeosPoint( points[0] ); 02913 break; 02914 02915 case QGis::Line: 02916 newPart = createGeosLineString( points.toVector() ); 02917 break; 02918 02919 case QGis::Polygon: 02920 { 02921 //create new polygon from ring 02922 GEOSGeometry *newRing = 0; 02923 02924 try 02925 { 02926 newRing = createGeosLinearRing( points.toVector() ); 02927 if ( !GEOSisValid( newRing ) ) 02928 throw GEOSException( "ring invalid" ); 02929 02930 newPart = createGeosPolygon( newRing ); 02931 } 02932 catch ( GEOSException &e ) 02933 { 02934 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02935 02936 if ( newRing ) 02937 GEOSGeom_destroy( newRing ); 02938 02939 return 2; 02940 } 02941 } 02942 break; 02943 02944 default: 02945 QgsDebugMsg( "unsupported type: " + QString::number( type() ) ); 02946 return 2; 02947 } 02948 02949 Q_ASSERT( newPart ); 02950 02951 try 02952 { 02953 if ( !GEOSisValid( newPart ) ) 02954 throw GEOSException( "new part geometry invalid" ); 02955 } 02956 catch ( GEOSException &e ) 02957 { 02958 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 02959 02960 if ( newPart ) 02961 GEOSGeom_destroy( newPart ); 02962 02963 QgsDebugMsg( "part invalid: " + e.what() ); 02964 return 2; 02965 } 02966 02967 QVector<GEOSGeometry*> parts; 02968 02969 //create new multipolygon 02970 int n = GEOSGetNumGeometries( mGeos ); 02971 int i; 02972 for ( i = 0; i < n; ++i ) 02973 { 02974 const GEOSGeometry *partN = GEOSGetGeometryN( mGeos, i ); 02975 02976 if ( geomType == QGis::Polygon && !GEOSDisjoint( partN, newPart ) ) 02977 //bail out if new polygon is not disjoint with existing ones 02978 break; 02979 02980 parts << GEOSGeom_clone( partN ); 02981 } 02982 02983 if ( i < n ) 02984 { 02985 // bailed out 02986 for ( int i = 0; i < parts.size(); i++ ) 02987 GEOSGeom_destroy( parts[i] ); 02988 02989 QgsDebugMsg( "new polygon part not disjoint" ); 02990 return 3; 02991 } 02992 02993 parts << newPart; 02994 02995 GEOSGeom_destroy( mGeos ); 02996 02997 mGeos = createGeosCollection( geosType, parts ); 02998 02999 mDirtyWkb = true; 03000 mDirtyGeos = false; 03001 03002 return 0; 03003 } 03004 03005 int QgsGeometry::translate( double dx, double dy ) 03006 { 03007 if ( mDirtyWkb ) 03008 { 03009 exportGeosToWkb(); 03010 } 03011 03012 if ( !mGeometry ) 03013 { 03014 QgsDebugMsg( "WKB geometry not available!" ); 03015 return 1; 03016 } 03017 03018 QGis::WkbType wkbType; 03019 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) ); 03020 bool hasZValue = false; 03021 int wkbPosition = 5; 03022 03023 switch ( wkbType ) 03024 { 03025 case QGis::WKBPoint25D: 03026 case QGis::WKBPoint: 03027 { 03028 translateVertex( wkbPosition, dx, dy, hasZValue ); 03029 } 03030 break; 03031 03032 case QGis::WKBLineString25D: 03033 hasZValue = true; 03034 case QGis::WKBLineString: 03035 { 03036 int* npoints = ( int* )( &mGeometry[wkbPosition] ); 03037 wkbPosition += sizeof( int ); 03038 for ( int index = 0; index < *npoints; ++index ) 03039 { 03040 translateVertex( wkbPosition, dx, dy, hasZValue ); 03041 } 03042 break; 03043 } 03044 03045 case QGis::WKBPolygon25D: 03046 hasZValue = true; 03047 case QGis::WKBPolygon: 03048 { 03049 int* nrings = ( int* )( &( mGeometry[wkbPosition] ) ); 03050 wkbPosition += sizeof( int ); 03051 int* npoints; 03052 03053 for ( int index = 0; index < *nrings; ++index ) 03054 { 03055 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03056 wkbPosition += sizeof( int ); 03057 for ( int index2 = 0; index2 < *npoints; ++index2 ) 03058 { 03059 translateVertex( wkbPosition, dx, dy, hasZValue ); 03060 } 03061 } 03062 break; 03063 } 03064 03065 case QGis::WKBMultiPoint25D: 03066 hasZValue = true; 03067 case QGis::WKBMultiPoint: 03068 { 03069 int* npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03070 wkbPosition += sizeof( int ); 03071 for ( int index = 0; index < *npoints; ++index ) 03072 { 03073 wkbPosition += ( sizeof( int ) + 1 ); 03074 translateVertex( wkbPosition, dx, dy, hasZValue ); 03075 } 03076 break; 03077 } 03078 03079 case QGis::WKBMultiLineString25D: 03080 hasZValue = true; 03081 case QGis::WKBMultiLineString: 03082 { 03083 int* nlines = ( int* )( &( mGeometry[wkbPosition] ) ); 03084 int* npoints = 0; 03085 wkbPosition += sizeof( int ); 03086 for ( int index = 0; index < *nlines; ++index ) 03087 { 03088 wkbPosition += ( sizeof( int ) + 1 ); 03089 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03090 wkbPosition += sizeof( int ); 03091 for ( int index2 = 0; index2 < *npoints; ++index2 ) 03092 { 03093 translateVertex( wkbPosition, dx, dy, hasZValue ); 03094 } 03095 } 03096 break; 03097 } 03098 03099 case QGis::WKBMultiPolygon25D: 03100 hasZValue = true; 03101 case QGis::WKBMultiPolygon: 03102 { 03103 int* npolys = ( int* )( &( mGeometry[wkbPosition] ) ); 03104 int* nrings; 03105 int* npoints; 03106 wkbPosition += sizeof( int ); 03107 for ( int index = 0; index < *npolys; ++index ) 03108 { 03109 wkbPosition += ( 1 + sizeof( int ) ); //skip endian and polygon type 03110 nrings = ( int* )( &( mGeometry[wkbPosition] ) ); 03111 wkbPosition += sizeof( int ); 03112 for ( int index2 = 0; index2 < *nrings; ++index2 ) 03113 { 03114 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03115 wkbPosition += sizeof( int ); 03116 for ( int index3 = 0; index3 < *npoints; ++index3 ) 03117 { 03118 translateVertex( wkbPosition, dx, dy, hasZValue ); 03119 } 03120 } 03121 } 03122 } 03123 03124 default: 03125 break; 03126 } 03127 mDirtyGeos = true; 03128 return 0; 03129 } 03130 03131 int QgsGeometry::transform( const QgsCoordinateTransform& ct ) 03132 { 03133 if ( mDirtyWkb ) 03134 { 03135 exportGeosToWkb(); 03136 } 03137 03138 if ( !mGeometry ) 03139 { 03140 QgsDebugMsg( "WKB geometry not available!" ); 03141 return 1; 03142 } 03143 03144 QGis::WkbType wkbType; 03145 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) ); 03146 bool hasZValue = false; 03147 int wkbPosition = 5; 03148 03149 switch ( wkbType ) 03150 { 03151 case QGis::WKBPoint25D: 03152 case QGis::WKBPoint: 03153 { 03154 transformVertex( wkbPosition, ct, hasZValue ); 03155 } 03156 break; 03157 03158 case QGis::WKBLineString25D: 03159 hasZValue = true; 03160 case QGis::WKBLineString: 03161 { 03162 int* npoints = ( int* )( &mGeometry[wkbPosition] ); 03163 wkbPosition += sizeof( int ); 03164 for ( int index = 0; index < *npoints; ++index ) 03165 { 03166 transformVertex( wkbPosition, ct, hasZValue ); 03167 } 03168 break; 03169 } 03170 03171 case QGis::WKBPolygon25D: 03172 hasZValue = true; 03173 case QGis::WKBPolygon: 03174 { 03175 int* nrings = ( int* )( &( mGeometry[wkbPosition] ) ); 03176 wkbPosition += sizeof( int ); 03177 int* npoints; 03178 03179 for ( int index = 0; index < *nrings; ++index ) 03180 { 03181 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03182 wkbPosition += sizeof( int ); 03183 for ( int index2 = 0; index2 < *npoints; ++index2 ) 03184 { 03185 transformVertex( wkbPosition, ct, hasZValue ); 03186 } 03187 } 03188 break; 03189 } 03190 03191 case QGis::WKBMultiPoint25D: 03192 hasZValue = true; 03193 case QGis::WKBMultiPoint: 03194 { 03195 int* npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03196 wkbPosition += sizeof( int ); 03197 for ( int index = 0; index < *npoints; ++index ) 03198 { 03199 wkbPosition += ( sizeof( int ) + 1 ); 03200 transformVertex( wkbPosition, ct, hasZValue ); 03201 } 03202 break; 03203 } 03204 03205 case QGis::WKBMultiLineString25D: 03206 hasZValue = true; 03207 case QGis::WKBMultiLineString: 03208 { 03209 int* nlines = ( int* )( &( mGeometry[wkbPosition] ) ); 03210 int* npoints = 0; 03211 wkbPosition += sizeof( int ); 03212 for ( int index = 0; index < *nlines; ++index ) 03213 { 03214 wkbPosition += ( sizeof( int ) + 1 ); 03215 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03216 wkbPosition += sizeof( int ); 03217 for ( int index2 = 0; index2 < *npoints; ++index2 ) 03218 { 03219 transformVertex( wkbPosition, ct, hasZValue ); 03220 } 03221 } 03222 break; 03223 } 03224 03225 case QGis::WKBMultiPolygon25D: 03226 hasZValue = true; 03227 case QGis::WKBMultiPolygon: 03228 { 03229 int* npolys = ( int* )( &( mGeometry[wkbPosition] ) ); 03230 int* nrings; 03231 int* npoints; 03232 wkbPosition += sizeof( int ); 03233 for ( int index = 0; index < *npolys; ++index ) 03234 { 03235 wkbPosition += ( 1 + sizeof( int ) ); //skip endian and polygon type 03236 nrings = ( int* )( &( mGeometry[wkbPosition] ) ); 03237 wkbPosition += sizeof( int ); 03238 for ( int index2 = 0; index2 < *nrings; ++index2 ) 03239 { 03240 npoints = ( int* )( &( mGeometry[wkbPosition] ) ); 03241 wkbPosition += sizeof( int ); 03242 for ( int index3 = 0; index3 < *npoints; ++index3 ) 03243 { 03244 transformVertex( wkbPosition, ct, hasZValue ); 03245 } 03246 } 03247 } 03248 } 03249 03250 default: 03251 break; 03252 } 03253 mDirtyGeos = true; 03254 return 0; 03255 } 03256 03257 int QgsGeometry::splitGeometry( const QList<QgsPoint>& splitLine, QList<QgsGeometry*>& newGeometries, bool topological, QList<QgsPoint> &topologyTestPoints ) 03258 { 03259 int returnCode = 0; 03260 03261 //return if this type is point/multipoint 03262 if ( type() == QGis::Point ) 03263 { 03264 return 1; //cannot split points 03265 } 03266 03267 //make sure, mGeos and mWkb are there and up-to-date 03268 if ( mDirtyWkb ) 03269 { 03270 exportGeosToWkb(); 03271 } 03272 03273 if ( !mGeos || mDirtyGeos ) 03274 { 03275 if ( !exportWkbToGeos() ) 03276 return 1; 03277 } 03278 03279 if ( !GEOSisValid( mGeos ) ) 03280 { 03281 return 7; 03282 } 03283 03284 //make sure splitLine is valid 03285 if ( splitLine.size() < 2 ) 03286 { 03287 return 1; 03288 } 03289 03290 newGeometries.clear(); 03291 03292 try 03293 { 03294 GEOSGeometry *splitLineGeos = createGeosLineString( splitLine.toVector() ); 03295 if ( !GEOSisValid( splitLineGeos ) || !GEOSisSimple( splitLineGeos ) ) 03296 { 03297 GEOSGeom_destroy( splitLineGeos ); 03298 return 1; 03299 } 03300 03301 if ( topological ) 03302 { 03303 //find out candidate points for topological corrections 03304 if ( topologicalTestPointsSplit( splitLineGeos, topologyTestPoints ) != 0 ) 03305 { 03306 return 1; 03307 } 03308 } 03309 03310 //call split function depending on geometry type 03311 if ( type() == QGis::Line ) 03312 { 03313 returnCode = splitLinearGeometry( splitLineGeos, newGeometries ); 03314 GEOSGeom_destroy( splitLineGeos ); 03315 } 03316 else if ( type() == QGis::Polygon ) 03317 { 03318 returnCode = splitPolygonGeometry( splitLineGeos, newGeometries ); 03319 GEOSGeom_destroy( splitLineGeos ); 03320 } 03321 else 03322 { 03323 return 1; 03324 } 03325 } 03326 CATCH_GEOS( 2 ) 03327 03328 return returnCode; 03329 } 03330 03332 int QgsGeometry::reshapeGeometry( const QList<QgsPoint>& reshapeWithLine ) 03333 { 03334 if ( reshapeWithLine.size() < 2 ) 03335 { 03336 return 1; 03337 } 03338 03339 if ( type() == QGis::Point ) 03340 { 03341 return 1; //cannot reshape points 03342 } 03343 03344 GEOSGeometry* reshapeLineGeos = createGeosLineString( reshapeWithLine.toVector() ); 03345 03346 //make sure this geos geometry is up-to-date 03347 if ( !mGeos || mDirtyGeos ) 03348 { 03349 exportWkbToGeos(); 03350 } 03351 03352 //single or multi? 03353 int numGeoms = GEOSGetNumGeometries( mGeos ); 03354 if ( numGeoms == -1 ) 03355 { 03356 return 1; 03357 } 03358 03359 bool isMultiGeom = false; 03360 int geosTypeId = GEOSGeomTypeId( mGeos ); 03361 if ( geosTypeId == GEOS_MULTILINESTRING || geosTypeId == GEOS_MULTIPOLYGON ) 03362 { 03363 isMultiGeom = true; 03364 } 03365 03366 bool isLine = ( type() == QGis::Line ); 03367 03368 //polygon or multipolygon? 03369 if ( !isMultiGeom ) 03370 { 03371 GEOSGeometry* reshapedGeometry; 03372 if ( isLine ) 03373 { 03374 reshapedGeometry = reshapeLine( mGeos, reshapeLineGeos ); 03375 } 03376 else 03377 { 03378 reshapedGeometry = reshapePolygon( mGeos, reshapeLineGeos ); 03379 } 03380 03381 GEOSGeom_destroy( reshapeLineGeos ); 03382 if ( reshapedGeometry ) 03383 { 03384 GEOSGeom_destroy( mGeos ); 03385 mGeos = reshapedGeometry; 03386 mDirtyWkb = true; 03387 return 0; 03388 } 03389 else 03390 { 03391 return 1; 03392 } 03393 } 03394 else 03395 { 03396 //call reshape for each geometry part and replace mGeos with new geometry if reshape took place 03397 bool reshapeTookPlace = false; 03398 03399 GEOSGeometry* currentReshapeGeometry = 0; 03400 GEOSGeometry** newGeoms = new GEOSGeometry*[numGeoms]; 03401 03402 for ( int i = 0; i < numGeoms; ++i ) 03403 { 03404 if ( isLine ) 03405 { 03406 currentReshapeGeometry = reshapeLine( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos ); 03407 } 03408 else 03409 { 03410 currentReshapeGeometry = reshapePolygon( GEOSGetGeometryN( mGeos, i ), reshapeLineGeos ); 03411 } 03412 03413 if ( currentReshapeGeometry ) 03414 { 03415 newGeoms[i] = currentReshapeGeometry; 03416 reshapeTookPlace = true; 03417 } 03418 else 03419 { 03420 newGeoms[i] = GEOSGeom_clone( GEOSGetGeometryN( mGeos, i ) ); 03421 } 03422 } 03423 GEOSGeom_destroy( reshapeLineGeos ); 03424 03425 GEOSGeometry* newMultiGeom = 0; 03426 if ( isLine ) 03427 { 03428 newMultiGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, newGeoms, numGeoms ); 03429 } 03430 else //multipolygon 03431 { 03432 newMultiGeom = GEOSGeom_createCollection( GEOS_MULTIPOLYGON, newGeoms, numGeoms ); 03433 } 03434 03435 delete[] newGeoms; 03436 if ( ! newMultiGeom ) 03437 { 03438 return 3; 03439 } 03440 03441 if ( reshapeTookPlace ) 03442 { 03443 GEOSGeom_destroy( mGeos ); 03444 mGeos = newMultiGeom; 03445 mDirtyWkb = true; 03446 return 0; 03447 } 03448 else 03449 { 03450 GEOSGeom_destroy( newMultiGeom ); 03451 return 1; 03452 } 03453 } 03454 } 03455 03456 int QgsGeometry::makeDifference( QgsGeometry* other ) 03457 { 03458 //make sure geos geometry is up to date 03459 if ( !mGeos || mDirtyGeos ) 03460 { 03461 exportWkbToGeos(); 03462 } 03463 03464 if ( !mGeos ) 03465 { 03466 return 1; 03467 } 03468 03469 if ( !GEOSisValid( mGeos ) ) 03470 { 03471 return 2; 03472 } 03473 03474 if ( !GEOSisSimple( mGeos ) ) 03475 { 03476 return 3; 03477 } 03478 03479 //convert other geometry to geos 03480 if ( !other->mGeos || other->mDirtyGeos ) 03481 { 03482 other->exportWkbToGeos(); 03483 } 03484 03485 if ( !other->mGeos ) 03486 { 03487 return 4; 03488 } 03489 03490 //make geometry::difference 03491 try 03492 { 03493 if ( GEOSIntersects( mGeos, other->mGeos ) ) 03494 { 03495 //check if multitype before and after 03496 bool multiType = isMultipart(); 03497 03498 mGeos = GEOSDifference( mGeos, other->mGeos ); 03499 mDirtyWkb = true; 03500 03501 if ( multiType && !isMultipart() ) 03502 { 03503 convertToMultiType(); 03504 exportWkbToGeos(); 03505 } 03506 } 03507 else 03508 { 03509 return 0; //nothing to do 03510 } 03511 } 03512 CATCH_GEOS( 5 ) 03513 03514 if ( !mGeos ) 03515 { 03516 mDirtyGeos = true; 03517 return 6; 03518 } 03519 03520 return 0; 03521 } 03522 03523 QgsRectangle QgsGeometry::boundingBox() 03524 { 03525 double xmin = std::numeric_limits<double>::max(); 03526 double ymin = std::numeric_limits<double>::max(); 03527 double xmax = -std::numeric_limits<double>::max(); 03528 double ymax = -std::numeric_limits<double>::max(); 03529 03530 double *x; 03531 double *y; 03532 int *nPoints; 03533 int *numRings; 03534 int *numPolygons; 03535 int numLineStrings; 03536 int idx, jdx, kdx; 03537 unsigned char *ptr; 03538 QgsPoint pt; 03539 QGis::WkbType wkbType; 03540 bool hasZValue = false; 03541 03542 // TODO: implement with GEOS 03543 if ( mDirtyWkb ) 03544 { 03545 exportGeosToWkb(); 03546 } 03547 03548 if ( !mGeometry ) 03549 { 03550 QgsDebugMsg( "WKB geometry not available!" ); 03551 return QgsRectangle( 0, 0, 0, 0 ); 03552 } 03553 // consider endian when fetching feature type 03554 //wkbType = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4]; //MH: Does not work for 25D geometries 03555 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) ); 03556 switch ( wkbType ) 03557 { 03558 case QGis::WKBPoint25D: 03559 case QGis::WKBPoint: 03560 x = ( double * )( mGeometry + 5 ); 03561 y = ( double * )( mGeometry + 5 + sizeof( double ) ); 03562 if ( *x < xmin ) 03563 { 03564 xmin = *x; 03565 } 03566 if ( *x > xmax ) 03567 { 03568 xmax = *x; 03569 } 03570 if ( *y < ymin ) 03571 { 03572 ymin = *y; 03573 } 03574 if ( *y > ymax ) 03575 { 03576 ymax = *y; 03577 } 03578 break; 03579 case QGis::WKBMultiPoint25D: 03580 hasZValue = true; 03581 case QGis::WKBMultiPoint: 03582 { 03583 ptr = mGeometry + 1 + sizeof( int ); 03584 nPoints = ( int * ) ptr; 03585 ptr += sizeof( int ); 03586 for ( idx = 0; idx < *nPoints; idx++ ) 03587 { 03588 ptr += ( 1 + sizeof( int ) ); 03589 x = ( double * ) ptr; 03590 ptr += sizeof( double ); 03591 y = ( double * ) ptr; 03592 ptr += sizeof( double ); 03593 if ( hasZValue ) 03594 { 03595 ptr += sizeof( double ); 03596 } 03597 if ( *x < xmin ) 03598 { 03599 xmin = *x; 03600 } 03601 if ( *x > xmax ) 03602 { 03603 xmax = *x; 03604 } 03605 if ( *y < ymin ) 03606 { 03607 ymin = *y; 03608 } 03609 if ( *y > ymax ) 03610 { 03611 ymax = *y; 03612 } 03613 } 03614 break; 03615 } 03616 case QGis::WKBLineString25D: 03617 hasZValue = true; 03618 case QGis::WKBLineString: 03619 { 03620 // get number of points in the line 03621 ptr = mGeometry + 5; 03622 nPoints = ( int * ) ptr; 03623 ptr = mGeometry + 1 + 2 * sizeof( int ); 03624 for ( idx = 0; idx < *nPoints; idx++ ) 03625 { 03626 x = ( double * ) ptr; 03627 ptr += sizeof( double ); 03628 y = ( double * ) ptr; 03629 ptr += sizeof( double ); 03630 if ( hasZValue ) 03631 { 03632 ptr += sizeof( double ); 03633 } 03634 if ( *x < xmin ) 03635 { 03636 xmin = *x; 03637 } 03638 if ( *x > xmax ) 03639 { 03640 xmax = *x; 03641 } 03642 if ( *y < ymin ) 03643 { 03644 ymin = *y; 03645 } 03646 if ( *y > ymax ) 03647 { 03648 ymax = *y; 03649 } 03650 } 03651 break; 03652 } 03653 case QGis::WKBMultiLineString25D: 03654 hasZValue = true; 03655 case QGis::WKBMultiLineString: 03656 { 03657 numLineStrings = ( int )( mGeometry[5] ); 03658 ptr = mGeometry + 9; 03659 for ( jdx = 0; jdx < numLineStrings; jdx++ ) 03660 { 03661 // each of these is a wbklinestring so must handle as such 03662 ptr += 5; // skip type since we know its 2 03663 nPoints = ( int * ) ptr; 03664 ptr += sizeof( int ); 03665 for ( idx = 0; idx < *nPoints; idx++ ) 03666 { 03667 x = ( double * ) ptr; 03668 ptr += sizeof( double ); 03669 y = ( double * ) ptr; 03670 ptr += sizeof( double ); 03671 if ( hasZValue ) 03672 { 03673 ptr += sizeof( double ); 03674 } 03675 if ( *x < xmin ) 03676 { 03677 xmin = *x; 03678 } 03679 if ( *x > xmax ) 03680 { 03681 xmax = *x; 03682 } 03683 if ( *y < ymin ) 03684 { 03685 ymin = *y; 03686 } 03687 if ( *y > ymax ) 03688 { 03689 ymax = *y; 03690 } 03691 } 03692 } 03693 break; 03694 } 03695 case QGis::WKBPolygon25D: 03696 hasZValue = true; 03697 case QGis::WKBPolygon: 03698 { 03699 // get number of rings in the polygon 03700 numRings = ( int * )( mGeometry + 1 + sizeof( int ) ); 03701 ptr = mGeometry + 1 + 2 * sizeof( int ); 03702 for ( idx = 0; idx < *numRings; idx++ ) 03703 { 03704 // get number of points in the ring 03705 nPoints = ( int * ) ptr; 03706 ptr += 4; 03707 for ( jdx = 0; jdx < *nPoints; jdx++ ) 03708 { 03709 // add points to a point array for drawing the polygon 03710 x = ( double * ) ptr; 03711 ptr += sizeof( double ); 03712 y = ( double * ) ptr; 03713 ptr += sizeof( double ); 03714 if ( hasZValue ) 03715 { 03716 ptr += sizeof( double ); 03717 } 03718 if ( *x < xmin ) 03719 { 03720 xmin = *x; 03721 } 03722 if ( *x > xmax ) 03723 { 03724 xmax = *x; 03725 } 03726 if ( *y < ymin ) 03727 { 03728 ymin = *y; 03729 } 03730 if ( *y > ymax ) 03731 { 03732 ymax = *y; 03733 } 03734 } 03735 } 03736 break; 03737 } 03738 case QGis::WKBMultiPolygon25D: 03739 hasZValue = true; 03740 case QGis::WKBMultiPolygon: 03741 { 03742 // get the number of polygons 03743 ptr = mGeometry + 5; 03744 numPolygons = ( int * ) ptr; 03745 ptr += 4; 03746 03747 for ( kdx = 0; kdx < *numPolygons; kdx++ ) 03748 { 03749 //skip the endian and mGeometry type info and 03750 // get number of rings in the polygon 03751 ptr += 5; 03752 numRings = ( int * ) ptr; 03753 ptr += 4; 03754 for ( idx = 0; idx < *numRings; idx++ ) 03755 { 03756 // get number of points in the ring 03757 nPoints = ( int * ) ptr; 03758 ptr += 4; 03759 for ( jdx = 0; jdx < *nPoints; jdx++ ) 03760 { 03761 // add points to a point array for drawing the polygon 03762 x = ( double * ) ptr; 03763 ptr += sizeof( double ); 03764 y = ( double * ) ptr; 03765 ptr += sizeof( double ); 03766 if ( hasZValue ) 03767 { 03768 ptr += sizeof( double ); 03769 } 03770 if ( *x < xmin ) 03771 { 03772 xmin = *x; 03773 } 03774 if ( *x > xmax ) 03775 { 03776 xmax = *x; 03777 } 03778 if ( *y < ymin ) 03779 { 03780 ymin = *y; 03781 } 03782 if ( *y > ymax ) 03783 { 03784 ymax = *y; 03785 } 03786 } 03787 } 03788 } 03789 break; 03790 } 03791 03792 default: 03793 QgsDebugMsg( "Unknown WkbType ENCOUNTERED" ); 03794 return QgsRectangle( 0, 0, 0, 0 ); 03795 break; 03796 03797 } 03798 return QgsRectangle( xmin, ymin, xmax, ymax ); 03799 } 03800 03801 bool QgsGeometry::intersects( const QgsRectangle& r ) 03802 { 03803 QgsGeometry* g = fromRect( r ); 03804 bool res = intersects( g ); 03805 delete g; 03806 return res; 03807 } 03808 03809 bool QgsGeometry::intersects( QgsGeometry* geometry ) 03810 { 03811 try // geos might throw exception on error 03812 { 03813 // ensure that both geometries have geos geometry 03814 exportWkbToGeos(); 03815 geometry->exportWkbToGeos(); 03816 03817 if ( !mGeos || !geometry->mGeos ) 03818 { 03819 QgsDebugMsg( "GEOS geometry not available!" ); 03820 return false; 03821 } 03822 03823 return GEOSIntersects( mGeos, geometry->mGeos ); 03824 } 03825 CATCH_GEOS( false ) 03826 } 03827 03828 03829 bool QgsGeometry::contains( QgsPoint* p ) 03830 { 03831 exportWkbToGeos(); 03832 03833 if ( !mGeos ) 03834 { 03835 QgsDebugMsg( "GEOS geometry not available!" ); 03836 return false; 03837 } 03838 03839 GEOSGeometry *geosPoint = 0; 03840 03841 bool returnval = false; 03842 03843 try 03844 { 03845 geosPoint = createGeosPoint( *p ); 03846 returnval = GEOSContains( mGeos, geosPoint ); 03847 } 03848 catch ( GEOSException &e ) 03849 { 03850 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 03851 returnval = false; 03852 } 03853 03854 if ( geosPoint ) 03855 GEOSGeom_destroy( geosPoint ); 03856 03857 return returnval; 03858 } 03859 03860 bool QgsGeometry::geosRelOp( 03861 char( *op )( const GEOSGeometry*, const GEOSGeometry * ), 03862 QgsGeometry *a, 03863 QgsGeometry *b ) 03864 { 03865 try // geos might throw exception on error 03866 { 03867 // ensure that both geometries have geos geometry 03868 a->exportWkbToGeos(); 03869 b->exportWkbToGeos(); 03870 03871 if ( !a->mGeos || !b->mGeos ) 03872 { 03873 QgsDebugMsg( "GEOS geometry not available!" ); 03874 return false; 03875 } 03876 return op( a->mGeos, b->mGeos ); 03877 } 03878 CATCH_GEOS( false ) 03879 } 03880 03881 bool QgsGeometry::contains( QgsGeometry* geometry ) 03882 { 03883 return geosRelOp( GEOSContains, this, geometry ); 03884 } 03885 03886 bool QgsGeometry::disjoint( QgsGeometry* geometry ) 03887 { 03888 return geosRelOp( GEOSDisjoint, this, geometry ); 03889 } 03890 03891 bool QgsGeometry::equals( QgsGeometry* geometry ) 03892 { 03893 return geosRelOp( GEOSEquals, this, geometry ); 03894 } 03895 03896 bool QgsGeometry::touches( QgsGeometry* geometry ) 03897 { 03898 return geosRelOp( GEOSTouches, this, geometry ); 03899 } 03900 03901 bool QgsGeometry::overlaps( QgsGeometry* geometry ) 03902 { 03903 return geosRelOp( GEOSOverlaps, this, geometry ); 03904 } 03905 03906 bool QgsGeometry::within( QgsGeometry* geometry ) 03907 { 03908 return geosRelOp( GEOSWithin, this, geometry ); 03909 } 03910 03911 bool QgsGeometry::crosses( QgsGeometry* geometry ) 03912 { 03913 return geosRelOp( GEOSCrosses, this, geometry ); 03914 } 03915 03916 QString QgsGeometry::exportToWkt() 03917 { 03918 QgsDebugMsg( "entered." ); 03919 03920 // TODO: implement with GEOS 03921 if ( mDirtyWkb ) 03922 { 03923 exportGeosToWkb(); 03924 } 03925 03926 if ( !mGeometry ) 03927 { 03928 QgsDebugMsg( "WKB geometry not available!" ); 03929 return QString::null; 03930 } 03931 03932 QGis::WkbType wkbType; 03933 bool hasZValue = false; 03934 double *x, *y; 03935 03936 QString mWkt; // TODO: rename 03937 03938 // Will this really work when mGeometry[0] == 0 ???? I (gavin) think not. 03939 //wkbType = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4]; 03940 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) ); 03941 03942 switch ( wkbType ) 03943 { 03944 case QGis::WKBPoint25D: 03945 case QGis::WKBPoint: 03946 { 03947 mWkt += "POINT("; 03948 x = ( double * )( mGeometry + 5 ); 03949 mWkt += QString::number( *x, 'f', 6 ); 03950 mWkt += " "; 03951 y = ( double * )( mGeometry + 5 + sizeof( double ) ); 03952 mWkt += QString::number( *y, 'f', 6 ); 03953 mWkt += ")"; 03954 return mWkt; 03955 } 03956 03957 case QGis::WKBLineString25D: 03958 hasZValue = true; 03959 case QGis::WKBLineString: 03960 { 03961 QgsDebugMsg( "LINESTRING found" ); 03962 unsigned char *ptr; 03963 int *nPoints; 03964 int idx; 03965 03966 mWkt += "LINESTRING("; 03967 // get number of points in the line 03968 ptr = mGeometry + 5; 03969 nPoints = ( int * ) ptr; 03970 ptr = mGeometry + 1 + 2 * sizeof( int ); 03971 for ( idx = 0; idx < *nPoints; ++idx ) 03972 { 03973 if ( idx != 0 ) 03974 { 03975 mWkt += ", "; 03976 } 03977 x = ( double * ) ptr; 03978 mWkt += QString::number( *x, 'f', 6 ); 03979 mWkt += " "; 03980 ptr += sizeof( double ); 03981 y = ( double * ) ptr; 03982 mWkt += QString::number( *y, 'f', 6 ); 03983 ptr += sizeof( double ); 03984 if ( hasZValue ) 03985 { 03986 ptr += sizeof( double ); 03987 } 03988 } 03989 mWkt += ")"; 03990 return mWkt; 03991 } 03992 03993 case QGis::WKBPolygon25D: 03994 hasZValue = true; 03995 case QGis::WKBPolygon: 03996 { 03997 QgsDebugMsg( "POLYGON found" ); 03998 unsigned char *ptr; 03999 int idx, jdx; 04000 int *numRings, *nPoints; 04001 04002 mWkt += "POLYGON("; 04003 // get number of rings in the polygon 04004 numRings = ( int * )( mGeometry + 1 + sizeof( int ) ); 04005 if ( !( *numRings ) ) // sanity check for zero rings in polygon 04006 { 04007 return QString(); 04008 } 04009 int *ringStart; // index of first point for each ring 04010 int *ringNumPoints; // number of points in each ring 04011 ringStart = new int[*numRings]; 04012 ringNumPoints = new int[*numRings]; 04013 ptr = mGeometry + 1 + 2 * sizeof( int ); // set pointer to the first ring 04014 for ( idx = 0; idx < *numRings; idx++ ) 04015 { 04016 if ( idx != 0 ) 04017 { 04018 mWkt += ","; 04019 } 04020 mWkt += "("; 04021 // get number of points in the ring 04022 nPoints = ( int * ) ptr; 04023 ringNumPoints[idx] = *nPoints; 04024 ptr += 4; 04025 04026 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04027 { 04028 if ( jdx != 0 ) 04029 { 04030 mWkt += ","; 04031 } 04032 x = ( double * ) ptr; 04033 mWkt += QString::number( *x, 'f', 6 ); 04034 mWkt += " "; 04035 ptr += sizeof( double ); 04036 y = ( double * ) ptr; 04037 mWkt += QString::number( *y, 'f', 6 ); 04038 ptr += sizeof( double ); 04039 if ( hasZValue ) 04040 { 04041 ptr += sizeof( double ); 04042 } 04043 } 04044 mWkt += ")"; 04045 } 04046 mWkt += ")"; 04047 delete [] ringStart; 04048 delete [] ringNumPoints; 04049 return mWkt; 04050 } 04051 04052 case QGis::WKBMultiPoint25D: 04053 hasZValue = true; 04054 case QGis::WKBMultiPoint: 04055 { 04056 unsigned char *ptr; 04057 int idx; 04058 int *nPoints; 04059 04060 mWkt += "MULTIPOINT("; 04061 nPoints = ( int* )( mGeometry + 5 ); 04062 ptr = mGeometry + 5 + sizeof( int ); 04063 for ( idx = 0; idx < *nPoints; ++idx ) 04064 { 04065 ptr += ( 1 + sizeof( int ) ); 04066 if ( idx != 0 ) 04067 { 04068 mWkt += ", "; 04069 } 04070 x = ( double * )( ptr ); 04071 mWkt += QString::number( *x, 'f', 6 ); 04072 mWkt += " "; 04073 ptr += sizeof( double ); 04074 y = ( double * )( ptr ); 04075 mWkt += QString::number( *y, 'f', 6 ); 04076 ptr += sizeof( double ); 04077 if ( hasZValue ) 04078 { 04079 ptr += sizeof( double ); 04080 } 04081 } 04082 mWkt += ")"; 04083 return mWkt; 04084 } 04085 04086 case QGis::WKBMultiLineString25D: 04087 hasZValue = true; 04088 case QGis::WKBMultiLineString: 04089 { 04090 QgsDebugMsg( "MULTILINESTRING found" ); 04091 unsigned char *ptr; 04092 int idx, jdx, numLineStrings; 04093 int *nPoints; 04094 04095 mWkt += "MULTILINESTRING("; 04096 numLineStrings = ( int )( mGeometry[5] ); 04097 ptr = mGeometry + 9; 04098 for ( jdx = 0; jdx < numLineStrings; jdx++ ) 04099 { 04100 if ( jdx != 0 ) 04101 { 04102 mWkt += ", "; 04103 } 04104 mWkt += "("; 04105 ptr += 5; // skip type since we know its 2 04106 nPoints = ( int * ) ptr; 04107 ptr += sizeof( int ); 04108 for ( idx = 0; idx < *nPoints; idx++ ) 04109 { 04110 if ( idx != 0 ) 04111 { 04112 mWkt += ", "; 04113 } 04114 x = ( double * ) ptr; 04115 mWkt += QString::number( *x, 'f', 6 ); 04116 ptr += sizeof( double ); 04117 mWkt += " "; 04118 y = ( double * ) ptr; 04119 mWkt += QString::number( *y, 'f', 6 ); 04120 ptr += sizeof( double ); 04121 if ( hasZValue ) 04122 { 04123 ptr += sizeof( double ); 04124 } 04125 } 04126 mWkt += ")"; 04127 } 04128 mWkt += ")"; 04129 return mWkt; 04130 } 04131 04132 case QGis::WKBMultiPolygon25D: 04133 hasZValue = true; 04134 case QGis::WKBMultiPolygon: 04135 { 04136 QgsDebugMsg( "MULTIPOLYGON found" ); 04137 unsigned char *ptr; 04138 int idx, jdx, kdx; 04139 int *numPolygons, *numRings, *nPoints; 04140 04141 mWkt += "MULTIPOLYGON("; 04142 ptr = mGeometry + 5; 04143 numPolygons = ( int * ) ptr; 04144 ptr = mGeometry + 9; 04145 for ( kdx = 0; kdx < *numPolygons; kdx++ ) 04146 { 04147 if ( kdx != 0 ) 04148 { 04149 mWkt += ","; 04150 } 04151 mWkt += "("; 04152 ptr += 5; 04153 numRings = ( int * ) ptr; 04154 ptr += 4; 04155 for ( idx = 0; idx < *numRings; idx++ ) 04156 { 04157 if ( idx != 0 ) 04158 { 04159 mWkt += ","; 04160 } 04161 mWkt += "("; 04162 nPoints = ( int * ) ptr; 04163 ptr += 4; 04164 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04165 { 04166 if ( jdx != 0 ) 04167 { 04168 mWkt += ","; 04169 } 04170 x = ( double * ) ptr; 04171 mWkt += QString::number( *x, 'f', 6 ); 04172 ptr += sizeof( double ); 04173 mWkt += " "; 04174 y = ( double * ) ptr; 04175 mWkt += QString::number( *y, 'f', 6 ); 04176 ptr += sizeof( double ); 04177 if ( hasZValue ) 04178 { 04179 ptr += sizeof( double ); 04180 } 04181 } 04182 mWkt += ")"; 04183 } 04184 mWkt += ")"; 04185 } 04186 mWkt += ")"; 04187 return mWkt; 04188 } 04189 04190 default: 04191 QgsDebugMsg( "error: mGeometry type not recognized" ); 04192 return QString::null; 04193 } 04194 } 04195 04196 QString QgsGeometry::exportToGeoJSON() 04197 { 04198 QgsDebugMsg( "entered." ); 04199 04200 // TODO: implement with GEOS 04201 if ( mDirtyWkb ) 04202 { 04203 exportGeosToWkb(); 04204 } 04205 04206 if ( !mGeometry ) 04207 { 04208 QgsDebugMsg( "WKB geometry not available!" ); 04209 return QString::null; 04210 } 04211 04212 QGis::WkbType wkbType; 04213 bool hasZValue = false; 04214 double *x, *y; 04215 04216 QString mWkt; // TODO: rename 04217 04218 // Will this really work when mGeometry[0] == 0 ???? I (gavin) think not. 04219 //wkbType = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4]; 04220 memcpy( &wkbType, &( mGeometry[1] ), sizeof( int ) ); 04221 04222 switch ( wkbType ) 04223 { 04224 case QGis::WKBPoint25D: 04225 case QGis::WKBPoint: 04226 { 04227 mWkt += "{ \"type\": \"Point\", \"coordinates\": ["; 04228 x = ( double * )( mGeometry + 5 ); 04229 mWkt += QString::number( *x, 'f', 6 ); 04230 mWkt += ", "; 04231 y = ( double * )( mGeometry + 5 + sizeof( double ) ); 04232 mWkt += QString::number( *y, 'f', 6 ); 04233 mWkt += "] }"; 04234 return mWkt; 04235 } 04236 04237 case QGis::WKBLineString25D: 04238 hasZValue = true; 04239 case QGis::WKBLineString: 04240 { 04241 QgsDebugMsg( "LINESTRING found" ); 04242 unsigned char *ptr; 04243 int *nPoints; 04244 int idx; 04245 04246 mWkt += "{ \"type\": \"LineString\", \"coordinates\": [ "; 04247 // get number of points in the line 04248 ptr = mGeometry + 5; 04249 nPoints = ( int * ) ptr; 04250 ptr = mGeometry + 1 + 2 * sizeof( int ); 04251 for ( idx = 0; idx < *nPoints; ++idx ) 04252 { 04253 if ( idx != 0 ) 04254 { 04255 mWkt += ", "; 04256 } 04257 mWkt += "["; 04258 x = ( double * ) ptr; 04259 mWkt += QString::number( *x, 'f', 6 ); 04260 mWkt += ", "; 04261 ptr += sizeof( double ); 04262 y = ( double * ) ptr; 04263 mWkt += QString::number( *y, 'f', 6 ); 04264 ptr += sizeof( double ); 04265 if ( hasZValue ) 04266 { 04267 ptr += sizeof( double ); 04268 } 04269 mWkt += "]"; 04270 } 04271 mWkt += " ] }"; 04272 return mWkt; 04273 } 04274 04275 case QGis::WKBPolygon25D: 04276 hasZValue = true; 04277 case QGis::WKBPolygon: 04278 { 04279 QgsDebugMsg( "POLYGON found" ); 04280 unsigned char *ptr; 04281 int idx, jdx; 04282 int *numRings, *nPoints; 04283 04284 mWkt += "{ \"type\": \"Polygon\", \"coordinates\": [ "; 04285 // get number of rings in the polygon 04286 numRings = ( int * )( mGeometry + 1 + sizeof( int ) ); 04287 if ( !( *numRings ) ) // sanity check for zero rings in polygon 04288 { 04289 return QString(); 04290 } 04291 int *ringStart; // index of first point for each ring 04292 int *ringNumPoints; // number of points in each ring 04293 ringStart = new int[*numRings]; 04294 ringNumPoints = new int[*numRings]; 04295 ptr = mGeometry + 1 + 2 * sizeof( int ); // set pointer to the first ring 04296 for ( idx = 0; idx < *numRings; idx++ ) 04297 { 04298 if ( idx != 0 ) 04299 { 04300 mWkt += ", "; 04301 } 04302 mWkt += "[ "; 04303 // get number of points in the ring 04304 nPoints = ( int * ) ptr; 04305 ringNumPoints[idx] = *nPoints; 04306 ptr += 4; 04307 04308 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04309 { 04310 if ( jdx != 0 ) 04311 { 04312 mWkt += ", "; 04313 } 04314 mWkt += "["; 04315 x = ( double * ) ptr; 04316 mWkt += QString::number( *x, 'f', 6 ); 04317 mWkt += ", "; 04318 ptr += sizeof( double ); 04319 y = ( double * ) ptr; 04320 mWkt += QString::number( *y, 'f', 6 ); 04321 ptr += sizeof( double ); 04322 if ( hasZValue ) 04323 { 04324 ptr += sizeof( double ); 04325 } 04326 mWkt += "]"; 04327 } 04328 mWkt += " ]"; 04329 } 04330 mWkt += " ] }"; 04331 delete [] ringStart; 04332 delete [] ringNumPoints; 04333 return mWkt; 04334 } 04335 04336 case QGis::WKBMultiPoint25D: 04337 hasZValue = true; 04338 case QGis::WKBMultiPoint: 04339 { 04340 unsigned char *ptr; 04341 int idx; 04342 int *nPoints; 04343 04344 mWkt += "{ \"type\": \"MultiPoint\", \"coordinates\": [ "; 04345 nPoints = ( int* )( mGeometry + 5 ); 04346 ptr = mGeometry + 5 + sizeof( int ); 04347 for ( idx = 0; idx < *nPoints; ++idx ) 04348 { 04349 ptr += ( 1 + sizeof( int ) ); 04350 if ( idx != 0 ) 04351 { 04352 mWkt += ", "; 04353 } 04354 mWkt += "["; 04355 x = ( double * )( ptr ); 04356 mWkt += QString::number( *x, 'f', 6 ); 04357 mWkt += ", "; 04358 ptr += sizeof( double ); 04359 y = ( double * )( ptr ); 04360 mWkt += QString::number( *y, 'f', 6 ); 04361 ptr += sizeof( double ); 04362 if ( hasZValue ) 04363 { 04364 ptr += sizeof( double ); 04365 } 04366 mWkt += "]"; 04367 } 04368 mWkt += " ] }"; 04369 return mWkt; 04370 } 04371 04372 case QGis::WKBMultiLineString25D: 04373 hasZValue = true; 04374 case QGis::WKBMultiLineString: 04375 { 04376 QgsDebugMsg( "MULTILINESTRING found" ); 04377 unsigned char *ptr; 04378 int idx, jdx, numLineStrings; 04379 int *nPoints; 04380 04381 mWkt += "{ \"type\": \"MultiLineString\", \"coordinates\": [ "; 04382 numLineStrings = ( int )( mGeometry[5] ); 04383 ptr = mGeometry + 9; 04384 for ( jdx = 0; jdx < numLineStrings; jdx++ ) 04385 { 04386 if ( jdx != 0 ) 04387 { 04388 mWkt += ", "; 04389 } 04390 mWkt += "[ "; 04391 ptr += 5; // skip type since we know its 2 04392 nPoints = ( int * ) ptr; 04393 ptr += sizeof( int ); 04394 for ( idx = 0; idx < *nPoints; idx++ ) 04395 { 04396 if ( idx != 0 ) 04397 { 04398 mWkt += ", "; 04399 } 04400 mWkt += "["; 04401 x = ( double * ) ptr; 04402 mWkt += QString::number( *x, 'f', 6 ); 04403 ptr += sizeof( double ); 04404 mWkt += ", "; 04405 y = ( double * ) ptr; 04406 mWkt += QString::number( *y, 'f', 6 ); 04407 ptr += sizeof( double ); 04408 if ( hasZValue ) 04409 { 04410 ptr += sizeof( double ); 04411 } 04412 mWkt += "]"; 04413 } 04414 mWkt += " ]"; 04415 } 04416 mWkt += " ] }"; 04417 return mWkt; 04418 } 04419 04420 case QGis::WKBMultiPolygon25D: 04421 hasZValue = true; 04422 case QGis::WKBMultiPolygon: 04423 { 04424 QgsDebugMsg( "MULTIPOLYGON found" ); 04425 unsigned char *ptr; 04426 int idx, jdx, kdx; 04427 int *numPolygons, *numRings, *nPoints; 04428 04429 mWkt += "{ \"type\": \"MultiPolygon\", \"coordinates\": [ "; 04430 ptr = mGeometry + 5; 04431 numPolygons = ( int * ) ptr; 04432 ptr = mGeometry + 9; 04433 for ( kdx = 0; kdx < *numPolygons; kdx++ ) 04434 { 04435 if ( kdx != 0 ) 04436 { 04437 mWkt += ", "; 04438 } 04439 mWkt += "[ "; 04440 ptr += 5; 04441 numRings = ( int * ) ptr; 04442 ptr += 4; 04443 for ( idx = 0; idx < *numRings; idx++ ) 04444 { 04445 if ( idx != 0 ) 04446 { 04447 mWkt += ", "; 04448 } 04449 mWkt += "[ "; 04450 nPoints = ( int * ) ptr; 04451 ptr += 4; 04452 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04453 { 04454 if ( jdx != 0 ) 04455 { 04456 mWkt += ", "; 04457 } 04458 mWkt += "["; 04459 x = ( double * ) ptr; 04460 mWkt += QString::number( *x, 'f', 6 ); 04461 ptr += sizeof( double ); 04462 mWkt += ", "; 04463 y = ( double * ) ptr; 04464 mWkt += QString::number( *y, 'f', 6 ); 04465 ptr += sizeof( double ); 04466 if ( hasZValue ) 04467 { 04468 ptr += sizeof( double ); 04469 } 04470 mWkt += "]"; 04471 } 04472 mWkt += " ]"; 04473 } 04474 mWkt += " ]"; 04475 } 04476 mWkt += " ] }"; 04477 return mWkt; 04478 } 04479 04480 default: 04481 QgsDebugMsg( "error: mGeometry type not recognized" ); 04482 return QString::null; 04483 } 04484 } 04485 04486 bool QgsGeometry::exportWkbToGeos() 04487 { 04488 QgsDebugMsgLevel( "entered.", 3 ); 04489 04490 if ( !mDirtyGeos ) 04491 { 04492 // No need to convert again 04493 return true; 04494 } 04495 04496 if ( mGeos ) 04497 { 04498 GEOSGeom_destroy( mGeos ); 04499 mGeos = 0; 04500 } 04501 04502 // this probably shouldn't return true 04503 if ( !mGeometry ) 04504 { 04505 // no WKB => no GEOS 04506 mDirtyGeos = false; 04507 return true; 04508 } 04509 04510 double *x; 04511 double *y; 04512 int *nPoints; 04513 int *numRings; 04514 int *numPolygons; 04515 int numLineStrings; 04516 int idx, jdx, kdx; 04517 unsigned char *ptr; 04518 QgsPoint pt; 04519 QGis::WkbType wkbtype; 04520 bool hasZValue = false; 04521 04522 //wkbtype = (mGeometry[0] == 1) ? mGeometry[1] : mGeometry[4]; 04523 memcpy( &wkbtype, &( mGeometry[1] ), sizeof( int ) ); 04524 04525 try 04526 { 04527 switch ( wkbtype ) 04528 { 04529 case QGis::WKBPoint25D: 04530 case QGis::WKBPoint: 04531 { 04532 x = ( double * )( mGeometry + 5 ); 04533 y = ( double * )( mGeometry + 5 + sizeof( double ) ); 04534 04535 mGeos = createGeosPoint( QgsPoint( *x, *y ) ); 04536 mDirtyGeos = false; 04537 break; 04538 } 04539 04540 case QGis::WKBMultiPoint25D: 04541 hasZValue = true; 04542 case QGis::WKBMultiPoint: 04543 { 04544 QVector<GEOSGeometry *> points; 04545 04546 ptr = mGeometry + 5; 04547 nPoints = ( int * ) ptr; 04548 ptr = mGeometry + 1 + 2 * sizeof( int ); 04549 for ( idx = 0; idx < *nPoints; idx++ ) 04550 { 04551 ptr += ( 1 + sizeof( int ) ); 04552 x = ( double * ) ptr; 04553 ptr += sizeof( double ); 04554 y = ( double * ) ptr; 04555 ptr += sizeof( double ); 04556 if ( hasZValue ) 04557 { 04558 ptr += sizeof( double ); 04559 } 04560 points << createGeosPoint( QgsPoint( *x, *y ) ); 04561 } 04562 mGeos = createGeosCollection( GEOS_MULTIPOINT, points ); 04563 mDirtyGeos = false; 04564 break; 04565 } 04566 04567 case QGis::WKBLineString25D: 04568 hasZValue = true; 04569 case QGis::WKBLineString: 04570 { 04571 QgsDebugMsgLevel( "Linestring found", 3 ); 04572 04573 QgsPolyline sequence; 04574 04575 ptr = mGeometry + 5; 04576 nPoints = ( int * ) ptr; 04577 ptr = mGeometry + 1 + 2 * sizeof( int ); 04578 for ( idx = 0; idx < *nPoints; idx++ ) 04579 { 04580 x = ( double * ) ptr; 04581 ptr += sizeof( double ); 04582 y = ( double * ) ptr; 04583 ptr += sizeof( double ); 04584 if ( hasZValue ) 04585 { 04586 ptr += sizeof( double ); 04587 } 04588 04589 sequence << QgsPoint( *x, *y ); 04590 } 04591 mDirtyGeos = false; 04592 mGeos = createGeosLineString( sequence ); 04593 break; 04594 } 04595 04596 case QGis::WKBMultiLineString25D: 04597 hasZValue = true; 04598 case QGis::WKBMultiLineString: 04599 { 04600 QVector<GEOSGeometry*> lines; 04601 numLineStrings = ( int )( mGeometry[5] ); 04602 ptr = ( mGeometry + 9 ); 04603 for ( jdx = 0; jdx < numLineStrings; jdx++ ) 04604 { 04605 QgsPolyline sequence; 04606 04607 // each of these is a wbklinestring so must handle as such 04608 ptr += 5; // skip type since we know its 2 04609 nPoints = ( int * ) ptr; 04610 ptr += sizeof( int ); 04611 for ( idx = 0; idx < *nPoints; idx++ ) 04612 { 04613 x = ( double * ) ptr; 04614 ptr += sizeof( double ); 04615 y = ( double * ) ptr; 04616 ptr += sizeof( double ); 04617 if ( hasZValue ) 04618 { 04619 ptr += sizeof( double ); 04620 } 04621 sequence << QgsPoint( *x, *y ); 04622 } 04623 lines << createGeosLineString( sequence ); 04624 } 04625 mGeos = createGeosCollection( GEOS_MULTILINESTRING, lines ); 04626 mDirtyGeos = false; 04627 break; 04628 } 04629 04630 case QGis::WKBPolygon25D: 04631 hasZValue = true; 04632 case QGis::WKBPolygon: 04633 { 04634 QgsDebugMsgLevel( "Polygon found", 3 ); 04635 04636 // get number of rings in the polygon 04637 numRings = ( int * )( mGeometry + 1 + sizeof( int ) ); 04638 ptr = mGeometry + 1 + 2 * sizeof( int ); 04639 04640 QVector<GEOSGeometry*> rings; 04641 04642 for ( idx = 0; idx < *numRings; idx++ ) 04643 { 04644 //QgsDebugMsg("Ring nr: "+QString::number(idx)); 04645 04646 QgsPolyline sequence; 04647 04648 // get number of points in the ring 04649 nPoints = ( int * ) ptr; 04650 ptr += 4; 04651 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04652 { 04653 // add points to a point array for drawing the polygon 04654 x = ( double * ) ptr; 04655 ptr += sizeof( double ); 04656 y = ( double * ) ptr; 04657 ptr += sizeof( double ); 04658 if ( hasZValue ) 04659 { 04660 ptr += sizeof( double ); 04661 } 04662 sequence << QgsPoint( *x, *y ); 04663 } 04664 04665 rings << createGeosLinearRing( sequence ); 04666 } 04667 mGeos = createGeosPolygon( rings ); 04668 mDirtyGeos = false; 04669 break; 04670 } 04671 04672 case QGis::WKBMultiPolygon25D: 04673 hasZValue = true; 04674 case QGis::WKBMultiPolygon: 04675 { 04676 QgsDebugMsgLevel( "Multipolygon found", 3 ); 04677 04678 QVector<GEOSGeometry*> polygons; 04679 04680 // get the number of polygons 04681 ptr = mGeometry + 5; 04682 numPolygons = ( int * ) ptr; 04683 ptr = mGeometry + 9; 04684 for ( kdx = 0; kdx < *numPolygons; kdx++ ) 04685 { 04686 //QgsDebugMsg("Polygon nr: "+QString::number(kdx)); 04687 QVector<GEOSGeometry*> rings; 04688 04689 //skip the endian and mGeometry type info and 04690 // get number of rings in the polygon 04691 ptr += 5; 04692 numRings = ( int * ) ptr; 04693 ptr += 4; 04694 for ( idx = 0; idx < *numRings; idx++ ) 04695 { 04696 //QgsDebugMsg("Ring nr: "+QString::number(idx)); 04697 04698 QgsPolyline sequence; 04699 04700 // get number of points in the ring 04701 nPoints = ( int * ) ptr; 04702 ptr += 4; 04703 for ( jdx = 0; jdx < *nPoints; jdx++ ) 04704 { 04705 // add points to a point array for drawing the polygon 04706 x = ( double * ) ptr; 04707 ptr += sizeof( double ); 04708 y = ( double * ) ptr; 04709 ptr += sizeof( double ); 04710 if ( hasZValue ) 04711 { 04712 ptr += sizeof( double ); 04713 } 04714 sequence << QgsPoint( *x, *y ); 04715 } 04716 04717 rings << createGeosLinearRing( sequence ); 04718 } 04719 04720 polygons << createGeosPolygon( rings ); 04721 } 04722 mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons ); 04723 mDirtyGeos = false; 04724 break; 04725 } 04726 04727 default: 04728 return false; 04729 } 04730 } 04731 CATCH_GEOS( false ) 04732 04733 return true; 04734 } 04735 04736 bool QgsGeometry::exportGeosToWkb() 04737 { 04738 //QgsDebugMsg("entered."); 04739 if ( !mDirtyWkb ) 04740 { 04741 // No need to convert again 04742 return true; 04743 } 04744 04745 // clear the WKB, ready to replace with the new one 04746 if ( mGeometry ) 04747 { 04748 delete [] mGeometry; 04749 mGeometry = 0; 04750 } 04751 04752 if ( !mGeos ) 04753 { 04754 // GEOS is null, therefore WKB is null. 04755 mDirtyWkb = false; 04756 return true; 04757 } 04758 04759 // set up byteOrder 04760 char byteOrder = QgsApplication::endian(); 04761 04762 switch ( GEOSGeomTypeId( mGeos ) ) 04763 { 04764 case GEOS_POINT: // a point 04765 { 04766 mGeometrySize = 1 + // sizeof(byte) 04767 4 + // sizeof(uint32) 04768 2 * sizeof( double ); 04769 mGeometry = new unsigned char[mGeometrySize]; 04770 04771 // assign byteOrder 04772 memcpy( mGeometry, &byteOrder, 1 ); 04773 04774 // assign wkbType 04775 int wkbType = QGis::WKBPoint; 04776 memcpy( mGeometry + 1, &wkbType, 4 ); 04777 04778 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos ); 04779 04780 double x, y; 04781 GEOSCoordSeq_getX( cs, 0, &x ); 04782 GEOSCoordSeq_getY( cs, 0, &y ); 04783 04784 memcpy( mGeometry + 5, &x, sizeof( double ) ); 04785 memcpy( mGeometry + 13, &y, sizeof( double ) ); 04786 04787 mDirtyWkb = false; 04788 return true; 04789 } // case GEOS_GEOM::GEOS_POINT 04790 04791 case GEOS_LINESTRING: // a linestring 04792 { 04793 //QgsDebugMsg("Got a geos::GEOS_LINESTRING."); 04794 04795 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( mGeos ); 04796 unsigned int numPoints; 04797 GEOSCoordSeq_getSize( cs, &numPoints ); 04798 04799 // allocate some space for the WKB 04800 mGeometrySize = 1 + // sizeof(byte) 04801 4 + // sizeof(uint32) 04802 4 + // sizeof(uint32) 04803 (( sizeof( double ) + 04804 sizeof( double ) ) * numPoints ); 04805 04806 mGeometry = new unsigned char[mGeometrySize]; 04807 04808 unsigned char* ptr = mGeometry; 04809 04810 // assign byteOrder 04811 memcpy( ptr, &byteOrder, 1 ); 04812 ptr += 1; 04813 04814 // assign wkbType 04815 int wkbType = QGis::WKBLineString; 04816 memcpy( ptr, &wkbType, 4 ); 04817 ptr += 4; 04818 04819 // assign numPoints 04820 memcpy( ptr, &numPoints, 4 ); 04821 ptr += 4; 04822 04823 const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq( mGeos ); 04824 04825 // assign points 04826 for ( unsigned int n = 0; n < numPoints; n++ ) 04827 { 04828 // assign x 04829 GEOSCoordSeq_getX( sequence, n, ( double * )ptr ); 04830 ptr += sizeof( double ); 04831 04832 // assign y 04833 GEOSCoordSeq_getY( sequence, n, ( double * )ptr ); 04834 ptr += sizeof( double ); 04835 } 04836 04837 mDirtyWkb = false; 04838 return true; 04839 04840 // TODO: Deal with endian-ness 04841 } // case GEOS_GEOM::GEOS_LINESTRING 04842 04843 case GEOS_LINEARRING: // a linear ring (linestring with 1st point == last point) 04844 { 04845 // TODO 04846 break; 04847 } // case GEOS_GEOM::GEOS_LINEARRING 04848 04849 case GEOS_POLYGON: // a polygon 04850 { 04851 int geometrySize; 04852 unsigned int nPointsInRing = 0; 04853 04854 //first calculate the geometry size 04855 geometrySize = 1 + 2 * sizeof( int ); //endian, type, number of rings 04856 const GEOSGeometry *theRing = GEOSGetExteriorRing( mGeos ); 04857 if ( theRing ) 04858 { 04859 geometrySize += sizeof( int ); 04860 geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double ); 04861 } 04862 for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); ++i ) 04863 { 04864 geometrySize += sizeof( int ); //number of points in ring 04865 theRing = GEOSGetInteriorRingN( mGeos, i ); 04866 if ( theRing ) 04867 { 04868 geometrySize += getNumGeosPoints( theRing ) * 2 * sizeof( double ); 04869 } 04870 } 04871 04872 mGeometry = new unsigned char[geometrySize]; 04873 mGeometrySize = geometrySize; 04874 04875 //then fill the geometry itself into the wkb 04876 int position = 0; 04877 // assign byteOrder 04878 memcpy( mGeometry, &byteOrder, 1 ); 04879 position += 1; 04880 int wkbtype = QGis::WKBPolygon; 04881 memcpy( &mGeometry[position], &wkbtype, sizeof( int ) ); 04882 position += sizeof( int ); 04883 int nRings = GEOSGetNumInteriorRings( mGeos ) + 1; 04884 memcpy( &mGeometry[position], &nRings, sizeof( int ) ); 04885 position += sizeof( int ); 04886 04887 //exterior ring first 04888 theRing = GEOSGetExteriorRing( mGeos ); 04889 if ( theRing ) 04890 { 04891 nPointsInRing = getNumGeosPoints( theRing ); 04892 memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) ); 04893 position += sizeof( int ); 04894 04895 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing ); 04896 unsigned int n; 04897 GEOSCoordSeq_getSize( cs, &n ); 04898 04899 for ( unsigned int j = 0; j < n; ++j ) 04900 { 04901 GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] ); 04902 position += sizeof( double ); 04903 GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] ); 04904 position += sizeof( double ); 04905 } 04906 } 04907 04908 //interior rings after 04909 for ( int i = 0; i < GEOSGetNumInteriorRings( mGeos ); i++ ) 04910 { 04911 theRing = GEOSGetInteriorRingN( mGeos, i ); 04912 04913 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing ); 04914 GEOSCoordSeq_getSize( cs, &nPointsInRing ); 04915 04916 memcpy( &mGeometry[position], &nPointsInRing, sizeof( int ) ); 04917 position += sizeof( int ); 04918 04919 for ( unsigned int j = 0; j < nPointsInRing; j++ ) 04920 { 04921 GEOSCoordSeq_getX( cs, j, ( double * )&mGeometry[position] ); 04922 position += sizeof( double ); 04923 GEOSCoordSeq_getY( cs, j, ( double * )&mGeometry[position] ); 04924 position += sizeof( double ); 04925 } 04926 } 04927 mDirtyWkb = false; 04928 return true; 04929 } // case GEOS_GEOM::GEOS_POLYGON 04930 break; 04931 04932 case GEOS_MULTIPOINT: // a collection of points 04933 { 04934 // determine size of geometry 04935 int geometrySize = 1 + 2 * sizeof( int ); 04936 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 04937 { 04938 geometrySize += 1 + sizeof( int ) + 2 * sizeof( double ); 04939 } 04940 04941 mGeometry = new unsigned char[geometrySize]; 04942 mGeometrySize = geometrySize; 04943 int wkbPosition = 0; //current position in the byte array 04944 04945 memcpy( mGeometry, &byteOrder, 1 ); 04946 wkbPosition += 1; 04947 int wkbtype = QGis::WKBMultiPoint; 04948 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) ); 04949 wkbPosition += sizeof( int ); 04950 int numPoints = GEOSGetNumGeometries( mGeos ); 04951 memcpy( &mGeometry[wkbPosition], &numPoints, sizeof( int ) ); 04952 wkbPosition += sizeof( int ); 04953 04954 int pointType = QGis::WKBPoint; 04955 const GEOSGeometry *currentPoint = 0; 04956 04957 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 04958 { 04959 //copy endian and point type 04960 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 ); 04961 wkbPosition += 1; 04962 memcpy( &mGeometry[wkbPosition], &pointType, sizeof( int ) ); 04963 wkbPosition += sizeof( int ); 04964 04965 currentPoint = GEOSGetGeometryN( mGeos, i ); 04966 04967 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( currentPoint ); 04968 04969 GEOSCoordSeq_getX( cs, 0, ( double* )&mGeometry[wkbPosition] ); 04970 wkbPosition += sizeof( double ); 04971 GEOSCoordSeq_getY( cs, 0, ( double* )&mGeometry[wkbPosition] ); 04972 wkbPosition += sizeof( double ); 04973 } 04974 mDirtyWkb = false; 04975 return true; 04976 } // case GEOS_GEOM::GEOS_MULTIPOINT 04977 04978 case GEOS_MULTILINESTRING: // a collection of linestrings 04979 { 04980 // determine size of geometry 04981 int geometrySize = 1 + 2 * sizeof( int ); 04982 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 04983 { 04984 geometrySize += 1 + 2 * sizeof( int ); 04985 geometrySize += getNumGeosPoints( GEOSGetGeometryN( mGeos, i ) ) * 2 * sizeof( double ); 04986 } 04987 04988 mGeometry = new unsigned char[geometrySize]; 04989 mGeometrySize = geometrySize; 04990 int wkbPosition = 0; //current position in the byte array 04991 04992 memcpy( mGeometry, &byteOrder, 1 ); 04993 wkbPosition += 1; 04994 int wkbtype = QGis::WKBMultiLineString; 04995 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) ); 04996 wkbPosition += sizeof( int ); 04997 int numLines = GEOSGetNumGeometries( mGeos ); 04998 memcpy( &mGeometry[wkbPosition], &numLines, sizeof( int ) ); 04999 wkbPosition += sizeof( int ); 05000 05001 //loop over lines 05002 int lineType = QGis::WKBLineString; 05003 const GEOSCoordSequence *cs = 0; 05004 unsigned int lineSize; 05005 05006 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 05007 { 05008 //endian and type WKBLineString 05009 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 ); 05010 wkbPosition += 1; 05011 memcpy( &mGeometry[wkbPosition], &lineType, sizeof( int ) ); 05012 wkbPosition += sizeof( int ); 05013 05014 cs = GEOSGeom_getCoordSeq( GEOSGetGeometryN( mGeos, i ) ); 05015 05016 //line size 05017 GEOSCoordSeq_getSize( cs, &lineSize ); 05018 memcpy( &mGeometry[wkbPosition], &lineSize, sizeof( int ) ); 05019 wkbPosition += sizeof( int ); 05020 05021 //vertex coordinates 05022 for ( unsigned int j = 0; j < lineSize; ++j ) 05023 { 05024 GEOSCoordSeq_getX( cs, j, ( double* )&mGeometry[wkbPosition] ); 05025 wkbPosition += sizeof( double ); 05026 GEOSCoordSeq_getY( cs, j, ( double* )&mGeometry[wkbPosition] ); 05027 wkbPosition += sizeof( double ); 05028 } 05029 } 05030 mDirtyWkb = false; 05031 return true; 05032 } // case GEOS_GEOM::GEOS_MULTILINESTRING 05033 05034 case GEOS_MULTIPOLYGON: // a collection of polygons 05035 { 05036 //first determine size of geometry 05037 int geometrySize = 1 + ( 2 * sizeof( int ) ); //endian, type, number of polygons 05038 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 05039 { 05040 const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i ); 05041 geometrySize += 1 + 2 * sizeof( int ); //endian, type, number of rings 05042 //exterior ring 05043 geometrySize += sizeof( int ); //number of points in exterior ring 05044 const GEOSGeometry *exRing = GEOSGetExteriorRing( thePoly ); 05045 geometrySize += 2 * sizeof( double ) * getNumGeosPoints( exRing ); 05046 05047 const GEOSGeometry *intRing = 0; 05048 for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ ) 05049 { 05050 geometrySize += sizeof( int ); //number of points in ring 05051 intRing = GEOSGetInteriorRingN( thePoly, j ); 05052 geometrySize += 2 * sizeof( double ) * getNumGeosPoints( intRing ); 05053 } 05054 } 05055 05056 mGeometry = new unsigned char[geometrySize]; 05057 mGeometrySize = geometrySize; 05058 int wkbPosition = 0; //current position in the byte array 05059 05060 memcpy( mGeometry, &byteOrder, 1 ); 05061 wkbPosition += 1; 05062 int wkbtype = QGis::WKBMultiPolygon; 05063 memcpy( &mGeometry[wkbPosition], &wkbtype, sizeof( int ) ); 05064 wkbPosition += sizeof( int ); 05065 int numPolygons = GEOSGetNumGeometries( mGeos ); 05066 memcpy( &mGeometry[wkbPosition], &numPolygons, sizeof( int ) ); 05067 wkbPosition += sizeof( int ); 05068 05069 //loop over polygons 05070 for ( int i = 0; i < GEOSGetNumGeometries( mGeos ); i++ ) 05071 { 05072 const GEOSGeometry *thePoly = GEOSGetGeometryN( mGeos, i ); 05073 memcpy( &mGeometry[wkbPosition], &byteOrder, 1 ); 05074 wkbPosition += 1; 05075 int polygonType = QGis::WKBPolygon; 05076 memcpy( &mGeometry[wkbPosition], &polygonType, sizeof( int ) ); 05077 wkbPosition += sizeof( int ); 05078 int numRings = GEOSGetNumInteriorRings( thePoly ) + 1; 05079 memcpy( &mGeometry[wkbPosition], &numRings, sizeof( int ) ); 05080 wkbPosition += sizeof( int ); 05081 05082 //exterior ring 05083 const GEOSGeometry *theRing = GEOSGetExteriorRing( thePoly ); 05084 int nPointsInRing = getNumGeosPoints( theRing ); 05085 memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) ); 05086 wkbPosition += sizeof( int ); 05087 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing ); 05088 05089 for ( int k = 0; k < nPointsInRing; ++k ) 05090 { 05091 GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] ); 05092 wkbPosition += sizeof( double ); 05093 GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] ); 05094 wkbPosition += sizeof( double ); 05095 } 05096 05097 //interior rings 05098 for ( int j = 0; j < GEOSGetNumInteriorRings( thePoly ); j++ ) 05099 { 05100 theRing = GEOSGetInteriorRingN( thePoly, j ); 05101 nPointsInRing = getNumGeosPoints( theRing ); 05102 memcpy( &mGeometry[wkbPosition], &nPointsInRing, sizeof( int ) ); 05103 wkbPosition += sizeof( int ); 05104 const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq( theRing ); 05105 05106 for ( int k = 0; k < nPointsInRing; ++k ) 05107 { 05108 GEOSCoordSeq_getX( cs, k, ( double * )&mGeometry[wkbPosition] ); 05109 wkbPosition += sizeof( double ); 05110 GEOSCoordSeq_getY( cs, k, ( double * )&mGeometry[wkbPosition] ); 05111 wkbPosition += sizeof( double ); 05112 } 05113 } 05114 } 05115 mDirtyWkb = false; 05116 return true; 05117 } // case GEOS_GEOM::GEOS_MULTIPOLYGON 05118 05119 case GEOS_GEOMETRYCOLLECTION: // a collection of heterogeneus geometries 05120 { 05121 // TODO 05122 QgsDebugMsg( "geometry collection - not supported" ); 05123 break; 05124 } // case GEOS_GEOM::GEOS_GEOMETRYCOLLECTION 05125 05126 } // switch (mGeos->getGeometryTypeId()) 05127 05128 return false; 05129 } 05130 05131 bool QgsGeometry::convertToMultiType() 05132 { 05133 if ( !mGeometry ) 05134 { 05135 return false; 05136 } 05137 05138 QGis::WkbType geomType = wkbType(); 05139 05140 if ( geomType == QGis::WKBMultiPoint || geomType == QGis::WKBMultiPoint25D || 05141 geomType == QGis::WKBMultiLineString || geomType == QGis::WKBMultiLineString25D || 05142 geomType == QGis::WKBMultiPolygon || geomType == QGis::WKBMultiPolygon25D || geomType == QGis::WKBUnknown ) 05143 { 05144 return false; //no need to convert 05145 } 05146 05147 int newGeomSize = mGeometrySize + 1 + 2 * sizeof( int ); //endian: 1, multitype: sizeof(int), number of geometries: sizeof(int) 05148 unsigned char* newGeometry = new unsigned char[newGeomSize]; 05149 05150 int currentWkbPosition = 0; 05151 //copy endian 05152 char byteOrder = QgsApplication::endian(); 05153 memcpy( &newGeometry[currentWkbPosition], &byteOrder, 1 ); 05154 currentWkbPosition += 1; 05155 05156 //copy wkbtype 05157 //todo 05158 QGis::WkbType newMultiType; 05159 switch ( geomType ) 05160 { 05161 case QGis::WKBPoint: 05162 newMultiType = QGis::WKBMultiPoint; 05163 break; 05164 case QGis::WKBPoint25D: 05165 newMultiType = QGis::WKBMultiPoint25D; 05166 break; 05167 case QGis::WKBLineString: 05168 newMultiType = QGis::WKBMultiLineString; 05169 break; 05170 case QGis::WKBLineString25D: 05171 newMultiType = QGis::WKBMultiLineString25D; 05172 break; 05173 case QGis::WKBPolygon: 05174 newMultiType = QGis::WKBMultiPolygon; 05175 break; 05176 case QGis::WKBPolygon25D: 05177 newMultiType = QGis::WKBMultiPolygon25D; 05178 break; 05179 default: 05180 delete newGeometry; 05181 return false; 05182 } 05183 memcpy( &newGeometry[currentWkbPosition], &newMultiType, sizeof( int ) ); 05184 currentWkbPosition += sizeof( int ); 05185 05186 //copy number of geometries 05187 int nGeometries = 1; 05188 memcpy( &newGeometry[currentWkbPosition], &nGeometries, sizeof( int ) ); 05189 currentWkbPosition += sizeof( int ); 05190 05191 //copy the existing single geometry 05192 memcpy( &newGeometry[currentWkbPosition], mGeometry, mGeometrySize ); 05193 05194 delete [] mGeometry; 05195 mGeometry = newGeometry; 05196 mGeometrySize = newGeomSize; 05197 mDirtyGeos = true; 05198 return true; 05199 } 05200 05201 void QgsGeometry::translateVertex( int& wkbPosition, double dx, double dy, bool hasZValue ) 05202 { 05203 double x, y, translated_x, translated_y; 05204 05205 //x-coordinate 05206 x = *(( double * )( &( mGeometry[wkbPosition] ) ) ); 05207 translated_x = x + dx; 05208 memcpy( &( mGeometry[wkbPosition] ), &translated_x, sizeof( double ) ); 05209 wkbPosition += sizeof( double ); 05210 05211 //y-coordinate 05212 y = *(( double * )( &( mGeometry[wkbPosition] ) ) ); 05213 translated_y = y + dy; 05214 memcpy( &( mGeometry[wkbPosition] ), &translated_y, sizeof( double ) ); 05215 wkbPosition += sizeof( double ); 05216 05217 if ( hasZValue ) 05218 { 05219 wkbPosition += sizeof( double ); 05220 } 05221 } 05222 05223 void QgsGeometry::transformVertex( int& wkbPosition, const QgsCoordinateTransform& ct, bool hasZValue ) 05224 { 05225 double x, y, z; 05226 05227 05228 x = *(( double * )( &( mGeometry[wkbPosition] ) ) ); 05229 y = *(( double * )( &( mGeometry[wkbPosition + sizeof( double )] ) ) ); 05230 z = 0.0; // Ignore Z for now. 05231 05232 ct.transformInPlace( x, y, z ); 05233 05234 // new x-coordinate 05235 memcpy( &( mGeometry[wkbPosition] ), &x, sizeof( double ) ); 05236 wkbPosition += sizeof( double ); 05237 05238 // new y-coordinate 05239 memcpy( &( mGeometry[wkbPosition] ), &y, sizeof( double ) ); 05240 wkbPosition += sizeof( double ); 05241 05242 if ( hasZValue ) 05243 { 05244 wkbPosition += sizeof( double ); 05245 } 05246 } 05247 05248 int QgsGeometry::splitLinearGeometry( GEOSGeometry *splitLine, QList<QgsGeometry*>& newGeometries ) 05249 { 05250 if ( !splitLine ) 05251 { 05252 return 2; 05253 } 05254 05255 if ( !mGeos || mDirtyGeos ) 05256 { 05257 if ( !exportWkbToGeos() ) 05258 return 5; 05259 } 05260 05261 //first test if linestring intersects geometry. If not, return straight away 05262 if ( !GEOSIntersects( splitLine, mGeos ) ) 05263 { 05264 return 1; 05265 } 05266 05267 //first union all the polygon rings together (to get them noded, see JTS developer guide) 05268 GEOSGeometry* nodedGeometry = nodeGeometries( splitLine, mGeos ); 05269 if ( !nodedGeometry ) 05270 { 05271 return 3; //an error occured during noding 05272 } 05273 05274 GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry ); 05275 if ( !mergedLines ) 05276 { 05277 GEOSGeom_destroy( nodedGeometry ); 05278 return 4; 05279 } 05280 05281 QVector<GEOSGeometry*> testedGeometries; 05282 05283 for ( int i = 0; i < GEOSGetNumGeometries( mergedLines ); i++ ) 05284 { 05285 const GEOSGeometry *testing = GEOSGetGeometryN( mergedLines, i ); 05286 if ( lineContainedInLine( testing, mGeos ) == 1 ) 05287 { 05288 testedGeometries << GEOSGeom_clone( testing ); 05289 } 05290 } 05291 05292 mergeGeometriesMultiTypeSplit( testedGeometries ); 05293 05294 if ( testedGeometries.size() > 0 ) 05295 { 05296 GEOSGeom_destroy( mGeos ); 05297 mGeos = testedGeometries[0]; 05298 mDirtyWkb = true; 05299 } 05300 05301 for ( int i = 1; i < testedGeometries.size(); ++i ) 05302 newGeometries << fromGeosGeom( testedGeometries[i] ); 05303 05304 GEOSGeom_destroy( nodedGeometry ); 05305 GEOSGeom_destroy( mergedLines ); 05306 return 0; 05307 } 05308 05309 int QgsGeometry::splitPolygonGeometry( GEOSGeometry* splitLine, QList<QgsGeometry*>& newGeometries ) 05310 { 05311 if ( !splitLine ) 05312 { 05313 return 2; 05314 } 05315 05316 if ( !mGeos || mDirtyGeos ) 05317 { 05318 if ( !exportWkbToGeos() ) 05319 return 5; 05320 } 05321 05322 //first test if linestring intersects geometry. If not, return straight away 05323 if ( !GEOSIntersects( splitLine, mGeos ) ) 05324 { 05325 return 1; 05326 } 05327 05328 //first union all the polygon rings together (to get them noded, see JTS developer guide) 05329 GEOSGeometry *nodedGeometry = nodeGeometries( splitLine, mGeos ); 05330 if ( !nodedGeometry ) 05331 { 05332 return 2; //an error occured during noding 05333 } 05334 05335 #if defined(GEOS_VERSION_MAJOR) && defined(GEOS_VERSION_MINOR) && \ 05336 ((GEOS_VERSION_MAJOR>3) || ((GEOS_VERSION_MAJOR==3) && (GEOS_VERSION_MINOR>=1))) 05337 GEOSGeometry *cutEdges = GEOSPolygonizer_getCutEdges( &nodedGeometry, 1 ); 05338 if ( cutEdges ) 05339 { 05340 if ( numberOfGeometries( cutEdges ) > 0 ) 05341 { 05342 GEOSGeom_destroy( cutEdges ); 05343 GEOSGeom_destroy( nodedGeometry ); 05344 return 3; 05345 } 05346 05347 GEOSGeom_destroy( cutEdges ); 05348 } 05349 #endif 05350 05351 GEOSGeometry *polygons = GEOSPolygonize( &nodedGeometry, 1 ); 05352 if ( !polygons || numberOfGeometries( polygons ) == 0 ) 05353 { 05354 if ( polygons ) 05355 GEOSGeom_destroy( polygons ); 05356 05357 GEOSGeom_destroy( nodedGeometry ); 05358 05359 return 4; 05360 } 05361 05362 GEOSGeom_destroy( nodedGeometry ); 05363 05364 //test every polygon if contained in original geometry 05365 //include in result if yes 05366 QVector<GEOSGeometry*> testedGeometries; 05367 GEOSGeometry *intersectGeometry = 0; 05368 05369 //ratio intersect geometry / geometry. This should be close to 1 05370 //if the polygon belongs to the input geometry 05371 05372 for ( int i = 0; i < numberOfGeometries( polygons ); i++ ) 05373 { 05374 const GEOSGeometry *polygon = GEOSGetGeometryN( polygons, i ); 05375 intersectGeometry = GEOSIntersection( mGeos, polygon ); 05376 if ( !intersectGeometry ) 05377 { 05378 QgsDebugMsg( "intersectGeometry is NULL" ); 05379 continue; 05380 } 05381 05382 double intersectionArea; 05383 GEOSArea( intersectGeometry, &intersectionArea ); 05384 05385 double polygonArea; 05386 GEOSArea( polygon, &polygonArea ); 05387 05388 const double areaRatio = intersectionArea / polygonArea; 05389 if ( areaRatio > 0.99 && areaRatio < 1.01 ) 05390 testedGeometries << GEOSGeom_clone( polygon ); 05391 05392 GEOSGeom_destroy( intersectGeometry ); 05393 } 05394 05395 bool splitDone = true; 05396 int nGeometriesThis = numberOfGeometries( mGeos ); //original number of geometries 05397 if ( testedGeometries.size() == nGeometriesThis ) 05398 { 05399 splitDone = false; 05400 } 05401 05402 mergeGeometriesMultiTypeSplit( testedGeometries ); 05403 05404 //no split done, preserve original geometry 05405 if ( !splitDone ) 05406 { 05407 for ( int i = 0; i < testedGeometries.size(); ++i ) 05408 { 05409 GEOSGeom_destroy( testedGeometries[i] ); 05410 } 05411 return 1; 05412 } 05413 else if ( testedGeometries.size() > 0 ) //split successfull 05414 { 05415 GEOSGeom_destroy( mGeos ); 05416 mGeos = testedGeometries[0]; 05417 mDirtyWkb = true; 05418 } 05419 05420 for ( int i = 1; i < testedGeometries.size(); ++i ) 05421 { 05422 newGeometries << fromGeosGeom( testedGeometries[i] ); 05423 } 05424 05425 GEOSGeom_destroy( polygons ); 05426 return 0; 05427 } 05428 05429 GEOSGeometry* QgsGeometry::reshapePolygon( const GEOSGeometry* polygon, const GEOSGeometry* reshapeLineGeos ) 05430 { 05431 //go through outer shell and all inner rings and check if there is exactly one intersection of a ring and the reshape line 05432 int nIntersections = 0; 05433 int lastIntersectingRing = -2; 05434 const GEOSGeometry* lastIntersectingGeom = 0; 05435 05436 int nRings = GEOSGetNumInteriorRings( polygon ); 05437 if ( nRings < 0 ) 05438 { 05439 return 0; 05440 } 05441 05442 //does outer ring intersect? 05443 const GEOSGeometry* outerRing = GEOSGetExteriorRing( polygon ); 05444 if ( GEOSIntersects( outerRing, reshapeLineGeos ) == 1 ) 05445 { 05446 ++nIntersections; 05447 lastIntersectingRing = -1; 05448 lastIntersectingGeom = outerRing; 05449 } 05450 05451 //do inner rings intersect? 05452 const GEOSGeometry **innerRings = new const GEOSGeometry*[nRings]; 05453 05454 try 05455 { 05456 for ( int i = 0; i < nRings; ++i ) 05457 { 05458 innerRings[i] = GEOSGetInteriorRingN( polygon, i ); 05459 if ( GEOSIntersects( innerRings[i], reshapeLineGeos ) == 1 ) 05460 { 05461 ++nIntersections; 05462 lastIntersectingRing = i; 05463 lastIntersectingGeom = innerRings[i]; 05464 } 05465 } 05466 } 05467 catch ( GEOSException &e ) 05468 { 05469 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 05470 nIntersections = 0; 05471 } 05472 05473 if ( nIntersections != 1 ) //reshape line is only allowed to intersect one ring 05474 { 05475 delete [] innerRings; 05476 return 0; 05477 } 05478 05479 //we have one intersecting ring, let's try to reshape it 05480 GEOSGeometry* reshapeResult = reshapeLine( lastIntersectingGeom, reshapeLineGeos ); 05481 if ( !reshapeResult ) 05482 { 05483 delete [] innerRings; 05484 return 0; 05485 } 05486 05487 //if reshaping took place, we need to reassemble the polygon and its rings 05488 GEOSGeometry* newRing = 0; 05489 const GEOSCoordSequence* reshapeSequence = GEOSGeom_getCoordSeq( reshapeResult ); 05490 GEOSCoordSequence* newCoordSequence = GEOSCoordSeq_clone( reshapeSequence ); 05491 05492 GEOSGeom_destroy( reshapeResult ); 05493 05494 newRing = GEOSGeom_createLinearRing( newCoordSequence ); 05495 if ( !newRing ) 05496 { 05497 delete [] innerRings; 05498 return 0; 05499 } 05500 05501 05502 GEOSGeometry* newOuterRing = 0; 05503 if ( lastIntersectingRing == -1 ) 05504 { 05505 newOuterRing = newRing; 05506 } 05507 else 05508 { 05509 newOuterRing = GEOSGeom_clone( outerRing ); 05510 } 05511 05512 //check if all the rings are still inside the outer boundary 05513 QList<GEOSGeometry*> ringList; 05514 if ( nRings > 0 ) 05515 { 05516 GEOSGeometry* outerRingPoly = GEOSGeom_createPolygon( GEOSGeom_clone( newOuterRing ), 0, 0 ); 05517 if ( outerRingPoly ) 05518 { 05519 GEOSGeometry* currentRing = 0; 05520 for ( int i = 0; i < nRings; ++i ) 05521 { 05522 if ( lastIntersectingRing == i ) 05523 { 05524 currentRing = newRing; 05525 } 05526 else 05527 { 05528 currentRing = GEOSGeom_clone( innerRings[i] ); 05529 } 05530 05531 //possibly a ring is no longer contained in the result polygon after reshape 05532 if ( GEOSContains( outerRingPoly, currentRing ) == 1 ) 05533 { 05534 ringList.push_back( currentRing ); 05535 } 05536 else 05537 { 05538 GEOSGeom_destroy( currentRing ); 05539 } 05540 } 05541 } 05542 GEOSGeom_destroy( outerRingPoly ); 05543 } 05544 05545 GEOSGeometry** newInnerRings = new GEOSGeometry*[ringList.size()]; 05546 for ( int i = 0; i < ringList.size(); ++i ) 05547 { 05548 newInnerRings[i] = ringList.at( i ); 05549 } 05550 05551 delete [] innerRings; 05552 05553 GEOSGeometry* reshapedPolygon = GEOSGeom_createPolygon( newOuterRing, newInnerRings, ringList.size() ); 05554 delete[] newInnerRings; 05555 if ( !reshapedPolygon ) 05556 { 05557 return 0; 05558 } 05559 return reshapedPolygon; 05560 } 05561 05562 GEOSGeometry* QgsGeometry::reshapeLine( const GEOSGeometry* line, const GEOSGeometry* reshapeLineGeos ) 05563 { 05564 if ( !line || !reshapeLineGeos ) 05565 { 05566 return 0; 05567 } 05568 05569 bool atLeastTwoIntersections; 05570 05571 try 05572 { 05573 //make sure there are at least two intersection between line and reshape geometry 05574 GEOSGeometry* intersectGeom = GEOSIntersection( line, reshapeLineGeos ); 05575 if ( intersectGeom ) 05576 { 05577 atLeastTwoIntersections = ( GEOSGeomTypeId( intersectGeom ) == GEOS_MULTIPOINT && GEOSGetNumGeometries( intersectGeom ) > 1 ); 05578 GEOSGeom_destroy( intersectGeom ); 05579 } 05580 } 05581 catch ( GEOSException &e ) 05582 { 05583 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 05584 atLeastTwoIntersections = false; 05585 } 05586 05587 if ( !atLeastTwoIntersections ) 05588 { 05589 return 0; 05590 } 05591 05592 //begin and end point of original line 05593 const GEOSCoordSequence* lineCoordSeq = GEOSGeom_getCoordSeq( line ); 05594 if ( !lineCoordSeq ) 05595 { 05596 return 0; 05597 } 05598 unsigned int lineCoordSeqSize; 05599 if ( GEOSCoordSeq_getSize( lineCoordSeq, &lineCoordSeqSize ) == 0 ) 05600 { 05601 return 0; 05602 } 05603 if ( lineCoordSeqSize < 2 ) 05604 { 05605 return 0; 05606 } 05607 //first and last vertex of line 05608 double x1, y1, x2, y2; 05609 GEOSCoordSeq_getX( lineCoordSeq, 0, &x1 ); 05610 GEOSCoordSeq_getY( lineCoordSeq, 0, &y1 ); 05611 GEOSCoordSeq_getX( lineCoordSeq, lineCoordSeqSize - 1, &x2 ); 05612 GEOSCoordSeq_getY( lineCoordSeq, lineCoordSeqSize - 1, &y2 ); 05613 GEOSGeometry* beginLineVertex = createGeosPoint( QgsPoint( x1, y1 ) ); 05614 GEOSGeometry* endLineVertex = createGeosPoint( QgsPoint( x2, y2 ) ); 05615 05616 bool isRing = false; 05617 if ( GEOSGeomTypeId( line ) == GEOS_LINEARRING || GEOSEquals( beginLineVertex, endLineVertex ) == 1 ) 05618 { 05619 isRing = true; 05620 } 05621 05622 //node line and reshape line 05623 GEOSGeometry* nodedGeometry = nodeGeometries( reshapeLineGeos, line ); 05624 if ( !nodedGeometry ) 05625 { 05626 GEOSGeom_destroy( beginLineVertex ); 05627 GEOSGeom_destroy( endLineVertex ); 05628 return 0; 05629 } 05630 05631 //and merge them together 05632 GEOSGeometry *mergedLines = GEOSLineMerge( nodedGeometry ); 05633 GEOSGeom_destroy( nodedGeometry ); 05634 if ( !mergedLines ) 05635 { 05636 GEOSGeom_destroy( beginLineVertex ); 05637 GEOSGeom_destroy( endLineVertex ); 05638 return 0; 05639 } 05640 05641 int numMergedLines = GEOSGetNumGeometries( mergedLines ); 05642 if ( numMergedLines < 2 ) //some special cases. Normally it is >2 05643 { 05644 GEOSGeom_destroy( beginLineVertex ); 05645 GEOSGeom_destroy( endLineVertex ); 05646 if ( numMergedLines == 1 ) //reshape line is from begin to endpoint. So we keep the reshapeline 05647 { 05648 return GEOSGeom_clone( reshapeLineGeos ); 05649 } 05650 else 05651 { 05652 return 0; 05653 } 05654 } 05655 05656 QList<GEOSGeometry*> resultLineParts; //collection with the line segments that will be contained in result 05657 QList<GEOSGeometry*> probableParts; //parts where we can decide on inclusion only after going through all the candidates 05658 05659 for ( int i = 0; i < numMergedLines; ++i ) 05660 { 05661 const GEOSGeometry* currentGeom; 05662 05663 currentGeom = GEOSGetGeometryN( mergedLines, i ); 05664 const GEOSCoordSequence* currentCoordSeq = GEOSGeom_getCoordSeq( currentGeom ); 05665 unsigned int currentCoordSeqSize; 05666 GEOSCoordSeq_getSize( currentCoordSeq, ¤tCoordSeqSize ); 05667 if ( currentCoordSeqSize < 2 ) 05668 { 05669 continue; 05670 } 05671 05672 //get the two endpoints of the current line merge result 05673 double xBegin, xEnd, yBegin, yEnd; 05674 GEOSCoordSeq_getX( currentCoordSeq, 0, &xBegin ); 05675 GEOSCoordSeq_getY( currentCoordSeq, 0, &yBegin ); 05676 GEOSCoordSeq_getX( currentCoordSeq, currentCoordSeqSize - 1, &xEnd ); 05677 GEOSCoordSeq_getY( currentCoordSeq, currentCoordSeqSize - 1, &yEnd ); 05678 GEOSGeometry* beginCurrentGeomVertex = createGeosPoint( QgsPoint( xBegin, yBegin ) ); 05679 GEOSGeometry* endCurrentGeomVertex = createGeosPoint( QgsPoint( xEnd, yEnd ) ); 05680 05681 //check how many endpoints of the line merge result are on the (original) line 05682 int nEndpointsOnOriginalLine = 0; 05683 if ( pointContainedInLine( beginCurrentGeomVertex, line ) == 1 ) 05684 { 05685 nEndpointsOnOriginalLine += 1; 05686 } 05687 05688 if ( pointContainedInLine( endCurrentGeomVertex, line ) == 1 ) 05689 { 05690 nEndpointsOnOriginalLine += 1; 05691 } 05692 05693 //check how many endpoints equal the endpoints of the original line 05694 int nEndpointsSameAsOriginalLine = 0; 05695 if ( GEOSEquals( beginCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( beginCurrentGeomVertex, endLineVertex ) == 1 ) 05696 { 05697 nEndpointsSameAsOriginalLine += 1; 05698 } 05699 if ( GEOSEquals( endCurrentGeomVertex, beginLineVertex ) == 1 || GEOSEquals( endCurrentGeomVertex, endLineVertex ) == 1 ) 05700 { 05701 nEndpointsSameAsOriginalLine += 1; 05702 } 05703 05704 //check if the current geometry overlaps the original geometry (GEOSOverlap does not seem to work with linestrings) 05705 bool currentGeomOverlapsOriginalGeom = false; 05706 bool currentGeomOverlapsReshapeLine = false; 05707 if ( QgsGeometry::lineContainedInLine( currentGeom, line ) == 1 ) 05708 { 05709 currentGeomOverlapsOriginalGeom = true; 05710 } 05711 if ( QgsGeometry::lineContainedInLine( currentGeom, reshapeLineGeos ) == 1 ) 05712 { 05713 currentGeomOverlapsReshapeLine = true; 05714 } 05715 05716 05717 //logic to decide if this part belongs to the result 05718 if ( nEndpointsSameAsOriginalLine == 1 && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom ) 05719 { 05720 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) ); 05721 } 05722 //for closed rings, we take one segment from the candidate list 05723 else if ( isRing && nEndpointsOnOriginalLine == 2 && currentGeomOverlapsOriginalGeom ) 05724 { 05725 probableParts.push_back( GEOSGeom_clone( currentGeom ) ); 05726 } 05727 else if ( nEndpointsOnOriginalLine == 2 && !currentGeomOverlapsOriginalGeom ) 05728 { 05729 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) ); 05730 } 05731 else if ( nEndpointsSameAsOriginalLine == 2 && !currentGeomOverlapsOriginalGeom ) 05732 { 05733 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) ); 05734 } 05735 else if ( currentGeomOverlapsOriginalGeom && currentGeomOverlapsReshapeLine ) 05736 { 05737 resultLineParts.push_back( GEOSGeom_clone( currentGeom ) ); 05738 } 05739 05740 GEOSGeom_destroy( beginCurrentGeomVertex ); 05741 GEOSGeom_destroy( endCurrentGeomVertex ); 05742 } 05743 05744 //add the longest segment from the probable list for rings (only used for polygon rings) 05745 if ( isRing && probableParts.size() > 0 ) 05746 { 05747 GEOSGeometry* maxGeom = 0; //the longest geometry in the probabla list 05748 GEOSGeometry* currentGeom = 0; 05749 double maxLength = -DBL_MAX; 05750 double currentLength = 0; 05751 for ( int i = 0; i < probableParts.size(); ++i ) 05752 { 05753 currentGeom = probableParts.at( i ); 05754 GEOSLength( currentGeom, ¤tLength ); 05755 if ( currentLength > maxLength ) 05756 { 05757 maxLength = currentLength; 05758 GEOSGeom_destroy( maxGeom ); 05759 maxGeom = currentGeom; 05760 } 05761 else 05762 { 05763 GEOSGeom_destroy( currentGeom ); 05764 } 05765 } 05766 resultLineParts.push_back( maxGeom ); 05767 } 05768 05769 GEOSGeom_destroy( beginLineVertex ); 05770 GEOSGeom_destroy( endLineVertex ); 05771 GEOSGeom_destroy( mergedLines ); 05772 05773 GEOSGeometry* result = 0; 05774 if ( resultLineParts.size() < 1 ) 05775 { 05776 return 0; 05777 } 05778 if ( resultLineParts.size() == 1 ) //the whole result was reshaped 05779 { 05780 result = resultLineParts[0]; 05781 } 05782 else //>1 05783 { 05784 GEOSGeometry **lineArray = new GEOSGeometry*[resultLineParts.size()]; 05785 for ( int i = 0; i < resultLineParts.size(); ++i ) 05786 { 05787 lineArray[i] = resultLineParts[i]; 05788 } 05789 05790 //create multiline from resultLineParts 05791 GEOSGeometry* multiLineGeom = GEOSGeom_createCollection( GEOS_MULTILINESTRING, lineArray, resultLineParts.size() ); 05792 delete [] lineArray; 05793 05794 //then do a linemerge with the newly combined partstrings 05795 result = GEOSLineMerge( multiLineGeom ); 05796 GEOSGeom_destroy( multiLineGeom ); 05797 } 05798 05799 //now test if the result is a linestring. Otherwise something went wrong 05800 if ( GEOSGeomTypeId( result ) != GEOS_LINESTRING ) 05801 { 05802 GEOSGeom_destroy( result ); 05803 return 0; 05804 } 05805 return result; 05806 } 05807 05808 int QgsGeometry::topologicalTestPointsSplit( const GEOSGeometry* splitLine, QList<QgsPoint>& testPoints ) const 05809 { 05810 //Find out the intersection points between splitLineGeos and this geometry. 05811 //These points need to be tested for topological correctness by the calling function 05812 //if topological editing is enabled 05813 05814 testPoints.clear(); 05815 GEOSGeometry* intersectionGeom = GEOSIntersection( mGeos, splitLine ); 05816 if ( !intersectionGeom ) 05817 { 05818 return 1; 05819 } 05820 05821 bool simple = false; 05822 int nIntersectGeoms = 1; 05823 if ( GEOSGeomTypeId( intersectionGeom ) == GEOS_LINESTRING || GEOSGeomTypeId( intersectionGeom ) == GEOS_POINT ) 05824 { 05825 simple = true; 05826 } 05827 05828 if ( !simple ) 05829 { 05830 nIntersectGeoms = GEOSGetNumGeometries( intersectionGeom ); 05831 } 05832 05833 for ( int i = 0; i < nIntersectGeoms; ++i ) 05834 { 05835 const GEOSGeometry* currentIntersectGeom; 05836 if ( simple ) 05837 { 05838 currentIntersectGeom = intersectionGeom; 05839 } 05840 else 05841 { 05842 currentIntersectGeom = GEOSGetGeometryN( intersectionGeom, i ); 05843 } 05844 05845 const GEOSCoordSequence* lineSequence = GEOSGeom_getCoordSeq( currentIntersectGeom ); 05846 unsigned int sequenceSize = 0; 05847 double x, y; 05848 if ( GEOSCoordSeq_getSize( lineSequence, &sequenceSize ) != 0 ) 05849 { 05850 for ( unsigned int i = 0; i < sequenceSize; ++i ) 05851 { 05852 if ( GEOSCoordSeq_getX( lineSequence, i, &x ) != 0 ) 05853 { 05854 if ( GEOSCoordSeq_getY( lineSequence, i, &y ) != 0 ) 05855 { 05856 testPoints.push_back( QgsPoint( x, y ) ); 05857 } 05858 } 05859 } 05860 } 05861 } 05862 GEOSGeom_destroy( intersectionGeom ); 05863 return 0; 05864 } 05865 05866 GEOSGeometry *QgsGeometry::nodeGeometries( const GEOSGeometry *splitLine, const GEOSGeometry *geom ) 05867 { 05868 if ( !splitLine || !geom ) 05869 { 05870 return 0; 05871 } 05872 05873 GEOSGeometry *geometryBoundary = 0; 05874 if ( GEOSGeomTypeId( geom ) == GEOS_POLYGON || GEOSGeomTypeId( geom ) == GEOS_MULTIPOLYGON ) 05875 { 05876 geometryBoundary = GEOSBoundary( geom ); 05877 } 05878 else 05879 { 05880 geometryBoundary = GEOSGeom_clone( geom ); 05881 } 05882 05883 GEOSGeometry *splitLineClone = GEOSGeom_clone( splitLine ); 05884 GEOSGeometry *unionGeometry = GEOSUnion( splitLineClone, geometryBoundary ); 05885 GEOSGeom_destroy( splitLineClone ); 05886 05887 GEOSGeom_destroy( geometryBoundary ); 05888 return unionGeometry; 05889 } 05890 05891 int QgsGeometry::lineContainedInLine( const GEOSGeometry* line1, const GEOSGeometry* line2 ) 05892 { 05893 if ( !line1 || !line2 ) 05894 { 05895 return -1; 05896 } 05897 05898 05899 double bufferDistance = 0.00001; 05900 if ( geomInDegrees( line2 ) ) //use more accurate tolerance for degrees 05901 { 05902 bufferDistance = 0.00000001; 05903 } 05904 GEOSGeometry* bufferGeom = GEOSBuffer( line2, bufferDistance, DEFAULT_QUADRANT_SEGMENTS ); 05905 if ( !bufferGeom ) 05906 { 05907 return -2; 05908 } 05909 05910 GEOSGeometry* intersectionGeom = GEOSIntersection( bufferGeom, line1 ); 05911 05912 //compare ratio between line1Length and intersectGeomLength (usually close to 1 if line1 is contained in line2) 05913 double intersectGeomLength; 05914 double line1Length; 05915 05916 GEOSLength( intersectionGeom, &intersectGeomLength ); 05917 GEOSLength( line1, &line1Length ); 05918 05919 GEOSGeom_destroy( bufferGeom ); 05920 GEOSGeom_destroy( intersectionGeom ); 05921 05922 double intersectRatio = line1Length / intersectGeomLength; 05923 if ( intersectRatio > 0.9 && intersectRatio < 1.1 ) 05924 { 05925 return 1; 05926 } 05927 return 0; 05928 } 05929 05930 int QgsGeometry::pointContainedInLine( const GEOSGeometry* point, const GEOSGeometry* line ) 05931 { 05932 if ( !point || !line ) 05933 { 05934 return -1; 05935 } 05936 05937 double bufferDistance = 0.000001; 05938 if ( geomInDegrees( line ) ) 05939 { 05940 bufferDistance = 0.00000001; 05941 } 05942 GEOSGeometry* lineBuffer = GEOSBuffer( line, bufferDistance, 8 ); 05943 if ( !lineBuffer ) 05944 { 05945 return -2; 05946 } 05947 05948 bool contained = false; 05949 if ( GEOSContains( lineBuffer, point ) == 1 ) 05950 { 05951 contained = true; 05952 } 05953 05954 GEOSGeom_destroy( lineBuffer ); 05955 return contained; 05956 } 05957 05958 bool QgsGeometry::geomInDegrees( const GEOSGeometry* geom ) 05959 { 05960 GEOSGeometry* bbox = GEOSEnvelope( geom ); 05961 if ( !bbox ) 05962 { 05963 return false; 05964 } 05965 05966 const GEOSGeometry* bBoxRing = GEOSGetExteriorRing( bbox ); 05967 if ( !bBoxRing ) 05968 { 05969 return false; 05970 } 05971 const GEOSCoordSequence* bBoxCoordSeq = GEOSGeom_getCoordSeq( bBoxRing ); 05972 05973 if ( !bBoxCoordSeq ) 05974 { 05975 return false; 05976 } 05977 05978 unsigned int nCoords = 0; 05979 if ( !GEOSCoordSeq_getSize( bBoxCoordSeq, &nCoords ) ) 05980 { 05981 return false; 05982 } 05983 05984 double x, y; 05985 for ( unsigned int i = 0; i < ( nCoords - 1 ); ++i ) 05986 { 05987 GEOSCoordSeq_getX( bBoxCoordSeq, i, &x ); 05988 if ( x > 180 || x < -180 ) 05989 { 05990 return false; 05991 } 05992 GEOSCoordSeq_getY( bBoxCoordSeq, i, &y ); 05993 if ( y > 90 || y < -90 ) 05994 { 05995 return false; 05996 } 05997 } 05998 05999 return true; 06000 } 06001 06002 int QgsGeometry::numberOfGeometries( GEOSGeometry* g ) const 06003 { 06004 if ( !g ) 06005 { 06006 return 0; 06007 } 06008 int geometryType = GEOSGeomTypeId( g ); 06009 if ( geometryType == GEOS_POINT || geometryType == GEOS_LINESTRING || geometryType == GEOS_LINEARRING 06010 || geometryType == GEOS_POLYGON ) 06011 { 06012 return 1; 06013 } 06014 06015 //calling GEOSGetNumGeometries is save for multi types and collections also in geos2 06016 return GEOSGetNumGeometries( g ); 06017 } 06018 06019 int QgsGeometry::mergeGeometriesMultiTypeSplit( QVector<GEOSGeometry*>& splitResult ) 06020 { 06021 if ( !mGeos || mDirtyGeos ) 06022 { 06023 if ( !exportWkbToGeos() ) 06024 return 1; 06025 } 06026 06027 //convert mGeos to geometry collection 06028 int type = GEOSGeomTypeId( mGeos ); 06029 if ( type != GEOS_GEOMETRYCOLLECTION && 06030 type != GEOS_MULTILINESTRING && 06031 type != GEOS_MULTIPOLYGON && 06032 type != GEOS_MULTIPOINT ) 06033 { 06034 return 0; 06035 } 06036 06037 QVector<GEOSGeometry*> copyList = splitResult; 06038 splitResult.clear(); 06039 06040 //collect all the geometries that belong to the initial multifeature 06041 QVector<GEOSGeometry*> unionGeom; 06042 06043 for ( int i = 0; i < copyList.size(); ++i ) 06044 { 06045 //is this geometry a part of the original multitype? 06046 bool isPart = false; 06047 for ( int j = 0; j < GEOSGetNumGeometries( mGeos ); j++ ) 06048 { 06049 if ( GEOSEquals( copyList[i], GEOSGetGeometryN( mGeos, j ) ) ) 06050 { 06051 isPart = true; 06052 break; 06053 } 06054 } 06055 06056 if ( isPart ) 06057 { 06058 unionGeom << copyList[i]; 06059 } 06060 else 06061 { 06062 QVector<GEOSGeometry*> geomVector; 06063 geomVector << copyList[i]; 06064 06065 if ( type == GEOS_MULTILINESTRING ) 06066 { 06067 splitResult << createGeosCollection( GEOS_MULTILINESTRING, geomVector ); 06068 } 06069 else if ( type == GEOS_MULTIPOLYGON ) 06070 { 06071 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, geomVector ); 06072 } 06073 else 06074 { 06075 GEOSGeom_destroy( copyList[i] ); 06076 } 06077 } 06078 } 06079 06080 //make multifeature out of unionGeom 06081 if ( unionGeom.size() > 0 ) 06082 { 06083 if ( type == GEOS_MULTILINESTRING ) 06084 { 06085 splitResult << createGeosCollection( GEOS_MULTILINESTRING, unionGeom ); 06086 } 06087 else if ( type == GEOS_MULTIPOLYGON ) 06088 { 06089 splitResult << createGeosCollection( GEOS_MULTIPOLYGON, unionGeom ); 06090 } 06091 } 06092 else 06093 { 06094 unionGeom.clear(); 06095 } 06096 06097 return 0; 06098 } 06099 06100 QgsPoint QgsGeometry::asPoint( unsigned char*& ptr, bool hasZValue ) 06101 { 06102 ptr += 5; 06103 double* x = ( double * )( ptr ); 06104 double* y = ( double * )( ptr + sizeof( double ) ); 06105 ptr += 2 * sizeof( double ); 06106 06107 if ( hasZValue ) 06108 ptr += sizeof( double ); 06109 06110 return QgsPoint( *x, *y ); 06111 } 06112 06113 06114 QgsPolyline QgsGeometry::asPolyline( unsigned char*& ptr, bool hasZValue ) 06115 { 06116 double x, y; 06117 ptr += 5; 06118 unsigned int nPoints = *(( int* )ptr ); 06119 ptr += 4; 06120 06121 QgsPolyline line( nPoints ); 06122 06123 // Extract the points from the WKB format into the x and y vectors. 06124 for ( uint i = 0; i < nPoints; ++i ) 06125 { 06126 x = *(( double * ) ptr ); 06127 y = *(( double * )( ptr + sizeof( double ) ) ); 06128 06129 ptr += 2 * sizeof( double ); 06130 06131 line[i] = QgsPoint( x, y ); 06132 06133 if ( hasZValue ) // ignore Z value 06134 ptr += sizeof( double ); 06135 } 06136 06137 return line; 06138 } 06139 06140 06141 QgsPolygon QgsGeometry::asPolygon( unsigned char*& ptr, bool hasZValue ) 06142 { 06143 double x, y; 06144 06145 ptr += 5; 06146 06147 // get number of rings in the polygon 06148 unsigned int numRings = *(( int* )ptr ); 06149 ptr += 4; 06150 06151 if ( numRings == 0 ) // sanity check for zero rings in polygon 06152 return QgsPolygon(); 06153 06154 QgsPolygon rings( numRings ); 06155 06156 for ( uint idx = 0; idx < numRings; idx++ ) 06157 { 06158 uint nPoints = *(( int* )ptr ); 06159 ptr += 4; 06160 06161 QgsPolyline ring( nPoints ); 06162 06163 for ( uint jdx = 0; jdx < nPoints; jdx++ ) 06164 { 06165 x = *(( double * ) ptr ); 06166 y = *(( double * )( ptr + sizeof( double ) ) ); 06167 06168 ptr += 2 * sizeof( double ); 06169 06170 if ( hasZValue ) 06171 ptr += sizeof( double ); 06172 06173 ring[jdx] = QgsPoint( x, y ); 06174 } 06175 06176 rings[idx] = ring; 06177 } 06178 06179 return rings; 06180 } 06181 06182 06183 QgsPoint QgsGeometry::asPoint() 06184 { 06185 QGis::WkbType type = wkbType(); 06186 if ( type != QGis::WKBPoint && type != QGis::WKBPoint25D ) 06187 return QgsPoint( 0, 0 ); 06188 06189 unsigned char* ptr = mGeometry; 06190 return asPoint( ptr, type == QGis::WKBPoint25D ); 06191 } 06192 06193 QgsPolyline QgsGeometry::asPolyline() 06194 { 06195 QGis::WkbType type = wkbType(); 06196 if ( type != QGis::WKBLineString && type != QGis::WKBLineString25D ) 06197 return QgsPolyline(); 06198 06199 unsigned char *ptr = mGeometry; 06200 return asPolyline( ptr, type == QGis::WKBLineString25D ); 06201 } 06202 06203 QgsPolygon QgsGeometry::asPolygon() 06204 { 06205 QGis::WkbType type = wkbType(); 06206 if ( type != QGis::WKBPolygon && type != QGis::WKBPolygon25D ) 06207 return QgsPolygon(); 06208 06209 unsigned char *ptr = mGeometry; 06210 return asPolygon( ptr, type == QGis::WKBPolygon25D ); 06211 } 06212 06213 QgsMultiPoint QgsGeometry::asMultiPoint() 06214 { 06215 QGis::WkbType type = wkbType(); 06216 if ( type != QGis::WKBMultiPoint && type != QGis::WKBMultiPoint25D ) 06217 return QgsMultiPoint(); 06218 06219 bool hasZValue = ( type == QGis::WKBMultiPoint25D ); 06220 06221 unsigned char* ptr = mGeometry + 5; 06222 unsigned int nPoints = *(( int* )ptr ); 06223 ptr += 4; 06224 06225 QgsMultiPoint points( nPoints ); 06226 for ( uint i = 0; i < nPoints; i++ ) 06227 { 06228 points[i] = asPoint( ptr, hasZValue ); 06229 } 06230 06231 return points; 06232 } 06233 06234 QgsMultiPolyline QgsGeometry::asMultiPolyline() 06235 { 06236 QGis::WkbType type = wkbType(); 06237 if ( type != QGis::WKBMultiLineString && type != QGis::WKBMultiLineString25D ) 06238 return QgsMultiPolyline(); 06239 06240 bool hasZValue = ( type == QGis::WKBMultiLineString25D ); 06241 06242 unsigned char* ptr = mGeometry + 5; 06243 unsigned int numLineStrings = *(( int* )ptr ); 06244 ptr += 4; 06245 06246 QgsMultiPolyline lines( numLineStrings ); 06247 06248 for ( uint i = 0; i < numLineStrings; i++ ) 06249 { 06250 lines[i] = asPolyline( ptr, hasZValue ); 06251 } 06252 06253 return lines; 06254 } 06255 06256 QgsMultiPolygon QgsGeometry::asMultiPolygon() 06257 { 06258 QGis::WkbType type = wkbType(); 06259 if ( type != QGis::WKBMultiPolygon && type != QGis::WKBMultiPolygon25D ) 06260 return QgsMultiPolygon(); 06261 06262 bool hasZValue = ( type == QGis::WKBMultiPolygon25D ); 06263 06264 unsigned char* ptr = mGeometry + 5; 06265 unsigned int numPolygons = *(( int* )ptr ); 06266 ptr += 4; 06267 06268 QgsMultiPolygon polygons( numPolygons ); 06269 06270 for ( uint i = 0; i < numPolygons; i++ ) 06271 { 06272 polygons[i] = asPolygon( ptr, hasZValue ); 06273 } 06274 06275 return polygons; 06276 } 06277 06278 double QgsGeometry::area() 06279 { 06280 if ( !mGeos ) 06281 { 06282 exportWkbToGeos(); 06283 } 06284 06285 double area; 06286 06287 try 06288 { 06289 if ( GEOSArea( mGeos, &area ) == 0 ) 06290 return -1.0; 06291 } 06292 CATCH_GEOS( -1.0 ) 06293 06294 return area; 06295 } 06296 06297 double QgsGeometry::length() 06298 { 06299 if ( !mGeos ) 06300 { 06301 exportWkbToGeos(); 06302 } 06303 06304 double length; 06305 06306 try 06307 { 06308 if ( GEOSLength( mGeos, &length ) == 0 ) 06309 return -1.0; 06310 } 06311 CATCH_GEOS( -1.0 ) 06312 06313 return length; 06314 } 06315 double QgsGeometry::distance( QgsGeometry& geom ) 06316 { 06317 if ( !mGeos ) 06318 { 06319 exportWkbToGeos(); 06320 } 06321 06322 if ( !geom.mGeos ) 06323 { 06324 geom.exportWkbToGeos(); 06325 } 06326 06327 if ( !mGeos || !geom.mGeos ) 06328 return -1.0; 06329 06330 double dist = -1.0; 06331 06332 try 06333 { 06334 GEOSDistance( mGeos, geom.mGeos, &dist ); 06335 } 06336 CATCH_GEOS( -1.0 ) 06337 06338 return dist; 06339 } 06340 06341 06342 QgsGeometry* QgsGeometry::buffer( double distance, int segments ) 06343 { 06344 if ( !mGeos ) 06345 { 06346 exportWkbToGeos(); 06347 } 06348 if ( !mGeos ) 06349 { 06350 return 0; 06351 } 06352 06353 try 06354 { 06355 return fromGeosGeom( GEOSBuffer( mGeos, distance, segments ) ); 06356 } 06357 CATCH_GEOS( 0 ) 06358 } 06359 06360 QgsGeometry* QgsGeometry::simplify( double tolerance ) 06361 { 06362 if ( !mGeos ) 06363 { 06364 exportWkbToGeos(); 06365 } 06366 if ( !mGeos ) 06367 { 06368 return 0; 06369 } 06370 try 06371 { 06372 return fromGeosGeom( GEOSSimplify( mGeos, tolerance ) ); 06373 } 06374 CATCH_GEOS( 0 ) 06375 } 06376 06377 QgsGeometry* QgsGeometry::centroid() 06378 { 06379 if ( !mGeos ) 06380 { 06381 exportWkbToGeos(); 06382 } 06383 if ( !mGeos ) 06384 { 06385 return 0; 06386 } 06387 try 06388 { 06389 return fromGeosGeom( GEOSGetCentroid( mGeos ) ); 06390 } 06391 CATCH_GEOS( 0 ) 06392 } 06393 06394 QgsGeometry* QgsGeometry::convexHull() 06395 { 06396 if ( !mGeos ) 06397 { 06398 exportWkbToGeos(); 06399 } 06400 if ( !mGeos ) 06401 { 06402 return 0; 06403 } 06404 06405 try 06406 { 06407 return fromGeosGeom( GEOSConvexHull( mGeos ) ); 06408 } 06409 CATCH_GEOS( 0 ) 06410 } 06411 06412 QgsGeometry* QgsGeometry::intersection( QgsGeometry* geometry ) 06413 { 06414 if ( !geometry ) 06415 { 06416 return NULL; 06417 } 06418 if ( !mGeos ) 06419 { 06420 exportWkbToGeos(); 06421 } 06422 if ( !geometry->mGeos ) 06423 { 06424 geometry->exportWkbToGeos(); 06425 } 06426 if ( !mGeos || !geometry->mGeos ) 06427 { 06428 return 0; 06429 } 06430 06431 try 06432 { 06433 return fromGeosGeom( GEOSIntersection( mGeos, geometry->mGeos ) ); 06434 } 06435 CATCH_GEOS( 0 ) 06436 } 06437 06438 QgsGeometry* QgsGeometry::combine( QgsGeometry* geometry ) 06439 { 06440 if ( !geometry ) 06441 { 06442 return NULL; 06443 } 06444 if ( !mGeos ) 06445 { 06446 exportWkbToGeos(); 06447 } 06448 if ( !geometry->mGeos ) 06449 { 06450 geometry->exportWkbToGeos(); 06451 } 06452 if ( !mGeos || !geometry->mGeos ) 06453 { 06454 return 0; 06455 } 06456 06457 try 06458 { 06459 GEOSGeometry* unionGeom = GEOSUnion( mGeos, geometry->mGeos ); 06460 if ( type() == QGis::Line ) 06461 { 06462 GEOSGeometry* mergedGeom = GEOSLineMerge( unionGeom ); 06463 if ( mergedGeom ) 06464 { 06465 GEOSGeom_destroy( unionGeom ); 06466 unionGeom = mergedGeom; 06467 } 06468 } 06469 return fromGeosGeom( unionGeom ); 06470 } 06471 CATCH_GEOS( new QgsGeometry( *this ) ) //return this geometry if union not possible 06472 } 06473 06474 QgsGeometry* QgsGeometry::difference( QgsGeometry* geometry ) 06475 { 06476 if ( !geometry ) 06477 { 06478 return NULL; 06479 } 06480 if ( !mGeos ) 06481 { 06482 exportWkbToGeos(); 06483 } 06484 if ( !geometry->mGeos ) 06485 { 06486 geometry->exportWkbToGeos(); 06487 } 06488 if ( !mGeos || !geometry->mGeos ) 06489 { 06490 return 0; 06491 } 06492 06493 try 06494 { 06495 return fromGeosGeom( GEOSDifference( mGeos, geometry->mGeos ) ); 06496 } 06497 CATCH_GEOS( 0 ) 06498 } 06499 06500 QgsGeometry* QgsGeometry::symDifference( QgsGeometry* geometry ) 06501 { 06502 if ( !geometry ) 06503 { 06504 return NULL; 06505 } 06506 if ( !mGeos ) 06507 { 06508 exportWkbToGeos(); 06509 } 06510 if ( !geometry->mGeos ) 06511 { 06512 geometry->exportWkbToGeos(); 06513 } 06514 if ( !mGeos || !geometry->mGeos ) 06515 { 06516 return 0; 06517 } 06518 06519 try 06520 { 06521 return fromGeosGeom( GEOSSymDifference( mGeos, geometry->mGeos ) ); 06522 } 06523 CATCH_GEOS( 0 ) 06524 } 06525 06526 QList<QgsGeometry*> QgsGeometry::asGeometryCollection() 06527 { 06528 if ( !mGeos ) 06529 { 06530 exportWkbToGeos(); 06531 if ( !mGeos ) 06532 return QList<QgsGeometry*>(); 06533 } 06534 06535 int type = GEOSGeomTypeId( mGeos ); 06536 QgsDebugMsg( "geom type: " + QString::number( type ) ); 06537 06538 QList<QgsGeometry*> geomCollection; 06539 06540 if ( type != GEOS_MULTIPOINT && 06541 type != GEOS_MULTILINESTRING && 06542 type != GEOS_MULTIPOLYGON && 06543 type != GEOS_GEOMETRYCOLLECTION ) 06544 { 06545 // we have a single-part geometry - put there a copy of this one 06546 geomCollection.append( new QgsGeometry( *this ) ); 06547 return geomCollection; 06548 } 06549 06550 int count = GEOSGetNumGeometries( mGeos ); 06551 QgsDebugMsg( "geom count: " + QString::number( count ) ); 06552 06553 for ( int i = 0; i < count; ++i ) 06554 { 06555 const GEOSGeometry * geometry = GEOSGetGeometryN( mGeos, i ); 06556 geomCollection.append( fromGeosGeom( GEOSGeom_clone( geometry ) ) ); 06557 } 06558 06559 return geomCollection; 06560 } 06561 06562 06563 bool QgsGeometry::deleteRing( int ringNum, int partNum ) 06564 { 06565 if ( ringNum <= 0 || partNum < 0 ) 06566 return false; 06567 06568 switch ( wkbType() ) 06569 { 06570 case QGis::WKBPolygon25D: 06571 case QGis::WKBPolygon: 06572 { 06573 if ( partNum != 0 ) 06574 return false; 06575 06576 QgsPolygon polygon = asPolygon(); 06577 if ( ringNum >= polygon.count() ) 06578 return false; 06579 06580 polygon.remove( ringNum ); 06581 06582 QgsGeometry* g2 = QgsGeometry::fromPolygon( polygon ); 06583 *this = *g2; 06584 delete g2; 06585 return true; 06586 } 06587 06588 case QGis::WKBMultiPolygon25D: 06589 case QGis::WKBMultiPolygon: 06590 { 06591 QgsMultiPolygon mpolygon = asMultiPolygon(); 06592 06593 if ( partNum >= mpolygon.count() ) 06594 return false; 06595 06596 if ( ringNum >= mpolygon[partNum].count() ) 06597 return false; 06598 06599 mpolygon[partNum].remove( ringNum ); 06600 06601 QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon ); 06602 *this = *g2; 06603 delete g2; 06604 return true; 06605 } 06606 06607 default: 06608 return false; // only makes sense with polygons and multipolygons 06609 } 06610 } 06611 06612 06613 bool QgsGeometry::deletePart( int partNum ) 06614 { 06615 if ( partNum < 0 ) 06616 return false; 06617 06618 switch ( wkbType() ) 06619 { 06620 case QGis::WKBMultiPoint25D: 06621 case QGis::WKBMultiPoint: 06622 { 06623 QgsMultiPoint mpoint = asMultiPoint(); 06624 06625 if ( partNum >= mpoint.size() || mpoint.size() == 1 ) 06626 return false; 06627 06628 mpoint.remove( partNum ); 06629 06630 QgsGeometry* g2 = QgsGeometry::fromMultiPoint( mpoint ); 06631 *this = *g2; 06632 delete g2; 06633 break; 06634 } 06635 06636 case QGis::WKBMultiLineString25D: 06637 case QGis::WKBMultiLineString: 06638 { 06639 QgsMultiPolyline mline = asMultiPolyline(); 06640 06641 if ( partNum >= mline.size() || mline.size() == 1 ) 06642 return false; 06643 06644 mline.remove( partNum ); 06645 06646 QgsGeometry* g2 = QgsGeometry::fromMultiPolyline( mline ); 06647 *this = *g2; 06648 delete g2; 06649 break; 06650 } 06651 06652 case QGis::WKBMultiPolygon25D: 06653 case QGis::WKBMultiPolygon: 06654 { 06655 QgsMultiPolygon mpolygon = asMultiPolygon(); 06656 06657 if ( partNum >= mpolygon.size() || mpolygon.size() == 1 ) 06658 return false; 06659 06660 mpolygon.remove( partNum ); 06661 06662 QgsGeometry* g2 = QgsGeometry::fromMultiPolygon( mpolygon ); 06663 *this = *g2; 06664 delete g2; 06665 break; 06666 } 06667 06668 default: 06669 // single part geometries are ignored 06670 return false; 06671 } 06672 06673 return true; 06674 } 06675 06676 int QgsGeometry::avoidIntersections() 06677 { 06678 int returnValue = 0; 06679 06680 //check if g has polygon type 06681 if ( type() != QGis::Polygon ) 06682 { 06683 return 1; 06684 } 06685 06686 QGis::WkbType geomTypeBeforeModification = wkbType(); 06687 06688 //read avoid intersections list from project properties 06689 bool listReadOk; 06690 QStringList avoidIntersectionsList = QgsProject::instance()->readListEntry( "Digitizing", "/AvoidIntersectionsList", &listReadOk ); 06691 if ( !listReadOk ) 06692 { 06693 return true; //no intersections stored in project does not mean error 06694 } 06695 06696 //go through list, convert each layer to vector layer and call QgsVectorLayer::removePolygonIntersections for each 06697 QgsVectorLayer* currentLayer = 0; 06698 QStringList::const_iterator aIt = avoidIntersectionsList.constBegin(); 06699 for ( ; aIt != avoidIntersectionsList.constEnd(); ++aIt ) 06700 { 06701 currentLayer = dynamic_cast<QgsVectorLayer*>( QgsMapLayerRegistry::instance()->mapLayer( *aIt ) ); 06702 if ( currentLayer ) 06703 { 06704 if ( currentLayer->removePolygonIntersections( this ) != 0 ) 06705 { 06706 returnValue = 3; 06707 } 06708 } 06709 } 06710 06711 //make sure the geometry still has the same type (e.g. no change from polygon to multipolygon) 06712 if ( wkbType() != geomTypeBeforeModification ) 06713 { 06714 return 2; 06715 } 06716 06717 return returnValue; 06718 } 06719 06720 void QgsGeometry::validateGeometry( QList<Error> &errors ) 06721 { 06722 QgsGeometryValidator::validateGeometry( this, errors ); 06723 } 06724 06725 bool QgsGeometry::isGeosValid() 06726 { 06727 try 06728 { 06729 GEOSGeometry *g = asGeos(); 06730 06731 if ( !g ) 06732 return false; 06733 06734 return GEOSisValid( g ); 06735 } 06736 catch ( GEOSException &e ) 06737 { 06738 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 06739 return false; 06740 } 06741 } 06742 06743 bool QgsGeometry::isGeosEqual( QgsGeometry &g ) 06744 { 06745 return geosRelOp( GEOSEquals, this, &g ); 06746 } 06747 06748 bool QgsGeometry::isGeosEmpty() 06749 { 06750 try 06751 { 06752 GEOSGeometry *g = asGeos(); 06753 06754 if ( !g ) 06755 return false; 06756 06757 return GEOSisEmpty( g ); 06758 } 06759 catch ( GEOSException &e ) 06760 { 06761 QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) ); 06762 return false; 06763 } 06764 } 06765 06766 double QgsGeometry::leftOf( double x, double y, double& x1, double& y1, double& x2, double& y2 ) 06767 { 06768 double f1 = x - x1; 06769 double f2 = y2 - y1; 06770 double f3 = y - y1; 06771 double f4 = x2 - x1; 06772 return f1*f2 - f3*f4; 06773 }