1 : /******************************************************************************
2 : * $Id: ogrgeometryfactory.cpp 18320 2009-12-17 02:22:03Z warmerdam $
3 : *
4 : * Project: OpenGIS Simple Features Reference Implementation
5 : * Purpose: Factory for converting geometry to and from well known binary
6 : * format.
7 : * Author: Frank Warmerdam, warmerdam@pobox.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 1999, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "ogr_geometry.h"
32 : #include "ogr_api.h"
33 : #include "ogr_p.h"
34 : #include "ogr_geos.h"
35 :
36 : CPL_CVSID("$Id: ogrgeometryfactory.cpp 18320 2009-12-17 02:22:03Z warmerdam $");
37 :
38 : #ifndef PI
39 : #define PI 3.14159265358979323846
40 : #endif
41 :
42 : /************************************************************************/
43 : /* createFromWkb() */
44 : /************************************************************************/
45 :
46 : /**
47 : * \brief Create a geometry object of the appropriate type from it's well known binary representation.
48 : *
49 : * Note that if nBytes is passed as zero, no checking can be done on whether
50 : * the pabyData is sufficient. This can result in a crash if the input
51 : * data is corrupt. This function returns no indication of the number of
52 : * bytes from the data source actually used to represent the returned
53 : * geometry object. Use OGRGeometry::WkbSize() on the returned geometry to
54 : * establish the number of bytes it required in WKB format.
55 : *
56 : * Also note that this is a static method, and that there
57 : * is no need to instantiate an OGRGeometryFactory object.
58 : *
59 : * The C function OGR_G_CreateFromWkb() is the same as this method.
60 : *
61 : * @param pabyData pointer to the input BLOB data.
62 : * @param poSR pointer to the spatial reference to be assigned to the
63 : * created geometry object. This may be NULL.
64 : * @param ppoReturn the newly created geometry object will be assigned to the
65 : * indicated pointer on return. This will be NULL in case
66 : * of failure.
67 : * @param nBytes the number of bytes available in pabyData, or -1 if it isn't
68 : * known.
69 : *
70 : * @return OGRERR_NONE if all goes well, otherwise any of
71 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
72 : * OGRERR_CORRUPT_DATA may be returned.
73 : */
74 :
75 1564 : OGRErr OGRGeometryFactory::createFromWkb(unsigned char *pabyData,
76 : OGRSpatialReference * poSR,
77 : OGRGeometry **ppoReturn,
78 : int nBytes )
79 :
80 : {
81 : OGRwkbGeometryType eGeometryType;
82 : OGRwkbByteOrder eByteOrder;
83 : OGRErr eErr;
84 : OGRGeometry *poGeom;
85 :
86 1564 : *ppoReturn = NULL;
87 :
88 1564 : if( nBytes < 9 && nBytes != -1 )
89 1014 : return OGRERR_NOT_ENOUGH_DATA;
90 :
91 : /* -------------------------------------------------------------------- */
92 : /* Get the byte order byte. The extra tests are to work around */
93 : /* bug sin the WKB of DB2 v7.2 as identified by Safe Software. */
94 : /* -------------------------------------------------------------------- */
95 550 : eByteOrder = DB2_V72_FIX_BYTE_ORDER((OGRwkbByteOrder) *pabyData);
96 :
97 :
98 550 : if( eByteOrder != wkbXDR && eByteOrder != wkbNDR )
99 : {
100 : CPLDebug( "OGR",
101 : "OGRGeometryFactory::createFromWkb() - got corrupt data.\n"
102 : "%02X%02X%02X%02X%02X%02X%02X%02X%02X\n",
103 0 : pabyData[0],
104 0 : pabyData[1],
105 0 : pabyData[2],
106 0 : pabyData[3],
107 0 : pabyData[4],
108 0 : pabyData[5],
109 0 : pabyData[6],
110 0 : pabyData[7],
111 0 : pabyData[8]);
112 0 : return OGRERR_CORRUPT_DATA;
113 : }
114 :
115 : /* -------------------------------------------------------------------- */
116 : /* Get the geometry feature type. For now we assume that */
117 : /* geometry type is between 0 and 255 so we only have to fetch */
118 : /* one byte. */
119 : /* -------------------------------------------------------------------- */
120 550 : if( eByteOrder == wkbNDR )
121 483 : eGeometryType = (OGRwkbGeometryType) pabyData[1];
122 : else
123 67 : eGeometryType = (OGRwkbGeometryType) pabyData[4];
124 :
125 : /* -------------------------------------------------------------------- */
126 : /* Instantiate a geometry of the appropriate type, and */
127 : /* initialize from the input stream. */
128 : /* -------------------------------------------------------------------- */
129 550 : poGeom = createGeometry( eGeometryType );
130 :
131 550 : if( poGeom == NULL )
132 6 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
133 :
134 : /* -------------------------------------------------------------------- */
135 : /* Import from binary. */
136 : /* -------------------------------------------------------------------- */
137 544 : eErr = poGeom->importFromWkb( pabyData, nBytes );
138 :
139 : /* -------------------------------------------------------------------- */
140 : /* Assign spatial reference system. */
141 : /* -------------------------------------------------------------------- */
142 544 : if( eErr == OGRERR_NONE )
143 : {
144 544 : poGeom->assignSpatialReference( poSR );
145 544 : *ppoReturn = poGeom;
146 : }
147 : else
148 : {
149 0 : delete poGeom;
150 : }
151 :
152 544 : return eErr;
153 : }
154 :
155 : /************************************************************************/
156 : /* OGR_G_CreateFromWkb() */
157 : /************************************************************************/
158 : /**
159 : * \brief Create a geometry object of the appropriate type from it's well known binary representation.
160 : *
161 : * Note that if nBytes is passed as zero, no checking can be done on whether
162 : * the pabyData is sufficient. This can result in a crash if the input
163 : * data is corrupt. This function returns no indication of the number of
164 : * bytes from the data source actually used to represent the returned
165 : * geometry object. Use OGR_G_WkbSize() on the returned geometry to
166 : * establish the number of bytes it required in WKB format.
167 : *
168 : * The OGRGeometryFactory::createFromWkb() CPP method is the same as this
169 : * function.
170 : *
171 : * @param pabyData pointer to the input BLOB data.
172 : * @param hSRS handle to the spatial reference to be assigned to the
173 : * created geometry object. This may be NULL.
174 : * @param phGeometry the newly created geometry object will
175 : * be assigned to the indicated handle on return. This will be NULL in case
176 : * of failure.
177 : * @param nBytes the number of bytes of data available in pabyData, or -1
178 : * if it is not known, but assumed to be sufficient.
179 : *
180 : * @return OGRERR_NONE if all goes well, otherwise any of
181 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
182 : * OGRERR_CORRUPT_DATA may be returned.
183 : */
184 :
185 62 : OGRErr CPL_DLL OGR_G_CreateFromWkb( unsigned char *pabyData,
186 : OGRSpatialReferenceH hSRS,
187 : OGRGeometryH *phGeometry,
188 : int nBytes )
189 :
190 : {
191 : return OGRGeometryFactory::createFromWkb( pabyData,
192 : (OGRSpatialReference *) hSRS,
193 : (OGRGeometry **) phGeometry,
194 62 : nBytes );
195 : }
196 :
197 : /************************************************************************/
198 : /* createFromWkt() */
199 : /************************************************************************/
200 :
201 : /**
202 : * \brief Create a geometry object of the appropriate type from it's well known text representation.
203 : *
204 : * The C function OGR_G_CreateFromWkt() is the same as this method.
205 : *
206 : * @param ppszData input zero terminated string containing well known text
207 : * representation of the geometry to be created. The pointer
208 : * is updated to point just beyond that last character consumed.
209 : * @param poSR pointer to the spatial reference to be assigned to the
210 : * created geometry object. This may be NULL.
211 : * @param ppoReturn the newly created geometry object will be assigned to the
212 : * indicated pointer on return. This will be NULL if the
213 : * method fails.
214 : *
215 : * <b>Example:</b>
216 : *
217 : * <pre>
218 : * const char* wkt= "POINT(0 0)";
219 : *
220 : * // cast because OGR_G_CreateFromWkt will move the pointer
221 : * char* pszWkt = (char*) wkt.c_str();
222 : * OGRSpatialReferenceH ref = OSRNewSpatialReference(NULL);
223 : * OGRGeometryH new_geom;
224 : * OGRErr err = OGR_G_CreateFromWkt(&pszWkt, ref, &new_geom);
225 : * </pre>
226 : *
227 : *
228 : *
229 : * @return OGRERR_NONE if all goes well, otherwise any of
230 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
231 : * OGRERR_CORRUPT_DATA may be returned.
232 : */
233 :
234 16351 : OGRErr OGRGeometryFactory::createFromWkt(char **ppszData,
235 : OGRSpatialReference * poSR,
236 : OGRGeometry **ppoReturn )
237 :
238 : {
239 : OGRErr eErr;
240 : char szToken[OGR_WKT_TOKEN_MAX];
241 16351 : char *pszInput = *ppszData;
242 : OGRGeometry *poGeom;
243 :
244 16351 : *ppoReturn = NULL;
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* Get the first token, which should be the geometry type. */
248 : /* -------------------------------------------------------------------- */
249 16351 : if( OGRWktReadToken( pszInput, szToken ) == NULL )
250 0 : return OGRERR_CORRUPT_DATA;
251 :
252 : /* -------------------------------------------------------------------- */
253 : /* Instantiate a geometry of the appropriate type. */
254 : /* -------------------------------------------------------------------- */
255 16351 : if( EQUAL(szToken,"POINT") )
256 : {
257 14881 : poGeom = new OGRPoint();
258 : }
259 :
260 1470 : else if( EQUAL(szToken,"LINESTRING") )
261 : {
262 128 : poGeom = new OGRLineString();
263 : }
264 :
265 1342 : else if( EQUAL(szToken,"POLYGON") )
266 : {
267 179 : poGeom = new OGRPolygon();
268 : }
269 :
270 1163 : else if( EQUAL(szToken,"GEOMETRYCOLLECTION") )
271 : {
272 21 : poGeom = new OGRGeometryCollection();
273 : }
274 :
275 1142 : else if( EQUAL(szToken,"MULTIPOLYGON") )
276 : {
277 54 : poGeom = new OGRMultiPolygon();
278 : }
279 :
280 1088 : else if( EQUAL(szToken,"MULTIPOINT") )
281 : {
282 36 : poGeom = new OGRMultiPoint();
283 : }
284 :
285 1052 : else if( EQUAL(szToken,"MULTILINESTRING") )
286 : {
287 42 : poGeom = new OGRMultiLineString();
288 : }
289 :
290 : else
291 : {
292 1010 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
293 : }
294 :
295 : /* -------------------------------------------------------------------- */
296 : /* Do the import. */
297 : /* -------------------------------------------------------------------- */
298 15341 : eErr = poGeom->importFromWkt( &pszInput );
299 :
300 : /* -------------------------------------------------------------------- */
301 : /* Assign spatial reference system. */
302 : /* -------------------------------------------------------------------- */
303 15341 : if( eErr == OGRERR_NONE )
304 : {
305 15341 : poGeom->assignSpatialReference( poSR );
306 15341 : *ppoReturn = poGeom;
307 15341 : *ppszData = pszInput;
308 : }
309 : else
310 : {
311 0 : delete poGeom;
312 : }
313 :
314 15341 : return eErr;
315 : }
316 :
317 : /************************************************************************/
318 : /* OGR_G_CreateFromWkt() */
319 : /************************************************************************/
320 : /**
321 : * \brief Create a geometry object of the appropriate type from it's well known text representation.
322 : *
323 : * The OGRGeometryFactory::createFromWkt CPP method is the same as this
324 : * function.
325 : *
326 : * @param ppszData input zero terminated string containing well known text
327 : * representation of the geometry to be created. The pointer
328 : * is updated to point just beyond that last character consumed.
329 : * @param hSRS handle to the spatial reference to be assigned to the
330 : * created geometry object. This may be NULL.
331 : * @param phGeometry the newly created geometry object will be assigned to the
332 : * indicated handle on return. This will be NULL if the
333 : * method fails.
334 : *
335 : * @return OGRERR_NONE if all goes well, otherwise any of
336 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
337 : * OGRERR_CORRUPT_DATA may be returned.
338 : */
339 :
340 15149 : OGRErr CPL_DLL OGR_G_CreateFromWkt( char **ppszData,
341 : OGRSpatialReferenceH hSRS,
342 : OGRGeometryH *phGeometry )
343 :
344 : {
345 : return OGRGeometryFactory::createFromWkt( ppszData,
346 : (OGRSpatialReference *) hSRS,
347 15149 : (OGRGeometry **) phGeometry );
348 : }
349 :
350 : /************************************************************************/
351 : /* createGeometry() */
352 : /************************************************************************/
353 :
354 : /**
355 : * \brief Create an empty geometry of desired type.
356 : *
357 : * This is equivelent to allocating the desired geometry with new, but
358 : * the allocation is guaranteed to take place in the context of the
359 : * GDAL/OGR heap.
360 : *
361 : * This method is the same as the C function OGR_G_CreateGeometry().
362 : *
363 : * @param eGeometryType the type code of the geometry class to be instantiated.
364 : *
365 : * @return the newly create geometry or NULL on failure.
366 : */
367 :
368 : OGRGeometry *
369 1111 : OGRGeometryFactory::createGeometry( OGRwkbGeometryType eGeometryType )
370 :
371 : {
372 1111 : switch( wkbFlatten(eGeometryType) )
373 : {
374 : case wkbPoint:
375 140 : return new OGRPoint();
376 :
377 : case wkbLineString:
378 265 : return new OGRLineString();
379 :
380 : case wkbPolygon:
381 403 : return new OGRPolygon();
382 :
383 : case wkbGeometryCollection:
384 18 : return new OGRGeometryCollection();
385 :
386 : case wkbMultiPolygon:
387 37 : return new OGRMultiPolygon();
388 :
389 : case wkbMultiPoint:
390 20 : return new OGRMultiPoint();
391 :
392 : case wkbMultiLineString:
393 26 : return new OGRMultiLineString();
394 :
395 : case wkbLinearRing:
396 196 : return new OGRLinearRing();
397 :
398 : default:
399 6 : return NULL;
400 : }
401 : }
402 :
403 : /************************************************************************/
404 : /* OGR_G_CreateGeometry() */
405 : /************************************************************************/
406 : /**
407 : * \brief Create an empty geometry of desired type.
408 : *
409 : * This is equivelent to allocating the desired geometry with new, but
410 : * the allocation is guaranteed to take place in the context of the
411 : * GDAL/OGR heap.
412 : *
413 : * This function is the same as the CPP method
414 : * OGRGeometryFactory::createGeometry.
415 : *
416 : * @param eGeometryType the type code of the geometry to be created.
417 : *
418 : * @return handle to the newly create geometry or NULL on failure.
419 : */
420 :
421 560 : OGRGeometryH OGR_G_CreateGeometry( OGRwkbGeometryType eGeometryType )
422 :
423 : {
424 560 : return (OGRGeometryH) OGRGeometryFactory::createGeometry( eGeometryType );
425 : }
426 :
427 :
428 : /************************************************************************/
429 : /* destroyGeometry() */
430 : /************************************************************************/
431 :
432 : /**
433 : * \brief Destroy geometry object.
434 : *
435 : * Equivalent to invoking delete on a geometry, but it guaranteed to take
436 : * place within the context of the GDAL/OGR heap.
437 : *
438 : * This method is the same as the C function OGR_G_DestroyGeometry().
439 : *
440 : * @param poGeom the geometry to deallocate.
441 : */
442 :
443 15974 : void OGRGeometryFactory::destroyGeometry( OGRGeometry *poGeom )
444 :
445 : {
446 15974 : delete poGeom;
447 15974 : }
448 :
449 :
450 : /************************************************************************/
451 : /* OGR_G_DestroyGeometry() */
452 : /************************************************************************/
453 : /**
454 : * \brief Destroy geometry object.
455 : *
456 : * Equivalent to invoking delete on a geometry, but it guaranteed to take
457 : * place within the context of the GDAL/OGR heap.
458 : *
459 : * This function is the same as the CPP method
460 : * OGRGeometryFactory::destroyGeometry.
461 : *
462 : * @param hGeom handle to the geometry to delete.
463 : */
464 :
465 15899 : void OGR_G_DestroyGeometry( OGRGeometryH hGeom )
466 :
467 : {
468 15899 : OGRGeometryFactory::destroyGeometry( (OGRGeometry *) hGeom );
469 15899 : }
470 :
471 : /************************************************************************/
472 : /* forceToPolygon() */
473 : /************************************************************************/
474 :
475 : /**
476 : * \brief Convert to polygon.
477 : *
478 : * Tries to force the provided geometry to be a polygon. Currently
479 : * this just effects a change on multipolygons. The passed in geometry is
480 : * consumed and a new one returned (or potentially the same one).
481 : *
482 : * @return new geometry.
483 : */
484 :
485 10 : OGRGeometry *OGRGeometryFactory::forceToPolygon( OGRGeometry *poGeom )
486 :
487 : {
488 10 : if( poGeom == NULL )
489 0 : return NULL;
490 :
491 10 : if( wkbFlatten(poGeom->getGeometryType()) != wkbGeometryCollection
492 0 : || wkbFlatten(poGeom->getGeometryType()) != wkbMultiPolygon )
493 10 : return poGeom;
494 :
495 : // build an aggregated polygon from all the polygon rings in the container.
496 0 : OGRPolygon *poPolygon = new OGRPolygon();
497 0 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
498 : int iGeom;
499 :
500 0 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
501 : {
502 0 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
503 : != wkbPolygon )
504 0 : continue;
505 :
506 0 : OGRPolygon *poOldPoly = (OGRPolygon *) poGC->getGeometryRef(iGeom);
507 : int iRing;
508 :
509 0 : poPolygon->addRing( poOldPoly->getExteriorRing() );
510 :
511 0 : for( iRing = 0; iRing < poOldPoly->getNumInteriorRings(); iRing++ )
512 0 : poPolygon->addRing( poOldPoly->getInteriorRing( iRing ) );
513 : }
514 :
515 0 : delete poGC;
516 :
517 0 : return poPolygon;
518 : }
519 :
520 : /************************************************************************/
521 : /* forceToMultiPolygon() */
522 : /************************************************************************/
523 :
524 : /**
525 : * \brief Convert to multipolygon.
526 : *
527 : * Tries to force the provided geometry to be a multipolygon. Currently
528 : * this just effects a change on polygons. The passed in geometry is
529 : * consumed and a new one returned (or potentially the same one).
530 : *
531 : * @return new geometry.
532 : */
533 :
534 0 : OGRGeometry *OGRGeometryFactory::forceToMultiPolygon( OGRGeometry *poGeom )
535 :
536 : {
537 0 : if( poGeom == NULL )
538 0 : return NULL;
539 :
540 : /* -------------------------------------------------------------------- */
541 : /* Check for the case of a geometrycollection that can be */
542 : /* promoted to MultiPolygon. */
543 : /* -------------------------------------------------------------------- */
544 0 : if( wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection )
545 : {
546 : int iGeom;
547 0 : int bAllPoly = TRUE;
548 0 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
549 :
550 0 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
551 : {
552 0 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
553 : != wkbPolygon )
554 0 : bAllPoly = FALSE;
555 : }
556 :
557 0 : if( !bAllPoly )
558 0 : return poGeom;
559 :
560 0 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
561 :
562 0 : while( poGC->getNumGeometries() > 0 )
563 : {
564 0 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
565 0 : poGC->removeGeometry( 0, FALSE );
566 : }
567 :
568 0 : delete poGC;
569 :
570 0 : return poMP;
571 : }
572 :
573 : /* -------------------------------------------------------------------- */
574 : /* Eventually we should try to split the polygon into component */
575 : /* island polygons. But thats alot of work and can be put off. */
576 : /* -------------------------------------------------------------------- */
577 0 : if( wkbFlatten(poGeom->getGeometryType()) != wkbPolygon )
578 0 : return poGeom;
579 :
580 0 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
581 0 : poMP->addGeometryDirectly( poGeom );
582 :
583 0 : return poMP;
584 : }
585 :
586 : /************************************************************************/
587 : /* forceToMultiPoint() */
588 : /************************************************************************/
589 :
590 : /**
591 : * \brief Convert to multipoint.
592 : *
593 : * Tries to force the provided geometry to be a multipoint. Currently
594 : * this just effects a change on points. The passed in geometry is
595 : * consumed and a new one returned (or potentially the same one).
596 : *
597 : * @return new geometry.
598 : */
599 :
600 0 : OGRGeometry *OGRGeometryFactory::forceToMultiPoint( OGRGeometry *poGeom )
601 :
602 : {
603 0 : if( poGeom == NULL )
604 0 : return NULL;
605 :
606 : /* -------------------------------------------------------------------- */
607 : /* Check for the case of a geometrycollection that can be */
608 : /* promoted to MultiPoint. */
609 : /* -------------------------------------------------------------------- */
610 0 : if( wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection )
611 : {
612 : int iGeom;
613 0 : int bAllPoint = TRUE;
614 0 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
615 :
616 0 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
617 : {
618 0 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
619 : != wkbPoint )
620 0 : bAllPoint = FALSE;
621 : }
622 :
623 0 : if( !bAllPoint )
624 0 : return poGeom;
625 :
626 0 : OGRMultiPoint *poMP = new OGRMultiPoint();
627 :
628 0 : while( poGC->getNumGeometries() > 0 )
629 : {
630 0 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
631 0 : poGC->removeGeometry( 0, FALSE );
632 : }
633 :
634 0 : delete poGC;
635 :
636 0 : return poMP;
637 : }
638 :
639 0 : if( wkbFlatten(poGeom->getGeometryType()) != wkbPoint )
640 0 : return poGeom;
641 :
642 0 : OGRMultiPoint *poMP = new OGRMultiPoint();
643 0 : poMP->addGeometryDirectly( poGeom );
644 :
645 0 : return poMP;
646 : }
647 :
648 : /************************************************************************/
649 : /* forceToMultiLinestring() */
650 : /************************************************************************/
651 :
652 : /**
653 : * \brief Convert to multilinestring.
654 : *
655 : * Tries to force the provided geometry to be a multilinestring.
656 : *
657 : * - linestrings are placed in a multilinestring.
658 : * - geometry collections will be converted to multilinestring if they only
659 : * contain linestrings.
660 : * - polygons will be changed to a collection of linestrings (one per ring).
661 : *
662 : * The passed in geometry is
663 : * consumed and a new one returned (or potentially the same one).
664 : *
665 : * @return new geometry.
666 : */
667 :
668 18 : OGRGeometry *OGRGeometryFactory::forceToMultiLineString( OGRGeometry *poGeom )
669 :
670 : {
671 18 : if( poGeom == NULL )
672 0 : return NULL;
673 :
674 : /* -------------------------------------------------------------------- */
675 : /* Check for the case of a geometrycollection that can be */
676 : /* promoted to MultiLineString. */
677 : /* -------------------------------------------------------------------- */
678 18 : if( wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection )
679 : {
680 : int iGeom;
681 0 : int bAllLines = TRUE;
682 0 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
683 :
684 0 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
685 : {
686 0 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
687 : != wkbLineString )
688 0 : bAllLines = FALSE;
689 : }
690 :
691 0 : if( !bAllLines )
692 0 : return poGeom;
693 :
694 0 : OGRMultiLineString *poMP = new OGRMultiLineString();
695 :
696 0 : while( poGC->getNumGeometries() > 0 )
697 : {
698 0 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
699 0 : poGC->removeGeometry( 0, FALSE );
700 : }
701 :
702 0 : delete poGC;
703 :
704 0 : return poMP;
705 : }
706 :
707 : /* -------------------------------------------------------------------- */
708 : /* Turn a linestring into a multilinestring. */
709 : /* -------------------------------------------------------------------- */
710 18 : if( wkbFlatten(poGeom->getGeometryType()) == wkbLineString )
711 : {
712 0 : OGRMultiLineString *poMP = new OGRMultiLineString();
713 0 : poMP->addGeometryDirectly( poGeom );
714 0 : return poMP;
715 : }
716 :
717 : /* -------------------------------------------------------------------- */
718 : /* Convert polygons into a multilinestring. */
719 : /* -------------------------------------------------------------------- */
720 18 : if( wkbFlatten(poGeom->getGeometryType()) == wkbPolygon )
721 : {
722 1 : OGRMultiLineString *poMP = new OGRMultiLineString();
723 1 : OGRPolygon *poPoly = (OGRPolygon *) poGeom;
724 : int iRing;
725 :
726 4 : for( iRing = 0; iRing < poPoly->getNumInteriorRings()+1; iRing++ )
727 : {
728 : OGRLineString *poNewLS, *poLR;
729 :
730 1 : if( iRing == 0 )
731 1 : poLR = poPoly->getExteriorRing();
732 : else
733 0 : poLR = poPoly->getInteriorRing(iRing-1);
734 :
735 1 : poNewLS = new OGRLineString();
736 1 : poNewLS->addSubLineString( poLR );
737 1 : poMP->addGeometryDirectly( poNewLS );
738 : }
739 :
740 1 : delete poPoly;
741 :
742 1 : return poMP;
743 : }
744 :
745 : /* -------------------------------------------------------------------- */
746 : /* Convert multi-polygons into a multilinestring. */
747 : /* -------------------------------------------------------------------- */
748 17 : if( wkbFlatten(poGeom->getGeometryType()) == wkbMultiPolygon )
749 : {
750 0 : OGRMultiLineString *poMP = new OGRMultiLineString();
751 0 : OGRMultiPolygon *poMPoly = (OGRMultiPolygon *) poGeom;
752 : int iPoly;
753 :
754 0 : for( iPoly = 0; iPoly < poMPoly->getNumGeometries(); iPoly++ )
755 : {
756 0 : OGRPolygon *poPoly = (OGRPolygon*) poMPoly->getGeometryRef(iPoly);
757 : int iRing;
758 :
759 0 : for( iRing = 0; iRing < poPoly->getNumInteriorRings()+1; iRing++ )
760 : {
761 : OGRLineString *poNewLS, *poLR;
762 :
763 0 : if( iRing == 0 )
764 0 : poLR = poPoly->getExteriorRing();
765 : else
766 0 : poLR = poPoly->getInteriorRing(iRing-1);
767 :
768 0 : poNewLS = new OGRLineString();
769 0 : poNewLS->addSubLineString( poLR );
770 0 : poMP->addGeometryDirectly( poNewLS );
771 : }
772 : }
773 0 : delete poMPoly;
774 :
775 0 : return poMP;
776 : }
777 :
778 17 : return poGeom;
779 : }
780 :
781 : /************************************************************************/
782 : /* organizePolygons() */
783 : /************************************************************************/
784 :
785 : typedef struct _sPolyExtended sPolyExtended;
786 :
787 : struct _sPolyExtended
788 566 : {
789 : OGRPolygon* poPolygon;
790 : OGREnvelope sEnvelope;
791 : OGRLinearRing* poExteriorRing;
792 : OGRPoint poAPoint;
793 : int nInitialIndex;
794 : int bIsTopLevel;
795 : OGRPolygon* poEnclosingPolygon;
796 : double dfArea;
797 : int bIsCW;
798 : };
799 :
800 337 : static int OGRGeometryFactoryCompareArea(const void* p1, const void* p2)
801 : {
802 337 : const sPolyExtended* psPoly1 = (const sPolyExtended*) p1;
803 337 : const sPolyExtended* psPoly2 = (const sPolyExtended*) p2;
804 337 : if (psPoly2->dfArea < psPoly1->dfArea)
805 188 : return -1;
806 149 : else if (psPoly2->dfArea > psPoly1->dfArea)
807 134 : return 1;
808 : else
809 15 : return 0;
810 : }
811 :
812 348 : static int OGRGeometryFactoryCompareByIndex(const void* p1, const void* p2)
813 : {
814 348 : const sPolyExtended* psPoly1 = (const sPolyExtended*) p1;
815 348 : const sPolyExtended* psPoly2 = (const sPolyExtended*) p2;
816 348 : if (psPoly1->nInitialIndex < psPoly2->nInitialIndex)
817 208 : return -1;
818 140 : else if (psPoly1->nInitialIndex > psPoly2->nInitialIndex)
819 140 : return 1;
820 : else
821 0 : return 0;
822 : }
823 :
824 : #define N_CRITICAL_PART_NUMBER 100
825 :
826 : typedef enum
827 : {
828 : METHOD_NORMAL,
829 : METHOD_SKIP,
830 : METHOD_ONLY_CCW
831 : } OrganizePolygonMethod;
832 :
833 : /**
834 : * \brief Organize polygons based on geometries.
835 : *
836 : * Analyse a set of rings (passed as simple polygons), and based on a
837 : * geometric analysis convert them into a polygon with inner rings,
838 : * or a MultiPolygon if dealing with more than one polygon.
839 : *
840 : * All the input geometries must be OGRPolygons with only a valid exterior
841 : * ring (at least 4 points) and no interior rings.
842 : *
843 : * The passed in geometries become the responsibility of the method, but the
844 : * papoPolygons "pointer array" remains owned by the caller.
845 : *
846 : * For faster computation, a polygon is considered to be inside
847 : * another one if a single point of its external ring is included into the other one.
848 : * (unless 'OGR_DEBUG_ORGANIZE_POLYGONS' configuration option is set to TRUE.
849 : * In that case, a slower algorithm that tests exact topological relationships
850 : * is used if GEOS is available.)
851 : *
852 : * In cases where a big number of polygons is passed to this function, the default processing
853 : * may be really slow. You can skip the processing by adding METHOD=SKIP
854 : * to the option list (the result of the function will be a multi-polygon with all polygons
855 : * as toplevel polygons) or only make it analyze counterclockwise polygons by adding
856 : * METHOD=ONLY_CCW to the option list if you can assume that the outline
857 : * of holes is counterclockwise defined (this is the convention for shapefiles e.g.)
858 : *
859 : * If the OGR_ORGANIZE_POLYGONS configuration option is defined, its value will override
860 : * the value of the METHOD option of papszOptions (usefull to modify the behaviour of the
861 : * shapefile driver)
862 : *
863 : * @param papoPolygons array of geometry pointers - should all be OGRPolygons.
864 : * Ownership of the geometries is passed, but not of the array itself.
865 : * @param nPolygonCount number of items in papoPolygons
866 : * @param pbIsValidGeometry value will be set TRUE if result is valid or
867 : * FALSE otherwise.
868 : * @param papszOptions a list of strings for passing options
869 : *
870 : * @return a single resulting geometry (either OGRPolygon or OGRMultiPolygon).
871 : */
872 :
873 70 : OGRGeometry* OGRGeometryFactory::organizePolygons( OGRGeometry **papoPolygons,
874 : int nPolygonCount,
875 : int *pbIsValidGeometry,
876 : const char** papszOptions )
877 : {
878 : int bUseFastVersion;
879 : int i, j;
880 70 : OGRGeometry* geom = NULL;
881 70 : OrganizePolygonMethod method = METHOD_NORMAL;
882 :
883 : /* -------------------------------------------------------------------- */
884 : /* Trivial case of a single polygon. */
885 : /* -------------------------------------------------------------------- */
886 70 : if (nPolygonCount == 1)
887 : {
888 0 : geom = papoPolygons[0];
889 0 : papoPolygons[0] = NULL;
890 :
891 0 : if( pbIsValidGeometry )
892 0 : *pbIsValidGeometry = TRUE;
893 :
894 0 : return geom;
895 : }
896 :
897 70 : if (CSLTestBoolean(CPLGetConfigOption("OGR_DEBUG_ORGANIZE_POLYGONS", "NO")))
898 : {
899 : /* -------------------------------------------------------------------- */
900 : /* A wee bit of a warning. */
901 : /* -------------------------------------------------------------------- */
902 : static int firstTime = 1;
903 0 : if (!haveGEOS() && firstTime)
904 : {
905 : CPLDebug("OGR",
906 : "In OGR_DEBUG_ORGANIZE_POLYGONS mode, GDAL should be built with GEOS support enabled in order "
907 0 : "OGRGeometryFactory::organizePolygons to provide reliable results on complex polygons.");
908 0 : firstTime = 0;
909 : }
910 0 : bUseFastVersion = !haveGEOS();
911 : }
912 : else
913 : {
914 70 : bUseFastVersion = TRUE;
915 : }
916 :
917 : /* -------------------------------------------------------------------- */
918 : /* Setup per polygon envelope and area information. */
919 : /* -------------------------------------------------------------------- */
920 70 : sPolyExtended* asPolyEx = new sPolyExtended[nPolygonCount];
921 :
922 70 : int go_on = TRUE;
923 70 : int bMixedUpGeometries = FALSE;
924 70 : int bNonPolygon = FALSE;
925 70 : int bFoundCCW = FALSE;
926 :
927 70 : const char* pszMethodValue = CSLFetchNameValue( (char**)papszOptions, "METHOD" );
928 70 : const char* pszMethodValueOption = CPLGetConfigOption("OGR_ORGANIZE_POLYGONS", NULL);
929 126 : if (pszMethodValueOption != NULL && pszMethodValueOption[0] != '\0')
930 1 : pszMethodValue = pszMethodValueOption;
931 :
932 70 : if (pszMethodValue != NULL)
933 : {
934 17 : if (EQUAL(pszMethodValue, "SKIP"))
935 : {
936 0 : method = METHOD_SKIP;
937 0 : bMixedUpGeometries = TRUE;
938 : }
939 17 : else if (EQUAL(pszMethodValue, "ONLY_CCW"))
940 : {
941 16 : method = METHOD_ONLY_CCW;
942 : }
943 1 : else if (!EQUAL(pszMethodValue, "DEFAULT"))
944 : {
945 : CPLError(CE_Warning, CPLE_AppDefined,
946 0 : "Unrecognized value for METHOD option : %s", pszMethodValue);
947 : }
948 : }
949 :
950 70 : int nCountCWPolygon = 0;
951 70 : int indexOfCWPolygon = -1;
952 :
953 353 : for(i=0;i<nPolygonCount;i++)
954 : {
955 283 : asPolyEx[i].nInitialIndex = i;
956 283 : asPolyEx[i].poPolygon = (OGRPolygon*)papoPolygons[i];
957 283 : papoPolygons[i]->getEnvelope(&asPolyEx[i].sEnvelope);
958 :
959 849 : if( wkbFlatten(papoPolygons[i]->getGeometryType()) == wkbPolygon
960 283 : && ((OGRPolygon *) papoPolygons[i])->getNumInteriorRings() == 0
961 283 : && ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints() >= 4)
962 : {
963 283 : asPolyEx[i].dfArea = asPolyEx[i].poPolygon->get_Area();
964 283 : asPolyEx[i].poExteriorRing = asPolyEx[i].poPolygon->getExteriorRing();
965 283 : asPolyEx[i].poExteriorRing->getPoint(0, &asPolyEx[i].poAPoint);
966 283 : asPolyEx[i].bIsCW = asPolyEx[i].poExteriorRing->isClockwise();
967 283 : if (asPolyEx[i].bIsCW)
968 : {
969 213 : indexOfCWPolygon = i;
970 213 : nCountCWPolygon ++;
971 : }
972 283 : if (!bFoundCCW)
973 163 : bFoundCCW = ! (asPolyEx[i].bIsCW);
974 : }
975 : else
976 : {
977 0 : if( !bMixedUpGeometries )
978 : {
979 : CPLError(
980 : CE_Warning, CPLE_AppDefined,
981 : "organizePolygons() received an unexpected geometry.\n"
982 : "Either a polygon with interior rings, or a polygon with less than 4 points,\n"
983 0 : "or a non-Polygon geometry. Return arguments as a collection." );
984 0 : bMixedUpGeometries = TRUE;
985 : }
986 0 : if( wkbFlatten(papoPolygons[i]->getGeometryType()) != wkbPolygon )
987 0 : bNonPolygon = TRUE;
988 : }
989 : }
990 :
991 : /* If we are in ONLY_CCW mode and that we have found that there is only one outer ring, */
992 : /* then it is pretty easy : we can assume that all other rings are inside */
993 70 : if (method == METHOD_ONLY_CCW && nCountCWPolygon == 1 && bUseFastVersion)
994 : {
995 3 : geom = asPolyEx[indexOfCWPolygon].poPolygon;
996 10 : for(i=0; i<nPolygonCount; i++)
997 : {
998 7 : if (i != indexOfCWPolygon)
999 : {
1000 4 : ((OGRPolygon*)geom)->addRing(asPolyEx[i].poPolygon->getExteriorRing());
1001 4 : delete asPolyEx[i].poPolygon;
1002 : }
1003 : }
1004 3 : delete [] asPolyEx;
1005 3 : if (pbIsValidGeometry)
1006 3 : *pbIsValidGeometry = TRUE;
1007 3 : return geom;
1008 : }
1009 :
1010 : /* Emits a warning if the number of parts is sufficiently big to anticipate for */
1011 : /* very long computation time, and the user didn't specify an explicit method */
1012 67 : if (nPolygonCount > N_CRITICAL_PART_NUMBER && method == METHOD_NORMAL && pszMethodValue == NULL)
1013 : {
1014 : static int firstTime = 1;
1015 0 : if (firstTime)
1016 : {
1017 0 : if (bFoundCCW)
1018 : {
1019 : CPLError( CE_Warning, CPLE_AppDefined,
1020 : "organizePolygons() received a polygon with more than %d parts. "
1021 : "The processing may be really slow.\n"
1022 : "You can skip the processing by setting METHOD=SKIP, "
1023 : "or only make it analyze counter-clock wise parts by setting "
1024 : "METHOD=ONLY_CCW if you can assume that the "
1025 0 : "outline of holes is counter-clock wise defined", N_CRITICAL_PART_NUMBER);
1026 : }
1027 : else
1028 : {
1029 : CPLError( CE_Warning, CPLE_AppDefined,
1030 : "organizePolygons() received a polygon with more than %d parts. "
1031 : "The processing may be really slow.\n"
1032 : "You can skip the processing by setting METHOD=SKIP.",
1033 0 : N_CRITICAL_PART_NUMBER);
1034 : }
1035 0 : firstTime = 0;
1036 : }
1037 : }
1038 :
1039 :
1040 : /* This a several steps algorithm :
1041 : 1) Sort polygons by descending areas
1042 : 2) For each polygon of rank i, find its smallest enclosing polygon
1043 : among the polygons of rank [i-1 ... 0]. If there are no such polygon,
1044 : this is a toplevel polygon. Otherwise, depending on if the enclosing
1045 : polygon is toplevel or not, we can decide if we are toplevel or not
1046 : 3) Re-sort the polygons to retrieve their inital order (nicer for some applications)
1047 : 4) For each non toplevel polygon (= inner ring), add it to its outer ring
1048 : 5) Add the toplevel polygons to the multipolygon
1049 :
1050 : Complexity : O(nPolygonCount^2)
1051 : */
1052 :
1053 : /* Compute how each polygon relate to the other ones
1054 : To save a bit of computation we always begin the computation by a test
1055 : on the enveloppe. We also take into account the areas to avoid some
1056 : useless tests. (A contains B implies envelop(A) contains envelop(B)
1057 : and area(A) > area(B)) In practise, we can hope that few full geometry
1058 : intersection of inclusion test is done:
1059 : * if the polygons are well separated geographically (a set of islands
1060 : for example), no full geometry intersection or inclusion test is done.
1061 : (the envelopes don't intersect each other)
1062 :
1063 : * if the polygons are 'lake inside an island inside a lake inside an
1064 : area' and that each polygon is much smaller than its enclosing one,
1065 : their bounding boxes are stricly contained into each oter, and thus,
1066 : no full geometry intersection or inclusion test is done
1067 : */
1068 :
1069 67 : if (!bMixedUpGeometries)
1070 : {
1071 : /* STEP 1 : Sort polygons by descending area */
1072 67 : qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareArea);
1073 : }
1074 67 : papoPolygons = NULL; /* just to use to avoid it afterwards */
1075 :
1076 : /* -------------------------------------------------------------------- */
1077 : /* Compute relationships, if things seem well structured. */
1078 : /* -------------------------------------------------------------------- */
1079 :
1080 : /* The first (largest) polygon is necessarily top-level */
1081 67 : asPolyEx[0].bIsTopLevel = TRUE;
1082 67 : asPolyEx[0].poEnclosingPolygon = NULL;
1083 :
1084 67 : int nCountTopLevel = 1;
1085 :
1086 : /* STEP 2 */
1087 276 : for(i=1; !bMixedUpGeometries && go_on && i<nPolygonCount; i++)
1088 : {
1089 209 : if (method == METHOD_ONLY_CCW && asPolyEx[i].bIsCW)
1090 : {
1091 20 : nCountTopLevel ++;
1092 20 : asPolyEx[i].bIsTopLevel = TRUE;
1093 20 : asPolyEx[i].poEnclosingPolygon = NULL;
1094 20 : continue;
1095 : }
1096 :
1097 409 : for(j=i-1; go_on && j>=0;j--)
1098 : {
1099 373 : int b_i_inside_j = FALSE;
1100 :
1101 373 : if (method == METHOD_ONLY_CCW && asPolyEx[j].bIsCW == FALSE)
1102 : {
1103 : /* In that mode, i which is CCW if we reach here can only be */
1104 : /* included in a CW polygon */
1105 5 : continue;
1106 : }
1107 :
1108 368 : if (asPolyEx[j].sEnvelope.Contains(asPolyEx[i].sEnvelope))
1109 : {
1110 154 : if (bUseFastVersion)
1111 : {
1112 : /* Note that isPointInRing only test strict inclusion in the ring */
1113 154 : if (asPolyEx[j].poExteriorRing->isPointInRing(&asPolyEx[i].poAPoint, FALSE))
1114 : {
1115 152 : b_i_inside_j = TRUE;
1116 : }
1117 2 : else if (asPolyEx[j].poExteriorRing->isPointOnRingBoundary(&asPolyEx[i].poAPoint, FALSE))
1118 : {
1119 : /* If the point of i is on the boundary of j, we will iterate over the other points of i */
1120 2 : int k, nPoints = asPolyEx[i].poExteriorRing->getNumPoints();
1121 2 : for(k=1;k<nPoints;k++)
1122 : {
1123 6 : OGRPoint point;
1124 6 : asPolyEx[i].poExteriorRing->getPoint(k, &point);
1125 6 : if (asPolyEx[j].poExteriorRing->isPointInRing(&point, FALSE))
1126 : {
1127 : /* If then point is strictly included in j, then i is considered inside j */
1128 1 : b_i_inside_j = TRUE;
1129 : break;
1130 : }
1131 5 : else if (asPolyEx[j].poExteriorRing->isPointOnRingBoundary(&point, FALSE))
1132 : {
1133 : /* If it is on the boundary of j, iterate again */
1134 : }
1135 : else
1136 : {
1137 : /* If it is outside, then i cannot be inside j */
1138 : break;
1139 : }
1140 : }
1141 : }
1142 : }
1143 0 : else if (asPolyEx[j].poPolygon->Contains(asPolyEx[i].poPolygon))
1144 : {
1145 0 : b_i_inside_j = TRUE;
1146 : }
1147 : }
1148 :
1149 :
1150 368 : if (b_i_inside_j)
1151 : {
1152 153 : if (asPolyEx[j].bIsTopLevel)
1153 : {
1154 : /* We are a lake */
1155 119 : asPolyEx[i].bIsTopLevel = FALSE;
1156 119 : asPolyEx[i].poEnclosingPolygon = asPolyEx[j].poPolygon;
1157 : }
1158 : else
1159 : {
1160 : /* We are included in a something not toplevel (a lake), */
1161 : /* so in OGCSF we are considered as toplevel too */
1162 34 : nCountTopLevel ++;
1163 34 : asPolyEx[i].bIsTopLevel = TRUE;
1164 34 : asPolyEx[i].poEnclosingPolygon = NULL;
1165 : }
1166 153 : break;
1167 : }
1168 : /* We use Overlaps instead of Intersects to be more
1169 : tolerant about touching polygons */
1170 215 : else if ( bUseFastVersion || !asPolyEx[i].sEnvelope.Intersects(asPolyEx[j].sEnvelope)
1171 0 : || !asPolyEx[i].poPolygon->Overlaps(asPolyEx[j].poPolygon) )
1172 : {
1173 :
1174 : }
1175 : else
1176 : {
1177 : /* Bad... The polygons are intersecting but no one is
1178 : contained inside the other one. This is a really broken
1179 : case. We just make a multipolygon with the whole set of
1180 : polygons */
1181 0 : go_on = FALSE;
1182 : #ifdef DEBUG
1183 : char* wkt1;
1184 : char* wkt2;
1185 : asPolyEx[i].poPolygon->exportToWkt(&wkt1);
1186 : asPolyEx[j].poPolygon->exportToWkt(&wkt2);
1187 : CPLDebug( "OGR",
1188 : "Bad intersection for polygons %d and %d\n"
1189 : "geom %d: %s\n"
1190 : "geom %d: %s",
1191 : i, j, i, wkt1, j, wkt2 );
1192 : CPLFree(wkt1);
1193 : CPLFree(wkt2);
1194 : #endif
1195 : }
1196 : }
1197 :
1198 189 : if (j < 0)
1199 : {
1200 : /* We come here because we are not included in anything */
1201 : /* We are toplevel */
1202 36 : nCountTopLevel ++;
1203 36 : asPolyEx[i].bIsTopLevel = TRUE;
1204 36 : asPolyEx[i].poEnclosingPolygon = NULL;
1205 : }
1206 : }
1207 :
1208 67 : if (pbIsValidGeometry)
1209 67 : *pbIsValidGeometry = go_on && !bMixedUpGeometries;
1210 :
1211 : /* -------------------------------------------------------------------- */
1212 : /* Things broke down - just turn everything into a multipolygon. */
1213 : /* -------------------------------------------------------------------- */
1214 67 : if ( !go_on || bMixedUpGeometries )
1215 : {
1216 0 : if( bNonPolygon )
1217 0 : geom = new OGRGeometryCollection();
1218 : else
1219 0 : geom = new OGRMultiPolygon();
1220 :
1221 0 : for( i=0; i < nPolygonCount; i++ )
1222 : {
1223 : ((OGRGeometryCollection*)geom)->
1224 0 : addGeometryDirectly( asPolyEx[i].poPolygon );
1225 : }
1226 : }
1227 :
1228 : /* -------------------------------------------------------------------- */
1229 : /* Try to turn into one or more polygons based on the ring */
1230 : /* relationships. */
1231 : /* -------------------------------------------------------------------- */
1232 : else
1233 : {
1234 : /* STEP 3: Resort in initial order */
1235 67 : qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareByIndex);
1236 :
1237 : /* STEP 4: Add holes as rings of their enclosing polygon */
1238 343 : for(i=0;i<nPolygonCount;i++)
1239 : {
1240 276 : if (asPolyEx[i].bIsTopLevel == FALSE)
1241 : {
1242 119 : asPolyEx[i].poEnclosingPolygon->addRing(
1243 238 : asPolyEx[i].poPolygon->getExteriorRing());
1244 119 : delete asPolyEx[i].poPolygon;
1245 : }
1246 157 : else if (nCountTopLevel == 1)
1247 : {
1248 19 : geom = asPolyEx[i].poPolygon;
1249 : }
1250 : }
1251 :
1252 : /* STEP 5: Add toplevel polygons */
1253 67 : if (nCountTopLevel > 1)
1254 : {
1255 284 : for(i=0;i<nPolygonCount;i++)
1256 : {
1257 236 : if (asPolyEx[i].bIsTopLevel)
1258 : {
1259 138 : if (geom == NULL)
1260 : {
1261 48 : geom = new OGRMultiPolygon();
1262 48 : ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon);
1263 : }
1264 : else
1265 : {
1266 90 : ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon);
1267 : }
1268 : }
1269 : }
1270 : }
1271 : }
1272 :
1273 67 : delete[] asPolyEx;
1274 :
1275 67 : return geom;
1276 : }
1277 :
1278 : /************************************************************************/
1279 : /* createFromGML() */
1280 : /************************************************************************/
1281 :
1282 : /**
1283 : * \brief Create geometry from GML.
1284 : *
1285 : * This method translates a fragment of GML containing only the geometry
1286 : * portion into a corresponding OGRGeometry. There are many limitations
1287 : * on the forms of GML geometries supported by this parser, but they are
1288 : * too numerous to list here.
1289 : *
1290 : * The C function OGR_G_CreateFromGML() is the same as this method.
1291 : *
1292 : * @param pszData The GML fragment for the geometry.
1293 : *
1294 : * @return a geometry on succes, or NULL on error.
1295 : */
1296 :
1297 52 : OGRGeometry *OGRGeometryFactory::createFromGML( const char *pszData )
1298 :
1299 : {
1300 : OGRGeometryH hGeom;
1301 :
1302 52 : hGeom = OGR_G_CreateFromGML( pszData );
1303 :
1304 52 : return (OGRGeometry *) hGeom;
1305 : }
1306 :
1307 : /************************************************************************/
1308 : /* createFromGEOS() */
1309 : /************************************************************************/
1310 :
1311 : OGRGeometry *
1312 23 : OGRGeometryFactory::createFromGEOS( GEOSGeom geosGeom )
1313 :
1314 : {
1315 : #ifndef HAVE_GEOS
1316 :
1317 : CPLError( CE_Failure, CPLE_NotSupported,
1318 : "GEOS support not enabled." );
1319 : return NULL;
1320 :
1321 : #else
1322 :
1323 23 : size_t nSize = 0;
1324 23 : unsigned char *pabyBuf = NULL;
1325 23 : OGRGeometry *poGeometry = NULL;
1326 :
1327 23 : pabyBuf = GEOSGeomToWKB_buf( geosGeom, &nSize );
1328 23 : if( pabyBuf == NULL || nSize == 0 )
1329 : {
1330 0 : return NULL;
1331 : }
1332 :
1333 23 : if( OGRGeometryFactory::createFromWkb( (unsigned char *) pabyBuf,
1334 : NULL, &poGeometry, (int) nSize )
1335 : != OGRERR_NONE )
1336 : {
1337 0 : poGeometry = NULL;
1338 : }
1339 :
1340 23 : if( pabyBuf != NULL )
1341 : {
1342 : #if GEOS_CAPI_VERSION_MAJOR >= 2 || (GEOS_CAPI_VERSION_MAJOR == 1 && GEOS_CAPI_VERSION_MINOR >= 6)
1343 23 : GEOSFree( pabyBuf );
1344 : #else
1345 : free( pabyBuf );
1346 : #endif
1347 : }
1348 :
1349 23 : return poGeometry;
1350 :
1351 : #endif /* HAVE_GEOS */
1352 : }
1353 :
1354 : /************************************************************************/
1355 : /* getGEOSGeometryFactory() */
1356 : /************************************************************************/
1357 :
1358 0 : void *OGRGeometryFactory::getGEOSGeometryFactory()
1359 :
1360 : {
1361 : // XXX - mloskot - What to do with this call
1362 : // after GEOS C++ API has been stripped?
1363 0 : return NULL;
1364 : }
1365 :
1366 : /************************************************************************/
1367 : /* haveGEOS() */
1368 : /************************************************************************/
1369 :
1370 : /**
1371 : * \brief Test if GEOS enabled.
1372 : *
1373 : * This static method returns TRUE if GEOS support is built into OGR,
1374 : * otherwise it returns FALSE.
1375 : *
1376 : * @return TRUE if available, otherwise FALSE.
1377 : */
1378 :
1379 34 : int OGRGeometryFactory::haveGEOS()
1380 :
1381 : {
1382 : #ifndef HAVE_GEOS
1383 : return FALSE;
1384 : #else
1385 34 : return TRUE;
1386 : #endif
1387 : }
1388 :
1389 : /************************************************************************/
1390 : /* createFromFgf() */
1391 : /************************************************************************/
1392 :
1393 : /**
1394 : * \brief Create a geometry object of the appropriate type from it's FGF (FDO Geometry Format) binary representation.
1395 : *
1396 : * Also note that this is a static method, and that there
1397 : * is no need to instantiate an OGRGeometryFactory object.
1398 : *
1399 : * The C function OGR_G_CreateFromFgf() is the same as this method.
1400 : *
1401 : * @param pabyData pointer to the input BLOB data.
1402 : * @param poSR pointer to the spatial reference to be assigned to the
1403 : * created geometry object. This may be NULL.
1404 : * @param ppoReturn the newly created geometry object will be assigned to the
1405 : * indicated pointer on return. This will be NULL in case
1406 : * of failure.
1407 : * @param nBytes the number of bytes available in pabyData.
1408 : * @param pnBytesConsumed if not NULL, it will be set to the number of bytes
1409 : * consumed (at most nBytes).
1410 : *
1411 : * @return OGRERR_NONE if all goes well, otherwise any of
1412 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
1413 : * OGRERR_CORRUPT_DATA may be returned.
1414 : */
1415 :
1416 10 : OGRErr OGRGeometryFactory::createFromFgf( unsigned char *pabyData,
1417 : OGRSpatialReference * poSR,
1418 : OGRGeometry **ppoReturn,
1419 : int nBytes,
1420 : int *pnBytesConsumed )
1421 :
1422 : {
1423 10 : OGRErr eErr = OGRERR_NONE;
1424 10 : OGRGeometry *poGeom = NULL;
1425 : GInt32 nGType, nGDim;
1426 10 : int nTupleSize = 0;
1427 10 : int iOrdinal = 0;
1428 :
1429 : (void) iOrdinal;
1430 :
1431 10 : *ppoReturn = NULL;
1432 :
1433 10 : if( nBytes < 4 )
1434 0 : return OGRERR_NOT_ENOUGH_DATA;
1435 :
1436 : /* -------------------------------------------------------------------- */
1437 : /* Decode the geometry type. */
1438 : /* -------------------------------------------------------------------- */
1439 10 : memcpy( &nGType, pabyData + 0, 4 );
1440 : CPL_LSBPTR32( &nGType );
1441 :
1442 10 : if( nGType < 0 || nGType > 13 )
1443 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
1444 :
1445 : /* -------------------------------------------------------------------- */
1446 : /* Decode the dimentionality if appropriate. */
1447 : /* -------------------------------------------------------------------- */
1448 10 : switch( nGType )
1449 : {
1450 : case 1: // Point
1451 : case 2: // LineString
1452 : case 3: // Polygon
1453 :
1454 8 : if( nBytes < 8 )
1455 0 : return OGRERR_NOT_ENOUGH_DATA;
1456 :
1457 8 : memcpy( &nGDim, pabyData + 4, 4 );
1458 : CPL_LSBPTR32( &nGDim );
1459 :
1460 8 : if( nGDim < 0 || nGDim > 3 )
1461 0 : return OGRERR_CORRUPT_DATA;
1462 :
1463 8 : nTupleSize = 2;
1464 8 : if( nGDim & 0x01 ) // Z
1465 1 : nTupleSize++;
1466 8 : if( nGDim & 0x02 ) // M
1467 0 : nTupleSize++;
1468 :
1469 : break;
1470 :
1471 : default:
1472 : break;
1473 : }
1474 :
1475 : /* -------------------------------------------------------------------- */
1476 : /* None */
1477 : /* -------------------------------------------------------------------- */
1478 10 : if( nGType == 0 )
1479 : {
1480 0 : if( pnBytesConsumed )
1481 0 : *pnBytesConsumed = 4;
1482 : }
1483 :
1484 : /* -------------------------------------------------------------------- */
1485 : /* Point */
1486 : /* -------------------------------------------------------------------- */
1487 10 : else if( nGType == 1 )
1488 : {
1489 : double adfTuple[4];
1490 :
1491 2 : if( nBytes < nTupleSize * 8 + 8 )
1492 0 : return OGRERR_NOT_ENOUGH_DATA;
1493 :
1494 2 : memcpy( adfTuple, pabyData + 8, nTupleSize*8 );
1495 : #ifdef CPL_MSB
1496 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1497 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1498 : #endif
1499 2 : if( nTupleSize > 2 )
1500 1 : poGeom = new OGRPoint( adfTuple[0], adfTuple[1], adfTuple[2] );
1501 : else
1502 1 : poGeom = new OGRPoint( adfTuple[0], adfTuple[1] );
1503 :
1504 2 : if( pnBytesConsumed )
1505 0 : *pnBytesConsumed = 8 + nTupleSize * 8;
1506 : }
1507 :
1508 : /* -------------------------------------------------------------------- */
1509 : /* LineString */
1510 : /* -------------------------------------------------------------------- */
1511 8 : else if( nGType == 2 )
1512 : {
1513 : double adfTuple[4];
1514 : GInt32 nPointCount;
1515 : int iPoint;
1516 : OGRLineString *poLS;
1517 :
1518 2 : if( nBytes < 12 )
1519 0 : return OGRERR_NOT_ENOUGH_DATA;
1520 :
1521 2 : memcpy( &nPointCount, pabyData + 8, 4 );
1522 : CPL_LSBPTR32( &nPointCount );
1523 :
1524 2 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
1525 0 : return OGRERR_CORRUPT_DATA;
1526 :
1527 2 : if( nBytes - 12 < nTupleSize * 8 * nPointCount )
1528 0 : return OGRERR_NOT_ENOUGH_DATA;
1529 :
1530 2 : poGeom = poLS = new OGRLineString();
1531 2 : poLS->setNumPoints( nPointCount );
1532 :
1533 4 : for( iPoint = 0; iPoint < nPointCount; iPoint++ )
1534 : {
1535 : memcpy( adfTuple, pabyData + 12 + 8*nTupleSize*iPoint,
1536 2 : nTupleSize*8 );
1537 : #ifdef CPL_MSB
1538 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1539 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1540 : #endif
1541 2 : if( nTupleSize > 2 )
1542 0 : poLS->setPoint( iPoint, adfTuple[0], adfTuple[1], adfTuple[2] );
1543 : else
1544 2 : poLS->setPoint( iPoint, adfTuple[0], adfTuple[1] );
1545 : }
1546 :
1547 2 : if( pnBytesConsumed )
1548 0 : *pnBytesConsumed = 12 + nTupleSize * 8 * nPointCount;
1549 : }
1550 :
1551 : /* -------------------------------------------------------------------- */
1552 : /* Polygon */
1553 : /* -------------------------------------------------------------------- */
1554 6 : else if( nGType == 3 )
1555 : {
1556 : double adfTuple[4];
1557 : GInt32 nPointCount;
1558 : GInt32 nRingCount;
1559 : int iPoint, iRing;
1560 : OGRLinearRing *poLR;
1561 : OGRPolygon *poPoly;
1562 : int nNextByte;
1563 :
1564 4 : if( nBytes < 12 )
1565 0 : return OGRERR_NOT_ENOUGH_DATA;
1566 :
1567 4 : memcpy( &nRingCount, pabyData + 8, 4 );
1568 : CPL_LSBPTR32( &nRingCount );
1569 :
1570 4 : if (nRingCount < 0 || nRingCount > INT_MAX / 4)
1571 0 : return OGRERR_CORRUPT_DATA;
1572 :
1573 : /* Each ring takes at least 4 bytes */
1574 4 : if (nBytes - 12 < nRingCount * 4)
1575 0 : return OGRERR_NOT_ENOUGH_DATA;
1576 :
1577 4 : nNextByte = 12;
1578 :
1579 4 : poGeom = poPoly = new OGRPolygon();
1580 :
1581 10 : for( iRing = 0; iRing < nRingCount; iRing++ )
1582 : {
1583 6 : if( nBytes - nNextByte < 4 )
1584 0 : return OGRERR_NOT_ENOUGH_DATA;
1585 :
1586 6 : memcpy( &nPointCount, pabyData + nNextByte, 4 );
1587 : CPL_LSBPTR32( &nPointCount );
1588 :
1589 6 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
1590 0 : return OGRERR_CORRUPT_DATA;
1591 :
1592 6 : nNextByte += 4;
1593 :
1594 6 : if( nBytes - nNextByte < nTupleSize * 8 * nPointCount )
1595 0 : return OGRERR_NOT_ENOUGH_DATA;
1596 :
1597 6 : poLR = new OGRLinearRing();
1598 6 : poLR->setNumPoints( nPointCount );
1599 :
1600 12 : for( iPoint = 0; iPoint < nPointCount; iPoint++ )
1601 : {
1602 6 : memcpy( adfTuple, pabyData + nNextByte, nTupleSize*8 );
1603 6 : nNextByte += nTupleSize * 8;
1604 :
1605 : #ifdef CPL_MSB
1606 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1607 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1608 : #endif
1609 6 : if( nTupleSize > 2 )
1610 0 : poLR->setPoint( iPoint, adfTuple[0], adfTuple[1], adfTuple[2] );
1611 : else
1612 6 : poLR->setPoint( iPoint, adfTuple[0], adfTuple[1] );
1613 : }
1614 :
1615 6 : poPoly->addRingDirectly( poLR );
1616 : }
1617 :
1618 4 : if( pnBytesConsumed )
1619 2 : *pnBytesConsumed = nNextByte;
1620 : }
1621 :
1622 : /* -------------------------------------------------------------------- */
1623 : /* GeometryCollections of various kinds. */
1624 : /* -------------------------------------------------------------------- */
1625 4 : else if( nGType == 4 // MultiPoint
1626 : || nGType == 5 // MultiLineString
1627 : || nGType == 6 // MultiPolygon
1628 : || nGType == 7 ) // MultiGeometry
1629 : {
1630 2 : OGRGeometryCollection *poGC = NULL;
1631 : GInt32 nGeomCount;
1632 : int iGeom, nBytesUsed;
1633 :
1634 2 : if( nGType == 4 )
1635 0 : poGC = new OGRMultiPoint();
1636 2 : else if( nGType == 5 )
1637 0 : poGC = new OGRMultiLineString();
1638 2 : else if( nGType == 6 )
1639 0 : poGC = new OGRMultiPolygon();
1640 2 : else if( nGType == 7 )
1641 2 : poGC = new OGRGeometryCollection();
1642 :
1643 2 : if( nBytes < 8 )
1644 0 : return OGRERR_NOT_ENOUGH_DATA;
1645 :
1646 2 : memcpy( &nGeomCount, pabyData + 4, 4 );
1647 : CPL_LSBPTR32( &nGeomCount );
1648 :
1649 2 : if (nGeomCount < 0 || nGeomCount > INT_MAX / 4)
1650 0 : return OGRERR_CORRUPT_DATA;
1651 :
1652 : /* Each geometry takes at least 4 bytes */
1653 2 : if (nBytes - 8 < 4 * nGeomCount)
1654 0 : return OGRERR_NOT_ENOUGH_DATA;
1655 :
1656 2 : nBytesUsed = 8;
1657 :
1658 4 : for( iGeom = 0; iGeom < nGeomCount; iGeom++ )
1659 : {
1660 : int nThisGeomSize;
1661 2 : OGRGeometry *poThisGeom = NULL;
1662 :
1663 : eErr = createFromFgf( pabyData + nBytesUsed, poSR, &poThisGeom,
1664 2 : nBytes - nBytesUsed, &nThisGeomSize);
1665 2 : if( eErr != OGRERR_NONE )
1666 : {
1667 0 : delete poGC;
1668 0 : return eErr;
1669 : }
1670 :
1671 2 : nBytesUsed += nThisGeomSize;
1672 2 : eErr = poGC->addGeometryDirectly( poThisGeom );
1673 2 : if( eErr != OGRERR_NONE )
1674 : {
1675 0 : delete poGC;
1676 0 : return eErr;
1677 : }
1678 : }
1679 :
1680 2 : poGeom = poGC;
1681 2 : if( pnBytesConsumed )
1682 0 : *pnBytesConsumed = nBytesUsed;
1683 : }
1684 :
1685 : /* -------------------------------------------------------------------- */
1686 : /* Currently unsupported geometry. */
1687 : /* */
1688 : /* We need to add 10/11/12/13 curve types in some fashion. */
1689 : /* -------------------------------------------------------------------- */
1690 : else
1691 : {
1692 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
1693 : }
1694 :
1695 : /* -------------------------------------------------------------------- */
1696 : /* Assign spatial reference system. */
1697 : /* -------------------------------------------------------------------- */
1698 10 : if( eErr == OGRERR_NONE )
1699 : {
1700 10 : if( poGeom != NULL && poSR )
1701 0 : poGeom->assignSpatialReference( poSR );
1702 10 : *ppoReturn = poGeom;
1703 : }
1704 : else
1705 : {
1706 0 : delete poGeom;
1707 : }
1708 :
1709 10 : return eErr;
1710 : }
1711 :
1712 : /************************************************************************/
1713 : /* OGR_G_CreateFromFgf() */
1714 : /************************************************************************/
1715 :
1716 0 : OGRErr CPL_DLL OGR_G_CreateFromFgf( unsigned char *pabyData,
1717 : OGRSpatialReferenceH hSRS,
1718 : OGRGeometryH *phGeometry,
1719 : int nBytes, int *pnBytesConsumed )
1720 :
1721 : {
1722 : return OGRGeometryFactory::createFromFgf( pabyData,
1723 : (OGRSpatialReference *) hSRS,
1724 : (OGRGeometry **) phGeometry,
1725 0 : nBytes, pnBytesConsumed );
1726 : }
1727 :
1728 :
1729 : /************************************************************************/
1730 : /* Add360ToNegLon() */
1731 : /************************************************************************/
1732 :
1733 2 : static void Add360ToNegLon( OGRGeometry* poGeom )
1734 : {
1735 2 : switch (wkbFlatten(poGeom->getGeometryType()))
1736 : {
1737 : case wkbPolygon:
1738 : case wkbMultiLineString:
1739 : case wkbMultiPolygon:
1740 : case wkbGeometryCollection:
1741 : {
1742 1 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
1743 2 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
1744 : {
1745 1 : Add360ToNegLon((OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom));
1746 : }
1747 :
1748 1 : break;
1749 : }
1750 :
1751 : case wkbLineString:
1752 : case wkbLinearRing:
1753 : {
1754 1 : OGRLineString* poLineString = (OGRLineString* )poGeom;
1755 1 : int nPointCount = poLineString->getNumPoints();
1756 1 : int nCoordDim = poLineString->getCoordinateDimension();
1757 6 : for( int iPoint = 0; iPoint < nPointCount; iPoint++)
1758 : {
1759 5 : if (poLineString->getX(iPoint) < -90)
1760 : {
1761 2 : if (nCoordDim == 2)
1762 : poLineString->setPoint(iPoint,
1763 : poLineString->getX(iPoint) + 360,
1764 2 : poLineString->getY(iPoint));
1765 : else
1766 : poLineString->setPoint(iPoint,
1767 : poLineString->getX(iPoint) + 360,
1768 : poLineString->getY(iPoint),
1769 0 : poLineString->getZ(iPoint));
1770 : }
1771 : }
1772 : break;
1773 : }
1774 :
1775 : default:
1776 : break;
1777 : }
1778 2 : }
1779 :
1780 : /************************************************************************/
1781 : /* Sub360ToLon() */
1782 : /************************************************************************/
1783 :
1784 2 : static void Sub360ToLon( OGRGeometry* poGeom )
1785 : {
1786 2 : switch (wkbFlatten(poGeom->getGeometryType()))
1787 : {
1788 : case wkbPolygon:
1789 : case wkbMultiLineString:
1790 : case wkbMultiPolygon:
1791 : case wkbGeometryCollection:
1792 : {
1793 1 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
1794 2 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
1795 : {
1796 1 : Sub360ToLon((OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom));
1797 : }
1798 :
1799 1 : break;
1800 : }
1801 :
1802 : case wkbLineString:
1803 : case wkbLinearRing:
1804 : {
1805 1 : OGRLineString* poLineString = (OGRLineString* )poGeom;
1806 1 : int nPointCount = poLineString->getNumPoints();
1807 1 : int nCoordDim = poLineString->getCoordinateDimension();
1808 6 : for( int iPoint = 0; iPoint < nPointCount; iPoint++)
1809 : {
1810 5 : if (nCoordDim == 2)
1811 : poLineString->setPoint(iPoint,
1812 : poLineString->getX(iPoint) - 360,
1813 5 : poLineString->getY(iPoint));
1814 : else
1815 : poLineString->setPoint(iPoint,
1816 : poLineString->getX(iPoint) - 360,
1817 : poLineString->getY(iPoint),
1818 0 : poLineString->getZ(iPoint));
1819 : }
1820 : break;
1821 : }
1822 :
1823 : default:
1824 : break;
1825 : }
1826 2 : }
1827 :
1828 : /************************************************************************/
1829 : /* AddSimpleGeomToMulti() */
1830 : /************************************************************************/
1831 :
1832 2 : static void AddSimpleGeomToMulti(OGRGeometryCollection* poMulti,
1833 : const OGRGeometry* poGeom)
1834 : {
1835 2 : switch (wkbFlatten(poGeom->getGeometryType()))
1836 : {
1837 : case wkbPolygon:
1838 : case wkbLineString:
1839 2 : poMulti->addGeometry(poGeom);
1840 2 : break;
1841 :
1842 : case wkbMultiLineString:
1843 : case wkbMultiPolygon:
1844 : case wkbGeometryCollection:
1845 : {
1846 0 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
1847 0 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
1848 : {
1849 : OGRGeometry* poSubGeom =
1850 0 : (OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom);
1851 0 : AddSimpleGeomToMulti(poMulti, poSubGeom);
1852 : }
1853 : break;
1854 : }
1855 :
1856 : default:
1857 : break;
1858 : }
1859 2 : }
1860 :
1861 : /************************************************************************/
1862 : /* CutGeometryOnDateLineAndAddToMulti() */
1863 : /************************************************************************/
1864 :
1865 1 : static void CutGeometryOnDateLineAndAddToMulti(OGRGeometryCollection* poMulti,
1866 : const OGRGeometry* poGeom)
1867 : {
1868 1 : switch (wkbFlatten(poGeom->getGeometryType()))
1869 : {
1870 : case wkbPolygon:
1871 : case wkbLineString:
1872 : {
1873 1 : int bWrapDateline = FALSE;
1874 1 : OGREnvelope oEnvelope;
1875 :
1876 1 : poGeom->getEnvelope(&oEnvelope);
1877 :
1878 : /* Naive heuristics... Place to improvement... */
1879 1 : OGRGeometry* poDupGeom = NULL;
1880 :
1881 2 : if (oEnvelope.MinX < -170 && oEnvelope.MaxX > 170)
1882 : {
1883 1 : bWrapDateline = TRUE;
1884 1 : poDupGeom = poGeom->clone();
1885 1 : Add360ToNegLon(poDupGeom);
1886 : }
1887 0 : else if (oEnvelope.MinX > 170 && oEnvelope.MaxX > 180)
1888 0 : bWrapDateline = TRUE;
1889 :
1890 1 : if (bWrapDateline)
1891 : {
1892 : #ifndef HAVE_GEOS
1893 : CPLError( CE_Failure, CPLE_NotSupported,
1894 : "GEOS support not enabled." );
1895 :
1896 : poMulti->addGeometry(poGeom);
1897 : #else
1898 1 : const OGRGeometry* poWorkGeom = (poDupGeom) ? poDupGeom : poGeom;
1899 1 : OGRGeometry* poRectangle1 = NULL;
1900 1 : OGRGeometry* poRectangle2 = NULL;
1901 1 : const char* pszWKT1 = "POLYGON((0 90,180 90,180 -90,0 -90,0 90))";
1902 1 : const char* pszWKT2 = "POLYGON((180 90,360 90,360 -90,180 -90,180 90))";
1903 1 : OGRGeometryFactory::createFromWkt((char**)&pszWKT1, NULL, &poRectangle1);
1904 1 : OGRGeometryFactory::createFromWkt((char**)&pszWKT2, NULL, &poRectangle2);
1905 1 : OGRGeometry* poGeom1 = poWorkGeom->Intersection(poRectangle1);
1906 1 : OGRGeometry* poGeom2 = poWorkGeom->Intersection(poRectangle2);
1907 1 : delete poRectangle1;
1908 1 : delete poRectangle2;
1909 :
1910 2 : if (poGeom1 != NULL && poGeom2 != NULL)
1911 : {
1912 1 : AddSimpleGeomToMulti(poMulti, poGeom1);
1913 1 : Sub360ToLon(poGeom2);
1914 1 : AddSimpleGeomToMulti(poMulti, poGeom2);
1915 : }
1916 : else
1917 : {
1918 0 : AddSimpleGeomToMulti(poMulti, poGeom);
1919 : }
1920 :
1921 1 : delete poGeom1;
1922 1 : delete poGeom2;
1923 1 : delete poDupGeom;
1924 : #endif
1925 : }
1926 : else
1927 : {
1928 0 : poMulti->addGeometry(poGeom);
1929 : }
1930 1 : break;
1931 : }
1932 :
1933 : case wkbMultiLineString:
1934 : case wkbMultiPolygon:
1935 : case wkbGeometryCollection:
1936 : {
1937 0 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
1938 0 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
1939 : {
1940 : OGRGeometry* poSubGeom =
1941 0 : (OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom);
1942 0 : CutGeometryOnDateLineAndAddToMulti(poMulti, poSubGeom);
1943 : }
1944 : break;
1945 : }
1946 :
1947 : default:
1948 : break;
1949 : }
1950 1 : }
1951 :
1952 : /************************************************************************/
1953 : /* transformWithOptions() */
1954 : /************************************************************************/
1955 :
1956 11 : OGRGeometry* OGRGeometryFactory::transformWithOptions( const OGRGeometry* poSrcGeom,
1957 : OGRCoordinateTransformation *poCT,
1958 : char** papszOptions )
1959 : {
1960 11 : OGRGeometry* poDstGeom = poSrcGeom->clone();
1961 11 : OGRErr eErr = poDstGeom->transform(poCT);
1962 11 : if (eErr != OGRERR_NONE)
1963 : {
1964 0 : delete poDstGeom;
1965 0 : return NULL;
1966 : }
1967 :
1968 11 : if (CSLTestBoolean(CSLFetchNameValueDef(papszOptions, "WRAPDATELINE", "NO")))
1969 : {
1970 1 : OGRwkbGeometryType eType = wkbFlatten(poSrcGeom->getGeometryType());
1971 : OGRwkbGeometryType eNewType;
1972 2 : if (eType == wkbPolygon || eType == wkbMultiPolygon)
1973 1 : eNewType = wkbMultiPolygon;
1974 0 : else if (eType == wkbLineString || eType == wkbMultiLineString)
1975 0 : eNewType = wkbMultiLineString;
1976 : else
1977 0 : eNewType = wkbGeometryCollection;
1978 :
1979 : OGRGeometryCollection* poMulti =
1980 1 : (OGRGeometryCollection* )createGeometry(eNewType);
1981 :
1982 1 : CutGeometryOnDateLineAndAddToMulti(poMulti, poDstGeom);
1983 :
1984 1 : if (poMulti->getNumGeometries() == 0)
1985 : {
1986 0 : delete poMulti;
1987 : }
1988 1 : else if (poMulti->getNumGeometries() == 1)
1989 : {
1990 0 : delete poDstGeom;
1991 0 : poDstGeom = poMulti->getGeometryRef(0)->clone();
1992 0 : delete poMulti;
1993 : }
1994 : else
1995 : {
1996 1 : delete poDstGeom;
1997 1 : poDstGeom = poMulti;
1998 : }
1999 : }
2000 :
2001 11 : return poDstGeom;
2002 : }
2003 :
2004 : /************************************************************************/
2005 : /* approximateArcAngles() */
2006 : /************************************************************************/
2007 :
2008 : /**
2009 : * Stroke arc to linestring.
2010 : *
2011 : * Stroke an arc of a circle to a linestring based on a center
2012 : * point, radius, start angle and end angle, all angles in degrees.
2013 : *
2014 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
2015 : * used. This is currently 4 degrees unless the user has overridden the
2016 : * value with the OGR_ARC_STEPSIZE configuration variable.
2017 : *
2018 : * @see CPLSetConfigOption()
2019 : *
2020 : * @param dfCenterX center X
2021 : * @param dfCenterY center Y
2022 : * @param dfZ center Z
2023 : * @param dfPrimaryRadius X radius of ellipse.
2024 : * @param dfSecondaryRadius Y radius of ellipse.
2025 : * @param dfRotation rotation of the ellipse clockwise.
2026 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
2027 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
2028 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
2029 : * arc, zero to use the default setting.
2030 : *
2031 : * @return OGRLineString geometry representing an approximation of the arc.
2032 : */
2033 :
2034 7 : OGRGeometry* OGRGeometryFactory::approximateArcAngles(
2035 : double dfCenterX, double dfCenterY, double dfZ,
2036 : double dfPrimaryRadius, double dfSecondaryRadius, double dfRotation,
2037 : double dfStartAngle, double dfEndAngle,
2038 : double dfMaxAngleStepSizeDegrees )
2039 :
2040 : {
2041 : double dfSlice;
2042 : int iPoint, nVertexCount;
2043 7 : OGRLineString *poLine = new OGRLineString();
2044 7 : double dfRotationRadians = dfRotation * PI / 180.0;
2045 :
2046 : // support default arc step setting.
2047 7 : if( dfMaxAngleStepSizeDegrees == 0 )
2048 : {
2049 : dfMaxAngleStepSizeDegrees =
2050 6 : atof(CPLGetConfigOption("OGR_ARC_STEPSIZE","4"));
2051 : }
2052 :
2053 : // switch direction
2054 7 : dfStartAngle *= -1;
2055 7 : dfEndAngle *= -1;
2056 :
2057 : // Figure out the number of slices to make this into.
2058 : nVertexCount = (int)
2059 7 : ceil(fabs(dfEndAngle - dfStartAngle)/dfMaxAngleStepSizeDegrees) + 1;
2060 7 : nVertexCount = MAX(2,nVertexCount);
2061 7 : dfSlice = (dfEndAngle-dfStartAngle)/(nVertexCount-1);
2062 :
2063 : /* -------------------------------------------------------------------- */
2064 : /* Compute the interpolated points. */
2065 : /* -------------------------------------------------------------------- */
2066 427 : for( iPoint=0; iPoint < nVertexCount; iPoint++ )
2067 : {
2068 : double dfAngleOnEllipse;
2069 : double dfArcX, dfArcY;
2070 : double dfEllipseX, dfEllipseY;
2071 :
2072 420 : dfAngleOnEllipse = (dfStartAngle + iPoint * dfSlice) * PI / 180.0;
2073 :
2074 : // Compute position on the unrotated ellipse.
2075 420 : dfEllipseX = cos(dfAngleOnEllipse) * dfPrimaryRadius;
2076 420 : dfEllipseY = sin(dfAngleOnEllipse) * dfSecondaryRadius;
2077 :
2078 : // Rotate this position around the center of the ellipse.
2079 : dfArcX = dfCenterX
2080 : + dfEllipseX * cos(dfRotationRadians)
2081 420 : + dfEllipseY * sin(dfRotationRadians);
2082 : dfArcY = dfCenterY
2083 : - dfEllipseX * sin(dfRotationRadians)
2084 420 : + dfEllipseY * cos(dfRotationRadians);
2085 :
2086 420 : poLine->setPoint( iPoint, dfArcX, dfArcY, dfZ );
2087 : }
2088 :
2089 7 : return poLine;
2090 : }
2091 :
2092 : /************************************************************************/
2093 : /* OGR_G_ApproximateArcAngles() */
2094 : /************************************************************************/
2095 :
2096 : OGRGeometryH CPL_DLL
2097 1 : OGR_G_ApproximateArcAngles(
2098 : double dfCenterX, double dfCenterY, double dfZ,
2099 : double dfPrimaryRadius, double dfSecondaryAxis, double dfRotation,
2100 : double dfStartAngle, double dfEndAngle,
2101 : double dfMaxAngleStepSizeDegrees )
2102 :
2103 : {
2104 : return (OGRGeometryH) OGRGeometryFactory::approximateArcAngles(
2105 : dfCenterX, dfCenterY, dfZ,
2106 : dfPrimaryRadius, dfSecondaryAxis, dfRotation,
2107 1 : dfStartAngle, dfEndAngle, dfMaxAngleStepSizeDegrees );
2108 : }
|
