1 : /******************************************************************************
2 : * $Id: ogrgeometryfactory.cpp 23589 2011-12-17 14:21:01Z rouault $
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 23589 2011-12-17 14:21:01Z rouault $");
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 6753 : 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 6753 : *ppoReturn = NULL;
87 :
88 6753 : if( nBytes < 9 && nBytes != -1 )
89 2506 : 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 4247 : eByteOrder = DB2_V72_FIX_BYTE_ORDER((OGRwkbByteOrder) *pabyData);
96 :
97 :
98 4247 : 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 4247 : if( eByteOrder == wkbNDR )
121 4033 : eGeometryType = (OGRwkbGeometryType) pabyData[1];
122 : else
123 214 : eGeometryType = (OGRwkbGeometryType) pabyData[4];
124 :
125 : /* -------------------------------------------------------------------- */
126 : /* Instantiate a geometry of the appropriate type, and */
127 : /* initialize from the input stream. */
128 : /* -------------------------------------------------------------------- */
129 4247 : poGeom = createGeometry( eGeometryType );
130 :
131 4247 : if( poGeom == NULL )
132 56 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
133 :
134 : /* -------------------------------------------------------------------- */
135 : /* Import from binary. */
136 : /* -------------------------------------------------------------------- */
137 4191 : eErr = poGeom->importFromWkb( pabyData, nBytes );
138 :
139 : /* -------------------------------------------------------------------- */
140 : /* Assign spatial reference system. */
141 : /* -------------------------------------------------------------------- */
142 4191 : if( eErr == OGRERR_NONE )
143 : {
144 4187 : poGeom->assignSpatialReference( poSR );
145 4187 : *ppoReturn = poGeom;
146 : }
147 : else
148 : {
149 4 : delete poGeom;
150 : }
151 :
152 4191 : 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 144 : 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 144 : 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;
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 63612 : OGRErr OGRGeometryFactory::createFromWkt(char **ppszData,
235 : OGRSpatialReference * poSR,
236 : OGRGeometry **ppoReturn )
237 :
238 : {
239 : OGRErr eErr;
240 : char szToken[OGR_WKT_TOKEN_MAX];
241 63612 : char *pszInput = *ppszData;
242 : OGRGeometry *poGeom;
243 :
244 63612 : *ppoReturn = NULL;
245 :
246 : /* -------------------------------------------------------------------- */
247 : /* Get the first token, which should be the geometry type. */
248 : /* -------------------------------------------------------------------- */
249 63612 : if( OGRWktReadToken( pszInput, szToken ) == NULL )
250 0 : return OGRERR_CORRUPT_DATA;
251 :
252 : /* -------------------------------------------------------------------- */
253 : /* Instantiate a geometry of the appropriate type. */
254 : /* -------------------------------------------------------------------- */
255 63612 : if( EQUAL(szToken,"POINT") )
256 : {
257 33376 : poGeom = new OGRPoint();
258 : }
259 :
260 30236 : else if( EQUAL(szToken,"LINESTRING") )
261 : {
262 652 : poGeom = new OGRLineString();
263 : }
264 :
265 29584 : else if( EQUAL(szToken,"POLYGON") )
266 : {
267 25288 : poGeom = new OGRPolygon();
268 : }
269 :
270 4296 : else if( EQUAL(szToken,"GEOMETRYCOLLECTION") )
271 : {
272 238 : poGeom = new OGRGeometryCollection();
273 : }
274 :
275 4058 : else if( EQUAL(szToken,"MULTIPOLYGON") )
276 : {
277 356 : poGeom = new OGRMultiPolygon();
278 : }
279 :
280 3702 : else if( EQUAL(szToken,"MULTIPOINT") )
281 : {
282 226 : poGeom = new OGRMultiPoint();
283 : }
284 :
285 3476 : else if( EQUAL(szToken,"MULTILINESTRING") )
286 : {
287 254 : poGeom = new OGRMultiLineString();
288 : }
289 :
290 : else
291 : {
292 3222 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
293 : }
294 :
295 : /* -------------------------------------------------------------------- */
296 : /* Do the import. */
297 : /* -------------------------------------------------------------------- */
298 60390 : eErr = poGeom->importFromWkt( &pszInput );
299 :
300 : /* -------------------------------------------------------------------- */
301 : /* Assign spatial reference system. */
302 : /* -------------------------------------------------------------------- */
303 60390 : if( eErr == OGRERR_NONE )
304 : {
305 60040 : poGeom->assignSpatialReference( poSR );
306 60040 : *ppoReturn = poGeom;
307 60040 : *ppszData = pszInput;
308 : }
309 : else
310 : {
311 350 : delete poGeom;
312 : }
313 :
314 60390 : 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 59784 : 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 59784 : (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 7061 : OGRGeometryFactory::createGeometry( OGRwkbGeometryType eGeometryType )
370 :
371 : {
372 7061 : switch( wkbFlatten(eGeometryType) )
373 : {
374 : case wkbPoint:
375 734 : return new OGRPoint();
376 :
377 : case wkbLineString:
378 2174 : return new OGRLineString();
379 :
380 : case wkbPolygon:
381 2672 : return new OGRPolygon();
382 :
383 : case wkbGeometryCollection:
384 478 : return new OGRGeometryCollection();
385 :
386 : case wkbMultiPolygon:
387 178 : return new OGRMultiPolygon();
388 :
389 : case wkbMultiPoint:
390 105 : return new OGRMultiPoint();
391 :
392 : case wkbMultiLineString:
393 122 : return new OGRMultiLineString();
394 :
395 : case wkbLinearRing:
396 542 : return new OGRLinearRing();
397 :
398 : default:
399 56 : 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 2664 : OGRGeometryH OGR_G_CreateGeometry( OGRwkbGeometryType eGeometryType )
422 :
423 : {
424 2664 : 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 89477 : void OGRGeometryFactory::destroyGeometry( OGRGeometry *poGeom )
444 :
445 : {
446 89477 : delete poGeom;
447 89477 : }
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 88901 : void OGR_G_DestroyGeometry( OGRGeometryH hGeom )
466 :
467 : {
468 88901 : OGRGeometryFactory::destroyGeometry( (OGRGeometry *) hGeom );
469 88901 : }
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 : * @param poGeom the input geometry - ownership is passed to the method.
483 : * @return new geometry.
484 : */
485 :
486 90 : OGRGeometry *OGRGeometryFactory::forceToPolygon( OGRGeometry *poGeom )
487 :
488 : {
489 90 : if( poGeom == NULL )
490 0 : return NULL;
491 :
492 90 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
493 :
494 90 : if( eGeomType != wkbGeometryCollection
495 : && eGeomType != wkbMultiPolygon )
496 62 : return poGeom;
497 :
498 : // build an aggregated polygon from all the polygon rings in the container.
499 28 : OGRPolygon *poPolygon = new OGRPolygon();
500 28 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
501 : int iGeom;
502 :
503 60 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
504 : {
505 32 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
506 : != wkbPolygon )
507 18 : continue;
508 :
509 14 : OGRPolygon *poOldPoly = (OGRPolygon *) poGC->getGeometryRef(iGeom);
510 : int iRing;
511 :
512 14 : if( poOldPoly->getExteriorRing() == NULL )
513 4 : continue;
514 :
515 10 : poPolygon->addRing( poOldPoly->getExteriorRing() );
516 :
517 18 : for( iRing = 0; iRing < poOldPoly->getNumInteriorRings(); iRing++ )
518 8 : poPolygon->addRing( poOldPoly->getInteriorRing( iRing ) );
519 : }
520 :
521 28 : delete poGC;
522 :
523 28 : return poPolygon;
524 : }
525 :
526 : /************************************************************************/
527 : /* OGR_G_ForceToPolygon() */
528 : /************************************************************************/
529 :
530 : /**
531 : * \brief Convert to polygon.
532 : *
533 : * This function is the same as the C++ method
534 : * OGRGeometryFactory::forceToPolygon().
535 : *
536 : * @param hGeom handle to the geometry to convert (ownership surrendered).
537 : * @return the converted geometry (ownership to caller).
538 : *
539 : * @since GDAL/OGR 1.8.0
540 : */
541 :
542 70 : OGRGeometryH OGR_G_ForceToPolygon( OGRGeometryH hGeom )
543 :
544 : {
545 : return (OGRGeometryH)
546 70 : OGRGeometryFactory::forceToPolygon( (OGRGeometry *) hGeom );
547 : }
548 :
549 : /************************************************************************/
550 : /* forceToMultiPolygon() */
551 : /************************************************************************/
552 :
553 : /**
554 : * \brief Convert to multipolygon.
555 : *
556 : * Tries to force the provided geometry to be a multipolygon. Currently
557 : * this just effects a change on polygons. The passed in geometry is
558 : * consumed and a new one returned (or potentially the same one).
559 : *
560 : * @return new geometry.
561 : */
562 :
563 690 : OGRGeometry *OGRGeometryFactory::forceToMultiPolygon( OGRGeometry *poGeom )
564 :
565 : {
566 690 : if( poGeom == NULL )
567 0 : return NULL;
568 :
569 690 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
570 :
571 : /* -------------------------------------------------------------------- */
572 : /* If this is already a MultiPolygon, nothing to do */
573 : /* -------------------------------------------------------------------- */
574 690 : if( eGeomType == wkbMultiPolygon )
575 : {
576 4 : return poGeom;
577 : }
578 :
579 : /* -------------------------------------------------------------------- */
580 : /* Check for the case of a geometrycollection that can be */
581 : /* promoted to MultiPolygon. */
582 : /* -------------------------------------------------------------------- */
583 686 : if( eGeomType == wkbGeometryCollection )
584 : {
585 : int iGeom;
586 24 : int bAllPoly = TRUE;
587 24 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
588 :
589 50 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
590 : {
591 26 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
592 : != wkbPolygon )
593 18 : bAllPoly = FALSE;
594 : }
595 :
596 24 : if( !bAllPoly )
597 16 : return poGeom;
598 :
599 8 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
600 :
601 24 : while( poGC->getNumGeometries() > 0 )
602 : {
603 8 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
604 8 : poGC->removeGeometry( 0, FALSE );
605 : }
606 :
607 8 : delete poGC;
608 :
609 8 : return poMP;
610 : }
611 :
612 : /* -------------------------------------------------------------------- */
613 : /* Eventually we should try to split the polygon into component */
614 : /* island polygons. But thats alot of work and can be put off. */
615 : /* -------------------------------------------------------------------- */
616 662 : if( eGeomType != wkbPolygon )
617 16 : return poGeom;
618 :
619 646 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
620 646 : poMP->addGeometryDirectly( poGeom );
621 :
622 646 : return poMP;
623 : }
624 :
625 : /************************************************************************/
626 : /* OGR_G_ForceToMultiPolygon() */
627 : /************************************************************************/
628 :
629 : /**
630 : * \brief Convert to multipolygon.
631 : *
632 : * This function is the same as the C++ method
633 : * OGRGeometryFactory::forceToMultiPolygon().
634 : *
635 : * @param hGeom handle to the geometry to convert (ownership surrendered).
636 : * @return the converted geometry (ownership to caller).
637 : *
638 : * @since GDAL/OGR 1.8.0
639 : */
640 :
641 60 : OGRGeometryH OGR_G_ForceToMultiPolygon( OGRGeometryH hGeom )
642 :
643 : {
644 : return (OGRGeometryH)
645 60 : OGRGeometryFactory::forceToMultiPolygon( (OGRGeometry *) hGeom );
646 : }
647 :
648 : /************************************************************************/
649 : /* forceToMultiPoint() */
650 : /************************************************************************/
651 :
652 : /**
653 : * \brief Convert to multipoint.
654 : *
655 : * Tries to force the provided geometry to be a multipoint. Currently
656 : * this just effects a change on points. The passed in geometry is
657 : * consumed and a new one returned (or potentially the same one).
658 : *
659 : * @return new geometry.
660 : */
661 :
662 52 : OGRGeometry *OGRGeometryFactory::forceToMultiPoint( OGRGeometry *poGeom )
663 :
664 : {
665 52 : if( poGeom == NULL )
666 0 : return NULL;
667 :
668 52 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* If this is already a MultiPoint, nothing to do */
672 : /* -------------------------------------------------------------------- */
673 52 : if( eGeomType == wkbMultiPoint )
674 : {
675 4 : return poGeom;
676 : }
677 :
678 : /* -------------------------------------------------------------------- */
679 : /* Check for the case of a geometrycollection that can be */
680 : /* promoted to MultiPoint. */
681 : /* -------------------------------------------------------------------- */
682 48 : if( eGeomType == wkbGeometryCollection )
683 : {
684 : int iGeom;
685 24 : int bAllPoint = TRUE;
686 24 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
687 :
688 52 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
689 : {
690 28 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
691 : != wkbPoint )
692 20 : bAllPoint = FALSE;
693 : }
694 :
695 24 : if( !bAllPoint )
696 16 : return poGeom;
697 :
698 8 : OGRMultiPoint *poMP = new OGRMultiPoint();
699 :
700 24 : while( poGC->getNumGeometries() > 0 )
701 : {
702 8 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
703 8 : poGC->removeGeometry( 0, FALSE );
704 : }
705 :
706 8 : delete poGC;
707 :
708 8 : return poMP;
709 : }
710 :
711 24 : if( eGeomType != wkbPoint )
712 18 : return poGeom;
713 :
714 6 : OGRMultiPoint *poMP = new OGRMultiPoint();
715 6 : poMP->addGeometryDirectly( poGeom );
716 :
717 6 : return poMP;
718 : }
719 :
720 : /************************************************************************/
721 : /* OGR_G_ForceToMultiPoint() */
722 : /************************************************************************/
723 :
724 : /**
725 : * \brief Convert to multipoint.
726 : *
727 : * This function is the same as the C++ method
728 : * OGRGeometryFactory::forceToMultiPoint().
729 : *
730 : * @param hGeom handle to the geometry to convert (ownership surrendered).
731 : * @return the converted geometry (ownership to caller).
732 : *
733 : * @since GDAL/OGR 1.8.0
734 : */
735 :
736 52 : OGRGeometryH OGR_G_ForceToMultiPoint( OGRGeometryH hGeom )
737 :
738 : {
739 : return (OGRGeometryH)
740 52 : OGRGeometryFactory::forceToMultiPoint( (OGRGeometry *) hGeom );
741 : }
742 :
743 : /************************************************************************/
744 : /* forceToMultiLinestring() */
745 : /************************************************************************/
746 :
747 : /**
748 : * \brief Convert to multilinestring.
749 : *
750 : * Tries to force the provided geometry to be a multilinestring.
751 : *
752 : * - linestrings are placed in a multilinestring.
753 : * - geometry collections will be converted to multilinestring if they only
754 : * contain linestrings.
755 : * - polygons will be changed to a collection of linestrings (one per ring).
756 : *
757 : * The passed in geometry is
758 : * consumed and a new one returned (or potentially the same one).
759 : *
760 : * @return new geometry.
761 : */
762 :
763 522 : OGRGeometry *OGRGeometryFactory::forceToMultiLineString( OGRGeometry *poGeom )
764 :
765 : {
766 522 : if( poGeom == NULL )
767 0 : return NULL;
768 :
769 522 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
770 :
771 : /* -------------------------------------------------------------------- */
772 : /* If this is already a MultiLineString, nothing to do */
773 : /* -------------------------------------------------------------------- */
774 522 : if( eGeomType == wkbMultiLineString )
775 : {
776 20 : return poGeom;
777 : }
778 :
779 : /* -------------------------------------------------------------------- */
780 : /* Check for the case of a geometrycollection that can be */
781 : /* promoted to MultiLineString. */
782 : /* -------------------------------------------------------------------- */
783 502 : if( eGeomType == wkbGeometryCollection )
784 : {
785 : int iGeom;
786 24 : int bAllLines = TRUE;
787 24 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
788 :
789 52 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
790 : {
791 28 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
792 : != wkbLineString )
793 18 : bAllLines = FALSE;
794 : }
795 :
796 24 : if( !bAllLines )
797 16 : return poGeom;
798 :
799 8 : OGRMultiLineString *poMP = new OGRMultiLineString();
800 :
801 26 : while( poGC->getNumGeometries() > 0 )
802 : {
803 10 : poMP->addGeometryDirectly( poGC->getGeometryRef(0) );
804 10 : poGC->removeGeometry( 0, FALSE );
805 : }
806 :
807 8 : delete poGC;
808 :
809 8 : return poMP;
810 : }
811 :
812 : /* -------------------------------------------------------------------- */
813 : /* Turn a linestring into a multilinestring. */
814 : /* -------------------------------------------------------------------- */
815 478 : if( eGeomType == wkbLineString )
816 : {
817 430 : OGRMultiLineString *poMP = new OGRMultiLineString();
818 430 : poMP->addGeometryDirectly( poGeom );
819 430 : return poMP;
820 : }
821 :
822 : /* -------------------------------------------------------------------- */
823 : /* Convert polygons into a multilinestring. */
824 : /* -------------------------------------------------------------------- */
825 48 : if( eGeomType == wkbPolygon )
826 : {
827 10 : OGRMultiLineString *poMP = new OGRMultiLineString();
828 10 : OGRPolygon *poPoly = (OGRPolygon *) poGeom;
829 : int iRing;
830 :
831 26 : for( iRing = 0; iRing < poPoly->getNumInteriorRings()+1; iRing++ )
832 : {
833 : OGRLineString *poNewLS, *poLR;
834 :
835 18 : if( iRing == 0 )
836 : {
837 10 : poLR = poPoly->getExteriorRing();
838 10 : if( poLR == NULL )
839 2 : break;
840 : }
841 : else
842 8 : poLR = poPoly->getInteriorRing(iRing-1);
843 :
844 16 : if (poLR == NULL || poLR->getNumPoints() == 0)
845 4 : continue;
846 :
847 12 : poNewLS = new OGRLineString();
848 12 : poNewLS->addSubLineString( poLR );
849 12 : poMP->addGeometryDirectly( poNewLS );
850 : }
851 :
852 10 : delete poPoly;
853 :
854 10 : return poMP;
855 : }
856 :
857 : /* -------------------------------------------------------------------- */
858 : /* Convert multi-polygons into a multilinestring. */
859 : /* -------------------------------------------------------------------- */
860 38 : if( eGeomType == wkbMultiPolygon )
861 : {
862 6 : OGRMultiLineString *poMP = new OGRMultiLineString();
863 6 : OGRMultiPolygon *poMPoly = (OGRMultiPolygon *) poGeom;
864 : int iPoly;
865 :
866 16 : for( iPoly = 0; iPoly < poMPoly->getNumGeometries(); iPoly++ )
867 : {
868 10 : OGRPolygon *poPoly = (OGRPolygon*) poMPoly->getGeometryRef(iPoly);
869 : int iRing;
870 :
871 24 : for( iRing = 0; iRing < poPoly->getNumInteriorRings()+1; iRing++ )
872 : {
873 : OGRLineString *poNewLS, *poLR;
874 :
875 16 : if( iRing == 0 )
876 : {
877 10 : poLR = poPoly->getExteriorRing();
878 10 : if( poLR == NULL )
879 2 : break;
880 : }
881 : else
882 6 : poLR = poPoly->getInteriorRing(iRing-1);
883 :
884 14 : if (poLR == NULL || poLR->getNumPoints() == 0)
885 2 : continue;
886 :
887 12 : poNewLS = new OGRLineString();
888 12 : poNewLS->addSubLineString( poLR );
889 12 : poMP->addGeometryDirectly( poNewLS );
890 : }
891 : }
892 6 : delete poMPoly;
893 :
894 6 : return poMP;
895 : }
896 :
897 32 : return poGeom;
898 : }
899 :
900 : /************************************************************************/
901 : /* OGR_G_ForceToMultiLineString() */
902 : /************************************************************************/
903 :
904 : /**
905 : * \brief Convert to multilinestring.
906 : *
907 : * This function is the same as the C++ method
908 : * OGRGeometryFactory::forceToMultiLineString().
909 : *
910 : * @param hGeom handle to the geometry to convert (ownership surrendered).
911 : * @return the converted geometry (ownership to caller).
912 : *
913 : * @since GDAL/OGR 1.8.0
914 : */
915 :
916 56 : OGRGeometryH OGR_G_ForceToMultiLineString( OGRGeometryH hGeom )
917 :
918 : {
919 : return (OGRGeometryH)
920 56 : OGRGeometryFactory::forceToMultiLineString( (OGRGeometry *) hGeom );
921 : }
922 :
923 : /************************************************************************/
924 : /* organizePolygons() */
925 : /************************************************************************/
926 :
927 : typedef struct _sPolyExtended sPolyExtended;
928 :
929 : struct _sPolyExtended
930 1260 : {
931 : OGRPolygon* poPolygon;
932 : OGREnvelope sEnvelope;
933 : OGRLinearRing* poExteriorRing;
934 : OGRPoint poAPoint;
935 : int nInitialIndex;
936 : int bIsTopLevel;
937 : OGRPolygon* poEnclosingPolygon;
938 : double dfArea;
939 : int bIsCW;
940 : };
941 :
942 692 : static int OGRGeometryFactoryCompareArea(const void* p1, const void* p2)
943 : {
944 692 : const sPolyExtended* psPoly1 = (const sPolyExtended*) p1;
945 692 : const sPolyExtended* psPoly2 = (const sPolyExtended*) p2;
946 692 : if (psPoly2->dfArea < psPoly1->dfArea)
947 384 : return -1;
948 308 : else if (psPoly2->dfArea > psPoly1->dfArea)
949 272 : return 1;
950 : else
951 36 : return 0;
952 : }
953 :
954 714 : static int OGRGeometryFactoryCompareByIndex(const void* p1, const void* p2)
955 : {
956 714 : const sPolyExtended* psPoly1 = (const sPolyExtended*) p1;
957 714 : const sPolyExtended* psPoly2 = (const sPolyExtended*) p2;
958 714 : if (psPoly1->nInitialIndex < psPoly2->nInitialIndex)
959 430 : return -1;
960 284 : else if (psPoly1->nInitialIndex > psPoly2->nInitialIndex)
961 284 : return 1;
962 : else
963 0 : return 0;
964 : }
965 :
966 : #define N_CRITICAL_PART_NUMBER 100
967 :
968 : typedef enum
969 : {
970 : METHOD_NORMAL,
971 : METHOD_SKIP,
972 : METHOD_ONLY_CCW
973 : } OrganizePolygonMethod;
974 :
975 : /**
976 : * \brief Organize polygons based on geometries.
977 : *
978 : * Analyse a set of rings (passed as simple polygons), and based on a
979 : * geometric analysis convert them into a polygon with inner rings,
980 : * or a MultiPolygon if dealing with more than one polygon.
981 : *
982 : * All the input geometries must be OGRPolygons with only a valid exterior
983 : * ring (at least 4 points) and no interior rings.
984 : *
985 : * The passed in geometries become the responsibility of the method, but the
986 : * papoPolygons "pointer array" remains owned by the caller.
987 : *
988 : * For faster computation, a polygon is considered to be inside
989 : * another one if a single point of its external ring is included into the other one.
990 : * (unless 'OGR_DEBUG_ORGANIZE_POLYGONS' configuration option is set to TRUE.
991 : * In that case, a slower algorithm that tests exact topological relationships
992 : * is used if GEOS is available.)
993 : *
994 : * In cases where a big number of polygons is passed to this function, the default processing
995 : * may be really slow. You can skip the processing by adding METHOD=SKIP
996 : * to the option list (the result of the function will be a multi-polygon with all polygons
997 : * as toplevel polygons) or only make it analyze counterclockwise polygons by adding
998 : * METHOD=ONLY_CCW to the option list if you can assume that the outline
999 : * of holes is counterclockwise defined (this is the convention for shapefiles e.g.)
1000 : *
1001 : * If the OGR_ORGANIZE_POLYGONS configuration option is defined, its value will override
1002 : * the value of the METHOD option of papszOptions (usefull to modify the behaviour of the
1003 : * shapefile driver)
1004 : *
1005 : * @param papoPolygons array of geometry pointers - should all be OGRPolygons.
1006 : * Ownership of the geometries is passed, but not of the array itself.
1007 : * @param nPolygonCount number of items in papoPolygons
1008 : * @param pbIsValidGeometry value will be set TRUE if result is valid or
1009 : * FALSE otherwise.
1010 : * @param papszOptions a list of strings for passing options
1011 : *
1012 : * @return a single resulting geometry (either OGRPolygon or OGRMultiPolygon).
1013 : */
1014 :
1015 980 : OGRGeometry* OGRGeometryFactory::organizePolygons( OGRGeometry **papoPolygons,
1016 : int nPolygonCount,
1017 : int *pbIsValidGeometry,
1018 : const char** papszOptions )
1019 : {
1020 : int bUseFastVersion;
1021 : int i, j;
1022 980 : OGRGeometry* geom = NULL;
1023 980 : OrganizePolygonMethod method = METHOD_NORMAL;
1024 :
1025 : /* -------------------------------------------------------------------- */
1026 : /* Trivial case of a single polygon. */
1027 : /* -------------------------------------------------------------------- */
1028 980 : if (nPolygonCount == 1)
1029 : {
1030 808 : geom = papoPolygons[0];
1031 808 : papoPolygons[0] = NULL;
1032 :
1033 808 : if( pbIsValidGeometry )
1034 808 : *pbIsValidGeometry = TRUE;
1035 :
1036 808 : return geom;
1037 : }
1038 :
1039 172 : if (CSLTestBoolean(CPLGetConfigOption("OGR_DEBUG_ORGANIZE_POLYGONS", "NO")))
1040 : {
1041 : /* -------------------------------------------------------------------- */
1042 : /* A wee bit of a warning. */
1043 : /* -------------------------------------------------------------------- */
1044 : static int firstTime = 1;
1045 0 : if (!haveGEOS() && firstTime)
1046 : {
1047 : CPLDebug("OGR",
1048 : "In OGR_DEBUG_ORGANIZE_POLYGONS mode, GDAL should be built with GEOS support enabled in order "
1049 0 : "OGRGeometryFactory::organizePolygons to provide reliable results on complex polygons.");
1050 0 : firstTime = 0;
1051 : }
1052 0 : bUseFastVersion = !haveGEOS();
1053 : }
1054 : else
1055 : {
1056 172 : bUseFastVersion = TRUE;
1057 : }
1058 :
1059 : /* -------------------------------------------------------------------- */
1060 : /* Setup per polygon envelope and area information. */
1061 : /* -------------------------------------------------------------------- */
1062 172 : sPolyExtended* asPolyEx = new sPolyExtended[nPolygonCount];
1063 :
1064 172 : int go_on = TRUE;
1065 172 : int bMixedUpGeometries = FALSE;
1066 172 : int bNonPolygon = FALSE;
1067 172 : int bFoundCCW = FALSE;
1068 :
1069 172 : const char* pszMethodValue = CSLFetchNameValue( (char**)papszOptions, "METHOD" );
1070 172 : const char* pszMethodValueOption = CPLGetConfigOption("OGR_ORGANIZE_POLYGONS", NULL);
1071 216 : if (pszMethodValueOption != NULL && pszMethodValueOption[0] != '\0')
1072 2 : pszMethodValue = pszMethodValueOption;
1073 :
1074 172 : if (pszMethodValue != NULL)
1075 : {
1076 64 : if (EQUAL(pszMethodValue, "SKIP"))
1077 : {
1078 0 : method = METHOD_SKIP;
1079 0 : bMixedUpGeometries = TRUE;
1080 : }
1081 64 : else if (EQUAL(pszMethodValue, "ONLY_CCW"))
1082 : {
1083 62 : method = METHOD_ONLY_CCW;
1084 : }
1085 2 : else if (!EQUAL(pszMethodValue, "DEFAULT"))
1086 : {
1087 : CPLError(CE_Warning, CPLE_AppDefined,
1088 0 : "Unrecognized value for METHOD option : %s", pszMethodValue);
1089 : }
1090 : }
1091 :
1092 172 : int nCountCWPolygon = 0;
1093 172 : int indexOfCWPolygon = -1;
1094 :
1095 802 : for(i=0;i<nPolygonCount;i++)
1096 : {
1097 630 : asPolyEx[i].nInitialIndex = i;
1098 630 : asPolyEx[i].poPolygon = (OGRPolygon*)papoPolygons[i];
1099 630 : papoPolygons[i]->getEnvelope(&asPolyEx[i].sEnvelope);
1100 :
1101 1890 : if( wkbFlatten(papoPolygons[i]->getGeometryType()) == wkbPolygon
1102 630 : && ((OGRPolygon *) papoPolygons[i])->getNumInteriorRings() == 0
1103 630 : && ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints() >= 4)
1104 : {
1105 630 : asPolyEx[i].dfArea = asPolyEx[i].poPolygon->get_Area();
1106 630 : asPolyEx[i].poExteriorRing = asPolyEx[i].poPolygon->getExteriorRing();
1107 630 : asPolyEx[i].poExteriorRing->getPoint(0, &asPolyEx[i].poAPoint);
1108 630 : asPolyEx[i].bIsCW = asPolyEx[i].poExteriorRing->isClockwise();
1109 630 : if (asPolyEx[i].bIsCW)
1110 : {
1111 476 : indexOfCWPolygon = i;
1112 476 : nCountCWPolygon ++;
1113 : }
1114 630 : if (!bFoundCCW)
1115 386 : bFoundCCW = ! (asPolyEx[i].bIsCW);
1116 : }
1117 : else
1118 : {
1119 0 : if( !bMixedUpGeometries )
1120 : {
1121 : CPLError(
1122 : CE_Warning, CPLE_AppDefined,
1123 : "organizePolygons() received an unexpected geometry.\n"
1124 : "Either a polygon with interior rings, or a polygon with less than 4 points,\n"
1125 0 : "or a non-Polygon geometry. Return arguments as a collection." );
1126 0 : bMixedUpGeometries = TRUE;
1127 : }
1128 0 : if( wkbFlatten(papoPolygons[i]->getGeometryType()) != wkbPolygon )
1129 0 : bNonPolygon = TRUE;
1130 : }
1131 : }
1132 :
1133 : /* If we are in ONLY_CCW mode and that we have found that there is only one outer ring, */
1134 : /* then it is pretty easy : we can assume that all other rings are inside */
1135 172 : if (method == METHOD_ONLY_CCW && nCountCWPolygon == 1 && bUseFastVersion)
1136 : {
1137 20 : geom = asPolyEx[indexOfCWPolygon].poPolygon;
1138 62 : for(i=0; i<nPolygonCount; i++)
1139 : {
1140 42 : if (i != indexOfCWPolygon)
1141 : {
1142 22 : ((OGRPolygon*)geom)->addRing(asPolyEx[i].poPolygon->getExteriorRing());
1143 22 : delete asPolyEx[i].poPolygon;
1144 : }
1145 : }
1146 20 : delete [] asPolyEx;
1147 20 : if (pbIsValidGeometry)
1148 20 : *pbIsValidGeometry = TRUE;
1149 20 : return geom;
1150 : }
1151 :
1152 : /* Emits a warning if the number of parts is sufficiently big to anticipate for */
1153 : /* very long computation time, and the user didn't specify an explicit method */
1154 152 : if (nPolygonCount > N_CRITICAL_PART_NUMBER && method == METHOD_NORMAL && pszMethodValue == NULL)
1155 : {
1156 : static int firstTime = 1;
1157 0 : if (firstTime)
1158 : {
1159 0 : if (bFoundCCW)
1160 : {
1161 : CPLError( CE_Warning, CPLE_AppDefined,
1162 : "organizePolygons() received a polygon with more than %d parts. "
1163 : "The processing may be really slow.\n"
1164 : "You can skip the processing by setting METHOD=SKIP, "
1165 : "or only make it analyze counter-clock wise parts by setting "
1166 : "METHOD=ONLY_CCW if you can assume that the "
1167 0 : "outline of holes is counter-clock wise defined", N_CRITICAL_PART_NUMBER);
1168 : }
1169 : else
1170 : {
1171 : CPLError( CE_Warning, CPLE_AppDefined,
1172 : "organizePolygons() received a polygon with more than %d parts. "
1173 : "The processing may be really slow.\n"
1174 : "You can skip the processing by setting METHOD=SKIP.",
1175 0 : N_CRITICAL_PART_NUMBER);
1176 : }
1177 0 : firstTime = 0;
1178 : }
1179 : }
1180 :
1181 :
1182 : /* This a several steps algorithm :
1183 : 1) Sort polygons by descending areas
1184 : 2) For each polygon of rank i, find its smallest enclosing polygon
1185 : among the polygons of rank [i-1 ... 0]. If there are no such polygon,
1186 : this is a toplevel polygon. Otherwise, depending on if the enclosing
1187 : polygon is toplevel or not, we can decide if we are toplevel or not
1188 : 3) Re-sort the polygons to retrieve their inital order (nicer for some applications)
1189 : 4) For each non toplevel polygon (= inner ring), add it to its outer ring
1190 : 5) Add the toplevel polygons to the multipolygon
1191 :
1192 : Complexity : O(nPolygonCount^2)
1193 : */
1194 :
1195 : /* Compute how each polygon relate to the other ones
1196 : To save a bit of computation we always begin the computation by a test
1197 : on the enveloppe. We also take into account the areas to avoid some
1198 : useless tests. (A contains B implies envelop(A) contains envelop(B)
1199 : and area(A) > area(B)) In practise, we can hope that few full geometry
1200 : intersection of inclusion test is done:
1201 : * if the polygons are well separated geographically (a set of islands
1202 : for example), no full geometry intersection or inclusion test is done.
1203 : (the envelopes don't intersect each other)
1204 :
1205 : * if the polygons are 'lake inside an island inside a lake inside an
1206 : area' and that each polygon is much smaller than its enclosing one,
1207 : their bounding boxes are stricly contained into each oter, and thus,
1208 : no full geometry intersection or inclusion test is done
1209 : */
1210 :
1211 152 : if (!bMixedUpGeometries)
1212 : {
1213 : /* STEP 1 : Sort polygons by descending area */
1214 152 : qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareArea);
1215 : }
1216 152 : papoPolygons = NULL; /* just to use to avoid it afterwards */
1217 :
1218 : /* -------------------------------------------------------------------- */
1219 : /* Compute relationships, if things seem well structured. */
1220 : /* -------------------------------------------------------------------- */
1221 :
1222 : /* The first (largest) polygon is necessarily top-level */
1223 152 : asPolyEx[0].bIsTopLevel = TRUE;
1224 152 : asPolyEx[0].poEnclosingPolygon = NULL;
1225 :
1226 152 : int nCountTopLevel = 1;
1227 :
1228 : /* STEP 2 */
1229 588 : for(i=1; !bMixedUpGeometries && go_on && i<nPolygonCount; i++)
1230 : {
1231 436 : if (method == METHOD_ONLY_CCW && asPolyEx[i].bIsCW)
1232 : {
1233 56 : nCountTopLevel ++;
1234 56 : asPolyEx[i].bIsTopLevel = TRUE;
1235 56 : asPolyEx[i].poEnclosingPolygon = NULL;
1236 56 : continue;
1237 : }
1238 :
1239 820 : for(j=i-1; go_on && j>=0;j--)
1240 : {
1241 748 : int b_i_inside_j = FALSE;
1242 :
1243 748 : if (method == METHOD_ONLY_CCW && asPolyEx[j].bIsCW == FALSE)
1244 : {
1245 : /* In that mode, i which is CCW if we reach here can only be */
1246 : /* included in a CW polygon */
1247 10 : continue;
1248 : }
1249 :
1250 738 : if (asPolyEx[j].sEnvelope.Contains(asPolyEx[i].sEnvelope))
1251 : {
1252 310 : if (bUseFastVersion)
1253 : {
1254 : /* Note that isPointInRing only test strict inclusion in the ring */
1255 310 : if (asPolyEx[j].poExteriorRing->isPointInRing(&asPolyEx[i].poAPoint, FALSE))
1256 : {
1257 306 : b_i_inside_j = TRUE;
1258 : }
1259 4 : else if (asPolyEx[j].poExteriorRing->isPointOnRingBoundary(&asPolyEx[i].poAPoint, FALSE))
1260 : {
1261 : /* If the point of i is on the boundary of j, we will iterate over the other points of i */
1262 4 : int k, nPoints = asPolyEx[i].poExteriorRing->getNumPoints();
1263 4 : for(k=1;k<nPoints;k++)
1264 : {
1265 12 : OGRPoint point;
1266 12 : asPolyEx[i].poExteriorRing->getPoint(k, &point);
1267 12 : if (asPolyEx[j].poExteriorRing->isPointInRing(&point, FALSE))
1268 : {
1269 : /* If then point is strictly included in j, then i is considered inside j */
1270 2 : b_i_inside_j = TRUE;
1271 : break;
1272 : }
1273 10 : else if (asPolyEx[j].poExteriorRing->isPointOnRingBoundary(&point, FALSE))
1274 : {
1275 : /* If it is on the boundary of j, iterate again */
1276 : }
1277 : else
1278 : {
1279 : /* If it is outside, then i cannot be inside j */
1280 : break;
1281 : }
1282 : }
1283 : }
1284 : }
1285 0 : else if (asPolyEx[j].poPolygon->Contains(asPolyEx[i].poPolygon))
1286 : {
1287 0 : b_i_inside_j = TRUE;
1288 : }
1289 : }
1290 :
1291 :
1292 738 : if (b_i_inside_j)
1293 : {
1294 308 : if (asPolyEx[j].bIsTopLevel)
1295 : {
1296 : /* We are a lake */
1297 240 : asPolyEx[i].bIsTopLevel = FALSE;
1298 240 : asPolyEx[i].poEnclosingPolygon = asPolyEx[j].poPolygon;
1299 : }
1300 : else
1301 : {
1302 : /* We are included in a something not toplevel (a lake), */
1303 : /* so in OGCSF we are considered as toplevel too */
1304 68 : nCountTopLevel ++;
1305 68 : asPolyEx[i].bIsTopLevel = TRUE;
1306 68 : asPolyEx[i].poEnclosingPolygon = NULL;
1307 : }
1308 308 : break;
1309 : }
1310 : /* We use Overlaps instead of Intersects to be more
1311 : tolerant about touching polygons */
1312 430 : else if ( bUseFastVersion || !asPolyEx[i].sEnvelope.Intersects(asPolyEx[j].sEnvelope)
1313 0 : || !asPolyEx[i].poPolygon->Overlaps(asPolyEx[j].poPolygon) )
1314 : {
1315 :
1316 : }
1317 : else
1318 : {
1319 : /* Bad... The polygons are intersecting but no one is
1320 : contained inside the other one. This is a really broken
1321 : case. We just make a multipolygon with the whole set of
1322 : polygons */
1323 0 : go_on = FALSE;
1324 : #ifdef DEBUG
1325 : char* wkt1;
1326 : char* wkt2;
1327 0 : asPolyEx[i].poPolygon->exportToWkt(&wkt1);
1328 0 : asPolyEx[j].poPolygon->exportToWkt(&wkt2);
1329 : CPLDebug( "OGR",
1330 : "Bad intersection for polygons %d and %d\n"
1331 : "geom %d: %s\n"
1332 : "geom %d: %s",
1333 0 : i, j, i, wkt1, j, wkt2 );
1334 0 : CPLFree(wkt1);
1335 0 : CPLFree(wkt2);
1336 : #endif
1337 : }
1338 : }
1339 :
1340 380 : if (j < 0)
1341 : {
1342 : /* We come here because we are not included in anything */
1343 : /* We are toplevel */
1344 72 : nCountTopLevel ++;
1345 72 : asPolyEx[i].bIsTopLevel = TRUE;
1346 72 : asPolyEx[i].poEnclosingPolygon = NULL;
1347 : }
1348 : }
1349 :
1350 152 : if (pbIsValidGeometry)
1351 152 : *pbIsValidGeometry = go_on && !bMixedUpGeometries;
1352 :
1353 : /* -------------------------------------------------------------------- */
1354 : /* Things broke down - just turn everything into a multipolygon. */
1355 : /* -------------------------------------------------------------------- */
1356 152 : if ( !go_on || bMixedUpGeometries )
1357 : {
1358 0 : if( bNonPolygon )
1359 0 : geom = new OGRGeometryCollection();
1360 : else
1361 0 : geom = new OGRMultiPolygon();
1362 :
1363 0 : for( i=0; i < nPolygonCount; i++ )
1364 : {
1365 : ((OGRGeometryCollection*)geom)->
1366 0 : addGeometryDirectly( asPolyEx[i].poPolygon );
1367 : }
1368 : }
1369 :
1370 : /* -------------------------------------------------------------------- */
1371 : /* Try to turn into one or more polygons based on the ring */
1372 : /* relationships. */
1373 : /* -------------------------------------------------------------------- */
1374 : else
1375 : {
1376 : /* STEP 3: Resort in initial order */
1377 152 : qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareByIndex);
1378 :
1379 : /* STEP 4: Add holes as rings of their enclosing polygon */
1380 740 : for(i=0;i<nPolygonCount;i++)
1381 : {
1382 588 : if (asPolyEx[i].bIsTopLevel == FALSE)
1383 : {
1384 240 : asPolyEx[i].poEnclosingPolygon->addRing(
1385 480 : asPolyEx[i].poPolygon->getExteriorRing());
1386 240 : delete asPolyEx[i].poPolygon;
1387 : }
1388 348 : else if (nCountTopLevel == 1)
1389 : {
1390 40 : geom = asPolyEx[i].poPolygon;
1391 : }
1392 : }
1393 :
1394 : /* STEP 5: Add toplevel polygons */
1395 152 : if (nCountTopLevel > 1)
1396 : {
1397 616 : for(i=0;i<nPolygonCount;i++)
1398 : {
1399 504 : if (asPolyEx[i].bIsTopLevel)
1400 : {
1401 308 : if (geom == NULL)
1402 : {
1403 112 : geom = new OGRMultiPolygon();
1404 112 : ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon);
1405 : }
1406 : else
1407 : {
1408 196 : ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon);
1409 : }
1410 : }
1411 : }
1412 : }
1413 : }
1414 :
1415 152 : delete[] asPolyEx;
1416 :
1417 152 : return geom;
1418 : }
1419 :
1420 : /************************************************************************/
1421 : /* createFromGML() */
1422 : /************************************************************************/
1423 :
1424 : /**
1425 : * \brief Create geometry from GML.
1426 : *
1427 : * This method translates a fragment of GML containing only the geometry
1428 : * portion into a corresponding OGRGeometry. There are many limitations
1429 : * on the forms of GML geometries supported by this parser, but they are
1430 : * too numerous to list here.
1431 : *
1432 : * The following GML2 elements are parsed : Point, LineString, Polygon,
1433 : * MultiPoint, MultiLineString, MultiPolygon, MultiGeometry.
1434 : *
1435 : * (OGR >= 1.8.0) The following GML3 elements are parsed : Surface, MultiSurface,
1436 : * PolygonPatch, Triangle, Rectangle, Curve, MultiCurve, LineStringSegment, Arc,
1437 : * Circle, CompositeSurface, OrientableSurface, Solid, Tin, TriangulatedSurface.
1438 : *
1439 : * Arc and Circle elements are stroked to linestring, by using a
1440 : * 4 degrees step, unless the user has overridden the value with the
1441 : * OGR_ARC_STEPSIZE configuration variable.
1442 : *
1443 : * The C function OGR_G_CreateFromGML() is the same as this method.
1444 : *
1445 : * @param pszData The GML fragment for the geometry.
1446 : *
1447 : * @return a geometry on succes, or NULL on error.
1448 : */
1449 :
1450 0 : OGRGeometry *OGRGeometryFactory::createFromGML( const char *pszData )
1451 :
1452 : {
1453 : OGRGeometryH hGeom;
1454 :
1455 0 : hGeom = OGR_G_CreateFromGML( pszData );
1456 :
1457 0 : return (OGRGeometry *) hGeom;
1458 : }
1459 :
1460 : /************************************************************************/
1461 : /* createFromGEOS() */
1462 : /************************************************************************/
1463 :
1464 : OGRGeometry *
1465 439 : OGRGeometryFactory::createFromGEOS( GEOSGeom geosGeom )
1466 :
1467 : {
1468 : #ifndef HAVE_GEOS
1469 :
1470 : CPLError( CE_Failure, CPLE_NotSupported,
1471 : "GEOS support not enabled." );
1472 : return NULL;
1473 :
1474 : #else
1475 :
1476 439 : size_t nSize = 0;
1477 439 : unsigned char *pabyBuf = NULL;
1478 439 : OGRGeometry *poGeometry = NULL;
1479 :
1480 : /* Special case as POINT EMPTY cannot be translated to WKB */
1481 439 : if (GEOSGeomTypeId(geosGeom) == GEOS_POINT &&
1482 : GEOSisEmpty(geosGeom))
1483 0 : return new OGRPoint();
1484 :
1485 : #if GEOS_VERSION_MAJOR > 3 || (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 3)
1486 : /* GEOSGeom_getCoordinateDimension only available in GEOS 3.3.0 (unreleased at time of writing) */
1487 : int nCoordDim = GEOSGeom_getCoordinateDimension(geosGeom);
1488 : GEOSWKBWriter* wkbwriter = GEOSWKBWriter_create();
1489 : GEOSWKBWriter_setOutputDimension(wkbwriter, nCoordDim);
1490 : pabyBuf = GEOSWKBWriter_write( wkbwriter, geosGeom, &nSize );
1491 : GEOSWKBWriter_destroy(wkbwriter);
1492 : #else
1493 439 : pabyBuf = GEOSGeomToWKB_buf( geosGeom, &nSize );
1494 : #endif
1495 439 : if( pabyBuf == NULL || nSize == 0 )
1496 : {
1497 0 : return NULL;
1498 : }
1499 :
1500 439 : if( OGRGeometryFactory::createFromWkb( (unsigned char *) pabyBuf,
1501 : NULL, &poGeometry, (int) nSize )
1502 : != OGRERR_NONE )
1503 : {
1504 0 : poGeometry = NULL;
1505 : }
1506 :
1507 439 : if( pabyBuf != NULL )
1508 : {
1509 : /* Since GEOS 3.1.1, so we test 3.2.0 */
1510 : #if GEOS_CAPI_VERSION_MAJOR >= 2 || (GEOS_CAPI_VERSION_MAJOR == 1 && GEOS_CAPI_VERSION_MINOR >= 6)
1511 439 : GEOSFree( pabyBuf );
1512 : #else
1513 : free( pabyBuf );
1514 : #endif
1515 : }
1516 :
1517 439 : return poGeometry;
1518 :
1519 : #endif /* HAVE_GEOS */
1520 : }
1521 :
1522 : /************************************************************************/
1523 : /* getGEOSGeometryFactory() */
1524 : /************************************************************************/
1525 :
1526 0 : void *OGRGeometryFactory::getGEOSGeometryFactory()
1527 :
1528 : {
1529 : // XXX - mloskot - What to do with this call
1530 : // after GEOS C++ API has been stripped?
1531 0 : return NULL;
1532 : }
1533 :
1534 : /************************************************************************/
1535 : /* haveGEOS() */
1536 : /************************************************************************/
1537 :
1538 : /**
1539 : * \brief Test if GEOS enabled.
1540 : *
1541 : * This static method returns TRUE if GEOS support is built into OGR,
1542 : * otherwise it returns FALSE.
1543 : *
1544 : * @return TRUE if available, otherwise FALSE.
1545 : */
1546 :
1547 1828 : int OGRGeometryFactory::haveGEOS()
1548 :
1549 : {
1550 : #ifndef HAVE_GEOS
1551 : return FALSE;
1552 : #else
1553 1828 : return TRUE;
1554 : #endif
1555 : }
1556 :
1557 : /************************************************************************/
1558 : /* createFromFgf() */
1559 : /************************************************************************/
1560 :
1561 : /**
1562 : * \brief Create a geometry object of the appropriate type from it's FGF (FDO Geometry Format) binary representation.
1563 : *
1564 : * Also note that this is a static method, and that there
1565 : * is no need to instantiate an OGRGeometryFactory object.
1566 : *
1567 : * The C function OGR_G_CreateFromFgf() is the same as this method.
1568 : *
1569 : * @param pabyData pointer to the input BLOB data.
1570 : * @param poSR pointer to the spatial reference to be assigned to the
1571 : * created geometry object. This may be NULL.
1572 : * @param ppoReturn the newly created geometry object will be assigned to the
1573 : * indicated pointer on return. This will be NULL in case
1574 : * of failure.
1575 : * @param nBytes the number of bytes available in pabyData.
1576 : * @param pnBytesConsumed if not NULL, it will be set to the number of bytes
1577 : * consumed (at most nBytes).
1578 : *
1579 : * @return OGRERR_NONE if all goes well, otherwise any of
1580 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
1581 : * OGRERR_CORRUPT_DATA may be returned.
1582 : */
1583 :
1584 22 : OGRErr OGRGeometryFactory::createFromFgf( unsigned char *pabyData,
1585 : OGRSpatialReference * poSR,
1586 : OGRGeometry **ppoReturn,
1587 : int nBytes,
1588 : int *pnBytesConsumed )
1589 :
1590 : {
1591 : return createFromFgfInternal(pabyData, poSR, ppoReturn, nBytes,
1592 22 : pnBytesConsumed, 0);
1593 : }
1594 :
1595 :
1596 : /************************************************************************/
1597 : /* createFromFgfInternal() */
1598 : /************************************************************************/
1599 :
1600 28 : OGRErr OGRGeometryFactory::createFromFgfInternal( unsigned char *pabyData,
1601 : OGRSpatialReference * poSR,
1602 : OGRGeometry **ppoReturn,
1603 : int nBytes,
1604 : int *pnBytesConsumed,
1605 : int nRecLevel )
1606 : {
1607 28 : OGRErr eErr = OGRERR_NONE;
1608 28 : OGRGeometry *poGeom = NULL;
1609 : GInt32 nGType, nGDim;
1610 28 : int nTupleSize = 0;
1611 28 : int iOrdinal = 0;
1612 :
1613 : (void) iOrdinal;
1614 :
1615 : /* Arbitrary value, but certainly large enough for reasonable usages ! */
1616 28 : if( nRecLevel == 32 )
1617 : {
1618 : CPLError( CE_Failure, CPLE_AppDefined,
1619 : "Too many recursiong level (%d) while parsing FGF geometry.",
1620 0 : nRecLevel );
1621 0 : return OGRERR_CORRUPT_DATA;
1622 : }
1623 :
1624 28 : *ppoReturn = NULL;
1625 :
1626 28 : if( nBytes < 4 )
1627 0 : return OGRERR_NOT_ENOUGH_DATA;
1628 :
1629 : /* -------------------------------------------------------------------- */
1630 : /* Decode the geometry type. */
1631 : /* -------------------------------------------------------------------- */
1632 28 : memcpy( &nGType, pabyData + 0, 4 );
1633 : CPL_LSBPTR32( &nGType );
1634 :
1635 28 : if( nGType < 0 || nGType > 13 )
1636 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
1637 :
1638 : /* -------------------------------------------------------------------- */
1639 : /* Decode the dimentionality if appropriate. */
1640 : /* -------------------------------------------------------------------- */
1641 28 : switch( nGType )
1642 : {
1643 : case 1: // Point
1644 : case 2: // LineString
1645 : case 3: // Polygon
1646 :
1647 18 : if( nBytes < 8 )
1648 0 : return OGRERR_NOT_ENOUGH_DATA;
1649 :
1650 18 : memcpy( &nGDim, pabyData + 4, 4 );
1651 : CPL_LSBPTR32( &nGDim );
1652 :
1653 18 : if( nGDim < 0 || nGDim > 3 )
1654 0 : return OGRERR_CORRUPT_DATA;
1655 :
1656 18 : nTupleSize = 2;
1657 18 : if( nGDim & 0x01 ) // Z
1658 2 : nTupleSize++;
1659 18 : if( nGDim & 0x02 ) // M
1660 0 : nTupleSize++;
1661 :
1662 : break;
1663 :
1664 : default:
1665 : break;
1666 : }
1667 :
1668 : /* -------------------------------------------------------------------- */
1669 : /* None */
1670 : /* -------------------------------------------------------------------- */
1671 28 : if( nGType == 0 )
1672 : {
1673 0 : if( pnBytesConsumed )
1674 0 : *pnBytesConsumed = 4;
1675 : }
1676 :
1677 : /* -------------------------------------------------------------------- */
1678 : /* Point */
1679 : /* -------------------------------------------------------------------- */
1680 28 : else if( nGType == 1 )
1681 : {
1682 : double adfTuple[4];
1683 :
1684 6 : if( nBytes < nTupleSize * 8 + 8 )
1685 0 : return OGRERR_NOT_ENOUGH_DATA;
1686 :
1687 6 : memcpy( adfTuple, pabyData + 8, nTupleSize*8 );
1688 : #ifdef CPL_MSB
1689 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1690 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1691 : #endif
1692 6 : if( nTupleSize > 2 )
1693 2 : poGeom = new OGRPoint( adfTuple[0], adfTuple[1], adfTuple[2] );
1694 : else
1695 4 : poGeom = new OGRPoint( adfTuple[0], adfTuple[1] );
1696 :
1697 6 : if( pnBytesConsumed )
1698 2 : *pnBytesConsumed = 8 + nTupleSize * 8;
1699 : }
1700 :
1701 : /* -------------------------------------------------------------------- */
1702 : /* LineString */
1703 : /* -------------------------------------------------------------------- */
1704 22 : else if( nGType == 2 )
1705 : {
1706 : double adfTuple[4];
1707 : GInt32 nPointCount;
1708 : int iPoint;
1709 : OGRLineString *poLS;
1710 :
1711 4 : if( nBytes < 12 )
1712 0 : return OGRERR_NOT_ENOUGH_DATA;
1713 :
1714 4 : memcpy( &nPointCount, pabyData + 8, 4 );
1715 : CPL_LSBPTR32( &nPointCount );
1716 :
1717 4 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
1718 0 : return OGRERR_CORRUPT_DATA;
1719 :
1720 4 : if( nBytes - 12 < nTupleSize * 8 * nPointCount )
1721 0 : return OGRERR_NOT_ENOUGH_DATA;
1722 :
1723 4 : poGeom = poLS = new OGRLineString();
1724 4 : poLS->setNumPoints( nPointCount );
1725 :
1726 8 : for( iPoint = 0; iPoint < nPointCount; iPoint++ )
1727 : {
1728 : memcpy( adfTuple, pabyData + 12 + 8*nTupleSize*iPoint,
1729 4 : nTupleSize*8 );
1730 : #ifdef CPL_MSB
1731 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1732 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1733 : #endif
1734 4 : if( nTupleSize > 2 )
1735 0 : poLS->setPoint( iPoint, adfTuple[0], adfTuple[1], adfTuple[2] );
1736 : else
1737 4 : poLS->setPoint( iPoint, adfTuple[0], adfTuple[1] );
1738 : }
1739 :
1740 4 : if( pnBytesConsumed )
1741 0 : *pnBytesConsumed = 12 + nTupleSize * 8 * nPointCount;
1742 : }
1743 :
1744 : /* -------------------------------------------------------------------- */
1745 : /* Polygon */
1746 : /* -------------------------------------------------------------------- */
1747 18 : else if( nGType == 3 )
1748 : {
1749 : double adfTuple[4];
1750 : GInt32 nPointCount;
1751 : GInt32 nRingCount;
1752 : int iPoint, iRing;
1753 : OGRLinearRing *poLR;
1754 : OGRPolygon *poPoly;
1755 : int nNextByte;
1756 :
1757 8 : if( nBytes < 12 )
1758 0 : return OGRERR_NOT_ENOUGH_DATA;
1759 :
1760 8 : memcpy( &nRingCount, pabyData + 8, 4 );
1761 : CPL_LSBPTR32( &nRingCount );
1762 :
1763 8 : if (nRingCount < 0 || nRingCount > INT_MAX / 4)
1764 0 : return OGRERR_CORRUPT_DATA;
1765 :
1766 : /* Each ring takes at least 4 bytes */
1767 8 : if (nBytes - 12 < nRingCount * 4)
1768 0 : return OGRERR_NOT_ENOUGH_DATA;
1769 :
1770 8 : nNextByte = 12;
1771 :
1772 8 : poGeom = poPoly = new OGRPolygon();
1773 :
1774 20 : for( iRing = 0; iRing < nRingCount; iRing++ )
1775 : {
1776 12 : if( nBytes - nNextByte < 4 )
1777 : {
1778 0 : delete poGeom;
1779 0 : return OGRERR_NOT_ENOUGH_DATA;
1780 : }
1781 :
1782 12 : memcpy( &nPointCount, pabyData + nNextByte, 4 );
1783 : CPL_LSBPTR32( &nPointCount );
1784 :
1785 12 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
1786 : {
1787 0 : delete poGeom;
1788 0 : return OGRERR_CORRUPT_DATA;
1789 : }
1790 :
1791 12 : nNextByte += 4;
1792 :
1793 12 : if( nBytes - nNextByte < nTupleSize * 8 * nPointCount )
1794 : {
1795 0 : delete poGeom;
1796 0 : return OGRERR_NOT_ENOUGH_DATA;
1797 : }
1798 :
1799 12 : poLR = new OGRLinearRing();
1800 12 : poLR->setNumPoints( nPointCount );
1801 :
1802 24 : for( iPoint = 0; iPoint < nPointCount; iPoint++ )
1803 : {
1804 12 : memcpy( adfTuple, pabyData + nNextByte, nTupleSize*8 );
1805 12 : nNextByte += nTupleSize * 8;
1806 :
1807 : #ifdef CPL_MSB
1808 : for( iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++ )
1809 : CPL_SWAP64PTR( adfTuple + iOrdinal );
1810 : #endif
1811 12 : if( nTupleSize > 2 )
1812 0 : poLR->setPoint( iPoint, adfTuple[0], adfTuple[1], adfTuple[2] );
1813 : else
1814 12 : poLR->setPoint( iPoint, adfTuple[0], adfTuple[1] );
1815 : }
1816 :
1817 12 : poPoly->addRingDirectly( poLR );
1818 : }
1819 :
1820 8 : if( pnBytesConsumed )
1821 4 : *pnBytesConsumed = nNextByte;
1822 : }
1823 :
1824 : /* -------------------------------------------------------------------- */
1825 : /* GeometryCollections of various kinds. */
1826 : /* -------------------------------------------------------------------- */
1827 14 : else if( nGType == 4 // MultiPoint
1828 : || nGType == 5 // MultiLineString
1829 : || nGType == 6 // MultiPolygon
1830 : || nGType == 7 ) // MultiGeometry
1831 : {
1832 10 : OGRGeometryCollection *poGC = NULL;
1833 : GInt32 nGeomCount;
1834 : int iGeom, nBytesUsed;
1835 :
1836 10 : if( nBytes < 8 )
1837 0 : return OGRERR_NOT_ENOUGH_DATA;
1838 :
1839 10 : memcpy( &nGeomCount, pabyData + 4, 4 );
1840 : CPL_LSBPTR32( &nGeomCount );
1841 :
1842 10 : if (nGeomCount < 0 || nGeomCount > INT_MAX / 4)
1843 0 : return OGRERR_CORRUPT_DATA;
1844 :
1845 : /* Each geometry takes at least 4 bytes */
1846 10 : if (nBytes - 8 < 4 * nGeomCount)
1847 4 : return OGRERR_NOT_ENOUGH_DATA;
1848 :
1849 6 : nBytesUsed = 8;
1850 :
1851 6 : if( nGType == 4 )
1852 0 : poGC = new OGRMultiPoint();
1853 6 : else if( nGType == 5 )
1854 0 : poGC = new OGRMultiLineString();
1855 6 : else if( nGType == 6 )
1856 2 : poGC = new OGRMultiPolygon();
1857 4 : else if( nGType == 7 )
1858 4 : poGC = new OGRGeometryCollection();
1859 :
1860 10 : for( iGeom = 0; iGeom < nGeomCount; iGeom++ )
1861 : {
1862 : int nThisGeomSize;
1863 6 : OGRGeometry *poThisGeom = NULL;
1864 :
1865 : eErr = createFromFgfInternal( pabyData + nBytesUsed, poSR, &poThisGeom,
1866 6 : nBytes - nBytesUsed, &nThisGeomSize, nRecLevel + 1);
1867 6 : if( eErr != OGRERR_NONE )
1868 : {
1869 0 : delete poGC;
1870 0 : return eErr;
1871 : }
1872 :
1873 6 : nBytesUsed += nThisGeomSize;
1874 6 : eErr = poGC->addGeometryDirectly( poThisGeom );
1875 6 : if( eErr != OGRERR_NONE )
1876 : {
1877 2 : delete poGC;
1878 2 : delete poThisGeom;
1879 2 : return eErr;
1880 : }
1881 : }
1882 :
1883 4 : poGeom = poGC;
1884 4 : if( pnBytesConsumed )
1885 0 : *pnBytesConsumed = nBytesUsed;
1886 : }
1887 :
1888 : /* -------------------------------------------------------------------- */
1889 : /* Currently unsupported geometry. */
1890 : /* */
1891 : /* We need to add 10/11/12/13 curve types in some fashion. */
1892 : /* -------------------------------------------------------------------- */
1893 : else
1894 : {
1895 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
1896 : }
1897 :
1898 : /* -------------------------------------------------------------------- */
1899 : /* Assign spatial reference system. */
1900 : /* -------------------------------------------------------------------- */
1901 22 : if( eErr == OGRERR_NONE )
1902 : {
1903 22 : if( poGeom != NULL && poSR )
1904 0 : poGeom->assignSpatialReference( poSR );
1905 22 : *ppoReturn = poGeom;
1906 : }
1907 : else
1908 : {
1909 0 : delete poGeom;
1910 : }
1911 :
1912 22 : return eErr;
1913 : }
1914 :
1915 : /************************************************************************/
1916 : /* OGR_G_CreateFromFgf() */
1917 : /************************************************************************/
1918 :
1919 0 : OGRErr CPL_DLL OGR_G_CreateFromFgf( unsigned char *pabyData,
1920 : OGRSpatialReferenceH hSRS,
1921 : OGRGeometryH *phGeometry,
1922 : int nBytes, int *pnBytesConsumed )
1923 :
1924 : {
1925 : return OGRGeometryFactory::createFromFgf( pabyData,
1926 : (OGRSpatialReference *) hSRS,
1927 : (OGRGeometry **) phGeometry,
1928 0 : nBytes, pnBytesConsumed );
1929 : }
1930 :
1931 : /************************************************************************/
1932 : /* SplitLineStringAtDateline() */
1933 : /************************************************************************/
1934 :
1935 : #define SWAP_DBL(a,b) do { double tmp = a; a = b; b = tmp; } while(0)
1936 :
1937 2 : static void SplitLineStringAtDateline(OGRGeometryCollection* poMulti,
1938 : const OGRLineString* poLS)
1939 : {
1940 : int i;
1941 2 : int bIs3D = poLS->getCoordinateDimension() == 3;
1942 2 : OGRLineString* poNewLS = new OGRLineString();
1943 2 : poMulti->addGeometryDirectly(poNewLS);
1944 24 : for(i=0;i<poLS->getNumPoints();i++)
1945 : {
1946 22 : double dfX = poLS->getX(i);
1947 22 : if (i > 0 && fabs(dfX - poLS->getX(i-1)) > 350)
1948 : {
1949 4 : double dfX1 = poLS->getX(i-1);
1950 4 : double dfY1 = poLS->getY(i-1);
1951 4 : double dfZ1 = poLS->getY(i-1);
1952 4 : double dfX2 = poLS->getX(i);
1953 4 : double dfY2 = poLS->getY(i);
1954 4 : double dfZ2 = poLS->getY(i);
1955 :
1956 4 : if (dfX1 > -180 && dfX1 < -170 && dfX2 == 180 &&
1957 : i+1 < poLS->getNumPoints() &&
1958 : poLS->getX(i+1) > -180 && poLS->getX(i+1) < -170)
1959 : {
1960 0 : if( bIs3D )
1961 0 : poNewLS->addPoint(-180, poLS->getY(i), poLS->getZ(i));
1962 : else
1963 0 : poNewLS->addPoint(-180, poLS->getY(i));
1964 :
1965 0 : i++;
1966 :
1967 0 : if( bIs3D )
1968 0 : poNewLS->addPoint(poLS->getX(i), poLS->getY(i), poLS->getZ(i));
1969 : else
1970 0 : poNewLS->addPoint(poLS->getX(i), poLS->getY(i));
1971 0 : continue;
1972 : }
1973 4 : else if (dfX1 > 170 && dfX1 < 180 && dfX2 == -180 &&
1974 : i+1 < poLS->getNumPoints() &&
1975 : poLS->getX(i+1) > 170 && poLS->getX(i+1) < 180)
1976 : {
1977 0 : if( bIs3D )
1978 0 : poNewLS->addPoint(180, poLS->getY(i), poLS->getZ(i));
1979 : else
1980 0 : poNewLS->addPoint(180, poLS->getY(i));
1981 :
1982 0 : i++;
1983 :
1984 0 : if( bIs3D )
1985 0 : poNewLS->addPoint(poLS->getX(i), poLS->getY(i), poLS->getZ(i));
1986 : else
1987 0 : poNewLS->addPoint(poLS->getX(i), poLS->getY(i));
1988 0 : continue;
1989 : }
1990 :
1991 4 : if (dfX1 < -170 && dfX2 > 170)
1992 : {
1993 2 : SWAP_DBL(dfX1, dfX2);
1994 2 : SWAP_DBL(dfY1, dfY2);
1995 2 : SWAP_DBL(dfZ1, dfZ2);
1996 : }
1997 4 : if (dfX1 > 170 && dfX2 < -170)
1998 4 : dfX2 += 360;
1999 :
2000 8 : if (dfX1 <= 180 && dfX2 >= 180 && dfX1 < dfX2)
2001 : {
2002 4 : double dfRatio = (180 - dfX1) / (dfX2 - dfX1);
2003 4 : double dfY = dfRatio * dfY2 + (1 - dfRatio) * dfY1;
2004 4 : double dfZ = dfRatio * dfZ2 + (1 - dfRatio) * dfZ1;
2005 4 : if( bIs3D )
2006 0 : poNewLS->addPoint(poLS->getX(i-1) > 170 ? 180 : -180, dfY, dfZ);
2007 : else
2008 4 : poNewLS->addPoint(poLS->getX(i-1) > 170 ? 180 : -180, dfY);
2009 4 : poNewLS = new OGRLineString();
2010 4 : if( bIs3D )
2011 0 : poNewLS->addPoint(poLS->getX(i-1) > 170 ? -180 : 180, dfY, dfZ);
2012 : else
2013 4 : poNewLS->addPoint(poLS->getX(i-1) > 170 ? -180 : 180, dfY);
2014 4 : poMulti->addGeometryDirectly(poNewLS);
2015 : }
2016 : else
2017 : {
2018 0 : poNewLS = new OGRLineString();
2019 0 : poMulti->addGeometryDirectly(poNewLS);
2020 : }
2021 : }
2022 22 : if( bIs3D )
2023 0 : poNewLS->addPoint(dfX, poLS->getY(i), poLS->getZ(i));
2024 : else
2025 22 : poNewLS->addPoint(dfX, poLS->getY(i));
2026 : }
2027 2 : }
2028 :
2029 : /************************************************************************/
2030 : /* FixPolygonCoordinatesAtDateLine() */
2031 : /************************************************************************/
2032 :
2033 : #ifdef HAVE_GEOS
2034 6 : static void FixPolygonCoordinatesAtDateLine(OGRPolygon* poPoly)
2035 : {
2036 : int i, iPart;
2037 12 : for(iPart = 0; iPart < 1 + poPoly->getNumInteriorRings(); iPart++)
2038 : {
2039 : OGRLineString* poLS = (iPart == 0) ? poPoly->getExteriorRing() :
2040 6 : poPoly->getInteriorRing(iPart-1);
2041 6 : int bGoEast = FALSE;
2042 6 : int bIs3D = poLS->getCoordinateDimension() == 3;
2043 62 : for(i=1;i<poLS->getNumPoints();i++)
2044 : {
2045 56 : double dfX = poLS->getX(i);
2046 56 : double dfPrevX = poLS->getX(i-1);
2047 56 : double dfDiffLong = fabs(dfX - dfPrevX);
2048 56 : if (dfDiffLong > 350)
2049 : {
2050 58 : if ((dfPrevX > 170 && dfX < -170) || (dfX < 0 && bGoEast))
2051 : {
2052 28 : dfX += 360;
2053 28 : bGoEast = TRUE;
2054 28 : if( bIs3D )
2055 0 : poLS->setPoint(i, dfX, poLS->getY(i), poLS->getZ(i));
2056 : else
2057 28 : poLS->setPoint(i, dfX, poLS->getY(i));
2058 : }
2059 4 : else if (dfPrevX < -170 && dfX > 170)
2060 : {
2061 : int j;
2062 12 : for(j=i-1;j>=0;j--)
2063 : {
2064 10 : dfX = poLS->getX(j);
2065 10 : if (dfX < 0)
2066 : {
2067 10 : if( bIs3D )
2068 0 : poLS->setPoint(j, dfX + 360, poLS->getY(j), poLS->getZ(j));
2069 : else
2070 10 : poLS->setPoint(j, dfX + 360, poLS->getY(j));
2071 : }
2072 : }
2073 2 : bGoEast = FALSE;
2074 : }
2075 : else
2076 : {
2077 0 : bGoEast = FALSE;
2078 : }
2079 : }
2080 : }
2081 : }
2082 6 : }
2083 : #endif
2084 :
2085 : /************************************************************************/
2086 : /* Sub360ToLon() */
2087 : /************************************************************************/
2088 :
2089 12 : static void Sub360ToLon( OGRGeometry* poGeom )
2090 : {
2091 12 : switch (wkbFlatten(poGeom->getGeometryType()))
2092 : {
2093 : case wkbPolygon:
2094 : case wkbMultiLineString:
2095 : case wkbMultiPolygon:
2096 : case wkbGeometryCollection:
2097 : {
2098 6 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
2099 12 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
2100 : {
2101 6 : Sub360ToLon((OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom));
2102 : }
2103 :
2104 6 : break;
2105 : }
2106 :
2107 : case wkbLineString:
2108 : case wkbLinearRing:
2109 : {
2110 6 : OGRLineString* poLineString = (OGRLineString* )poGeom;
2111 6 : int nPointCount = poLineString->getNumPoints();
2112 6 : int nCoordDim = poLineString->getCoordinateDimension();
2113 60 : for( int iPoint = 0; iPoint < nPointCount; iPoint++)
2114 : {
2115 54 : if (nCoordDim == 2)
2116 : poLineString->setPoint(iPoint,
2117 : poLineString->getX(iPoint) - 360,
2118 54 : poLineString->getY(iPoint));
2119 : else
2120 : poLineString->setPoint(iPoint,
2121 : poLineString->getX(iPoint) - 360,
2122 : poLineString->getY(iPoint),
2123 0 : poLineString->getZ(iPoint));
2124 : }
2125 : break;
2126 : }
2127 :
2128 : default:
2129 : break;
2130 : }
2131 12 : }
2132 :
2133 : /************************************************************************/
2134 : /* AddSimpleGeomToMulti() */
2135 : /************************************************************************/
2136 :
2137 12 : static void AddSimpleGeomToMulti(OGRGeometryCollection* poMulti,
2138 : const OGRGeometry* poGeom)
2139 : {
2140 12 : switch (wkbFlatten(poGeom->getGeometryType()))
2141 : {
2142 : case wkbPolygon:
2143 : case wkbLineString:
2144 12 : poMulti->addGeometry(poGeom);
2145 12 : break;
2146 :
2147 : case wkbMultiLineString:
2148 : case wkbMultiPolygon:
2149 : case wkbGeometryCollection:
2150 : {
2151 0 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
2152 0 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
2153 : {
2154 : OGRGeometry* poSubGeom =
2155 0 : (OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom);
2156 0 : AddSimpleGeomToMulti(poMulti, poSubGeom);
2157 : }
2158 : break;
2159 : }
2160 :
2161 : default:
2162 : break;
2163 : }
2164 12 : }
2165 :
2166 : /************************************************************************/
2167 : /* CutGeometryOnDateLineAndAddToMulti() */
2168 : /************************************************************************/
2169 :
2170 8 : static void CutGeometryOnDateLineAndAddToMulti(OGRGeometryCollection* poMulti,
2171 : const OGRGeometry* poGeom)
2172 : {
2173 8 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
2174 8 : switch (eGeomType)
2175 : {
2176 : case wkbPolygon:
2177 : case wkbLineString:
2178 : {
2179 8 : int bWrapDateline = FALSE;
2180 8 : int bSplitLineStringAtDateline = FALSE;
2181 8 : OGREnvelope oEnvelope;
2182 :
2183 8 : poGeom->getEnvelope(&oEnvelope);
2184 :
2185 : /* Naive heuristics... Place to improvement... */
2186 8 : OGRGeometry* poDupGeom = NULL;
2187 :
2188 8 : if (oEnvelope.MinX > 170 && oEnvelope.MaxX > 180)
2189 : {
2190 : #ifndef HAVE_GEOS
2191 : CPLError( CE_Failure, CPLE_NotSupported,
2192 : "GEOS support not enabled." );
2193 : #else
2194 0 : bWrapDateline = TRUE;
2195 : #endif
2196 : }
2197 : else
2198 : {
2199 : OGRLineString* poLS;
2200 8 : if (eGeomType == wkbPolygon)
2201 6 : poLS = ((OGRPolygon*)poGeom)->getExteriorRing();
2202 : else
2203 2 : poLS = (OGRLineString*)poGeom;
2204 8 : if (poLS)
2205 : {
2206 : int i;
2207 8 : double dfMaxSmallDiffLong = 0;
2208 8 : int bHasBigDiff = FALSE;
2209 : /* Detect big gaps in longitude */
2210 84 : for(i=1;i<poLS->getNumPoints();i++)
2211 : {
2212 76 : double dfPrevX = poLS->getX(i-1);
2213 76 : double dfX = poLS->getX(i);
2214 76 : double dfDiffLong = fabs(dfX - dfPrevX);
2215 92 : if (dfDiffLong > 350 &&
2216 : ((dfX > 170 && dfPrevX < -170) || (dfPrevX > 170 && dfX < -170)))
2217 16 : bHasBigDiff = TRUE;
2218 60 : else if (dfDiffLong > dfMaxSmallDiffLong)
2219 14 : dfMaxSmallDiffLong = dfDiffLong;
2220 : }
2221 8 : if (bHasBigDiff && dfMaxSmallDiffLong < 10)
2222 : {
2223 8 : if (eGeomType == wkbLineString)
2224 2 : bSplitLineStringAtDateline = TRUE;
2225 : else
2226 : {
2227 : #ifndef HAVE_GEOS
2228 : CPLError( CE_Failure, CPLE_NotSupported,
2229 : "GEOS support not enabled." );
2230 : #else
2231 6 : bWrapDateline = TRUE;
2232 6 : poDupGeom = poGeom->clone();
2233 6 : FixPolygonCoordinatesAtDateLine((OGRPolygon*)poDupGeom);
2234 : #endif
2235 : }
2236 : }
2237 : }
2238 : }
2239 :
2240 8 : if (bSplitLineStringAtDateline)
2241 : {
2242 2 : SplitLineStringAtDateline(poMulti, (OGRLineString*)poGeom);
2243 : }
2244 6 : else if (bWrapDateline)
2245 : {
2246 6 : const OGRGeometry* poWorkGeom = (poDupGeom) ? poDupGeom : poGeom;
2247 6 : OGRGeometry* poRectangle1 = NULL;
2248 6 : OGRGeometry* poRectangle2 = NULL;
2249 6 : const char* pszWKT1 = "POLYGON((0 90,180 90,180 -90,0 -90,0 90))";
2250 6 : const char* pszWKT2 = "POLYGON((180 90,360 90,360 -90,180 -90,180 90))";
2251 6 : OGRGeometryFactory::createFromWkt((char**)&pszWKT1, NULL, &poRectangle1);
2252 6 : OGRGeometryFactory::createFromWkt((char**)&pszWKT2, NULL, &poRectangle2);
2253 6 : OGRGeometry* poGeom1 = poWorkGeom->Intersection(poRectangle1);
2254 6 : OGRGeometry* poGeom2 = poWorkGeom->Intersection(poRectangle2);
2255 6 : delete poRectangle1;
2256 6 : delete poRectangle2;
2257 :
2258 12 : if (poGeom1 != NULL && poGeom2 != NULL)
2259 : {
2260 6 : AddSimpleGeomToMulti(poMulti, poGeom1);
2261 6 : Sub360ToLon(poGeom2);
2262 6 : AddSimpleGeomToMulti(poMulti, poGeom2);
2263 : }
2264 : else
2265 : {
2266 0 : AddSimpleGeomToMulti(poMulti, poGeom);
2267 : }
2268 :
2269 6 : delete poGeom1;
2270 6 : delete poGeom2;
2271 6 : delete poDupGeom;
2272 : }
2273 : else
2274 : {
2275 0 : poMulti->addGeometry(poGeom);
2276 : }
2277 8 : break;
2278 : }
2279 :
2280 : case wkbMultiLineString:
2281 : case wkbMultiPolygon:
2282 : case wkbGeometryCollection:
2283 : {
2284 0 : int nSubGeomCount = OGR_G_GetGeometryCount((OGRGeometryH)poGeom);
2285 0 : for( int iGeom = 0; iGeom < nSubGeomCount; iGeom++ )
2286 : {
2287 : OGRGeometry* poSubGeom =
2288 0 : (OGRGeometry*)OGR_G_GetGeometryRef((OGRGeometryH)poGeom, iGeom);
2289 0 : CutGeometryOnDateLineAndAddToMulti(poMulti, poSubGeom);
2290 : }
2291 : break;
2292 : }
2293 :
2294 : default:
2295 : break;
2296 : }
2297 8 : }
2298 :
2299 : /************************************************************************/
2300 : /* transformWithOptions() */
2301 : /************************************************************************/
2302 :
2303 28 : OGRGeometry* OGRGeometryFactory::transformWithOptions( const OGRGeometry* poSrcGeom,
2304 : OGRCoordinateTransformation *poCT,
2305 : char** papszOptions )
2306 : {
2307 28 : OGRGeometry* poDstGeom = poSrcGeom->clone();
2308 28 : if (poCT != NULL)
2309 : {
2310 22 : OGRErr eErr = poDstGeom->transform(poCT);
2311 22 : if (eErr != OGRERR_NONE)
2312 : {
2313 0 : delete poDstGeom;
2314 0 : return NULL;
2315 : }
2316 : }
2317 :
2318 28 : if (CSLTestBoolean(CSLFetchNameValueDef(papszOptions, "WRAPDATELINE", "NO")))
2319 : {
2320 8 : OGRwkbGeometryType eType = wkbFlatten(poSrcGeom->getGeometryType());
2321 : OGRwkbGeometryType eNewType;
2322 14 : if (eType == wkbPolygon || eType == wkbMultiPolygon)
2323 6 : eNewType = wkbMultiPolygon;
2324 4 : else if (eType == wkbLineString || eType == wkbMultiLineString)
2325 2 : eNewType = wkbMultiLineString;
2326 : else
2327 0 : eNewType = wkbGeometryCollection;
2328 :
2329 : OGRGeometryCollection* poMulti =
2330 8 : (OGRGeometryCollection* )createGeometry(eNewType);
2331 :
2332 8 : CutGeometryOnDateLineAndAddToMulti(poMulti, poDstGeom);
2333 :
2334 8 : if (poMulti->getNumGeometries() == 0)
2335 : {
2336 0 : delete poMulti;
2337 : }
2338 8 : else if (poMulti->getNumGeometries() == 1)
2339 : {
2340 0 : delete poDstGeom;
2341 0 : poDstGeom = poMulti->getGeometryRef(0)->clone();
2342 0 : delete poMulti;
2343 : }
2344 : else
2345 : {
2346 8 : delete poDstGeom;
2347 8 : poDstGeom = poMulti;
2348 : }
2349 : }
2350 :
2351 28 : return poDstGeom;
2352 : }
2353 :
2354 : /************************************************************************/
2355 : /* approximateArcAngles() */
2356 : /************************************************************************/
2357 :
2358 : /**
2359 : * Stroke arc to linestring.
2360 : *
2361 : * Stroke an arc of a circle to a linestring based on a center
2362 : * point, radius, start angle and end angle, all angles in degrees.
2363 : *
2364 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
2365 : * used. This is currently 4 degrees unless the user has overridden the
2366 : * value with the OGR_ARC_STEPSIZE configuration variable.
2367 : *
2368 : * @see CPLSetConfigOption()
2369 : *
2370 : * @param dfCenterX center X
2371 : * @param dfCenterY center Y
2372 : * @param dfZ center Z
2373 : * @param dfPrimaryRadius X radius of ellipse.
2374 : * @param dfSecondaryRadius Y radius of ellipse.
2375 : * @param dfRotation rotation of the ellipse clockwise.
2376 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
2377 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
2378 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
2379 : * arc, zero to use the default setting.
2380 : *
2381 : * @return OGRLineString geometry representing an approximation of the arc.
2382 : *
2383 : * @since OGR 1.8.0
2384 : */
2385 :
2386 36 : OGRGeometry* OGRGeometryFactory::approximateArcAngles(
2387 : double dfCenterX, double dfCenterY, double dfZ,
2388 : double dfPrimaryRadius, double dfSecondaryRadius, double dfRotation,
2389 : double dfStartAngle, double dfEndAngle,
2390 : double dfMaxAngleStepSizeDegrees )
2391 :
2392 : {
2393 : double dfSlice;
2394 : int iPoint, nVertexCount;
2395 36 : OGRLineString *poLine = new OGRLineString();
2396 36 : double dfRotationRadians = dfRotation * PI / 180.0;
2397 :
2398 : // support default arc step setting.
2399 36 : if( dfMaxAngleStepSizeDegrees == 0 )
2400 : {
2401 : dfMaxAngleStepSizeDegrees =
2402 34 : atof(CPLGetConfigOption("OGR_ARC_STEPSIZE","4"));
2403 : }
2404 :
2405 : // switch direction
2406 36 : dfStartAngle *= -1;
2407 36 : dfEndAngle *= -1;
2408 :
2409 : // Figure out the number of slices to make this into.
2410 : nVertexCount = (int)
2411 36 : ceil(fabs(dfEndAngle - dfStartAngle)/dfMaxAngleStepSizeDegrees) + 1;
2412 36 : nVertexCount = MAX(2,nVertexCount);
2413 36 : dfSlice = (dfEndAngle-dfStartAngle)/(nVertexCount-1);
2414 :
2415 : /* -------------------------------------------------------------------- */
2416 : /* Compute the interpolated points. */
2417 : /* -------------------------------------------------------------------- */
2418 1856 : for( iPoint=0; iPoint < nVertexCount; iPoint++ )
2419 : {
2420 : double dfAngleOnEllipse;
2421 : double dfArcX, dfArcY;
2422 : double dfEllipseX, dfEllipseY;
2423 :
2424 1820 : dfAngleOnEllipse = (dfStartAngle + iPoint * dfSlice) * PI / 180.0;
2425 :
2426 : // Compute position on the unrotated ellipse.
2427 1820 : dfEllipseX = cos(dfAngleOnEllipse) * dfPrimaryRadius;
2428 1820 : dfEllipseY = sin(dfAngleOnEllipse) * dfSecondaryRadius;
2429 :
2430 : // Rotate this position around the center of the ellipse.
2431 : dfArcX = dfCenterX
2432 : + dfEllipseX * cos(dfRotationRadians)
2433 1820 : + dfEllipseY * sin(dfRotationRadians);
2434 : dfArcY = dfCenterY
2435 : - dfEllipseX * sin(dfRotationRadians)
2436 1820 : + dfEllipseY * cos(dfRotationRadians);
2437 :
2438 1820 : poLine->setPoint( iPoint, dfArcX, dfArcY, dfZ );
2439 : }
2440 :
2441 36 : return poLine;
2442 : }
2443 :
2444 : /************************************************************************/
2445 : /* OGR_G_ApproximateArcAngles() */
2446 : /************************************************************************/
2447 :
2448 : /**
2449 : * Stroke arc to linestring.
2450 : *
2451 : * Stroke an arc of a circle to a linestring based on a center
2452 : * point, radius, start angle and end angle, all angles in degrees.
2453 : *
2454 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
2455 : * used. This is currently 4 degrees unless the user has overridden the
2456 : * value with the OGR_ARC_STEPSIZE configuration variable.
2457 : *
2458 : * @see CPLSetConfigOption()
2459 : *
2460 : * @param dfCenterX center X
2461 : * @param dfCenterY center Y
2462 : * @param dfZ center Z
2463 : * @param dfPrimaryRadius X radius of ellipse.
2464 : * @param dfSecondaryRadius Y radius of ellipse.
2465 : * @param dfRotation rotation of the ellipse clockwise.
2466 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
2467 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
2468 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
2469 : * arc, zero to use the default setting.
2470 : *
2471 : * @return OGRLineString geometry representing an approximation of the arc.
2472 : *
2473 : * @since OGR 1.8.0
2474 : */
2475 :
2476 : OGRGeometryH CPL_DLL
2477 2 : OGR_G_ApproximateArcAngles(
2478 : double dfCenterX, double dfCenterY, double dfZ,
2479 : double dfPrimaryRadius, double dfSecondaryRadius, double dfRotation,
2480 : double dfStartAngle, double dfEndAngle,
2481 : double dfMaxAngleStepSizeDegrees )
2482 :
2483 : {
2484 : return (OGRGeometryH) OGRGeometryFactory::approximateArcAngles(
2485 : dfCenterX, dfCenterY, dfZ,
2486 : dfPrimaryRadius, dfSecondaryRadius, dfRotation,
2487 2 : dfStartAngle, dfEndAngle, dfMaxAngleStepSizeDegrees );
2488 : }
|
