qgis-geos.diff

update patch (support for GEOS<3, addRing & addIsland fixed) - Jürgen Fischer, 2008-08-07 09:33 AM

Download (105 KB)

View differences:

src/core/qgsgeometry.cpp (Arbeitskopie)
15 15
/* $Id$ */
16 16

  
17 17
#include <limits>
18
#include <cstdarg>
19
#include <cstdio>
18 20

  
19 21
#include "qgis.h"
20 22
#include "qgsgeometry.h"
......
23 25
#include "qgspoint.h"
24 26
#include "qgsrect.h"
25 27

  
26
#if GEOS_VERSION_MAJOR < 3
27
#include <geos/io.h>
28
#include <geos/opLinemerge.h>
29
#include <geos/opPolygonize.h>
30
#include <geos/util.h>
31
#define GEOS_IO geos
32
#define GEOS_LINEMERGE geos
33
#define GEOS_POLYGONIZE geos
34
#define GEOS_UTIL geos
35
#define GEOS_SIZE_T int
36
#define COORD_SEQ_FACTORY DefaultCoordinateSequenceFactory
37
#define GEOS_EXCEPTION GEOS_UTIL::GEOSException* 
28
#include <geos_c.h>
38 29

  
30
#define DEFAULT_QUADRANT_SEGMENTS 8
31

  
39 32
#define CATCH_GEOS(r) \
40
  catch (GEOS_EXCEPTION e) \
33
  catch (GEOSException &e) \
41 34
  { \
42
    QString error = e->toString().c_str(); \
43
    delete e; \
44
    QgsDebugMsg("GEOS: " + error); \
35
    QgsDebugMsg("GEOS: " + QString( e.what() ) ); \
45 36
    return r; \
46 37
  }
47
#else
48
#include <geos/geom/CoordinateArraySequence.h>
49
#include <geos/geom/CoordinateArraySequenceFactory.h>
50
#include <geos/geom/LineString.h>
51
#include <geos/geom/LinearRing.h>
52
#include <geos/geom/MultiLineString.h>
53
#include <geos/geom/MultiPoint.h>
54
#include <geos/geom/MultiPolygon.h>
55
#include <geos/geom/Point.h>
56
#include <geos/geom/Polygon.h>
57
#include <geos/io/WKTReader.h>
58
#include <geos/operation/linemerge/LineMerger.h>
59
#include <geos/operation/polygonize/Polygonizer.h>
60
#include <geos/util/IllegalArgumentException.h>
61
#define GEOS_IO geos::io
62
#define GEOS_LINEMERGE geos::operation::linemerge
63
#define GEOS_POLYGONIZE geos::operation::polygonize
64
#define GEOS_UTIL geos::util
65
#define GEOS_SIZE_T size_t
66
#define COORD_SEQ_FACTORY CoordinateArraySequenceFactory
67
#define GEOS_EXCEPTION GEOS_UTIL::GEOSException&
68 38

  
69
#define CATCH_GEOS(r) \
70
  catch (GEOS_EXCEPTION e) \
39
#define RELAY_GEOS() \
40
  catch (GEOSException &e) \
71 41
  { \
72
    QgsDebugMsg("GEOS: " + QString( e.what() ) ); \
73
    return r; \
42
    Q_UNUSED(e); \
43
    throw; \
74 44
  }
45

  
46
class GEOSException
47
{
48
public:
49
  GEOSException(const char *theMsg)
50
  {
51
    msg  = theMsg;
52
  }
53

  
54
  ~GEOSException()
55
  {
56
    delete [] msg;
57
  }
58

  
59

  
60
  const char *what()
61
  {
62
    return msg;
63
  }
64

  
65
private:
66
  const char *msg;
67
};
68

  
69
extern "C" void throwGEOSException(const char *fmt, ...)
70
#ifdef _MSC_VER
71
throw()
75 72
#endif
73
{
74
  va_list ap;
75
  va_start(ap, fmt);
76
  size_t buflen = vsnprintf(NULL, 0, fmt, ap);
77
  char *msg=new char[buflen+1];
78
  vsnprintf(msg, buflen+1, fmt, ap);
79
  va_end(ap);
76 80

  
77
// Set up static GEOS geometry factory
78
static GEOS_GEOM::GeometryFactory* geosGeometryFactory = new GEOS_GEOM::GeometryFactory();
81
  throw GEOSException(msg);
82
}
79 83

  
84
extern "C" void printGEOSNotice(const char *fmt, ...)
85
{
86
#if defined(QGISDEBUG)
87
  va_list ap;
88
  va_start(ap, fmt);
89
  size_t buflen = vsnprintf(NULL, 0, fmt, ap);
90
  char *msg=new char[buflen+1];
91
  vsnprintf(msg, buflen+1, fmt, ap);
92
  va_end(ap);
93
  QgsDebugMsg( QString("GEOS notice: ").arg(msg) );
94
  delete [] msg;
95
#endif
96
}
80 97

  
98

  
99
int QgsGeometry::refcount=0;
100

  
81 101
QgsGeometry::QgsGeometry()
82 102
  : mGeometry(0),
83
  mGeometrySize(0),
84
  mGeos(0),
85

  
86
  mDirtyWkb(FALSE),
87
mDirtyGeos(FALSE)
103
    mGeometrySize(0),
104
    mGeos(0),
105
    mDirtyWkb(FALSE),
106
    mDirtyGeos(FALSE)
88 107
{
108
  if(refcount++==0) {
109
    initGEOS(printGEOSNotice, throwGEOSException);
110
  }
89 111
  // NOOP
90 112
}    
91 113

  
92 114

  
93 115
QgsGeometry::QgsGeometry( QgsGeometry const & rhs )
94 116
  : mGeometry(0),
95
  mGeometrySize( rhs.mGeometrySize ),
96
  mDirtyWkb( rhs.mDirtyWkb ),
97
  mDirtyGeos( rhs.mDirtyGeos )
117
    mGeometrySize( rhs.mGeometrySize ),
118
    mDirtyWkb( rhs.mDirtyWkb ),
119
    mDirtyGeos( rhs.mDirtyGeos )
98 120
{
121
  if(refcount++==0) {
122
    initGEOS(printGEOSNotice, throwGEOSException);
123
  }
124

  
99 125
  if ( mGeometrySize && rhs.mGeometry )
100 126
  {
101 127
    mGeometry = new unsigned char[mGeometrySize];
......
105 131
  // deep-copy the GEOS Geometry if appropriate
106 132
  if (rhs.mGeos)
107 133
  {  
108
    if(rhs.mGeos->getGeometryTypeId() == GEOS_GEOM::GEOS_MULTIPOLYGON)//MH:problems with cloning for multipolygons in geos 2
109
    {
110
      GEOS_GEOM::MultiPolygon* multiPoly = dynamic_cast<GEOS_GEOM::MultiPolygon*>(rhs.mGeos);
111
      if(multiPoly)
112
      {
113
        std::vector<GEOS_GEOM::Geometry*> polygonVector;
114
        for(GEOS_SIZE_T i = 0; i < multiPoly->getNumGeometries(); ++i)
115
        {
116
          polygonVector.push_back((GEOS_GEOM::Geometry*)(multiPoly->getGeometryN(i)));
117
        }
118
        mGeos = geosGeometryFactory->createMultiPolygon(polygonVector);
119
      }
120
    }
121
    else if(rhs.mGeos->getGeometryTypeId() == GEOS_GEOM::GEOS_MULTILINESTRING) //MH: and also for cloning multilines
122
    {
123
      GEOS_GEOM::MultiLineString* multiLine = dynamic_cast<GEOS_GEOM::MultiLineString*>(rhs.mGeos);
124
      if(multiLine)
125
      {
126
        std::vector<GEOS_GEOM::Geometry*> lineVector;
127
        for(GEOS_SIZE_T i = 0; i < multiLine->getNumGeometries(); ++i)
128
        {
129
          lineVector.push_back((GEOS_GEOM::Geometry*)(multiLine->getGeometryN(i)));
130
        }
131
        mGeos = geosGeometryFactory->createMultiLineString(lineVector);
132
      }
133
    }
134
    else
135
    {
136
      mGeos = rhs.mGeos->clone();
137
    }
134
    mGeos = GEOSGeom_clone(rhs.mGeos);
138 135
  }
139 136
  else
140 137
  {
......
142 139
  }
143 140
}
144 141

  
145
QgsGeometry* QgsGeometry::fromWkt(QString wkt)
142
//! Destructor
143
QgsGeometry::~QgsGeometry()
146 144
{
147
  GEOS_IO::WKTReader reader(geosGeometryFactory);
148
  GEOS_GEOM::Geometry* geom = reader.read(wkt.toLocal8Bit().data());
149
  QgsGeometry* g = new QgsGeometry;
150
  g->setGeos(geom);
151
  return g;
152
}
145
  if (mGeometry)
146
  {
147
    delete [] mGeometry;
148
  }
153 149

  
154
QgsGeometry* QgsGeometry::fromPoint(const QgsPoint& point)
155
{
156
  GEOS_GEOM::Coordinate coord = GEOS_GEOM::Coordinate(point.x(), point.y());
157
  GEOS_GEOM::Geometry* geom = 0;
158
  try
150
  if (mGeos)
159 151
  {
160
    geom = geosGeometryFactory->createPoint(coord);
152
    GEOSGeom_destroy(mGeos);
161 153
  }
162
  CATCH_GEOS(0)
163 154

  
164
  QgsGeometry* g = new QgsGeometry;
165
  g->setGeos(geom);
166
  return g;
155
  if(--refcount==0)
156
    finishGEOS();
167 157
}
168 158

  
169
QgsGeometry* QgsGeometry::fromMultiPoint(const QgsMultiPoint& multipoint)
159
static unsigned int getNumGeosPoints(const GEOSGeometry *geom)
170 160
{
171
  std::vector<GEOS_GEOM::Geometry*>* pointVector = new std::vector<GEOS_GEOM::Geometry*>(multipoint.size());
172
  GEOS_GEOM::Coordinate currentCoord;
173
  std::auto_ptr< std::vector<GEOS_GEOM::Geometry*> > owner(pointVector);
161
  unsigned int n;
162
  const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq(geom);
163
  GEOSCoordSeq_getSize(cs, &n);
164
  return n;
165
}
174 166

  
175
  for(int i = 0; i < multipoint.size(); ++i)
176
  {
177
    currentCoord.x = multipoint.at(i).x();
178
    currentCoord.y = multipoint.at(i).y(); 
179
    try
180
    { 
181
      (*pointVector)[i] = geosGeometryFactory->createPoint(currentCoord);
167
static GEOSGeometry *createGeosPoint(const QgsPoint &point)
168
{
169
  GEOSCoordSequence *coord = GEOSCoordSeq_create(1, 2);
170
  GEOSCoordSeq_setX(coord, 0, point.x());
171
  GEOSCoordSeq_setY(coord, 0, point.y());
172
  return GEOSGeom_createPoint(coord);
173
}
174

  
175
static GEOSCoordSequence *createGeosCoordSequence(const QgsPolyline& points)
176
{
177
  GEOSCoordSequence *coord;
178

  
179
  try {
180
    coord = GEOSCoordSeq_create(points.count(), 2);
181
    int i;
182
    for(i=0; i<points.count(); i++) {
183
      GEOSCoordSeq_setX(coord, i, points[i].x());
184
      GEOSCoordSeq_setY(coord, i, points[i].y());
182 185
    }
183
    CATCH_GEOS(0)
186
    return coord;
187
  } catch(GEOSException &e) {
188
    Q_UNUSED(e);
189
    if(coord)
190
      GEOSCoordSeq_destroy(coord);
191
    throw;
184 192
  }
185

  
186
  GEOS_GEOM::Geometry* geom = 0;
187
  try
188
  {
189
    geom = geosGeometryFactory->createMultiPoint(pointVector);
190
    owner.release();
191
  }
192
  CATCH_GEOS(0)
193
  
194
  QgsGeometry* g = new QgsGeometry;
195
  g->setGeos(geom);
196
  return g;
197 193
}
198 194

  
199
QgsGeometry* QgsGeometry::fromPolyline(const QgsPolyline& polyline)
195
static GEOSGeometry *createGeosCollection(int typeId, QVector<GEOSGeometry*> geoms)
200 196
{
201
  const GEOS_GEOM::CoordinateSequenceFactory* seqFactory = GEOS_GEOM::COORD_SEQ_FACTORY::instance();
202
  GEOS_GEOM::CoordinateSequence* seq = seqFactory->create(polyline.count(), 2);
203
  std::auto_ptr< GEOS_GEOM::CoordinateSequence > owner(seq);
197
  GEOSGeometry **geomarr = new GEOSGeometry*[ geoms.size() ];
198
  if(!geomarr)
199
    return 0;
204 200

  
205
  QgsPolyline::const_iterator it;
206
  int i = 0;
207
  for (it = polyline.begin(); it != polyline.end(); ++it)
208
  {
209
    seq->setAt(GEOS_GEOM::Coordinate(it->x(), it->y()), i++);
201
  for(int i=0; i<geoms.size(); i++)
202
    geomarr[i] = geoms[i];
203

  
204
  GEOSGeometry *geom = 0;
205

  
206
  try {
207
     geom = GEOSGeom_createCollection(typeId, geomarr, geoms.size());
208
  } catch(GEOSException &e) {
209
    Q_UNUSED(e);
210 210
  }
211 211

  
212
  // new geometry takes ownership of the sequence
213
  GEOS_GEOM::Geometry* geom = 0;
214
  try
215
    {
216
      geom = geosGeometryFactory->createLineString(seq);
217
      owner.release();
218
    }
219
  CATCH_GEOS(0)
212
  delete [] geomarr;
220 213

  
221
  QgsGeometry* g = new QgsGeometry;
222
  g->setGeos(geom);
223
  return g;
214
  return geom;
224 215
}
225 216

  
226
QgsGeometry* QgsGeometry::fromMultiPolyline(const QgsMultiPolyline& multiline)
217
static GEOSGeometry *createGeosLineString(const QgsPolyline& polyline)
227 218
{
228
  const GEOS_GEOM::CoordinateSequenceFactory* seqFactory = GEOS_GEOM::COORD_SEQ_FACTORY::instance();
229
  std::vector<GEOS_GEOM::Geometry*>* lineVector = new std::vector<GEOS_GEOM::Geometry*>(multiline.count());
230
  GEOS_GEOM::LineString* currentLineString = 0;
231
  std::auto_ptr< std::vector<GEOS_GEOM::Geometry*> > owner(lineVector);
219
  GEOSCoordSequence *coord = 0;
232 220

  
233
  for(int i = 0; i < multiline.count(); ++i)
221
  try {
222
    coord = createGeosCoordSequence(polyline);
223
    return GEOSGeom_createLineString(coord);
224
  }
225
  catch(GEOSException &e)
234 226
  {
235
    QgsPolyline currentLine = multiline.at(i);
236
    GEOS_GEOM::CoordinateSequence* seq = seqFactory->create(currentLine.count(), 2);
237
    std::auto_ptr< GEOS_GEOM::CoordinateSequence > owner(seq);
227
    Q_UNUSED(e);
228
    if(coord)
229
      GEOSCoordSeq_destroy(coord);
230
    return 0;
231
  }
232
}
238 233

  
239
    for(int j = 0; j < currentLine.count(); ++j)
240
    {
241
      seq->setAt(GEOS_GEOM::Coordinate(currentLine.at(j).x(), currentLine.at(j).y()), j);
234
static GEOSGeometry *createGeosLinearRing(const QgsPolyline& polyline)
235
{
236
  GEOSCoordSequence *coord = 0;
237

  
238
  try {
239
    if( polyline[0]!=polyline[polyline.size()-1] ) {
240
      // Ring not closed
241
      QgsPolyline closed(polyline);
242
      closed << closed[0];
243
      coord = createGeosCoordSequence(closed);
244
    } else {
245
      coord = createGeosCoordSequence(polyline);
242 246
    }
243
    try
244
    {
245
      currentLineString = geosGeometryFactory->createLineString(seq);
246
      owner.release();
247
    }
248
    CATCH_GEOS(0)
249 247

  
250
    (*lineVector)[i] = currentLineString;
248
    return GEOSGeom_createLinearRing(coord);
251 249
  }
252

  
253
  GEOS_GEOM::Geometry* geom = 0;
254
  try
250
  catch(GEOSException &e)
255 251
  {
256
    geom = geosGeometryFactory->createMultiLineString(lineVector);
257
    owner.release();
252
    Q_UNUSED(e);
253
    if(coord)
254
      GEOSCoordSeq_destroy(coord);
255
    return 0;
258 256
  }
259
  CATCH_GEOS(0)
260

  
261
  QgsGeometry* g = new QgsGeometry;
262
  g->setGeos(geom);
263
  return g; 
264 257
}
265 258

  
266
static GEOS_GEOM::LinearRing* _createGeosLinearRing(const QgsPolyline& ring)
259
static GEOSGeometry *createGeosPolygon(const QVector<GEOSGeometry*> &rings)
267 260
{
268
  // LinearRing in GEOS must have the first point the same as the last one
269
  bool needRepeatLastPnt = (ring[0] != ring[ring.count()-1]);
261
  GEOSGeometry *shell = rings[0];
262
  GEOSGeometry **holes = NULL;
270 263

  
271
  const GEOS_GEOM::CoordinateSequenceFactory* seqFactory = GEOS_GEOM::COORD_SEQ_FACTORY::instance();
264
  if(rings.size()>1) {
265
    holes = new GEOSGeometry*[ rings.size()-1 ];
266
    if(!holes)
267
      return 0;
272 268

  
273
  GEOS_GEOM::CoordinateSequence* seq = seqFactory->create(ring.count() + (needRepeatLastPnt ? 1 : 0), 2);
274
  QgsPolyline::const_iterator it;
275
  int i = 0;
276
  for (it = ring.begin(); it != ring.end(); ++it)
277
  {
278
    seq->setAt(GEOS_GEOM::Coordinate(it->x(), it->y()), i++);
269
    for(int i=0; i<rings.size()-1; i++)
270
      holes[i] = rings[i+1];
279 271
  }
280 272

  
281
  // add the first point to close the ring if needed
282
  if (needRepeatLastPnt)
283
    seq->setAt(GEOS_GEOM::Coordinate(ring[0].x(), ring[0].y()), ring.count());
273
  GEOSGeometry *geom = GEOSGeom_createPolygon(shell, holes, rings.size()-1);
284 274

  
285
  // ring takes ownership of the sequence
286
  GEOS_GEOM::LinearRing* linRing = 0;
287
  try
288
  {
289
    linRing = geosGeometryFactory->createLinearRing(seq);
290
  }
291
  CATCH_GEOS(0)
275
  if(holes)
276
    delete [] holes;
277
 
278
  return geom;
279
}
292 280

  
293
  return linRing;
281
static GEOSGeometry *createGeosPolygon(GEOSGeometry *shell)
282
{
283
  return createGeosPolygon( QVector<GEOSGeometry*>() << shell );
294 284
}
295 285

  
296
QgsGeometry* QgsGeometry::fromPolygon(const QgsPolygon& polygon)
286
static GEOSGeometry *createGeosPolygon(const QgsPolygon& polygon)
297 287
{
298 288
  if (polygon.count() == 0)
299
    return NULL;
289
    return 0;
300 290

  
301
  const QgsPolyline& ring0 = polygon[0];
291
  QVector<GEOSGeometry *> geoms;
302 292

  
303
  // outer ring
304
  GEOS_GEOM::LinearRing* outerRing = _createGeosLinearRing(ring0);
305

  
306
  // holes
307
  std::vector<GEOS_GEOM::Geometry*>* holes = new std::vector<GEOS_GEOM::Geometry*> (polygon.count()-1);
308
  for (int i = 1; i < polygon.count(); i++)
309
  {
310
    (*holes)[i-1] = _createGeosLinearRing(polygon[i]);
293
  try {
294
    for(int i=0; i<polygon.count(); i++)
295
      geoms << createGeosLinearRing( polygon[i] );
296
 
297
    return createGeosPolygon( geoms );
298
  } catch(GEOSException &e) {
299
    Q_UNUSED(e);
300
    for(int i=0; i<geoms.count(); i++)
301
      GEOSGeom_destroy( geoms[i] );
302
    return 0;
311 303
  }
304
}
312 305

  
313
  // new geometry takes ownership of outerRing and vector of holes
314
  GEOS_GEOM::Geometry* geom = 0;
315
  try
316
  {
317
    geom = geosGeometryFactory->createPolygon(outerRing, holes);
318
  }
319
  CATCH_GEOS(0)
306
static QgsGeometry *fromGeosGeom(GEOSGeometry *geom)
307
{
308
  if(!geom)
309
    return 0;
320 310

  
321 311
  QgsGeometry* g = new QgsGeometry;
322 312
  g->setGeos(geom);
323 313
  return g;
324 314
}
325 315

  
326
QgsGeometry* QgsGeometry::fromMultiPolygon(const QgsMultiPolygon& multipoly)
316
QgsGeometry* QgsGeometry::fromWkt(QString wkt)
327 317
{
328
  if(multipoly.count() == 0)
329
  {
330
    return 0;
331
  }
318
  GEOSWKTReader *reader = GEOSWKTReader_create();;
319
  QgsGeometry *g = fromGeosGeom( GEOSWKTReader_read(reader, wkt.toLocal8Bit().data()) );
320
  GEOSWKTReader_destroy(reader);
321
  return g;
322
}
332 323

  
333
  std::vector<GEOS_GEOM::Geometry*>* polygons = new std::vector<GEOS_GEOM::Geometry*>(multipoly.count());
334
  std::auto_ptr< std::vector<GEOS_GEOM::Geometry*> > owner(polygons);
335
  GEOS_GEOM::Polygon* currentPolygon = 0;
336
  GEOS_GEOM::LinearRing* currentOuterRing = 0;
337
  std::vector<GEOS_GEOM::Geometry*>* currentHoles = 0;
324
QgsGeometry* QgsGeometry::fromPoint(const QgsPoint& point)
325
{
326
  return fromGeosGeom( createGeosPoint(point) );
327
}
338 328

  
339
  for(int i = 0; i < multipoly.count(); ++i)
340
  {
341
    currentOuterRing = _createGeosLinearRing(multipoly[i].at(0));
342
    currentHoles = new std::vector<GEOS_GEOM::Geometry*>(multipoly[i].count() - 1);
343
    for(int j = 1; j < multipoly[i].count(); ++j)
344
    {
345
      (*currentHoles)[j-1] =  _createGeosLinearRing(multipoly[i].at(j));
329
QgsGeometry* QgsGeometry::fromPolyline(const QgsPolyline& polyline)
330
{
331
  return fromGeosGeom( createGeosLineString(polyline) );
332
}
333

  
334
QgsGeometry* QgsGeometry::fromPolygon(const QgsPolygon& polygon)
335
{
336
  return fromGeosGeom( createGeosPolygon( polygon ) );
337
}
338

  
339
QgsGeometry* QgsGeometry::fromMultiPoint(const QgsMultiPoint& multipoint)
340
{
341
  QVector<GEOSGeometry *> geoms;
342

  
343
  try {
344
    for(int i = 0; i < multipoint.size(); ++i) {
345
      geoms << createGeosPoint( multipoint[i] );
346 346
    }
347
    try
348
    {
349
      currentPolygon = geosGeometryFactory->createPolygon(currentOuterRing, currentHoles);
350
    }
351
    CATCH_GEOS(0)
347
   
348
    return fromGeosGeom( createGeosCollection(GEOS_MULTIPOINT, geoms) );
349
  } catch(GEOSException &e) {
350
    Q_UNUSED(e);
351
    for(int i=0; i < geoms.size(); ++i)
352
      GEOSGeom_destroy(geoms[i]);
352 353

  
353
    (*polygons)[i] = currentPolygon;
354
    return 0;
354 355
  }
356
}
355 357

  
356
  GEOS_GEOM::Geometry* geom = 0;
357
  try
358
  {
359
    geom = geosGeometryFactory->createMultiPolygon(polygons);
360
    owner.release();
358
QgsGeometry* QgsGeometry::fromMultiPolyline(const QgsMultiPolyline& multiline)
359
{
360
  QVector<GEOSGeometry *> geoms;
361

  
362
  try {
363
    for(int i=0; i<multiline.count(); i++)
364
      geoms << createGeosLineString( multiline[i] );
365

  
366
    return fromGeosGeom( createGeosCollection(GEOS_MULTILINESTRING, geoms) );
367
  } catch(GEOSException &e) {
368
    Q_UNUSED(e);
369
    for(int i=0; i<geoms.count(); i++)
370
      GEOSGeom_destroy(geoms[i]);
371

  
372
    return 0;
361 373
  }
362
  CATCH_GEOS(0)
374
}
363 375

  
364
  QgsGeometry* g = new QgsGeometry;
365
  g->setGeos(geom);
366
  return g;
376
QgsGeometry* QgsGeometry::fromMultiPolygon(const QgsMultiPolygon& multipoly)
377
{
378
  if(multipoly.count() == 0)
379
    return 0;
380

  
381
  QVector<GEOSGeometry *> geoms;
382

  
383
  try {
384
    for(int i=0; i<multipoly.count(); i++)
385
      geoms << createGeosPolygon( multipoly[i] );
386

  
387
    return fromGeosGeom( createGeosCollection(GEOS_MULTIPOLYGON, geoms) );
388
  } catch(GEOSException &e) {
389
    Q_UNUSED(e);
390
    for(int i=0; i<geoms.count(); i++)
391
      GEOSGeom_destroy(geoms[i]);
392

  
393
    return 0;
394
  }
367 395
}
368 396

  
369 397
QgsGeometry* QgsGeometry::fromRect(const QgsRect& rect)
......
398 426
  mGeometrySize    = rhs.mGeometrySize;
399 427

  
400 428
  // deep-copy the GEOS Geometry if appropriate
401
  if (rhs.mGeos)
402
  {  
403
    if(rhs.mGeos->getGeometryTypeId() == GEOS_GEOM::GEOS_MULTIPOLYGON)//MH:problems with cloning for multipolygons in geos 2
404
    {
405
      GEOS_GEOM::MultiPolygon* multiPoly = dynamic_cast<GEOS_GEOM::MultiPolygon*>(rhs.mGeos);
406
      if(multiPoly)
407
      {
408
        std::vector<GEOS_GEOM::Geometry*> polygonVector;
409
        for(GEOS_SIZE_T i = 0; i < multiPoly->getNumGeometries(); ++i)
410
        {
411
          polygonVector.push_back((GEOS_GEOM::Geometry*)(multiPoly->getGeometryN(i)));
412
        }
413
        mGeos = geosGeometryFactory->createMultiPolygon(polygonVector);
414
      }
415
    }
416
    else if(rhs.mGeos->getGeometryTypeId() == GEOS_GEOM::GEOS_MULTILINESTRING) //MH: and also for cloning multilines
417
    {
418
      GEOS_GEOM::MultiLineString* multiLine = dynamic_cast<GEOS_GEOM::MultiLineString*>(rhs.mGeos);
419
      if(multiLine)
420
      {
421
        std::vector<GEOS_GEOM::Geometry*> lineVector;
422
        for(GEOS_SIZE_T i = 0; i < multiLine->getNumGeometries(); ++i)
423
        {
424
          lineVector.push_back((GEOS_GEOM::Geometry*)(multiLine->getGeometryN(i)));
425
        }
426
        mGeos = geosGeometryFactory->createMultiLineString(lineVector);
427
      }
428
    }
429
    else
430
    {
431
      mGeos = rhs.mGeos->clone();
432
    }
433
  }
434
  else
435
  {
436
    mGeos = 0;
437
  }    
429
  mGeos = rhs.mGeos ? GEOSGeom_clone(rhs.mGeos) : 0;
438 430

  
439 431
  mDirtyGeos = rhs.mDirtyGeos;
440 432
  mDirtyWkb  = rhs.mDirtyWkb;
......
449 441
} // QgsGeometry::operator=( QgsGeometry const & rhs )
450 442

  
451 443

  
452
//! Destructor
453
QgsGeometry::~QgsGeometry()
454
{
455
  if (mGeometry)
456
  {
457
    delete [] mGeometry;
458
  }
459

  
460
  if (mGeos)
461
  {
462
    delete mGeos;
463
  }
464
}
465

  
466

  
467 444
void QgsGeometry::setWkbAndOwnership(unsigned char * wkb, size_t length)
468 445
{
469 446
  // delete any existing WKB geometry before assigning new one
......
474 451
  }
475 452
  if (mGeos)
476 453
  {
477
    delete mGeos;
454
    GEOSGeom_destroy(mGeos);
478 455
    mGeos = 0;
479 456
  }
480 457

  
......
566 543
}
567 544

  
568 545

  
569
void QgsGeometry::setGeos(GEOS_GEOM::Geometry* geos)
546
void QgsGeometry::setGeos(GEOSGeometry* geos)
570 547
{
571 548
  // TODO - make this more heap-friendly
572 549

  
573 550
  if (mGeos)
574 551
  {
575
    delete mGeos;
552
    GEOSGeom_destroy(mGeos);
576 553
    mGeos = 0;
577 554
  }
578 555
  if (mGeometry)
......
1093 1070

  
1094 1071
bool QgsGeometry::insertVertexBefore(double x, double y,
1095 1072
    int beforeVertex,
1096
    const GEOS_GEOM::CoordinateSequence*  old_sequence,
1097
    GEOS_GEOM::CoordinateSequence** new_sequence)
1073
    const GEOSCoordSequence*  old_sequence,
1074
    GEOSCoordSequence** new_sequence)
1098 1075

  
1099 1076
{
1100 1077
  // Bounds checking
......
1104 1081
    return FALSE;
1105 1082
  }
1106 1083

  
1107
  int numPoints = old_sequence->getSize();
1084
  unsigned int numPoints; 
1085
  GEOSCoordSeq_getSize(old_sequence, &numPoints);
1108 1086

  
1109
  // Copy to the new sequence, including the new vertex
1110
  (*new_sequence) = new GEOS_GEOM::DefaultCoordinateSequence();
1087
  *new_sequence = GEOSCoordSeq_create(numPoints+1, 2);
1088
  if(!*new_sequence)
1089
    return false;
1111 1090

  
1112
  bool inserted = FALSE;
1113
  for (int i = 0; i < numPoints; i++)
1091
  bool inserted = false;
1092
  for (unsigned int i=0, j=0; i<numPoints; i++, j++)
1114 1093
  {
1115 1094
    // Do we insert the new vertex here?
1116 1095
    if (beforeVertex == i)
1117 1096
    {
1118
      (*new_sequence)->add( GEOS_GEOM::Coordinate(x, y) );
1119
      inserted = TRUE;
1097
      GEOSCoordSeq_setX(*new_sequence, j, x);
1098
      GEOSCoordSeq_setY(*new_sequence, j, y);
1099
      j++;
1100
      inserted = true;
1120 1101
    }
1121
    (*new_sequence)->add( old_sequence->getAt(i) );
1102

  
1103
    double aX, aY;
1104
    GEOSCoordSeq_getX(old_sequence, i, &aX);
1105
    GEOSCoordSeq_getY(old_sequence, i, &aY);
1106

  
1107
    GEOSCoordSeq_setX(*new_sequence, j, aX);
1108
    GEOSCoordSeq_setY(*new_sequence, j, aY);
1122 1109
  }
1123 1110

  
1124 1111
  if (!inserted)
1125 1112
  {
1126 1113
    // The beforeVertex is greater than the actual number of vertices
1127 1114
    // in the geometry - append it.
1128
    (*new_sequence)->add( GEOS_GEOM::Coordinate(x, y) );
1115
    GEOSCoordSeq_setX(*new_sequence, numPoints, x);
1116
    GEOSCoordSeq_setY(*new_sequence, numPoints, y);
1129 1117
  }
1130 1118
  // TODO: Check that the sequence is still simple, e.g. with GEOS_GEOM::Geometry->isSimple()
1131 1119

  
......
2018 2006

  
2019 2007
  if (mGeos) //try to find the vertex from the Geos geometry (it present)
2020 2008
  {
2021
    GEOS_GEOM::CoordinateSequence* cs = mGeos->getCoordinates();
2022
    if(cs)
2023
    {
2024
      const GEOS_GEOM::Coordinate& coord = cs->getAt(atVertex);
2025
      x = coord.x;
2026
      y = coord.y;
2027
      delete cs;
2009
    try {
2010
      const GEOSCoordSequence *cs = GEOSGeom_getCoordSeq(mGeos);
2011
      GEOSCoordSeq_getX(cs, atVertex, &x);
2012
      GEOSCoordSeq_getY(cs, atVertex, &y);
2028 2013
      return QgsPoint(x,y);
2014
    } catch(GEOSException &e) {
2015
      Q_UNUSED(e);
2016
      return QgsPoint(0,0);
2029 2017
    }
2030 2018
  }
2031 2019
  else if(mGeometry)
......
2244 2232

  
2245 2233
double QgsGeometry::closestVertexWithContext(const QgsPoint& point, int& atVertex)
2246 2234
{
2247
  // Initialise some stuff
2248 2235
  double sqrDist = std::numeric_limits<double>::max();
2249
  int closestVertexIndex = 0;
2250 2236

  
2251
  // set up the GEOS geometry
2252
  exportWkbToGeos();
2237
  try {
2238
    // Initialise some stuff
2239
    int closestVertexIndex = 0;
2253 2240

  
2254
  if (!mGeos)
2255
  {
2256
    QgsDebugMsg("GEOS geometry not available!");
2257
    return -1;
2258
  }
2241
    // set up the GEOS geometry
2242
    exportWkbToGeos();
2259 2243

  
2260
  GEOS_GEOM::CoordinateSequence* sequence = mGeos->getCoordinates();
2261
  if(sequence)
2262
  {
2263
    for(GEOS_SIZE_T i = 0; i < sequence->getSize(); ++i)
2244
    const GEOSCoordSequence *sequence = GEOSGeom_getCoordSeq(mGeos);
2245

  
2246
    unsigned int n;
2247
    GEOSCoordSeq_getSize(sequence, &n);
2248

  
2249
    for(unsigned int i=0; i<n; i++)
2264 2250
    {
2265
      double testDist = point.sqrDist(sequence->getAt(i).x, sequence->getAt(i).y);
2251
      double x, y;
2252
      GEOSCoordSeq_getX(sequence, i, &x);
2253
      GEOSCoordSeq_getY(sequence, i, &y);
2254

  
2255
      double testDist = point.sqrDist(x,y);
2266 2256
      if(testDist < sqrDist)
2267 2257
      {
2268 2258
        closestVertexIndex = i;
2269 2259
        sqrDist = testDist;
2270 2260
      }
2271 2261
    }
2262

  
2263
    atVertex = closestVertexIndex;
2264
  } catch(GEOSException &e) {
2265
    Q_UNUSED(e);
2266
    return -1;
2272 2267
  }
2273
  atVertex = closestVertexIndex;
2274 2268

  
2275 2269
  return sqrDist;
2276 2270
}
......
2505 2499
{
2506 2500
  //bail out if this geometry is not polygon/multipolygon
2507 2501
  if(vectorType() != QGis::Polygon)
2508
  {
2509 2502
    return 1;
2510
  }
2511 2503

  
2512 2504
  //test for invalid geometries
2513 2505
  if(ring.size() < 4)
2514
  {
2515 2506
    return 3;
2516
  }
2517 2507

  
2518 2508
  //ring must be closed
2519 2509
  if(ring.first() != ring.last())
2520
  {
2521 2510
    return 2;
2522
  }
2523 2511

  
2524 2512
  //create geos geometry from wkb if not already there
2525 2513
  if(!mGeos || mDirtyGeos)
2526
  {
2527 2514
    exportWkbToGeos();
2528
  }
2529 2515

  
2516
  int type = GEOSGeomTypeId(mGeos);
2517

  
2530 2518
  //Fill GEOS Polygons of the feature into list
2531
  std::list<GEOS_GEOM::Polygon*> polygonList; //list of polygon pointers (only one for polygon geometries)
2532
  GEOS_GEOM::Polygon* thisPolygon = 0;
2533
  GEOS_GEOM::MultiPolygon* thisMultiPolygon = 0;
2519
  QVector<const GEOSGeometry*> polygonList;
2534 2520

  
2535 2521
  if(this->wkbType() == QGis::WKBPolygon)
2536 2522
  {
2537
    thisPolygon = dynamic_cast<GEOS_GEOM::Polygon*>(mGeos);
2538
    if(!thisPolygon)
2539
    {
2523
    if(type!=GEOS_POLYGON)
2540 2524
      return 1;
2541
    }
2542
    polygonList.push_back(thisPolygon);
2525

  
2526
    polygonList << mGeos;
2543 2527
  }
2544 2528
  else if(this->wkbType() == QGis::WKBMultiPolygon)
2545 2529
  {
2546
    thisMultiPolygon = dynamic_cast<GEOS_GEOM::MultiPolygon*>(mGeos);
2547
    if(!thisMultiPolygon)
2548
    {
2530
    if(type!=GEOS_MULTIPOLYGON)
2549 2531
      return 1;
2550
    }
2551
    int numPolys = thisMultiPolygon->getNumGeometries();
2552
    for(int i = 0; i < numPolys; ++i)
2553
    {
2554
      polygonList.push_back((GEOS_GEOM::Polygon*)(thisMultiPolygon->getGeometryN(i)));
2555
    }
2556
  }
2557 2532

  
2558
  //create new ring
2559
  GEOS_GEOM::DefaultCoordinateSequence* newSequence=new GEOS_GEOM::DefaultCoordinateSequence();
2560
  std::auto_ptr< GEOS_GEOM::DefaultCoordinateSequence > owner(newSequence);
2561
  for(QList<QgsPoint>::const_iterator it = ring.begin(); it != ring.end(); ++it)
2562
  {
2563
    newSequence->add(GEOS_GEOM::Coordinate(it->x(),it->y()));
2533
    for(int i=0; i < GEOSGetNumGeometries(mGeos); ++i)
2534
      polygonList << GEOSGetGeometryN(mGeos, i);
2564 2535
  }
2565 2536

  
2566 2537
  //create new ring
2567
  GEOS_GEOM::LinearRing* newRing = 0;
2568
  try
2569
  {
2570
    newRing = geosGeometryFactory->createLinearRing(newSequence);
2571
    owner.release();
2572
  }
2573
  CATCH_GEOS(3)
2538
  GEOSGeometry *newRing = 0;
2539
  GEOSGeometry *newRingPolygon = 0;
2574 2540

  
2575
  std::vector<GEOS_GEOM::Geometry*> dummyVector;
2541
  try {
2542
    newRing = createGeosLinearRing(ring.toVector());
2543
    if( !GEOSisValid(newRing) )
2544
      throw GEOSException("ring is invalid");
2576 2545

  
2577
  //create polygon from new ring because there is a problem with geos operations and linear rings
2578
  GEOS_GEOM::Polygon* newRingPolygon = geosGeometryFactory->createPolygon(*newRing, dummyVector);
2579
  if(!newRing || !newRingPolygon || !newRing->isValid() || !newRingPolygon->isValid())
2580
  {
2581
    QgsDebugMsg("ring is not valid");
2582
    delete newRing;
2583
    delete newRingPolygon;
2546
    newRingPolygon = createGeosPolygon( newRing );
2547
    if( !GEOSisValid(newRingPolygon) )
2548
      throw GEOSException("ring polygon is invalid");
2549
  } catch(GEOSException &e) {
2550
    Q_UNUSED(e);
2551
    if(newRingPolygon)
2552
      GEOSGeom_destroy(newRingPolygon);
2553
    else if(newRing)
2554
      GEOSGeom_destroy(newRing);
2555
      
2584 2556
    return 3;
2585 2557
  }
2586 2558

  
2587
  GEOS_GEOM::LinearRing* outerRing = 0; //outer ring of already existing feature
2588
  std::vector<GEOS_GEOM::Geometry*>* inner = 0; //vector of inner rings. The existing rings and the new one will be added
2589
  int numberOfPolyContainingRing = 0; //for multipolygons: store index of the polygon where the ring is
2590
  bool foundPoly = false; //set to true as soon we found a polygon containing the ring
2559
  QVector<GEOSGeometry*> rings;
2591 2560

  
2592
  for(std::list<GEOS_GEOM::Polygon*>::const_iterator it = polygonList.begin(); it != polygonList.end(); ++it)
2593
  {
2594
    /***********inner rings*****************************************/
2595
    //consider already existing rings
2596
    inner=new std::vector<GEOS_GEOM::Geometry*>();
2597
    int numExistingRings = (*it)->getNumInteriorRing();
2598
    inner->resize(numExistingRings + 1);
2599
    for(int i = 0; i < numExistingRings; ++i)
2600
    {
2601
      GEOS_GEOM::LinearRing* existingRing = geosGeometryFactory->createLinearRing((*it)->getInteriorRingN(i)->getCoordinates());
2602
      //create polygon from new ring because there is a problem with geos operations and linear rings
2603
      GEOS_GEOM::Polygon* existingRingPolygon = geosGeometryFactory->createPolygon(*existingRing, dummyVector);
2561
  int i;
2562
  for(i=0; i<polygonList.size(); i++) {
2563
    for(int j=0; j<rings.size(); j++)
2564
      GEOSGeom_destroy(rings[j]);
2565
    rings.clear();
2604 2566

  
2605
      //check, if the new ring intersects the existing one and bail out if yes
2606
      //if(existingRing->disjoint(newRing))
2607
      if(!existingRingPolygon->disjoint(newRingPolygon)) //does only work with polygons, not linear rings
2608
      {
2609
        QgsDebugMsg("new ring not disjoint with existing ring");
2610
        //delete objects wich are no longer needed
2611
        delete existingRing;
2612
        delete existingRingPolygon;
2613
        for(std::vector<GEOS_GEOM::Geometry*>::iterator it = inner->begin(); it != inner->end(); ++it)
2614
        {
2615
          delete *it;
2616
        }
2617
        delete inner;
2618
        delete newRing;
2619
        delete newRingPolygon;
2620
        return 4; //error: ring not disjoint with existing rings
2567
    GEOSGeometry *shellRing = 0;
2568
    GEOSGeometry *shell = 0;
2569
    try {
2570
      shellRing = GEOSGeom_clone( GEOSGetExteriorRing( polygonList[i] ) );
2571
      shell = createGeosPolygon(shellRing);
2572

  
2573
      if( !GEOSWithin(newRingPolygon, shell) ) {
2574
        GEOSGeom_destroy(shell);
2575
        continue;
2621 2576
      }
2622
      (*inner)[i] = existingRing;
2623
      delete existingRingPolygon; //delete since this polygon has only be created for disjoint() test
2577
    } catch(GEOSException &e) {
2578
      Q_UNUSED(e);
2579
      if(shell)
2580
        GEOSGeom_destroy(shell);
2581
      else if(shellRing)
2582
        GEOSGeom_destroy(shellRing);
2583
      return 4;
2624 2584
    }
2625 2585

  
2626
    //also add new ring to the vector
2627
    if(newRing)
2628
    {
2629
      (*inner)[numExistingRings] = newRing;
2630
    }
2586
    // add outer ring
2587
    rings << GEOSGeom_clone(shellRing);
2631 2588

  
2632
    /*****************outer ring****************/
2633
    outerRing = geosGeometryFactory->createLinearRing((*it)->getExteriorRing()->getCoordinates());
2634
    //create polygon from new ring because there is a problem with geos operations and linear rings
2635
    GEOS_GEOM::Polygon* outerRingPolygon = geosGeometryFactory->createPolygon(*outerRing, dummyVector);
2589
    // check inner rings
2590
    int n = GEOSGetNumInteriorRings(polygonList[i]);
2591
     
2592
    int j;
2593
    for(j=0; j<n; j++) {
2594
      GEOSGeometry *holeRing = 0;
2595
      GEOSGeometry *hole = 0;
2596
      try {
2597
        holeRing = GEOSGeom_clone( GEOSGetInteriorRingN( polygonList[i], j ) );
2598
        hole = createGeosPolygon(holeRing);
2636 2599

  
2637
    //check if the new ring is within the outer shell of the polygon and bail out if not
2638
    //if(newRing->within(outerRing))
2639
    if(newRingPolygon->within(outerRingPolygon)) //does only work for polygons, not linear rings
2640
    {
2641
      QgsDebugMsg("new ring within outer ring");
2642
      foundPoly = true;
2643
      delete outerRingPolygon;
2644
      break; //ring is in geometry and does not intersect existing rings -> proceed with adding ring to feature
2645
    }
2646

  
2647
    //we need to search in other polygons...
2648
    for(std::vector<GEOS_GEOM::Geometry*>::iterator it = inner->begin(); it != inner->end(); ++it)
2649
    {
2650
      if( (*it) != newRing) //we need newRing for later polygons
2651
      {
2652
        delete *it;
2600
        if( GEOSDisjoint(hole, newRingPolygon) )
2601
          continue;
2602
      } catch(GEOSException &e) {
2603
        Q_UNUSED(e);
2604
        if(hole)
2605
          GEOSGeom_destroy(hole);
2606
        else if(holeRing)
2607
          GEOSGeom_destroy(holeRing);
2608
        throw;
2653 2609
      }
2610

  
2611
      rings << GEOSGeom_clone(holeRing);
2612

  
2613
      GEOSGeom_destroy(hole);
2654 2614
    }
2655
    delete inner;
2656
    delete outerRing;
2657
    delete outerRingPolygon;
2658 2615

  
2659
    ++numberOfPolyContainingRing;
2616
    if(j==n)
2617
      // this is it...
2618
      break;
2660 2619
  }
2661 2620

  
2662
  delete newRingPolygon;
2621
  if( i==polygonList.size() ) {
2622
    // clear rings
2623
    for(int j=0; j<rings.size(); j++)
2624
      GEOSGeom_destroy(rings[j]);
2625
    rings.clear();
2663 2626

  
2664
  if(foundPoly)
2665
  {
2666
    GEOS_GEOM::Polygon* newPolygon = geosGeometryFactory->createPolygon(outerRing,inner);
2667
    if(this->wkbType() == QGis::WKBPolygon)
2668
    {
2669
      delete mGeos;
2670
      mGeos = newPolygon;
2671
    }
2672
    else if(this->wkbType() == QGis::WKBMultiPolygon)
2673
    {
2674
      //remember in which polygon the ring is and replace only this polygon
2675
      std::vector<GEOS_GEOM::Geometry*>* polygons = new std::vector<GEOS_GEOM::Geometry*>();
2676
      int numPolys = thisMultiPolygon->getNumGeometries();
2677
      for(int i = 0; i < numPolys; ++i)
2678
      {
2679
        if(i == numberOfPolyContainingRing)
2680
        {
2681
          polygons->push_back(newPolygon);
2682
        }
2683
        else
2684
        {
2685
          GEOS_GEOM::Polygon* p = (GEOS_GEOM::Polygon*)(thisMultiPolygon->getGeometryN(i)->clone());
2686
          polygons->push_back(p);
2687
        }
2688
      }
2689
      delete mGeos;
2690
      mGeos = geosGeometryFactory->createMultiPolygon(polygons);
2691
    }
2692
    mDirtyWkb = true;
2693
    mDirtyGeos = false;
2694
    return 0;
2695
  }
2696
  else
2697
  {
2698
    delete newRing;
2627
    // no containing polygon found
2699 2628
    return 5;
2700 2629
  }
2630

  
2631
  rings << newRing;
2632
 
2633
  GEOSGeometry *newPolygon = createGeosPolygon(rings);
2634

  
2635
  if(this->wkbType() == QGis::WKBPolygon) {
2636
    GEOSGeom_destroy(mGeos);
2637
    mGeos = newPolygon;
2638
  } else if(this->wkbType() == QGis::WKBMultiPolygon) {
2639
    QVector<GEOSGeometry*> newPolygons;
2640

  
2641
    for(int j=0; j<polygonList.size(); j++) {
2642
      newPolygons << (i==j ? newPolygon : GEOSGeom_clone(polygonList[j]) );
2643
    }
2644

  
2645
    GEOSGeom_destroy(mGeos);
2646
    mGeos = createGeosCollection(GEOS_MULTIPOLYGON, newPolygons);
2647
  }
2648

  
2649
  mDirtyWkb = true;
2650
  mDirtyGeos = false;
2651
  return 0;
2701 2652
}
2702 2653

  
2703 2654
int QgsGeometry::addIsland(const QList<QgsPoint>& ring)
......
2734 2685
    exportWkbToGeos();
2735 2686
  }
2736 2687

  
2737
  //this multipolygon
2738
  GEOS_GEOM::MultiPolygon* thisMultiPolygon = dynamic_cast<GEOS_GEOM::MultiPolygon*>(mGeos);
2739
  if(!thisMultiPolygon)
2740
  {
2688
  if( GEOSGeomTypeId(mGeos)!=GEOS_MULTIPOLYGON) 
2741 2689
    return 1;
2742
  }
2743 2690

  
2744 2691
  //create new polygon from ring
2745 2692

  
2746 2693
  //coordinate sequence first
2747
  GEOS_GEOM::DefaultCoordinateSequence* newSequence=new GEOS_GEOM::DefaultCoordinateSequence();
2748
  std::auto_ptr< GEOS_GEOM::DefaultCoordinateSequence > owner(newSequence);
2749
  for(QList<QgsPoint>::const_iterator it = ring.begin(); it != ring.end(); ++it)
2750
  {
2751
    newSequence->add(GEOS_GEOM::Coordinate(it->x(),it->y()));
2752
  }
2694
  GEOSGeometry *newRing = 0;
2695
  GEOSGeometry *newPolygon = 0;
2753 2696

  
2754
  //then linear ring
2755
  GEOS_GEOM::LinearRing* newRing = 0;
2756
  try
2757
  {
2758
    newRing = geosGeometryFactory->createLinearRing(newSequence);
2759
    owner.release();
2760
  }
2761
  CATCH_GEOS(2)
2762

  
2763
  //finally the polygon
2764
  std::vector<GEOS_GEOM::Geometry*> dummyVector;
2765
  GEOS_GEOM::Polygon* newPolygon = geosGeometryFactory->createPolygon(*newRing, dummyVector);
2766
  delete newRing;
2767

  
2768
  if(!newPolygon || !newPolygon->isValid())
2769
  {
2770
    delete newPolygon;
2697
  try {
2698
    newRing = createGeosLinearRing(ring.toVector());
2699
    if( !GEOSisValid(newRing) )
2700
      throw GEOSException("ring invalid");
2701
    newPolygon = createGeosPolygon(newRing);
2702
    if( !GEOSisValid(newPolygon) )
2703
      throw GEOSException("polygon invalid");
2704
  } catch(GEOSException &e) {
2705
    Q_UNUSED(e);
2706
    if(newPolygon)
2707
      GEOSGeom_destroy(newPolygon);
2708
    else if(newRing)
2709
      GEOSGeom_destroy(newRing);
2771 2710
    return 2;
2772 2711
  }
2773 2712

  
2713
  QVector<GEOSGeometry*> polygons;
2714

  
2774 2715
  //create new multipolygon
2775
  std::vector<GEOS_GEOM::Geometry*>* newMultiPolygonVector = new std::vector<GEOS_GEOM::Geometry*>();
2776
  for(GEOS_SIZE_T i = 0; i < thisMultiPolygon->getNumGeometries(); ++i)
2716
  int n = GEOSGetNumGeometries(mGeos);
2717
  int i;
2718
  for(i=0; i<n; ++i)
2777 2719
  {
2778
    const GEOS_GEOM::Geometry* polygonN = thisMultiPolygon->getGeometryN(i);
2720
    const GEOSGeometry *polygonN = GEOSGetGeometryN(mGeos, i);
2779 2721

  
2780
    //bail out if new polygon is not disjoint with existing ones
2781
    if(!polygonN->disjoint(newPolygon))
2782
    {
2783
      delete newPolygon;
2784
      for(std::vector<GEOS_GEOM::Geometry*>::iterator it =newMultiPolygonVector->begin(); it != newMultiPolygonVector->end(); ++it)
2785
      {
2786
        delete *it;
2787
      }
2788
      delete newMultiPolygonVector;
2789
      return 3;
2790
    }
2791
    newMultiPolygonVector->push_back(polygonN->clone());
2722
    if( !GEOSDisjoint(polygonN, newPolygon) )
2723
      //bail out if new polygon is not disjoint with existing ones
2724
      break;
2725

  
2726
    polygons << GEOSGeom_clone(polygonN);
2792 2727
  }
2793
  newMultiPolygonVector->push_back(newPolygon);
2794
  GEOS_GEOM::MultiPolygon* newMultiPolygon = geosGeometryFactory->createMultiPolygon(newMultiPolygonVector);
2795 2728

  
2796
  delete mGeos;
2797
  mGeos = newMultiPolygon;
2729
  if(i<n) {
2730
    // bailed out
2731
    for(int i=0; i<polygons.size(); i++)
2732
      GEOSGeom_destroy(polygons[i]);
2733
    return 3;
2734
  }
2798 2735

  
2736
  polygons << newPolygon;
2737

  
2738
  GEOSGeom_destroy(mGeos);
2739
  mGeos = createGeosCollection( GEOS_MULTIPOLYGON, polygons );
2799 2740
  mDirtyWkb = true;
2800 2741
  mDirtyGeos = false;
2801 2742
  return 0;
......
2957 2898

  
2958 2899
  try
2959 2900
  {
2960
    GEOS_GEOM::DefaultCoordinateSequence* splitLineCoords = new GEOS_GEOM::DefaultCoordinateSequence();
2961
    QList<QgsPoint>::const_iterator lineIt;
2962
    for(lineIt = splitLine.constBegin(); lineIt != splitLine.constEnd(); ++lineIt)
2901
    GEOSGeometry *splitLineGeos = createGeosLineString(splitLine.toVector());
2902
    if(!GEOSisValid(splitLineGeos) || !GEOSisSimple(splitLineGeos))
2963 2903
    {
2964
      splitLineCoords->add(GEOS_GEOM::Coordinate(lineIt->x(),lineIt->y()));
2965
    }
2966
    GEOS_GEOM::LineString* splitLineGeos = geosGeometryFactory->createLineString(splitLineCoords);
2967
    if(!splitLineGeos)
2968
    {
2969
      delete splitLineCoords;
2904
      GEOSGeom_destroy(splitLineGeos);
2970 2905
      return 1;
2971 2906
    }
2972 2907

  
2973
    if(!splitLineGeos->isValid() || !splitLineGeos->isSimple())
2974
    {
2975
      delete splitLineGeos;
2976
      return 1;
2977
    }
2978

  
2979 2908
    //for line/multiline: call splitLinearGeometry
2980 2909
    if(vectorType() == QGis::Line)
2981 2910
    {
2982 2911
      returnCode = splitLinearGeometry(splitLineGeos, newGeometries);
2983
      delete splitLineGeos;
2912
      GEOSGeom_destroy(splitLineGeos);
2984 2913
    }
2985 2914
    else if(vectorType() == QGis::Polygon)
2986 2915
    {
2987 2916
      returnCode = splitPolygonGeometry(splitLineGeos, newGeometries);
2988
      delete splitLineGeos;
2917
      GEOSGeom_destroy(splitLineGeos);
2989 2918
    }
2990 2919
    else
2991 2920
    {
......
3010 2939
    return 1;
3011 2940
  }
3012 2941

  
3013
  if(!mGeos->isValid())
2942
  if(!GEOSisValid(mGeos))
3014 2943
  {
3015 2944
    return 2;
3016 2945
  }
3017 2946

  
3018
  if(!mGeos->isSimple())
2947
  if(!GEOSisSimple(mGeos))
3019 2948
  {
3020 2949
    return 3;
3021 2950
  }
......
3034 2963
  //make geometry::difference
3035 2964
  try
3036 2965
  {
3037
    if(mGeos->intersects(other->mGeos))
2966
    if(GEOSIntersects(mGeos, other->mGeos))
3038 2967
    {
3039 2968
      //check if multitype before and after
3040 2969
      bool multiType = isMultipart();
3041 2970
 
3042
      mGeos = mGeos->difference(other->mGeos);
2971
      mGeos = GEOSDifference(mGeos, other->mGeos);
3043 2972
      mDirtyWkb = true;
3044 2973

  
3045 2974
      if(multiType && !isMultipart())
......
3365 3294
      return false;
3366 3295
    }
3367 3296

  
3368
     return mGeos->intersects(geometry->mGeos);
3297
     return GEOSIntersects(mGeos, geometry->mGeos);
3369 3298
  }
3370 3299
  CATCH_GEOS(false)
3371 3300
}
......
3381 3310
    return false;
3382 3311
  }
3383 3312

  
3384
  GEOS_GEOM::Point* geosPoint = geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(p->x(), p->y()));
3385
  std::auto_ptr< GEOS_GEOM::Point > owner(geosPoint);
3313
  GEOSGeometry *geosPoint = 0;
3386 3314

  
3387
  bool returnval;
3388
  try
3389
  {
3390
    returnval = mGeos->contains(geosPoint);
3391
    owner.release();  
3315
  bool returnval = false;
3316

  
3317
  try {
3318
    geosPoint = createGeosPoint( *p );
3319
    returnval = GEOSContains(mGeos, geosPoint);
3320
  } catch(GEOSException &e) {
3321
    Q_UNUSED(e);
3322
    returnval = false;
3392 3323
  }
3393
  CATCH_GEOS(false)
3394 3324

  
3325
  if(geosPoint)
3326
    GEOSGeom_destroy(geosPoint);
3327

  
3395 3328
  return returnval;
3396 3329
}
3397 3330

  
......
3690 3623

  
3691 3624
  if (mGeos)
3692 3625
  {
3693
    delete mGeos;
3626
    GEOSGeom_destroy(mGeos);
3694 3627
    mGeos = 0;
3695 3628
  }
3696 3629

  
......
3727 3660
        x = (double *) (mGeometry + 5);
3728 3661
        y = (double *) (mGeometry + 5 + sizeof(double));
3729 3662

  
3730
        mGeos = geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(*x,*y));
3663
        mGeos = createGeosPoint( QgsPoint(*x,*y) );
3731 3664
        mDirtyGeos = FALSE;
3732 3665
        break;
3733 3666
      }
......
3736 3669
      hasZValue = true;
3737 3670
    case QGis::WKBMultiPoint:
3738 3671
      {
3739
        std::vector<GEOS_GEOM::Geometry*>* points=new std::vector<GEOS_GEOM::Geometry*>;
3672
        QVector<GEOSGeometry *> points;
3673

  
3740 3674
        ptr = mGeometry + 5;
3741 3675
        nPoints = (int *) ptr;
3742 3676
        ptr = mGeometry + 1 + 2 * sizeof(int);
......
3751 3685
          {
3752 3686
            ptr += sizeof(double);
3753 3687
          }
3754
          points->push_back(geosGeometryFactory->createPoint(GEOS_GEOM::Coordinate(*x,*y)));
3688
          points << createGeosPoint( QgsPoint(*x,*y) );
3755 3689
        }
3756
        mGeos = geosGeometryFactory->createMultiPoint(points);
3690
        mGeos = createGeosCollection(GEOS_MULTIPOINT, points);
3757 3691
        mDirtyGeos = FALSE;
3758 3692
        break;
3759 3693
      }
......
3764 3698
      {
3765 3699
        QgsDebugMsg("QgsGeometry::geosGeometry: Linestring found");
3766 3700

  
3767
        GEOS_GEOM::DefaultCoordinateSequence* sequence=new GEOS_GEOM::DefaultCoordinateSequence();
3701
        QgsPolyline sequence;
3702

  
3768 3703
        ptr = mGeometry + 5;
3769 3704
        nPoints = (int *) ptr;
3770 3705
        ptr = mGeometry + 1 + 2 * sizeof(int);
......
3778 3713
          {
3779 3714
            ptr += sizeof(double);
3780 3715
          }
3781
          sequence->add(GEOS_GEOM::Coordinate(*x,*y));
3716

  
3717
          sequence << QgsPoint(*x,*y);
3782 3718
        }
3783 3719
        mDirtyGeos = FALSE;
3784
        mGeos = geosGeometryFactory->createLineString(sequence); 
3720
        mGeos = createGeosLineString(sequence);
3785 3721
        break;
3786 3722
      }
3787 3723

  
......
3789 3725
      hasZValue = true;
3790 3726
    case QGis::WKBMultiLineString:
3791 3727
      {
3792
        std::vector<GEOS_GEOM::Geometry*>* lines=new std::vector<GEOS_GEOM::Geometry*>;
3728
        QVector<GEOSGeometry*> lines;
3793 3729
        numLineStrings = (int) (mGeometry[5]);
3794 3730
        ptr = (mGeometry + 9);
3795 3731
        for (jdx = 0; jdx < numLineStrings; jdx++)
3796 3732
        {
3797
          GEOS_GEOM::DefaultCoordinateSequence* sequence=new GEOS_GEOM::DefaultCoordinateSequence();
3733
          QgsPolyline sequence;
3734

  
3798 3735
          // each of these is a wbklinestring so must handle as such
3799 3736
          lsb = *ptr;
3800 3737
          ptr += 5;   // skip type since we know its 2
......
3810 3747
            {
3811 3748
              ptr += sizeof(double);
3812 3749
            }
3813
            sequence->add(GEOS_GEOM::Coordinate(*x,*y));
3750
            sequence << QgsPoint(*x,*y);
3814 3751
          }
3815
          lines->push_back(geosGeometryFactory->createLineString(sequence));
3752
          lines << createGeosLineString(sequence);
3816 3753
        }
3817
        mGeos = geosGeometryFactory->createMultiLineString(lines);
3754
        mGeos = createGeosCollection(GEOS_MULTILINESTRING, lines);
3818 3755
        mDirtyGeos = FALSE;
3819 3756
        break;
3820 3757
      }
......
3829 3766
        numRings = (int *) (mGeometry + 1 + sizeof(int));
3830 3767
        ptr = mGeometry + 1 + 2 * sizeof(int);
3831 3768

  
3832
        GEOS_GEOM::LinearRing* outer=0;
3833
        std::vector<GEOS_GEOM::Geometry*>* inner=new std::vector<GEOS_GEOM::Geometry*>;
3769
        QVector<GEOSGeometry*> rings;
3834 3770

  
3835 3771
        for (idx = 0; idx < *numRings; idx++)
3836 3772
        {
3837

  
3838 3773
          //QgsDebugMsg("Ring nr: "+QString::number(idx));
3839 3774

  
3840
          GEOS_GEOM::DefaultCoordinateSequence* sequence=new GEOS_GEOM::DefaultCoordinateSequence();
3775
          QgsPolyline sequence;
3776

  
3841 3777
          // get number of points in the ring
3842 3778
          nPoints = (int *) ptr;
3843 3779
          ptr += 4;
......
3852 3788
            {
3853 3789
              ptr += sizeof(double);
3854 3790
            }
3855
            sequence->add(GEOS_GEOM::Coordinate(*x,*y));
3791
            sequence << QgsPoint(*x,*y);
3856 3792
          }
3857
          GEOS_GEOM::LinearRing* ring=geosGeometryFactory->createLinearRing(sequence);
3858
          if(idx==0)
3859
          {
3860
            outer=ring;
3861
          }
... This diff was truncated because it exceeds the maximum size that can be displayed.