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