Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: Factory for converting geometry to and from well known binary
5 : * format.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, Frank Warmerdam
10 : * Copyright (c) 2008-2014, Even Rouault <even dot rouault at spatialys dot com>
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 "cpl_port.h"
32 :
33 : #include "cpl_conv.h"
34 : #include "cpl_error.h"
35 : #include "cpl_string.h"
36 : #include "ogr_geometry.h"
37 : #include "ogr_api.h"
38 : #include "ogr_core.h"
39 : #include "ogr_geos.h"
40 : #include "ogr_sfcgal.h"
41 : #include "ogr_p.h"
42 : #include "ogr_spatialref.h"
43 : #include "ogr_srs_api.h"
44 : #ifdef HAVE_GEOS
45 : #include "geos_c.h"
46 : #endif
47 : #include "ogrgeojsonreader.h"
48 :
49 : #include <cassert>
50 : #include <climits>
51 : #include <cmath>
52 : #include <cstdlib>
53 : #include <cstring>
54 : #include <cstddef>
55 :
56 : #include <algorithm>
57 : #include <limits>
58 : #include <new>
59 : #include <utility>
60 : #include <vector>
61 :
62 : #ifndef HAVE_GEOS
63 : #define UNUSED_IF_NO_GEOS CPL_UNUSED
64 : #else
65 : #define UNUSED_IF_NO_GEOS
66 : #endif
67 :
68 : /************************************************************************/
69 : /* createFromWkb() */
70 : /************************************************************************/
71 :
72 : /**
73 : * \brief Create a geometry object of the appropriate type from its
74 : * well known binary representation.
75 : *
76 : * Note that if nBytes is passed as zero, no checking can be done on whether
77 : * the pabyData is sufficient. This can result in a crash if the input
78 : * data is corrupt. This function returns no indication of the number of
79 : * bytes from the data source actually used to represent the returned
80 : * geometry object. Use OGRGeometry::WkbSize() on the returned geometry to
81 : * establish the number of bytes it required in WKB format.
82 : *
83 : * Also note that this is a static method, and that there
84 : * is no need to instantiate an OGRGeometryFactory object.
85 : *
86 : * The C function OGR_G_CreateFromWkb() is the same as this method.
87 : *
88 : * @param pabyData pointer to the input BLOB data.
89 : * @param poSR pointer to the spatial reference to be assigned to the
90 : * created geometry object. This may be NULL.
91 : * @param ppoReturn the newly created geometry object will be assigned to the
92 : * indicated pointer on return. This will be NULL in case
93 : * of failure. If not NULL, *ppoReturn should be freed with
94 : * OGRGeometryFactory::destroyGeometry() after use.
95 : * @param nBytes the number of bytes available in pabyData, or -1 if it isn't
96 : * known
97 : * @param eWkbVariant WKB variant.
98 : *
99 : * @return OGRERR_NONE if all goes well, otherwise any of
100 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
101 : * OGRERR_CORRUPT_DATA may be returned.
102 : */
103 :
104 65558 : OGRErr OGRGeometryFactory::createFromWkb(const void *pabyData,
105 : const OGRSpatialReference *poSR,
106 : OGRGeometry **ppoReturn, size_t nBytes,
107 : OGRwkbVariant eWkbVariant)
108 :
109 : {
110 65558 : size_t nBytesConsumedOutIgnored = 0;
111 65558 : return createFromWkb(pabyData, poSR, ppoReturn, nBytes, eWkbVariant,
112 131114 : nBytesConsumedOutIgnored);
113 : }
114 :
115 : /**
116 : * \brief Create a geometry object of the appropriate type from its
117 : * well known binary representation.
118 : *
119 : * Note that if nBytes is passed as zero, no checking can be done on whether
120 : * the pabyData is sufficient. This can result in a crash if the input
121 : * data is corrupt. This function returns no indication of the number of
122 : * bytes from the data source actually used to represent the returned
123 : * geometry object. Use OGRGeometry::WkbSize() on the returned geometry to
124 : * establish the number of bytes it required in WKB format.
125 : *
126 : * Also note that this is a static method, and that there
127 : * is no need to instantiate an OGRGeometryFactory object.
128 : *
129 : * The C function OGR_G_CreateFromWkb() is the same as this method.
130 : *
131 : * @param pabyData pointer to the input BLOB data.
132 : * @param poSR pointer to the spatial reference to be assigned to the
133 : * created geometry object. This may be NULL.
134 : * @param ppoReturn the newly created geometry object will be assigned to the
135 : * indicated pointer on return. This will be NULL in case
136 : * of failure. If not NULL, *ppoReturn should be freed with
137 : * OGRGeometryFactory::destroyGeometry() after use.
138 : * @param nBytes the number of bytes available in pabyData, or -1 if it isn't
139 : * known
140 : * @param eWkbVariant WKB variant.
141 : * @param nBytesConsumedOut output parameter. Number of bytes consumed.
142 : *
143 : * @return OGRERR_NONE if all goes well, otherwise any of
144 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
145 : * OGRERR_CORRUPT_DATA may be returned.
146 : * @since GDAL 2.3
147 : */
148 :
149 103065 : OGRErr OGRGeometryFactory::createFromWkb(const void *pabyData,
150 : const OGRSpatialReference *poSR,
151 : OGRGeometry **ppoReturn, size_t nBytes,
152 : OGRwkbVariant eWkbVariant,
153 : size_t &nBytesConsumedOut)
154 :
155 : {
156 103065 : const GByte *l_pabyData = static_cast<const GByte *>(pabyData);
157 103065 : nBytesConsumedOut = 0;
158 103065 : *ppoReturn = nullptr;
159 :
160 103065 : if (nBytes < 9 && nBytes != static_cast<size_t>(-1))
161 1379 : return OGRERR_NOT_ENOUGH_DATA;
162 :
163 : /* -------------------------------------------------------------------- */
164 : /* Get the byte order byte. The extra tests are to work around */
165 : /* bug sin the WKB of DB2 v7.2 as identified by Safe Software. */
166 : /* -------------------------------------------------------------------- */
167 101686 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(*l_pabyData);
168 101686 : if (nByteOrder != wkbXDR && nByteOrder != wkbNDR)
169 : {
170 354 : CPLDebug("OGR",
171 : "OGRGeometryFactory::createFromWkb() - got corrupt data.\n"
172 : "%02X%02X%02X%02X%02X%02X%02X%02X%02X",
173 354 : l_pabyData[0], l_pabyData[1], l_pabyData[2], l_pabyData[3],
174 354 : l_pabyData[4], l_pabyData[5], l_pabyData[6], l_pabyData[7],
175 354 : l_pabyData[8]);
176 354 : return OGRERR_CORRUPT_DATA;
177 : }
178 :
179 : /* -------------------------------------------------------------------- */
180 : /* Get the geometry feature type. For now we assume that */
181 : /* geometry type is between 0 and 255 so we only have to fetch */
182 : /* one byte. */
183 : /* -------------------------------------------------------------------- */
184 :
185 101332 : OGRwkbGeometryType eGeometryType = wkbUnknown;
186 : const OGRErr err =
187 101332 : OGRReadWKBGeometryType(l_pabyData, eWkbVariant, &eGeometryType);
188 :
189 101333 : if (err != OGRERR_NONE)
190 563 : return err;
191 :
192 : /* -------------------------------------------------------------------- */
193 : /* Instantiate a geometry of the appropriate type, and */
194 : /* initialize from the input stream. */
195 : /* -------------------------------------------------------------------- */
196 100770 : OGRGeometry *poGeom = createGeometry(eGeometryType);
197 :
198 100769 : if (poGeom == nullptr)
199 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
200 :
201 : /* -------------------------------------------------------------------- */
202 : /* Import from binary. */
203 : /* -------------------------------------------------------------------- */
204 201539 : const OGRErr eErr = poGeom->importFromWkb(l_pabyData, nBytes, eWkbVariant,
205 100769 : nBytesConsumedOut);
206 100770 : if (eErr != OGRERR_NONE)
207 : {
208 7315 : delete poGeom;
209 7315 : return eErr;
210 : }
211 :
212 : /* -------------------------------------------------------------------- */
213 : /* Assign spatial reference system. */
214 : /* -------------------------------------------------------------------- */
215 :
216 97011 : if (poGeom->hasCurveGeometry() &&
217 3556 : CPLTestBool(CPLGetConfigOption("OGR_STROKE_CURVE", "FALSE")))
218 : {
219 5 : OGRGeometry *poNewGeom = poGeom->getLinearGeometry();
220 5 : delete poGeom;
221 5 : poGeom = poNewGeom;
222 : }
223 93455 : poGeom->assignSpatialReference(poSR);
224 93454 : *ppoReturn = poGeom;
225 :
226 93454 : return OGRERR_NONE;
227 : }
228 :
229 : /************************************************************************/
230 : /* OGR_G_CreateFromWkb() */
231 : /************************************************************************/
232 : /**
233 : * \brief Create a geometry object of the appropriate type from its
234 : * well known binary representation.
235 : *
236 : * Note that if nBytes is passed as zero, no checking can be done on whether
237 : * the pabyData is sufficient. This can result in a crash if the input
238 : * data is corrupt. This function returns no indication of the number of
239 : * bytes from the data source actually used to represent the returned
240 : * geometry object. Use OGR_G_WkbSize() on the returned geometry to
241 : * establish the number of bytes it required in WKB format.
242 : *
243 : * The OGRGeometryFactory::createFromWkb() CPP method is the same as this
244 : * function.
245 : *
246 : * @param pabyData pointer to the input BLOB data.
247 : * @param hSRS handle to the spatial reference to be assigned to the
248 : * created geometry object. This may be NULL.
249 : * @param phGeometry the newly created geometry object will
250 : * be assigned to the indicated handle on return. This will be NULL in case
251 : * of failure. If not NULL, *phGeometry should be freed with
252 : * OGR_G_DestroyGeometry() after use.
253 : * @param nBytes the number of bytes of data available in pabyData, or -1
254 : * if it is not known, but assumed to be sufficient.
255 : *
256 : * @return OGRERR_NONE if all goes well, otherwise any of
257 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
258 : * OGRERR_CORRUPT_DATA may be returned.
259 : */
260 :
261 2 : OGRErr CPL_DLL OGR_G_CreateFromWkb(const void *pabyData,
262 : OGRSpatialReferenceH hSRS,
263 : OGRGeometryH *phGeometry, int nBytes)
264 :
265 : {
266 2 : return OGRGeometryFactory::createFromWkb(
267 2 : pabyData, OGRSpatialReference::FromHandle(hSRS),
268 2 : reinterpret_cast<OGRGeometry **>(phGeometry), nBytes);
269 : }
270 :
271 : /************************************************************************/
272 : /* OGR_G_CreateFromWkbEx() */
273 : /************************************************************************/
274 : /**
275 : * \brief Create a geometry object of the appropriate type from its
276 : * well known binary representation.
277 : *
278 : * Note that if nBytes is passed as zero, no checking can be done on whether
279 : * the pabyData is sufficient. This can result in a crash if the input
280 : * data is corrupt. This function returns no indication of the number of
281 : * bytes from the data source actually used to represent the returned
282 : * geometry object. Use OGR_G_WkbSizeEx() on the returned geometry to
283 : * establish the number of bytes it required in WKB format.
284 : *
285 : * The OGRGeometryFactory::createFromWkb() CPP method is the same as this
286 : * function.
287 : *
288 : * @param pabyData pointer to the input BLOB data.
289 : * @param hSRS handle to the spatial reference to be assigned to the
290 : * created geometry object. This may be NULL.
291 : * @param phGeometry the newly created geometry object will
292 : * be assigned to the indicated handle on return. This will be NULL in case
293 : * of failure. If not NULL, *phGeometry should be freed with
294 : * OGR_G_DestroyGeometry() after use.
295 : * @param nBytes the number of bytes of data available in pabyData, or -1
296 : * if it is not known, but assumed to be sufficient.
297 : *
298 : * @return OGRERR_NONE if all goes well, otherwise any of
299 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
300 : * OGRERR_CORRUPT_DATA may be returned.
301 : * @since GDAL 3.3
302 : */
303 :
304 31028 : OGRErr CPL_DLL OGR_G_CreateFromWkbEx(const void *pabyData,
305 : OGRSpatialReferenceH hSRS,
306 : OGRGeometryH *phGeometry, size_t nBytes)
307 :
308 : {
309 31028 : return OGRGeometryFactory::createFromWkb(
310 31028 : pabyData, OGRSpatialReference::FromHandle(hSRS),
311 31028 : reinterpret_cast<OGRGeometry **>(phGeometry), nBytes);
312 : }
313 :
314 : /************************************************************************/
315 : /* createFromWkt() */
316 : /************************************************************************/
317 :
318 : /**
319 : * \brief Create a geometry object of the appropriate type from its
320 : * well known text representation.
321 : *
322 : * The C function OGR_G_CreateFromWkt() is the same as this method.
323 : *
324 : * @param ppszData input zero terminated string containing well known text
325 : * representation of the geometry to be created. The pointer
326 : * is updated to point just beyond that last character consumed.
327 : * @param poSR pointer to the spatial reference to be assigned to the
328 : * created geometry object. This may be NULL.
329 : * @param ppoReturn the newly created geometry object will be assigned to the
330 : * indicated pointer on return. This will be NULL if the
331 : * method fails. If not NULL, *ppoReturn should be freed with
332 : * OGRGeometryFactory::destroyGeometry() after use.
333 : *
334 : * <b>Example:</b>
335 : *
336 : * \code{.cpp}
337 : * const char* wkt= "POINT(0 0)";
338 : *
339 : * // cast because OGR_G_CreateFromWkt will move the pointer
340 : * char* pszWkt = (char*) wkt;
341 : * OGRSpatialReferenceH ref = OSRNewSpatialReference(NULL);
342 : * OGRGeometryH new_geom;
343 : * OSRSetAxisMappingStrategy(poSR, OAMS_TRADITIONAL_GIS_ORDER);
344 : * OGRErr err = OGR_G_CreateFromWkt(&pszWkt, ref, &new_geom);
345 : * \endcode
346 : *
347 : *
348 : *
349 : * @return OGRERR_NONE if all goes well, otherwise any of
350 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
351 : * OGRERR_CORRUPT_DATA may be returned.
352 : */
353 :
354 123304 : OGRErr OGRGeometryFactory::createFromWkt(const char **ppszData,
355 : const OGRSpatialReference *poSR,
356 : OGRGeometry **ppoReturn)
357 :
358 : {
359 123304 : const char *pszInput = *ppszData;
360 123304 : *ppoReturn = nullptr;
361 :
362 : /* -------------------------------------------------------------------- */
363 : /* Get the first token, which should be the geometry type. */
364 : /* -------------------------------------------------------------------- */
365 123304 : char szToken[OGR_WKT_TOKEN_MAX] = {};
366 123304 : if (OGRWktReadToken(pszInput, szToken) == nullptr)
367 0 : return OGRERR_CORRUPT_DATA;
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Instantiate a geometry of the appropriate type. */
371 : /* -------------------------------------------------------------------- */
372 123304 : OGRGeometry *poGeom = nullptr;
373 123304 : if (STARTS_WITH_CI(szToken, "POINT"))
374 : {
375 97491 : poGeom = new OGRPoint();
376 : }
377 25813 : else if (STARTS_WITH_CI(szToken, "LINESTRING"))
378 : {
379 1602 : poGeom = new OGRLineString();
380 : }
381 24211 : else if (STARTS_WITH_CI(szToken, "POLYGON"))
382 : {
383 15739 : poGeom = new OGRPolygon();
384 : }
385 8472 : else if (STARTS_WITH_CI(szToken, "TRIANGLE"))
386 : {
387 61 : poGeom = new OGRTriangle();
388 : }
389 8411 : else if (STARTS_WITH_CI(szToken, "GEOMETRYCOLLECTION"))
390 : {
391 523 : poGeom = new OGRGeometryCollection();
392 : }
393 7888 : else if (STARTS_WITH_CI(szToken, "MULTIPOLYGON"))
394 : {
395 937 : poGeom = new OGRMultiPolygon();
396 : }
397 6951 : else if (STARTS_WITH_CI(szToken, "MULTIPOINT"))
398 : {
399 578 : poGeom = new OGRMultiPoint();
400 : }
401 6373 : else if (STARTS_WITH_CI(szToken, "MULTILINESTRING"))
402 : {
403 637 : poGeom = new OGRMultiLineString();
404 : }
405 5736 : else if (STARTS_WITH_CI(szToken, "CIRCULARSTRING"))
406 : {
407 3501 : poGeom = new OGRCircularString();
408 : }
409 2235 : else if (STARTS_WITH_CI(szToken, "COMPOUNDCURVE"))
410 : {
411 289 : poGeom = new OGRCompoundCurve();
412 : }
413 1946 : else if (STARTS_WITH_CI(szToken, "CURVEPOLYGON"))
414 : {
415 318 : poGeom = new OGRCurvePolygon();
416 : }
417 1628 : else if (STARTS_WITH_CI(szToken, "MULTICURVE"))
418 : {
419 135 : poGeom = new OGRMultiCurve();
420 : }
421 1493 : else if (STARTS_WITH_CI(szToken, "MULTISURFACE"))
422 : {
423 152 : poGeom = new OGRMultiSurface();
424 : }
425 :
426 1341 : else if (STARTS_WITH_CI(szToken, "POLYHEDRALSURFACE"))
427 : {
428 63 : poGeom = new OGRPolyhedralSurface();
429 : }
430 :
431 1278 : else if (STARTS_WITH_CI(szToken, "TIN"))
432 : {
433 200 : poGeom = new OGRTriangulatedSurface();
434 : }
435 :
436 : else
437 : {
438 1078 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
439 : }
440 :
441 : /* -------------------------------------------------------------------- */
442 : /* Do the import. */
443 : /* -------------------------------------------------------------------- */
444 122226 : const OGRErr eErr = poGeom->importFromWkt(&pszInput);
445 :
446 : /* -------------------------------------------------------------------- */
447 : /* Assign spatial reference system. */
448 : /* -------------------------------------------------------------------- */
449 122226 : if (eErr == OGRERR_NONE)
450 : {
451 126356 : if (poGeom->hasCurveGeometry() &&
452 4367 : CPLTestBool(CPLGetConfigOption("OGR_STROKE_CURVE", "FALSE")))
453 : {
454 9 : OGRGeometry *poNewGeom = poGeom->getLinearGeometry();
455 9 : delete poGeom;
456 9 : poGeom = poNewGeom;
457 : }
458 121989 : poGeom->assignSpatialReference(poSR);
459 121989 : *ppoReturn = poGeom;
460 121989 : *ppszData = pszInput;
461 : }
462 : else
463 : {
464 237 : delete poGeom;
465 : }
466 :
467 122226 : return eErr;
468 : }
469 :
470 : /**
471 : * \brief Create a geometry object of the appropriate type from its
472 : * well known text representation.
473 : *
474 : * The C function OGR_G_CreateFromWkt() is the same as this method.
475 : *
476 : * @param pszData input zero terminated string containing well known text
477 : * representation of the geometry to be created.
478 : * @param poSR pointer to the spatial reference to be assigned to the
479 : * created geometry object. This may be NULL.
480 : * @param ppoReturn the newly created geometry object will be assigned to the
481 : * indicated pointer on return. This will be NULL if the
482 : * method fails. If not NULL, *ppoReturn should be freed with
483 : * OGRGeometryFactory::destroyGeometry() after use.
484 :
485 : * @return OGRERR_NONE if all goes well, otherwise any of
486 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
487 : * OGRERR_CORRUPT_DATA may be returned.
488 : * @since GDAL 2.3
489 : */
490 :
491 5652 : OGRErr OGRGeometryFactory::createFromWkt(const char *pszData,
492 : const OGRSpatialReference *poSR,
493 : OGRGeometry **ppoReturn)
494 :
495 : {
496 5652 : return createFromWkt(&pszData, poSR, ppoReturn);
497 : }
498 :
499 : /************************************************************************/
500 : /* OGR_G_CreateFromWkt() */
501 : /************************************************************************/
502 : /**
503 : * \brief Create a geometry object of the appropriate type from its well known
504 : * text representation.
505 : *
506 : * The OGRGeometryFactory::createFromWkt CPP method is the same as this
507 : * function.
508 : *
509 : * @param ppszData input zero terminated string containing well known text
510 : * representation of the geometry to be created. The pointer
511 : * is updated to point just beyond that last character consumed.
512 : * @param hSRS handle to the spatial reference to be assigned to the
513 : * created geometry object. This may be NULL.
514 : * @param phGeometry the newly created geometry object will be assigned to the
515 : * indicated handle on return. This will be NULL if the
516 : * method fails. If not NULL, *phGeometry should be freed with
517 : * OGR_G_DestroyGeometry() after use.
518 : *
519 : * @return OGRERR_NONE if all goes well, otherwise any of
520 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
521 : * OGRERR_CORRUPT_DATA may be returned.
522 : */
523 :
524 115929 : OGRErr CPL_DLL OGR_G_CreateFromWkt(char **ppszData, OGRSpatialReferenceH hSRS,
525 : OGRGeometryH *phGeometry)
526 :
527 : {
528 115929 : return OGRGeometryFactory::createFromWkt(
529 : const_cast<const char **>(ppszData),
530 115929 : OGRSpatialReference::FromHandle(hSRS),
531 115929 : reinterpret_cast<OGRGeometry **>(phGeometry));
532 : }
533 :
534 : /************************************************************************/
535 : /* createGeometry() */
536 : /************************************************************************/
537 :
538 : /**
539 : * \brief Create an empty geometry of desired type.
540 : *
541 : * This is equivalent to allocating the desired geometry with new, but
542 : * the allocation is guaranteed to take place in the context of the
543 : * GDAL/OGR heap.
544 : *
545 : * This method is the same as the C function OGR_G_CreateGeometry().
546 : *
547 : * @param eGeometryType the type code of the geometry class to be instantiated.
548 : *
549 : * @return the newly create geometry or NULL on failure. Should be freed with
550 : * OGRGeometryFactory::destroyGeometry() after use.
551 : */
552 :
553 : OGRGeometry *
554 289759 : OGRGeometryFactory::createGeometry(OGRwkbGeometryType eGeometryType)
555 :
556 : {
557 289759 : OGRGeometry *poGeom = nullptr;
558 289759 : switch (wkbFlatten(eGeometryType))
559 : {
560 192274 : case wkbPoint:
561 384548 : poGeom = new (std::nothrow) OGRPoint();
562 192274 : break;
563 :
564 23395 : case wkbLineString:
565 46790 : poGeom = new (std::nothrow) OGRLineString();
566 23395 : break;
567 :
568 30066 : case wkbPolygon:
569 60133 : poGeom = new (std::nothrow) OGRPolygon();
570 30067 : break;
571 :
572 1953 : case wkbGeometryCollection:
573 3906 : poGeom = new (std::nothrow) OGRGeometryCollection();
574 1953 : break;
575 :
576 2720 : case wkbMultiPolygon:
577 5440 : poGeom = new (std::nothrow) OGRMultiPolygon();
578 2720 : break;
579 :
580 1356 : case wkbMultiPoint:
581 2712 : poGeom = new (std::nothrow) OGRMultiPoint();
582 1356 : break;
583 :
584 1655 : case wkbMultiLineString:
585 3310 : poGeom = new (std::nothrow) OGRMultiLineString();
586 1655 : break;
587 :
588 3706 : case wkbLinearRing:
589 7412 : poGeom = new (std::nothrow) OGRLinearRing();
590 3706 : break;
591 :
592 68 : case wkbCircularString:
593 136 : poGeom = new (std::nothrow) OGRCircularString();
594 68 : break;
595 :
596 1980 : case wkbCompoundCurve:
597 3960 : poGeom = new (std::nothrow) OGRCompoundCurve();
598 1980 : break;
599 :
600 44 : case wkbCurvePolygon:
601 88 : poGeom = new (std::nothrow) OGRCurvePolygon();
602 44 : break;
603 :
604 1120 : case wkbMultiCurve:
605 2240 : poGeom = new (std::nothrow) OGRMultiCurve();
606 1120 : break;
607 :
608 1182 : case wkbMultiSurface:
609 2364 : poGeom = new (std::nothrow) OGRMultiSurface();
610 1182 : break;
611 :
612 14638 : case wkbTriangle:
613 29276 : poGeom = new (std::nothrow) OGRTriangle();
614 14638 : break;
615 :
616 7378 : case wkbPolyhedralSurface:
617 14756 : poGeom = new (std::nothrow) OGRPolyhedralSurface();
618 7378 : break;
619 :
620 6223 : case wkbTIN:
621 12446 : poGeom = new (std::nothrow) OGRTriangulatedSurface();
622 6223 : break;
623 :
624 1 : case wkbUnknown:
625 1 : break;
626 :
627 0 : default:
628 0 : CPLAssert(false);
629 : break;
630 : }
631 289760 : if (poGeom)
632 : {
633 289759 : if (OGR_GT_HasZ(eGeometryType))
634 63967 : poGeom->set3D(true);
635 289758 : if (OGR_GT_HasM(eGeometryType))
636 59799 : poGeom->setMeasured(true);
637 : }
638 289760 : return poGeom;
639 : }
640 :
641 : /************************************************************************/
642 : /* OGR_G_CreateGeometry() */
643 : /************************************************************************/
644 : /**
645 : * \brief Create an empty geometry of desired type.
646 : *
647 : * This is equivalent to allocating the desired geometry with new, but
648 : * the allocation is guaranteed to take place in the context of the
649 : * GDAL/OGR heap.
650 : *
651 : * This function is the same as the CPP method
652 : * OGRGeometryFactory::createGeometry.
653 : *
654 : * @param eGeometryType the type code of the geometry to be created.
655 : *
656 : * @return handle to the newly create geometry or NULL on failure. Should be
657 : * freed with OGR_G_DestroyGeometry() after use.
658 : */
659 :
660 185971 : OGRGeometryH OGR_G_CreateGeometry(OGRwkbGeometryType eGeometryType)
661 :
662 : {
663 185971 : return OGRGeometry::ToHandle(
664 185971 : OGRGeometryFactory::createGeometry(eGeometryType));
665 : }
666 :
667 : /************************************************************************/
668 : /* destroyGeometry() */
669 : /************************************************************************/
670 :
671 : /**
672 : * \brief Destroy geometry object.
673 : *
674 : * Equivalent to invoking delete on a geometry, but it guaranteed to take
675 : * place within the context of the GDAL/OGR heap.
676 : *
677 : * This method is the same as the C function OGR_G_DestroyGeometry().
678 : *
679 : * @param poGeom the geometry to deallocate.
680 : */
681 :
682 18 : void OGRGeometryFactory::destroyGeometry(OGRGeometry *poGeom)
683 :
684 : {
685 18 : delete poGeom;
686 18 : }
687 :
688 : /************************************************************************/
689 : /* OGR_G_DestroyGeometry() */
690 : /************************************************************************/
691 : /**
692 : * \brief Destroy geometry object.
693 : *
694 : * Equivalent to invoking delete on a geometry, but it guaranteed to take
695 : * place within the context of the GDAL/OGR heap.
696 : *
697 : * This function is the same as the CPP method
698 : * OGRGeometryFactory::destroyGeometry.
699 : *
700 : * @param hGeom handle to the geometry to delete.
701 : */
702 :
703 290562 : void OGR_G_DestroyGeometry(OGRGeometryH hGeom)
704 :
705 : {
706 290562 : delete OGRGeometry::FromHandle(hGeom);
707 290562 : }
708 :
709 : /************************************************************************/
710 : /* forceToPolygon() */
711 : /************************************************************************/
712 :
713 : /**
714 : * \brief Convert to polygon.
715 : *
716 : * Tries to force the provided geometry to be a polygon. This effects a change
717 : * on multipolygons.
718 : * Starting with GDAL 2.0, curve polygons or closed curves will be changed to
719 : * polygons. The passed in geometry is consumed and a new one returned (or
720 : * potentially the same one).
721 : *
722 : * Note: the resulting polygon may break the Simple Features rules for polygons,
723 : * for example when converting from a multi-part multipolygon.
724 : *
725 : * @param poGeom the input geometry - ownership is passed to the method.
726 : * @return new geometry.
727 : */
728 :
729 148 : OGRGeometry *OGRGeometryFactory::forceToPolygon(OGRGeometry *poGeom)
730 :
731 : {
732 148 : if (poGeom == nullptr)
733 0 : return nullptr;
734 :
735 148 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
736 :
737 148 : if (eGeomType == wkbCurvePolygon)
738 : {
739 34 : OGRCurvePolygon *poCurve = poGeom->toCurvePolygon();
740 :
741 34 : if (!poGeom->hasCurveGeometry(TRUE))
742 14 : return OGRSurface::CastToPolygon(poCurve);
743 :
744 20 : OGRPolygon *poPoly = poCurve->CurvePolyToPoly();
745 20 : delete poGeom;
746 20 : return poPoly;
747 : }
748 :
749 : // base polygon or triangle
750 114 : if (OGR_GT_IsSubClassOf(eGeomType, wkbPolygon))
751 : {
752 7 : return OGRSurface::CastToPolygon(poGeom->toSurface());
753 : }
754 :
755 107 : if (OGR_GT_IsCurve(eGeomType))
756 : {
757 60 : OGRCurve *poCurve = poGeom->toCurve();
758 60 : if (poCurve->getNumPoints() >= 3 && poCurve->get_IsClosed())
759 : {
760 40 : OGRPolygon *poPolygon = new OGRPolygon();
761 40 : poPolygon->assignSpatialReference(poGeom->getSpatialReference());
762 :
763 40 : if (!poGeom->hasCurveGeometry(TRUE))
764 : {
765 26 : poPolygon->addRingDirectly(OGRCurve::CastToLinearRing(poCurve));
766 : }
767 : else
768 : {
769 14 : OGRLineString *poLS = poCurve->CurveToLine();
770 14 : poPolygon->addRingDirectly(OGRCurve::CastToLinearRing(poLS));
771 14 : delete poGeom;
772 : }
773 40 : return poPolygon;
774 : }
775 : }
776 :
777 67 : if (OGR_GT_IsSubClassOf(eGeomType, wkbPolyhedralSurface))
778 : {
779 6 : OGRPolyhedralSurface *poPS = poGeom->toPolyhedralSurface();
780 6 : if (poPS->getNumGeometries() == 1)
781 : {
782 5 : poGeom = OGRSurface::CastToPolygon(
783 5 : poPS->getGeometryRef(0)->clone()->toSurface());
784 5 : delete poPS;
785 5 : return poGeom;
786 : }
787 : }
788 :
789 62 : if (eGeomType != wkbGeometryCollection && eGeomType != wkbMultiPolygon &&
790 : eGeomType != wkbMultiSurface)
791 38 : return poGeom;
792 :
793 : // Build an aggregated polygon from all the polygon rings in the container.
794 24 : OGRPolygon *poPolygon = new OGRPolygon();
795 24 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
796 24 : if (poGeom->hasCurveGeometry())
797 : {
798 : OGRGeometryCollection *poNewGC =
799 5 : poGC->getLinearGeometry()->toGeometryCollection();
800 5 : delete poGC;
801 5 : poGeom = poNewGC;
802 5 : poGC = poNewGC;
803 : }
804 :
805 24 : poPolygon->assignSpatialReference(poGeom->getSpatialReference());
806 :
807 53 : for (int iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++)
808 : {
809 29 : if (wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType()) !=
810 : wkbPolygon)
811 12 : continue;
812 :
813 17 : OGRPolygon *poOldPoly = poGC->getGeometryRef(iGeom)->toPolygon();
814 :
815 17 : if (poOldPoly->getExteriorRing() == nullptr)
816 3 : continue;
817 :
818 14 : poPolygon->addRingDirectly(poOldPoly->stealExteriorRing());
819 :
820 22 : for (int iRing = 0; iRing < poOldPoly->getNumInteriorRings(); iRing++)
821 8 : poPolygon->addRingDirectly(poOldPoly->stealInteriorRing(iRing));
822 : }
823 :
824 24 : delete poGC;
825 :
826 24 : return poPolygon;
827 : }
828 :
829 : /************************************************************************/
830 : /* OGR_G_ForceToPolygon() */
831 : /************************************************************************/
832 :
833 : /**
834 : * \brief Convert to polygon.
835 : *
836 : * This function is the same as the C++ method
837 : * OGRGeometryFactory::forceToPolygon().
838 : *
839 : * @param hGeom handle to the geometry to convert (ownership surrendered).
840 : * @return the converted geometry (ownership to caller).
841 : *
842 : * @since GDAL/OGR 1.8.0
843 : */
844 :
845 46 : OGRGeometryH OGR_G_ForceToPolygon(OGRGeometryH hGeom)
846 :
847 : {
848 46 : return OGRGeometry::ToHandle(
849 46 : OGRGeometryFactory::forceToPolygon(OGRGeometry::FromHandle(hGeom)));
850 : }
851 :
852 : /************************************************************************/
853 : /* forceToMultiPolygon() */
854 : /************************************************************************/
855 :
856 : /**
857 : * \brief Convert to multipolygon.
858 : *
859 : * Tries to force the provided geometry to be a multipolygon. Currently
860 : * this just effects a change on polygons. The passed in geometry is
861 : * consumed and a new one returned (or potentially the same one).
862 : *
863 : * @return new geometry.
864 : */
865 :
866 2941 : OGRGeometry *OGRGeometryFactory::forceToMultiPolygon(OGRGeometry *poGeom)
867 :
868 : {
869 2941 : if (poGeom == nullptr)
870 0 : return nullptr;
871 :
872 2941 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
873 :
874 : /* -------------------------------------------------------------------- */
875 : /* If this is already a MultiPolygon, nothing to do */
876 : /* -------------------------------------------------------------------- */
877 2941 : if (eGeomType == wkbMultiPolygon)
878 : {
879 60 : return poGeom;
880 : }
881 :
882 : /* -------------------------------------------------------------------- */
883 : /* If this is already a MultiSurface with compatible content, */
884 : /* just cast */
885 : /* -------------------------------------------------------------------- */
886 2881 : if (eGeomType == wkbMultiSurface)
887 : {
888 9 : OGRMultiSurface *poMS = poGeom->toMultiSurface();
889 9 : if (!poMS->hasCurveGeometry(TRUE))
890 : {
891 4 : return OGRMultiSurface::CastToMultiPolygon(poMS);
892 : }
893 : }
894 :
895 : /* -------------------------------------------------------------------- */
896 : /* Check for the case of a geometrycollection that can be */
897 : /* promoted to MultiPolygon. */
898 : /* -------------------------------------------------------------------- */
899 2877 : if (eGeomType == wkbGeometryCollection || eGeomType == wkbMultiSurface)
900 : {
901 91 : bool bAllPoly = true;
902 91 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
903 91 : if (poGeom->hasCurveGeometry())
904 : {
905 : OGRGeometryCollection *poNewGC =
906 6 : poGC->getLinearGeometry()->toGeometryCollection();
907 6 : delete poGC;
908 6 : poGeom = poNewGC;
909 6 : poGC = poNewGC;
910 : }
911 :
912 91 : bool bCanConvertToMultiPoly = true;
913 372 : for (int iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++)
914 : {
915 : OGRwkbGeometryType eSubGeomType =
916 281 : wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType());
917 281 : if (eSubGeomType != wkbPolygon)
918 252 : bAllPoly = false;
919 281 : if (eSubGeomType != wkbMultiPolygon && eSubGeomType != wkbPolygon &&
920 194 : eSubGeomType != wkbPolyhedralSurface && eSubGeomType != wkbTIN)
921 : {
922 16 : bCanConvertToMultiPoly = false;
923 : }
924 : }
925 :
926 91 : if (!bCanConvertToMultiPoly)
927 12 : return poGeom;
928 :
929 79 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
930 79 : poMP->assignSpatialReference(poGeom->getSpatialReference());
931 :
932 342 : while (poGC->getNumGeometries() > 0)
933 : {
934 263 : OGRGeometry *poSubGeom = poGC->getGeometryRef(0);
935 263 : poGC->removeGeometry(0, FALSE);
936 263 : if (bAllPoly)
937 : {
938 27 : poMP->addGeometryDirectly(poSubGeom);
939 : }
940 : else
941 : {
942 236 : poSubGeom = forceToMultiPolygon(poSubGeom);
943 236 : OGRMultiPolygon *poSubMP = poSubGeom->toMultiPolygon();
944 614 : while (poSubMP != nullptr && poSubMP->getNumGeometries() > 0)
945 : {
946 378 : poMP->addGeometryDirectly(poSubMP->getGeometryRef(0));
947 378 : poSubMP->removeGeometry(0, FALSE);
948 : }
949 236 : delete poSubMP;
950 : }
951 : }
952 :
953 79 : delete poGC;
954 :
955 79 : return poMP;
956 : }
957 :
958 2786 : if (eGeomType == wkbCurvePolygon)
959 : {
960 5 : OGRPolygon *poPoly = poGeom->toCurvePolygon()->CurvePolyToPoly();
961 5 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
962 5 : poMP->assignSpatialReference(poGeom->getSpatialReference());
963 5 : poMP->addGeometryDirectly(poPoly);
964 5 : delete poGeom;
965 5 : return poMP;
966 : }
967 :
968 : /* -------------------------------------------------------------------- */
969 : /* If it is PolyhedralSurface or TIN, then pretend it is a */
970 : /* multipolygon. */
971 : /* -------------------------------------------------------------------- */
972 2781 : if (OGR_GT_IsSubClassOf(eGeomType, wkbPolyhedralSurface))
973 : {
974 274 : return OGRPolyhedralSurface::CastToMultiPolygon(
975 274 : poGeom->toPolyhedralSurface());
976 : }
977 :
978 2507 : if (eGeomType == wkbTriangle)
979 : {
980 2 : return forceToMultiPolygon(forceToPolygon(poGeom));
981 : }
982 :
983 : /* -------------------------------------------------------------------- */
984 : /* Eventually we should try to split the polygon into component */
985 : /* island polygons. But that is a lot of work and can be put off. */
986 : /* -------------------------------------------------------------------- */
987 2505 : if (eGeomType != wkbPolygon)
988 30 : return poGeom;
989 :
990 2475 : OGRMultiPolygon *poMP = new OGRMultiPolygon();
991 2475 : poMP->assignSpatialReference(poGeom->getSpatialReference());
992 2475 : poMP->addGeometryDirectly(poGeom);
993 :
994 2475 : return poMP;
995 : }
996 :
997 : /************************************************************************/
998 : /* OGR_G_ForceToMultiPolygon() */
999 : /************************************************************************/
1000 :
1001 : /**
1002 : * \brief Convert to multipolygon.
1003 : *
1004 : * This function is the same as the C++ method
1005 : * OGRGeometryFactory::forceToMultiPolygon().
1006 : *
1007 : * @param hGeom handle to the geometry to convert (ownership surrendered).
1008 : * @return the converted geometry (ownership to caller).
1009 : *
1010 : * @since GDAL/OGR 1.8.0
1011 : */
1012 :
1013 47 : OGRGeometryH OGR_G_ForceToMultiPolygon(OGRGeometryH hGeom)
1014 :
1015 : {
1016 47 : return OGRGeometry::ToHandle(OGRGeometryFactory::forceToMultiPolygon(
1017 47 : OGRGeometry::FromHandle(hGeom)));
1018 : }
1019 :
1020 : /************************************************************************/
1021 : /* forceToMultiPoint() */
1022 : /************************************************************************/
1023 :
1024 : /**
1025 : * \brief Convert to multipoint.
1026 : *
1027 : * Tries to force the provided geometry to be a multipoint. Currently
1028 : * this just effects a change on points or collection of points.
1029 : * The passed in geometry is
1030 : * consumed and a new one returned (or potentially the same one).
1031 : *
1032 : * @return new geometry.
1033 : */
1034 :
1035 63 : OGRGeometry *OGRGeometryFactory::forceToMultiPoint(OGRGeometry *poGeom)
1036 :
1037 : {
1038 63 : if (poGeom == nullptr)
1039 0 : return nullptr;
1040 :
1041 63 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
1042 :
1043 : /* -------------------------------------------------------------------- */
1044 : /* If this is already a MultiPoint, nothing to do */
1045 : /* -------------------------------------------------------------------- */
1046 63 : if (eGeomType == wkbMultiPoint)
1047 : {
1048 2 : return poGeom;
1049 : }
1050 :
1051 : /* -------------------------------------------------------------------- */
1052 : /* Check for the case of a geometrycollection that can be */
1053 : /* promoted to MultiPoint. */
1054 : /* -------------------------------------------------------------------- */
1055 61 : if (eGeomType == wkbGeometryCollection)
1056 : {
1057 14 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
1058 18 : for (const auto &poMember : poGC)
1059 : {
1060 14 : if (wkbFlatten(poMember->getGeometryType()) != wkbPoint)
1061 10 : return poGeom;
1062 : }
1063 :
1064 4 : OGRMultiPoint *poMP = new OGRMultiPoint();
1065 4 : poMP->assignSpatialReference(poGeom->getSpatialReference());
1066 :
1067 8 : while (poGC->getNumGeometries() > 0)
1068 : {
1069 4 : poMP->addGeometryDirectly(poGC->getGeometryRef(0));
1070 4 : poGC->removeGeometry(0, FALSE);
1071 : }
1072 :
1073 4 : delete poGC;
1074 :
1075 4 : return poMP;
1076 : }
1077 :
1078 47 : if (eGeomType != wkbPoint)
1079 44 : return poGeom;
1080 :
1081 3 : OGRMultiPoint *poMP = new OGRMultiPoint();
1082 3 : poMP->assignSpatialReference(poGeom->getSpatialReference());
1083 3 : poMP->addGeometryDirectly(poGeom);
1084 :
1085 3 : return poMP;
1086 : }
1087 :
1088 : /************************************************************************/
1089 : /* OGR_G_ForceToMultiPoint() */
1090 : /************************************************************************/
1091 :
1092 : /**
1093 : * \brief Convert to multipoint.
1094 : *
1095 : * This function is the same as the C++ method
1096 : * OGRGeometryFactory::forceToMultiPoint().
1097 : *
1098 : * @param hGeom handle to the geometry to convert (ownership surrendered).
1099 : * @return the converted geometry (ownership to caller).
1100 : *
1101 : * @since GDAL/OGR 1.8.0
1102 : */
1103 :
1104 41 : OGRGeometryH OGR_G_ForceToMultiPoint(OGRGeometryH hGeom)
1105 :
1106 : {
1107 41 : return OGRGeometry::ToHandle(
1108 41 : OGRGeometryFactory::forceToMultiPoint(OGRGeometry::FromHandle(hGeom)));
1109 : }
1110 :
1111 : /************************************************************************/
1112 : /* forceToMultiLinestring() */
1113 : /************************************************************************/
1114 :
1115 : /**
1116 : * \brief Convert to multilinestring.
1117 : *
1118 : * Tries to force the provided geometry to be a multilinestring.
1119 : *
1120 : * - linestrings are placed in a multilinestring.
1121 : * - circularstrings and compoundcurves will be approximated and placed in a
1122 : * multilinestring.
1123 : * - geometry collections will be converted to multilinestring if they only
1124 : * contain linestrings.
1125 : * - polygons will be changed to a collection of linestrings (one per ring).
1126 : * - curvepolygons will be approximated and changed to a collection of
1127 : ( linestrings (one per ring).
1128 : *
1129 : * The passed in geometry is
1130 : * consumed and a new one returned (or potentially the same one).
1131 : *
1132 : * @return new geometry.
1133 : */
1134 :
1135 2130 : OGRGeometry *OGRGeometryFactory::forceToMultiLineString(OGRGeometry *poGeom)
1136 :
1137 : {
1138 2130 : if (poGeom == nullptr)
1139 0 : return nullptr;
1140 :
1141 2130 : OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
1142 :
1143 : /* -------------------------------------------------------------------- */
1144 : /* If this is already a MultiLineString, nothing to do */
1145 : /* -------------------------------------------------------------------- */
1146 2130 : if (eGeomType == wkbMultiLineString)
1147 : {
1148 2 : return poGeom;
1149 : }
1150 :
1151 : /* -------------------------------------------------------------------- */
1152 : /* Check for the case of a geometrycollection that can be */
1153 : /* promoted to MultiLineString. */
1154 : /* -------------------------------------------------------------------- */
1155 2128 : if (eGeomType == wkbGeometryCollection)
1156 : {
1157 16 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
1158 16 : if (poGeom->hasCurveGeometry())
1159 : {
1160 : OGRGeometryCollection *poNewGC =
1161 1 : poGC->getLinearGeometry()->toGeometryCollection();
1162 1 : delete poGC;
1163 1 : poGeom = poNewGC;
1164 1 : poGC = poNewGC;
1165 : }
1166 :
1167 24 : for (auto &&poMember : poGC)
1168 : {
1169 18 : if (wkbFlatten(poMember->getGeometryType()) != wkbLineString)
1170 : {
1171 10 : return poGeom;
1172 : }
1173 : }
1174 :
1175 6 : OGRMultiLineString *poMP = new OGRMultiLineString();
1176 6 : poMP->assignSpatialReference(poGeom->getSpatialReference());
1177 :
1178 14 : while (poGC->getNumGeometries() > 0)
1179 : {
1180 8 : poMP->addGeometryDirectly(poGC->getGeometryRef(0));
1181 8 : poGC->removeGeometry(0, FALSE);
1182 : }
1183 :
1184 6 : delete poGC;
1185 :
1186 6 : return poMP;
1187 : }
1188 :
1189 : /* -------------------------------------------------------------------- */
1190 : /* Turn a linestring into a multilinestring. */
1191 : /* -------------------------------------------------------------------- */
1192 2112 : if (eGeomType == wkbLineString)
1193 : {
1194 2029 : OGRMultiLineString *poMP = new OGRMultiLineString();
1195 2029 : poMP->assignSpatialReference(poGeom->getSpatialReference());
1196 2029 : poMP->addGeometryDirectly(poGeom);
1197 2029 : return poMP;
1198 : }
1199 :
1200 : /* -------------------------------------------------------------------- */
1201 : /* Convert polygons into a multilinestring. */
1202 : /* -------------------------------------------------------------------- */
1203 83 : if (OGR_GT_IsSubClassOf(eGeomType, wkbCurvePolygon))
1204 : {
1205 25 : OGRMultiLineString *poMLS = new OGRMultiLineString();
1206 25 : poMLS->assignSpatialReference(poGeom->getSpatialReference());
1207 :
1208 51 : const auto AddRingFromSrcPoly = [poMLS](const OGRPolygon *poPoly)
1209 : {
1210 55 : for (int iRing = 0; iRing < poPoly->getNumInteriorRings() + 1;
1211 : iRing++)
1212 : {
1213 : const OGRLineString *poLR;
1214 :
1215 32 : if (iRing == 0)
1216 : {
1217 25 : poLR = poPoly->getExteriorRing();
1218 25 : if (poLR == nullptr)
1219 2 : break;
1220 : }
1221 : else
1222 7 : poLR = poPoly->getInteriorRing(iRing - 1);
1223 :
1224 30 : if (poLR == nullptr || poLR->getNumPoints() == 0)
1225 4 : continue;
1226 :
1227 26 : auto poNewLS = new OGRLineString();
1228 26 : poNewLS->addSubLineString(poLR);
1229 26 : poMLS->addGeometryDirectly(poNewLS);
1230 : }
1231 25 : };
1232 :
1233 25 : if (OGR_GT_IsSubClassOf(eGeomType, wkbPolygon))
1234 : {
1235 21 : AddRingFromSrcPoly(poGeom->toPolygon());
1236 : }
1237 : else
1238 : {
1239 : auto poTmpPoly = std::unique_ptr<OGRPolygon>(
1240 8 : poGeom->toCurvePolygon()->CurvePolyToPoly());
1241 4 : AddRingFromSrcPoly(poTmpPoly.get());
1242 : }
1243 :
1244 25 : delete poGeom;
1245 :
1246 25 : return poMLS;
1247 : }
1248 :
1249 : /* -------------------------------------------------------------------- */
1250 : /* If it is PolyhedralSurface or TIN, then pretend it is a */
1251 : /* multipolygon. */
1252 : /* -------------------------------------------------------------------- */
1253 58 : if (OGR_GT_IsSubClassOf(eGeomType, wkbPolyhedralSurface))
1254 : {
1255 0 : poGeom = CPLAssertNotNull(forceToMultiPolygon(poGeom));
1256 0 : eGeomType = wkbMultiPolygon;
1257 : }
1258 :
1259 : /* -------------------------------------------------------------------- */
1260 : /* Convert multi-polygons into a multilinestring. */
1261 : /* -------------------------------------------------------------------- */
1262 58 : if (eGeomType == wkbMultiPolygon || eGeomType == wkbMultiSurface)
1263 : {
1264 9 : OGRMultiLineString *poMLS = new OGRMultiLineString();
1265 9 : poMLS->assignSpatialReference(poGeom->getSpatialReference());
1266 :
1267 22 : const auto AddRingFromSrcMP = [poMLS](const OGRMultiPolygon *poSrcMP)
1268 : {
1269 21 : for (auto &&poPoly : poSrcMP)
1270 : {
1271 27 : for (auto &&poLR : poPoly)
1272 : {
1273 15 : if (poLR->IsEmpty())
1274 2 : continue;
1275 :
1276 13 : OGRLineString *poNewLS = new OGRLineString();
1277 13 : poNewLS->addSubLineString(poLR);
1278 13 : poMLS->addGeometryDirectly(poNewLS);
1279 : }
1280 : }
1281 9 : };
1282 :
1283 9 : if (eGeomType == wkbMultiPolygon)
1284 : {
1285 6 : AddRingFromSrcMP(poGeom->toMultiPolygon());
1286 : }
1287 : else
1288 : {
1289 : auto poTmpMPoly = std::unique_ptr<OGRMultiPolygon>(
1290 6 : poGeom->getLinearGeometry()->toMultiPolygon());
1291 3 : AddRingFromSrcMP(poTmpMPoly.get());
1292 : }
1293 :
1294 9 : delete poGeom;
1295 9 : return poMLS;
1296 : }
1297 :
1298 : /* -------------------------------------------------------------------- */
1299 : /* If it is a curve line, approximate it and wrap in a multilinestring
1300 : */
1301 : /* -------------------------------------------------------------------- */
1302 49 : if (eGeomType == wkbCircularString || eGeomType == wkbCompoundCurve)
1303 : {
1304 20 : OGRMultiLineString *poMP = new OGRMultiLineString();
1305 20 : poMP->assignSpatialReference(poGeom->getSpatialReference());
1306 20 : poMP->addGeometryDirectly(poGeom->toCurve()->CurveToLine());
1307 20 : delete poGeom;
1308 20 : return poMP;
1309 : }
1310 :
1311 : /* -------------------------------------------------------------------- */
1312 : /* If this is already a MultiCurve with compatible content, */
1313 : /* just cast */
1314 : /* -------------------------------------------------------------------- */
1315 38 : if (eGeomType == wkbMultiCurve &&
1316 9 : !poGeom->toMultiCurve()->hasCurveGeometry(TRUE))
1317 : {
1318 3 : return OGRMultiCurve::CastToMultiLineString(poGeom->toMultiCurve());
1319 : }
1320 :
1321 : /* -------------------------------------------------------------------- */
1322 : /* If it is a multicurve, call getLinearGeometry() */
1323 : /* -------------------------------------------------------------------- */
1324 26 : if (eGeomType == wkbMultiCurve)
1325 : {
1326 6 : OGRGeometry *poNewGeom = poGeom->getLinearGeometry();
1327 6 : CPLAssert(wkbFlatten(poNewGeom->getGeometryType()) ==
1328 : wkbMultiLineString);
1329 6 : delete poGeom;
1330 6 : return poNewGeom->toMultiLineString();
1331 : }
1332 :
1333 20 : return poGeom;
1334 : }
1335 :
1336 : /************************************************************************/
1337 : /* OGR_G_ForceToMultiLineString() */
1338 : /************************************************************************/
1339 :
1340 : /**
1341 : * \brief Convert to multilinestring.
1342 : *
1343 : * This function is the same as the C++ method
1344 : * OGRGeometryFactory::forceToMultiLineString().
1345 : *
1346 : * @param hGeom handle to the geometry to convert (ownership surrendered).
1347 : * @return the converted geometry (ownership to caller).
1348 : *
1349 : * @since GDAL/OGR 1.8.0
1350 : */
1351 :
1352 50 : OGRGeometryH OGR_G_ForceToMultiLineString(OGRGeometryH hGeom)
1353 :
1354 : {
1355 50 : return OGRGeometry::ToHandle(OGRGeometryFactory::forceToMultiLineString(
1356 50 : OGRGeometry::FromHandle(hGeom)));
1357 : }
1358 :
1359 : /************************************************************************/
1360 : /* removeLowerDimensionSubGeoms() */
1361 : /************************************************************************/
1362 :
1363 : /** \brief Remove sub-geometries from a geometry collection that do not have
1364 : * the maximum topological dimensionality of the collection.
1365 : *
1366 : * This is typically to be used as a cleanup phase after running
1367 : * OGRGeometry::MakeValid()
1368 : *
1369 : * For example, MakeValid() on a polygon can return a geometry collection of
1370 : * polygons and linestrings. Calling this method will return either a polygon
1371 : * or multipolygon by dropping those linestrings.
1372 : *
1373 : * On a non-geometry collection, this will return a clone of the passed
1374 : * geometry.
1375 : *
1376 : * @param poGeom input geometry
1377 : * @return a new geometry.
1378 : *
1379 : * @since GDAL 3.1.0
1380 : */
1381 : OGRGeometry *
1382 20 : OGRGeometryFactory::removeLowerDimensionSubGeoms(const OGRGeometry *poGeom)
1383 : {
1384 20 : if (poGeom == nullptr)
1385 0 : return nullptr;
1386 34 : if (wkbFlatten(poGeom->getGeometryType()) != wkbGeometryCollection ||
1387 14 : poGeom->IsEmpty())
1388 : {
1389 7 : return poGeom->clone();
1390 : }
1391 13 : const OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
1392 13 : int nMaxDim = 0;
1393 13 : OGRBoolean bHasCurve = FALSE;
1394 36 : for (const auto poSubGeom : *poGC)
1395 : {
1396 23 : nMaxDim = std::max(nMaxDim, poSubGeom->getDimension());
1397 23 : bHasCurve |= poSubGeom->hasCurveGeometry();
1398 : }
1399 13 : int nCountAtMaxDim = 0;
1400 13 : const OGRGeometry *poGeomAtMaxDim = nullptr;
1401 36 : for (const auto poSubGeom : *poGC)
1402 : {
1403 23 : if (poSubGeom->getDimension() == nMaxDim)
1404 : {
1405 18 : poGeomAtMaxDim = poSubGeom;
1406 18 : nCountAtMaxDim++;
1407 : }
1408 : }
1409 13 : if (nCountAtMaxDim == 1 && poGeomAtMaxDim != nullptr)
1410 : {
1411 8 : return poGeomAtMaxDim->clone();
1412 : }
1413 : OGRGeometryCollection *poRet =
1414 5 : (nMaxDim == 0)
1415 10 : ? static_cast<OGRGeometryCollection *>(new OGRMultiPoint())
1416 5 : : (nMaxDim == 1)
1417 10 : ? (!bHasCurve
1418 4 : ? static_cast<OGRGeometryCollection *>(
1419 1 : new OGRMultiLineString())
1420 1 : : static_cast<OGRGeometryCollection *>(new OGRMultiCurve()))
1421 3 : : (nMaxDim == 2 && !bHasCurve)
1422 6 : ? static_cast<OGRGeometryCollection *>(new OGRMultiPolygon())
1423 1 : : static_cast<OGRGeometryCollection *>(new OGRMultiSurface());
1424 15 : for (const auto poSubGeom : *poGC)
1425 : {
1426 10 : if (poSubGeom->getDimension() == nMaxDim)
1427 : {
1428 10 : if (OGR_GT_IsSubClassOf(poSubGeom->getGeometryType(),
1429 10 : wkbGeometryCollection))
1430 : {
1431 : const OGRGeometryCollection *poSubGeomAsGC =
1432 1 : poSubGeom->toGeometryCollection();
1433 2 : for (const auto poSubSubGeom : *poSubGeomAsGC)
1434 : {
1435 1 : if (poSubSubGeom->getDimension() == nMaxDim)
1436 : {
1437 1 : poRet->addGeometryDirectly(poSubSubGeom->clone());
1438 : }
1439 : }
1440 : }
1441 : else
1442 : {
1443 9 : poRet->addGeometryDirectly(poSubGeom->clone());
1444 : }
1445 : }
1446 : }
1447 5 : return poRet;
1448 : }
1449 :
1450 : /************************************************************************/
1451 : /* OGR_G_RemoveLowerDimensionSubGeoms() */
1452 : /************************************************************************/
1453 :
1454 : /** \brief Remove sub-geometries from a geometry collection that do not have
1455 : * the maximum topological dimensionality of the collection.
1456 : *
1457 : * This function is the same as the C++ method
1458 : * OGRGeometryFactory::removeLowerDimensionSubGeoms().
1459 : *
1460 : * @param hGeom handle to the geometry to convert
1461 : * @return a new geometry.
1462 : *
1463 : * @since GDAL 3.1.0
1464 : */
1465 :
1466 18 : OGRGeometryH OGR_G_RemoveLowerDimensionSubGeoms(const OGRGeometryH hGeom)
1467 :
1468 : {
1469 18 : return OGRGeometry::ToHandle(
1470 : OGRGeometryFactory::removeLowerDimensionSubGeoms(
1471 36 : OGRGeometry::FromHandle(hGeom)));
1472 : }
1473 :
1474 : /************************************************************************/
1475 : /* organizePolygons() */
1476 : /************************************************************************/
1477 :
1478 4025 : struct sPolyExtended
1479 : {
1480 : CPL_DISALLOW_COPY_ASSIGN(sPolyExtended)
1481 2264 : sPolyExtended() = default;
1482 4971 : sPolyExtended(sPolyExtended &&) = default;
1483 : sPolyExtended &operator=(sPolyExtended &&) = default;
1484 :
1485 : OGRGeometry *poGeometry = nullptr;
1486 : OGRCurvePolygon *poPolygon = nullptr;
1487 : OGREnvelope sEnvelope{};
1488 : OGRCurve *poExteriorRing = nullptr;
1489 : OGRPoint poAPoint{};
1490 : int nInitialIndex = 0;
1491 : OGRCurvePolygon *poEnclosingPolygon = nullptr;
1492 : double dfArea = 0.0;
1493 : bool bIsTopLevel = false;
1494 : bool bIsCW = false;
1495 : bool bIsPolygon = false;
1496 : };
1497 :
1498 3313 : static bool OGRGeometryFactoryCompareArea(const sPolyExtended &sPoly1,
1499 : const sPolyExtended &sPoly2)
1500 : {
1501 3313 : return sPoly2.dfArea < sPoly1.dfArea;
1502 : }
1503 :
1504 3327 : static bool OGRGeometryFactoryCompareByIndex(const sPolyExtended &sPoly1,
1505 : const sPolyExtended &sPoly2)
1506 : {
1507 3327 : return sPoly1.nInitialIndex < sPoly2.nInitialIndex;
1508 : }
1509 :
1510 : constexpr int N_CRITICAL_PART_NUMBER = 100;
1511 :
1512 : enum OrganizePolygonMethod
1513 : {
1514 : METHOD_NORMAL,
1515 : METHOD_SKIP,
1516 : METHOD_ONLY_CCW,
1517 : METHOD_CCW_INNER_JUST_AFTER_CW_OUTER
1518 : };
1519 :
1520 : /**
1521 : * \brief Organize polygons based on geometries.
1522 : *
1523 : * Analyse a set of rings (passed as simple polygons), and based on a
1524 : * geometric analysis convert them into a polygon with inner rings,
1525 : * (or a MultiPolygon if dealing with more than one polygon) that follow the
1526 : * OGC Simple Feature specification.
1527 : *
1528 : * All the input geometries must be OGRPolygon/OGRCurvePolygon with only a valid
1529 : * exterior ring (at least 4 points) and no interior rings.
1530 : *
1531 : * The passed in geometries become the responsibility of the method, but the
1532 : * papoPolygons "pointer array" remains owned by the caller.
1533 : *
1534 : * For faster computation, a polygon is considered to be inside
1535 : * another one if a single point of its external ring is included into the other
1536 : * one. (unless 'OGR_DEBUG_ORGANIZE_POLYGONS' configuration option is set to
1537 : * TRUE. In that case, a slower algorithm that tests exact topological
1538 : * relationships is used if GEOS is available.)
1539 : *
1540 : * In cases where a big number of polygons is passed to this function, the
1541 : * default processing may be really slow. You can skip the processing by adding
1542 : * METHOD=SKIP to the option list (the result of the function will be a
1543 : * multi-polygon with all polygons as toplevel polygons) or only make it analyze
1544 : * counterclockwise polygons by adding METHOD=ONLY_CCW to the option list if you
1545 : * can assume that the outline of holes is counterclockwise defined (this is the
1546 : * convention for example in shapefiles, Personal Geodatabases or File
1547 : * Geodatabases).
1548 : *
1549 : * For FileGDB, in most cases, but not always, a faster method than ONLY_CCW can
1550 : * be used. It is CCW_INNER_JUST_AFTER_CW_OUTER. When using it, inner rings are
1551 : * assumed to be counterclockwise oriented, and following immediately the outer
1552 : * ring (clockwise oriented) that they belong to. If that assumption is not met,
1553 : * an inner ring could be attached to the wrong outer ring, so this method must
1554 : * be used with care.
1555 : *
1556 : * If the OGR_ORGANIZE_POLYGONS configuration option is defined, its value will
1557 : * override the value of the METHOD option of papszOptions (useful to modify the
1558 : * behavior of the shapefile driver)
1559 : *
1560 : * @param papoPolygons array of geometry pointers - should all be OGRPolygons
1561 : * or OGRCurvePolygons. Ownership of the geometries is passed, but not of the
1562 : * array itself.
1563 : * @param nPolygonCount number of items in papoPolygons
1564 : * @param pbIsValidGeometry value may be set to FALSE if an invalid result is
1565 : * detected. Validity checks vary according to the method used and are are limited
1566 : * to what is needed to link inner rings to outer rings, so a result of TRUE
1567 : * does not mean that OGRGeometry::IsValid() returns TRUE.
1568 : * @param papszOptions a list of strings for passing options
1569 : *
1570 : * @return a single resulting geometry (either OGRPolygon, OGRCurvePolygon,
1571 : * OGRMultiPolygon, OGRMultiSurface or OGRGeometryCollection). Returns a
1572 : * POLYGON EMPTY in the case of nPolygonCount being 0.
1573 : */
1574 :
1575 811 : OGRGeometry *OGRGeometryFactory::organizePolygons(OGRGeometry **papoPolygons,
1576 : int nPolygonCount,
1577 : int *pbIsValidGeometry,
1578 : const char **papszOptions)
1579 : {
1580 811 : if (nPolygonCount == 0)
1581 : {
1582 4 : if (pbIsValidGeometry)
1583 0 : *pbIsValidGeometry = TRUE;
1584 :
1585 4 : return new OGRPolygon();
1586 : }
1587 :
1588 807 : OGRGeometry *geom = nullptr;
1589 807 : OrganizePolygonMethod method = METHOD_NORMAL;
1590 807 : bool bHasCurves = false;
1591 :
1592 : /* -------------------------------------------------------------------- */
1593 : /* Trivial case of a single polygon. */
1594 : /* -------------------------------------------------------------------- */
1595 807 : if (nPolygonCount == 1)
1596 : {
1597 : OGRwkbGeometryType eType =
1598 20 : wkbFlatten(papoPolygons[0]->getGeometryType());
1599 :
1600 20 : bool bIsValid = true;
1601 :
1602 20 : if (eType != wkbPolygon && eType != wkbCurvePolygon)
1603 : {
1604 3 : CPLError(CE_Warning, CPLE_AppDefined,
1605 : "organizePolygons() received a non-Polygon geometry.");
1606 3 : bIsValid = false;
1607 3 : delete papoPolygons[0];
1608 3 : geom = new OGRPolygon();
1609 : }
1610 : else
1611 : {
1612 17 : geom = papoPolygons[0];
1613 : }
1614 :
1615 20 : papoPolygons[0] = nullptr;
1616 :
1617 20 : if (pbIsValidGeometry)
1618 9 : *pbIsValidGeometry = bIsValid;
1619 :
1620 20 : return geom;
1621 : }
1622 :
1623 787 : bool bUseFastVersion = true;
1624 787 : if (CPLTestBool(CPLGetConfigOption("OGR_DEBUG_ORGANIZE_POLYGONS", "NO")))
1625 : {
1626 : /* ------------------------------------------------------------------ */
1627 : /* A wee bit of a warning. */
1628 : /* ------------------------------------------------------------------ */
1629 : static int firstTime = 1;
1630 : // cppcheck-suppress knownConditionTrueFalse
1631 0 : if (!haveGEOS() && firstTime)
1632 : {
1633 0 : CPLDebug(
1634 : "OGR",
1635 : "In OGR_DEBUG_ORGANIZE_POLYGONS mode, GDAL should be built "
1636 : "with GEOS support enabled in order "
1637 : "OGRGeometryFactory::organizePolygons to provide reliable "
1638 : "results on complex polygons.");
1639 0 : firstTime = 0;
1640 : }
1641 : // cppcheck-suppress knownConditionTrueFalse
1642 0 : bUseFastVersion = !haveGEOS();
1643 : }
1644 :
1645 : /* -------------------------------------------------------------------- */
1646 : /* Setup per polygon envelope and area information. */
1647 : /* -------------------------------------------------------------------- */
1648 1574 : std::vector<sPolyExtended> asPolyEx;
1649 787 : asPolyEx.reserve(nPolygonCount);
1650 :
1651 787 : bool bValidTopology = true;
1652 787 : bool bMixedUpGeometries = false;
1653 787 : bool bFoundCCW = false;
1654 :
1655 787 : const char *pszMethodValue = CSLFetchNameValue(papszOptions, "METHOD");
1656 : const char *pszMethodValueOption =
1657 787 : CPLGetConfigOption("OGR_ORGANIZE_POLYGONS", nullptr);
1658 787 : if (pszMethodValueOption != nullptr && pszMethodValueOption[0] != '\0')
1659 3 : pszMethodValue = pszMethodValueOption;
1660 :
1661 787 : if (pszMethodValue != nullptr)
1662 : {
1663 307 : if (EQUAL(pszMethodValue, "SKIP"))
1664 : {
1665 7 : method = METHOD_SKIP;
1666 7 : bMixedUpGeometries = true;
1667 : }
1668 300 : else if (EQUAL(pszMethodValue, "ONLY_CCW"))
1669 : {
1670 246 : method = METHOD_ONLY_CCW;
1671 : }
1672 54 : else if (EQUAL(pszMethodValue, "CCW_INNER_JUST_AFTER_CW_OUTER"))
1673 : {
1674 0 : method = METHOD_CCW_INNER_JUST_AFTER_CW_OUTER;
1675 : }
1676 54 : else if (!EQUAL(pszMethodValue, "DEFAULT"))
1677 : {
1678 0 : CPLError(CE_Warning, CPLE_AppDefined,
1679 : "Unrecognized value for METHOD option : %s",
1680 : pszMethodValue);
1681 : }
1682 : }
1683 :
1684 787 : int nCountCWPolygon = 0;
1685 787 : int indexOfCWPolygon = -1;
1686 :
1687 3054 : for (int i = 0; i < nPolygonCount; i++)
1688 : {
1689 : OGRwkbGeometryType eType =
1690 2267 : wkbFlatten(papoPolygons[i]->getGeometryType());
1691 :
1692 2267 : if (eType != wkbPolygon && eType != wkbCurvePolygon)
1693 : {
1694 : // Ignore any points or lines that find their way in here.
1695 3 : CPLError(CE_Warning, CPLE_AppDefined,
1696 : "organizePolygons() received a non-Polygon geometry.");
1697 3 : delete papoPolygons[i];
1698 3 : continue;
1699 : }
1700 :
1701 4528 : sPolyExtended sPolyEx;
1702 :
1703 2264 : sPolyEx.nInitialIndex = i;
1704 2264 : sPolyEx.poGeometry = papoPolygons[i];
1705 2264 : sPolyEx.poPolygon = papoPolygons[i]->toCurvePolygon();
1706 :
1707 2264 : papoPolygons[i]->getEnvelope(&sPolyEx.sEnvelope);
1708 :
1709 2264 : if (eType == wkbCurvePolygon)
1710 33 : bHasCurves = true;
1711 2264 : if (!sPolyEx.poPolygon->IsEmpty() &&
1712 4528 : sPolyEx.poPolygon->getNumInteriorRings() == 0 &&
1713 2264 : sPolyEx.poPolygon->getExteriorRingCurve()->getNumPoints() >= 4)
1714 : {
1715 2264 : if (method != METHOD_CCW_INNER_JUST_AFTER_CW_OUTER)
1716 2264 : sPolyEx.dfArea = sPolyEx.poPolygon->get_Area();
1717 2264 : sPolyEx.poExteriorRing = sPolyEx.poPolygon->getExteriorRingCurve();
1718 2264 : sPolyEx.poExteriorRing->StartPoint(&sPolyEx.poAPoint);
1719 2264 : if (eType == wkbPolygon)
1720 : {
1721 2231 : sPolyEx.bIsCW = CPL_TO_BOOL(
1722 2231 : sPolyEx.poExteriorRing->toLinearRing()->isClockwise());
1723 2231 : sPolyEx.bIsPolygon = true;
1724 : }
1725 : else
1726 : {
1727 33 : OGRLineString *poLS = sPolyEx.poExteriorRing->CurveToLine();
1728 66 : OGRLinearRing oLR;
1729 33 : oLR.addSubLineString(poLS);
1730 33 : sPolyEx.bIsCW = CPL_TO_BOOL(oLR.isClockwise());
1731 33 : sPolyEx.bIsPolygon = false;
1732 33 : delete poLS;
1733 : }
1734 2264 : if (sPolyEx.bIsCW)
1735 : {
1736 1426 : indexOfCWPolygon = i;
1737 1426 : nCountCWPolygon++;
1738 : }
1739 2264 : if (!bFoundCCW)
1740 1549 : bFoundCCW = !(sPolyEx.bIsCW);
1741 : }
1742 : else
1743 : {
1744 0 : if (!bMixedUpGeometries)
1745 : {
1746 0 : CPLError(CE_Warning, CPLE_AppDefined,
1747 : "organizePolygons() received an unexpected geometry. "
1748 : "Either a polygon with interior rings, or a polygon "
1749 : "with less than 4 points, or a non-Polygon geometry. "
1750 : "Return arguments as a collection.");
1751 0 : bMixedUpGeometries = true;
1752 : }
1753 : }
1754 :
1755 2264 : asPolyEx.push_back(std::move(sPolyEx));
1756 : }
1757 :
1758 : // If we are in ONLY_CCW mode and that we have found that there is only one
1759 : // outer ring, then it is pretty easy : we can assume that all other rings
1760 : // are inside.
1761 787 : if ((method == METHOD_ONLY_CCW ||
1762 246 : method == METHOD_CCW_INNER_JUST_AFTER_CW_OUTER) &&
1763 77 : nCountCWPolygon == 1 && bUseFastVersion)
1764 : {
1765 77 : OGRCurvePolygon *poCP = asPolyEx[indexOfCWPolygon].poPolygon;
1766 279 : for (int i = 0; i < static_cast<int>(asPolyEx.size()); i++)
1767 : {
1768 202 : if (i != indexOfCWPolygon)
1769 : {
1770 125 : poCP->addRingDirectly(
1771 125 : asPolyEx[i].poPolygon->stealExteriorRingCurve());
1772 125 : delete asPolyEx[i].poPolygon;
1773 : }
1774 : }
1775 :
1776 77 : if (pbIsValidGeometry)
1777 75 : *pbIsValidGeometry = TRUE;
1778 77 : return poCP;
1779 : }
1780 :
1781 710 : if (method == METHOD_CCW_INNER_JUST_AFTER_CW_OUTER && asPolyEx[0].bIsCW)
1782 : {
1783 : // Inner rings are CCW oriented and follow immediately the outer
1784 : // ring (that is CW oriented) in which they are included.
1785 0 : OGRMultiSurface *poMulti = nullptr;
1786 0 : OGRCurvePolygon *poCur = asPolyEx[0].poPolygon;
1787 0 : OGRGeometry *poRet = poCur;
1788 : // We have already checked that the first ring is CW.
1789 0 : OGREnvelope *psEnvelope = &(asPolyEx[0].sEnvelope);
1790 0 : for (std::size_t i = 1; i < asPolyEx.size(); i++)
1791 : {
1792 0 : if (asPolyEx[i].bIsCW)
1793 : {
1794 0 : if (poMulti == nullptr)
1795 : {
1796 0 : if (bHasCurves)
1797 0 : poMulti = new OGRMultiSurface();
1798 : else
1799 0 : poMulti = new OGRMultiPolygon();
1800 0 : poRet = poMulti;
1801 0 : poMulti->addGeometryDirectly(poCur);
1802 : }
1803 0 : poCur = asPolyEx[i].poPolygon;
1804 0 : poMulti->addGeometryDirectly(poCur);
1805 0 : psEnvelope = &(asPolyEx[i].sEnvelope);
1806 : }
1807 : else
1808 : {
1809 0 : poCur->addRingDirectly(
1810 0 : asPolyEx[i].poPolygon->stealExteriorRingCurve());
1811 0 : if (!(asPolyEx[i].poAPoint.getX() >= psEnvelope->MinX &&
1812 0 : asPolyEx[i].poAPoint.getX() <= psEnvelope->MaxX &&
1813 0 : asPolyEx[i].poAPoint.getY() >= psEnvelope->MinY &&
1814 0 : asPolyEx[i].poAPoint.getY() <= psEnvelope->MaxY))
1815 : {
1816 0 : CPLError(CE_Warning, CPLE_AppDefined,
1817 : "Part %d does not respect "
1818 : "CCW_INNER_JUST_AFTER_CW_OUTER rule",
1819 : static_cast<int>(i));
1820 : }
1821 0 : delete asPolyEx[i].poPolygon;
1822 : }
1823 : }
1824 :
1825 0 : if (pbIsValidGeometry)
1826 0 : *pbIsValidGeometry = TRUE;
1827 0 : return poRet;
1828 : }
1829 710 : else if (method == METHOD_CCW_INNER_JUST_AFTER_CW_OUTER)
1830 : {
1831 0 : method = METHOD_ONLY_CCW;
1832 0 : for (std::size_t i = 0; i < asPolyEx.size(); i++)
1833 0 : asPolyEx[i].dfArea = asPolyEx[i].poPolygon->get_Area();
1834 : }
1835 :
1836 : // Emits a warning if the number of parts is sufficiently big to anticipate
1837 : // for very long computation time, and the user didn't specify an explicit
1838 : // method.
1839 710 : if (nPolygonCount > N_CRITICAL_PART_NUMBER && method == METHOD_NORMAL &&
1840 : pszMethodValue == nullptr)
1841 : {
1842 : static int firstTime = 1;
1843 0 : if (firstTime)
1844 : {
1845 0 : if (bFoundCCW)
1846 : {
1847 0 : CPLError(
1848 : CE_Warning, CPLE_AppDefined,
1849 : "organizePolygons() received a polygon with more than %d "
1850 : "parts. The processing may be really slow. "
1851 : "You can skip the processing by setting METHOD=SKIP, "
1852 : "or only make it analyze counter-clock wise parts by "
1853 : "setting METHOD=ONLY_CCW if you can assume that the "
1854 : "outline of holes is counter-clock wise defined",
1855 : N_CRITICAL_PART_NUMBER);
1856 : }
1857 : else
1858 : {
1859 0 : CPLError(
1860 : CE_Warning, CPLE_AppDefined,
1861 : "organizePolygons() received a polygon with more than %d "
1862 : "parts. The processing may be really slow. "
1863 : "You can skip the processing by setting METHOD=SKIP.",
1864 : N_CRITICAL_PART_NUMBER);
1865 : }
1866 0 : firstTime = 0;
1867 : }
1868 : }
1869 :
1870 : /* This a nulti-step algorithm :
1871 : 1) Sort polygons by descending areas
1872 : 2) For each polygon of rank i, find its smallest enclosing polygon
1873 : among the polygons of rank [i-1 ... 0]. If there are no such polygon,
1874 : this is a top-level polygon. Otherwise, depending on if the enclosing
1875 : polygon is top-level or not, we can decide if we are top-level or not
1876 : 3) Re-sort the polygons to retrieve their initial order (nicer for
1877 : some applications)
1878 : 4) For each non top-level polygon (= inner ring), add it to its
1879 : outer ring
1880 : 5) Add the top-level polygons to the multipolygon
1881 :
1882 : Complexity : O(nPolygonCount^2)
1883 : */
1884 :
1885 : /* Compute how each polygon relate to the other ones
1886 : To save a bit of computation we always begin the computation by a test
1887 : on the envelope. We also take into account the areas to avoid some
1888 : useless tests. (A contains B implies envelop(A) contains envelop(B)
1889 : and area(A) > area(B)) In practice, we can hope that few full geometry
1890 : intersection of inclusion test is done:
1891 : * if the polygons are well separated geographically (a set of islands
1892 : for example), no full geometry intersection or inclusion test is done.
1893 : (the envelopes don't intersect each other)
1894 :
1895 : * if the polygons are 'lake inside an island inside a lake inside an
1896 : area' and that each polygon is much smaller than its enclosing one,
1897 : their bounding boxes are strictly contained into each other, and thus,
1898 : no full geometry intersection or inclusion test is done
1899 : */
1900 :
1901 710 : if (!bMixedUpGeometries)
1902 : {
1903 : // STEP 1: Sort polygons by descending area.
1904 703 : std::sort(asPolyEx.begin(), asPolyEx.end(),
1905 : OGRGeometryFactoryCompareArea);
1906 : }
1907 710 : papoPolygons = nullptr; // Just to use to avoid it afterwards.
1908 :
1909 : /* -------------------------------------------------------------------- */
1910 : /* Compute relationships, if things seem well structured. */
1911 : /* -------------------------------------------------------------------- */
1912 :
1913 : // The first (largest) polygon is necessarily top-level.
1914 710 : asPolyEx[0].bIsTopLevel = true;
1915 710 : asPolyEx[0].poEnclosingPolygon = nullptr;
1916 :
1917 710 : int nCountTopLevel = 1;
1918 :
1919 : // STEP 2.
1920 4103 : for (int i = 1; !bMixedUpGeometries && bValidTopology &&
1921 2048 : i < static_cast<int>(asPolyEx.size());
1922 : i++)
1923 : {
1924 1345 : if (method == METHOD_ONLY_CCW && asPolyEx[i].bIsCW)
1925 : {
1926 196 : nCountTopLevel++;
1927 196 : asPolyEx[i].bIsTopLevel = true;
1928 196 : asPolyEx[i].poEnclosingPolygon = nullptr;
1929 196 : continue;
1930 : }
1931 :
1932 1149 : int j = i - 1; // Used after for.
1933 2439 : for (; bValidTopology && j >= 0; j--)
1934 : {
1935 1968 : bool b_i_inside_j = false;
1936 :
1937 1968 : if (method == METHOD_ONLY_CCW && asPolyEx[j].bIsCW == false)
1938 : {
1939 : // In that mode, i which is CCW if we reach here can only be
1940 : // included in a CW polygon.
1941 8 : continue;
1942 : }
1943 :
1944 1960 : if (asPolyEx[j].sEnvelope.Contains(asPolyEx[i].sEnvelope))
1945 : {
1946 684 : if (bUseFastVersion)
1947 : {
1948 684 : if (method == METHOD_ONLY_CCW && j == 0)
1949 : {
1950 : // We are testing if a CCW ring is in the biggest CW
1951 : // ring It *must* be inside as this is the last
1952 : // candidate, otherwise the winding order rules is
1953 : // broken.
1954 123 : b_i_inside_j = true;
1955 : }
1956 1122 : else if (asPolyEx[i].bIsPolygon && asPolyEx[j].bIsPolygon &&
1957 561 : asPolyEx[j]
1958 561 : .poExteriorRing->toLinearRing()
1959 561 : ->isPointOnRingBoundary(&asPolyEx[i].poAPoint,
1960 : FALSE))
1961 : {
1962 : OGRLinearRing *poLR_i =
1963 16 : asPolyEx[i].poExteriorRing->toLinearRing();
1964 : OGRLinearRing *poLR_j =
1965 16 : asPolyEx[j].poExteriorRing->toLinearRing();
1966 :
1967 : // If the point of i is on the boundary of j, we will
1968 : // iterate over the other points of i.
1969 16 : const int nPoints = poLR_i->getNumPoints();
1970 16 : int k = 1; // Used after for.
1971 32 : OGRPoint previousPoint = asPolyEx[i].poAPoint;
1972 31 : for (; k < nPoints; k++)
1973 : {
1974 30 : OGRPoint point;
1975 30 : poLR_i->getPoint(k, &point);
1976 32 : if (point.getX() == previousPoint.getX() &&
1977 2 : point.getY() == previousPoint.getY())
1978 : {
1979 0 : continue;
1980 : }
1981 30 : if (poLR_j->isPointOnRingBoundary(&point, FALSE))
1982 : {
1983 : // If it is on the boundary of j, iterate again.
1984 : }
1985 15 : else if (poLR_j->isPointInRing(&point, FALSE))
1986 : {
1987 : // If then point is strictly included in j, then
1988 : // i is considered inside j.
1989 13 : b_i_inside_j = true;
1990 13 : break;
1991 : }
1992 : else
1993 : {
1994 : // If it is outside, then i cannot be inside j.
1995 2 : break;
1996 : }
1997 15 : previousPoint = point;
1998 : }
1999 16 : if (!b_i_inside_j && k == nPoints && nPoints > 2)
2000 : {
2001 : // All points of i are on the boundary of j.
2002 : // Take a point in the middle of a segment of i and
2003 : // test it against j.
2004 1 : poLR_i->getPoint(0, &previousPoint);
2005 2 : for (k = 1; k < nPoints; k++)
2006 : {
2007 2 : OGRPoint point;
2008 2 : poLR_i->getPoint(k, &point);
2009 2 : if (point.getX() == previousPoint.getX() &&
2010 0 : point.getY() == previousPoint.getY())
2011 : {
2012 0 : continue;
2013 : }
2014 2 : OGRPoint pointMiddle;
2015 2 : pointMiddle.setX(
2016 2 : (point.getX() + previousPoint.getX()) / 2);
2017 2 : pointMiddle.setY(
2018 2 : (point.getY() + previousPoint.getY()) / 2);
2019 2 : if (poLR_j->isPointOnRingBoundary(&pointMiddle,
2020 2 : FALSE))
2021 : {
2022 : // If it is on the boundary of j, iterate
2023 : // again.
2024 : }
2025 1 : else if (poLR_j->isPointInRing(&pointMiddle,
2026 1 : FALSE))
2027 : {
2028 : // If then point is strictly included in j,
2029 : // then i is considered inside j.
2030 1 : b_i_inside_j = true;
2031 1 : break;
2032 : }
2033 : else
2034 : {
2035 : // If it is outside, then i cannot be inside
2036 : // j.
2037 0 : break;
2038 : }
2039 1 : previousPoint = point;
2040 : }
2041 : }
2042 : }
2043 : // Note that isPointInRing only test strict inclusion in the
2044 : // ring.
2045 1090 : else if (asPolyEx[i].bIsPolygon && asPolyEx[j].bIsPolygon &&
2046 545 : asPolyEx[j]
2047 545 : .poExteriorRing->toLinearRing()
2048 545 : ->isPointInRing(&asPolyEx[i].poAPoint, FALSE))
2049 : {
2050 541 : b_i_inside_j = true;
2051 : }
2052 : }
2053 0 : else if (asPolyEx[j].poPolygon->Contains(asPolyEx[i].poPolygon))
2054 : {
2055 0 : b_i_inside_j = true;
2056 : }
2057 : }
2058 :
2059 1960 : if (b_i_inside_j)
2060 : {
2061 678 : if (asPolyEx[j].bIsTopLevel)
2062 : {
2063 : // We are a lake.
2064 677 : asPolyEx[i].bIsTopLevel = false;
2065 677 : asPolyEx[i].poEnclosingPolygon = asPolyEx[j].poPolygon;
2066 : }
2067 : else
2068 : {
2069 : // We are included in a something not toplevel (a lake),
2070 : // so in OGCSF we are considered as toplevel too.
2071 1 : nCountTopLevel++;
2072 1 : asPolyEx[i].bIsTopLevel = true;
2073 1 : asPolyEx[i].poEnclosingPolygon = nullptr;
2074 : }
2075 678 : break;
2076 : }
2077 : // Use Overlaps instead of Intersects to be more
2078 : // tolerant about touching polygons.
2079 0 : else if (bUseFastVersion ||
2080 1282 : !asPolyEx[i].sEnvelope.Intersects(asPolyEx[j].sEnvelope) ||
2081 0 : !asPolyEx[i].poPolygon->Overlaps(asPolyEx[j].poPolygon))
2082 : {
2083 : }
2084 : else
2085 : {
2086 : // Bad... The polygons are intersecting but no one is
2087 : // contained inside the other one. This is a really broken
2088 : // case. We just make a multipolygon with the whole set of
2089 : // polygons.
2090 0 : bValidTopology = false;
2091 : #ifdef DEBUG
2092 0 : char *wkt1 = nullptr;
2093 0 : char *wkt2 = nullptr;
2094 0 : asPolyEx[i].poPolygon->exportToWkt(&wkt1);
2095 0 : asPolyEx[j].poPolygon->exportToWkt(&wkt2);
2096 0 : CPLDebug("OGR",
2097 : "Bad intersection for polygons %d and %d\n"
2098 : "geom %d: %s\n"
2099 : "geom %d: %s",
2100 : static_cast<int>(i), j, static_cast<int>(i), wkt1, j,
2101 : wkt2);
2102 0 : CPLFree(wkt1);
2103 0 : CPLFree(wkt2);
2104 : #endif
2105 : }
2106 : }
2107 :
2108 1149 : if (j < 0)
2109 : {
2110 : // We come here because we are not included in anything.
2111 : // We are toplevel.
2112 471 : nCountTopLevel++;
2113 471 : asPolyEx[i].bIsTopLevel = true;
2114 471 : asPolyEx[i].poEnclosingPolygon = nullptr;
2115 : }
2116 : }
2117 :
2118 710 : if (pbIsValidGeometry)
2119 220 : *pbIsValidGeometry = bValidTopology && !bMixedUpGeometries;
2120 :
2121 : /* --------------------------------------------------------------------- */
2122 : /* Things broke down - just mark everything as top-level so it gets */
2123 : /* turned into a multipolygon. */
2124 : /* --------------------------------------------------------------------- */
2125 710 : if (!bValidTopology || bMixedUpGeometries)
2126 : {
2127 21 : for (auto &sPolyEx : asPolyEx)
2128 : {
2129 14 : sPolyEx.bIsTopLevel = true;
2130 : }
2131 7 : nCountTopLevel = static_cast<int>(asPolyEx.size());
2132 : }
2133 :
2134 : /* -------------------------------------------------------------------- */
2135 : /* Try to turn into one or more polygons based on the ring */
2136 : /* relationships. */
2137 : /* -------------------------------------------------------------------- */
2138 : // STEP 3: Sort again in initial order.
2139 710 : std::sort(asPolyEx.begin(), asPolyEx.end(),
2140 : OGRGeometryFactoryCompareByIndex);
2141 :
2142 : // STEP 4: Add holes as rings of their enclosing polygon.
2143 2772 : for (auto &sPolyEx : asPolyEx)
2144 : {
2145 2062 : if (sPolyEx.bIsTopLevel == false)
2146 : {
2147 677 : sPolyEx.poEnclosingPolygon->addRingDirectly(
2148 677 : sPolyEx.poPolygon->stealExteriorRingCurve());
2149 677 : delete sPolyEx.poPolygon;
2150 : }
2151 1385 : else if (nCountTopLevel == 1)
2152 : {
2153 68 : geom = sPolyEx.poPolygon;
2154 : }
2155 : }
2156 :
2157 : // STEP 5: Add toplevel polygons.
2158 710 : if (nCountTopLevel > 1)
2159 : {
2160 : OGRGeometryCollection *poGC =
2161 642 : bHasCurves ? new OGRMultiSurface() : new OGRMultiPolygon();
2162 2541 : for (auto &sPolyEx : asPolyEx)
2163 : {
2164 1899 : if (sPolyEx.bIsTopLevel)
2165 : {
2166 1317 : poGC->addGeometryDirectly(sPolyEx.poPolygon);
2167 : }
2168 : }
2169 642 : geom = poGC;
2170 : }
2171 :
2172 710 : return geom;
2173 : }
2174 :
2175 : /************************************************************************/
2176 : /* createFromGML() */
2177 : /************************************************************************/
2178 :
2179 : /**
2180 : * \brief Create geometry from GML.
2181 : *
2182 : * This method translates a fragment of GML containing only the geometry
2183 : * portion into a corresponding OGRGeometry. There are many limitations
2184 : * on the forms of GML geometries supported by this parser, but they are
2185 : * too numerous to list here.
2186 : *
2187 : * The following GML2 elements are parsed : Point, LineString, Polygon,
2188 : * MultiPoint, MultiLineString, MultiPolygon, MultiGeometry.
2189 : *
2190 : * The following GML3 elements are parsed : Surface,
2191 : * MultiSurface, PolygonPatch, Triangle, Rectangle, Curve, MultiCurve,
2192 : * LineStringSegment, Arc, Circle, CompositeSurface, OrientableSurface, Solid,
2193 : * Tin, TriangulatedSurface.
2194 : *
2195 : * Arc and Circle elements are returned as curves by default. Stroking to
2196 : * linestrings can be done with
2197 : * OGR_G_ForceTo(hGeom, OGR_GT_GetLinear(OGR_G_GetGeometryType(hGeom)), NULL).
2198 : * A 4 degrees step is used by default, unless the user
2199 : * has overridden the value with the OGR_ARC_STEPSIZE configuration variable.
2200 : *
2201 : * The C function OGR_G_CreateFromGML() is the same as this method.
2202 : *
2203 : * @param pszData The GML fragment for the geometry.
2204 : *
2205 : * @return a geometry on success, or NULL on error.
2206 : *
2207 : * @see OGR_G_ForceTo()
2208 : * @see OGR_GT_GetLinear()
2209 : * @see OGR_G_GetGeometryType()
2210 : */
2211 :
2212 0 : OGRGeometry *OGRGeometryFactory::createFromGML(const char *pszData)
2213 :
2214 : {
2215 : OGRGeometryH hGeom;
2216 :
2217 0 : hGeom = OGR_G_CreateFromGML(pszData);
2218 :
2219 0 : return OGRGeometry::FromHandle(hGeom);
2220 : }
2221 :
2222 : /************************************************************************/
2223 : /* createFromGEOS() */
2224 : /************************************************************************/
2225 :
2226 : /** Builds a OGRGeometry* from a GEOSGeom.
2227 : * @param hGEOSCtxt GEOS context
2228 : * @param geosGeom GEOS geometry
2229 : * @return a OGRGeometry*
2230 : */
2231 1311 : OGRGeometry *OGRGeometryFactory::createFromGEOS(
2232 : UNUSED_IF_NO_GEOS GEOSContextHandle_t hGEOSCtxt,
2233 : UNUSED_IF_NO_GEOS GEOSGeom geosGeom)
2234 :
2235 : {
2236 : #ifndef HAVE_GEOS
2237 :
2238 : CPLError(CE_Failure, CPLE_NotSupported, "GEOS support not enabled.");
2239 : return nullptr;
2240 :
2241 : #else
2242 :
2243 1311 : size_t nSize = 0;
2244 1311 : unsigned char *pabyBuf = nullptr;
2245 1311 : OGRGeometry *poGeometry = nullptr;
2246 :
2247 : // Special case as POINT EMPTY cannot be translated to WKB.
2248 1366 : if (GEOSGeomTypeId_r(hGEOSCtxt, geosGeom) == GEOS_POINT &&
2249 55 : GEOSisEmpty_r(hGEOSCtxt, geosGeom))
2250 4 : return new OGRPoint();
2251 :
2252 : const int nCoordDim =
2253 1306 : GEOSGeom_getCoordinateDimension_r(hGEOSCtxt, geosGeom);
2254 1307 : GEOSWKBWriter *wkbwriter = GEOSWKBWriter_create_r(hGEOSCtxt);
2255 1306 : GEOSWKBWriter_setOutputDimension_r(hGEOSCtxt, wkbwriter, nCoordDim);
2256 1306 : pabyBuf = GEOSWKBWriter_write_r(hGEOSCtxt, wkbwriter, geosGeom, &nSize);
2257 1306 : GEOSWKBWriter_destroy_r(hGEOSCtxt, wkbwriter);
2258 :
2259 1307 : if (pabyBuf == nullptr || nSize == 0)
2260 : {
2261 0 : return nullptr;
2262 : }
2263 :
2264 1307 : if (OGRGeometryFactory::createFromWkb(pabyBuf, nullptr, &poGeometry,
2265 1307 : static_cast<int>(nSize)) !=
2266 : OGRERR_NONE)
2267 : {
2268 0 : poGeometry = nullptr;
2269 : }
2270 :
2271 1307 : GEOSFree_r(hGEOSCtxt, pabyBuf);
2272 :
2273 1305 : return poGeometry;
2274 :
2275 : #endif // HAVE_GEOS
2276 : }
2277 :
2278 : /************************************************************************/
2279 : /* haveGEOS() */
2280 : /************************************************************************/
2281 :
2282 : /**
2283 : * \brief Test if GEOS enabled.
2284 : *
2285 : * This static method returns TRUE if GEOS support is built into OGR,
2286 : * otherwise it returns FALSE.
2287 : *
2288 : * @return TRUE if available, otherwise FALSE.
2289 : */
2290 :
2291 39988 : bool OGRGeometryFactory::haveGEOS()
2292 :
2293 : {
2294 : #ifndef HAVE_GEOS
2295 : return false;
2296 : #else
2297 39988 : return true;
2298 : #endif
2299 : }
2300 :
2301 : /************************************************************************/
2302 : /* createFromFgf() */
2303 : /************************************************************************/
2304 :
2305 : /**
2306 : * \brief Create a geometry object of the appropriate type from its FGF (FDO
2307 : * Geometry Format) binary representation.
2308 : *
2309 : * Also note that this is a static method, and that there
2310 : * is no need to instantiate an OGRGeometryFactory object.
2311 : *
2312 : * The C function OGR_G_CreateFromFgf() is the same as this method.
2313 : *
2314 : * @param pabyData pointer to the input BLOB data.
2315 : * @param poSR pointer to the spatial reference to be assigned to the
2316 : * created geometry object. This may be NULL.
2317 : * @param ppoReturn the newly created geometry object will be assigned to the
2318 : * indicated pointer on return. This will be NULL in case
2319 : * of failure, but NULL might be a valid return for a NULL
2320 : * shape.
2321 : * @param nBytes the number of bytes available in pabyData.
2322 : * @param pnBytesConsumed if not NULL, it will be set to the number of bytes
2323 : * consumed (at most nBytes).
2324 : *
2325 : * @return OGRERR_NONE if all goes well, otherwise any of
2326 : * OGRERR_NOT_ENOUGH_DATA, OGRERR_UNSUPPORTED_GEOMETRY_TYPE, or
2327 : * OGRERR_CORRUPT_DATA may be returned.
2328 : */
2329 :
2330 352 : OGRErr OGRGeometryFactory::createFromFgf(const void *pabyData,
2331 : OGRSpatialReference *poSR,
2332 : OGRGeometry **ppoReturn, int nBytes,
2333 : int *pnBytesConsumed)
2334 :
2335 : {
2336 352 : return createFromFgfInternal(static_cast<const GByte *>(pabyData), poSR,
2337 352 : ppoReturn, nBytes, pnBytesConsumed, 0);
2338 : }
2339 :
2340 : /************************************************************************/
2341 : /* createFromFgfInternal() */
2342 : /************************************************************************/
2343 :
2344 355 : OGRErr OGRGeometryFactory::createFromFgfInternal(
2345 : const unsigned char *pabyData, OGRSpatialReference *poSR,
2346 : OGRGeometry **ppoReturn, int nBytes, int *pnBytesConsumed, int nRecLevel)
2347 : {
2348 : // Arbitrary value, but certainly large enough for reasonable usages.
2349 355 : if (nRecLevel == 32)
2350 : {
2351 0 : CPLError(CE_Failure, CPLE_AppDefined,
2352 : "Too many recursion levels (%d) while parsing FGF geometry.",
2353 : nRecLevel);
2354 0 : return OGRERR_CORRUPT_DATA;
2355 : }
2356 :
2357 355 : *ppoReturn = nullptr;
2358 :
2359 355 : if (nBytes < 4)
2360 109 : return OGRERR_NOT_ENOUGH_DATA;
2361 :
2362 : /* -------------------------------------------------------------------- */
2363 : /* Decode the geometry type. */
2364 : /* -------------------------------------------------------------------- */
2365 246 : GInt32 nGType = 0;
2366 246 : memcpy(&nGType, pabyData + 0, 4);
2367 246 : CPL_LSBPTR32(&nGType);
2368 :
2369 246 : if (nGType < 0 || nGType > 13)
2370 232 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
2371 :
2372 : /* -------------------------------------------------------------------- */
2373 : /* Decode the dimensionality if appropriate. */
2374 : /* -------------------------------------------------------------------- */
2375 14 : int nTupleSize = 0;
2376 14 : GInt32 nGDim = 0;
2377 :
2378 : // TODO: Why is this a switch?
2379 14 : switch (nGType)
2380 : {
2381 9 : case 1: // Point
2382 : case 2: // LineString
2383 : case 3: // Polygon
2384 9 : if (nBytes < 8)
2385 0 : return OGRERR_NOT_ENOUGH_DATA;
2386 :
2387 9 : memcpy(&nGDim, pabyData + 4, 4);
2388 9 : CPL_LSBPTR32(&nGDim);
2389 :
2390 9 : if (nGDim < 0 || nGDim > 3)
2391 0 : return OGRERR_CORRUPT_DATA;
2392 :
2393 9 : nTupleSize = 2;
2394 9 : if (nGDim & 0x01) // Z
2395 1 : nTupleSize++;
2396 9 : if (nGDim & 0x02) // M
2397 0 : nTupleSize++;
2398 :
2399 9 : break;
2400 :
2401 5 : default:
2402 5 : break;
2403 : }
2404 :
2405 14 : OGRGeometry *poGeom = nullptr;
2406 :
2407 : /* -------------------------------------------------------------------- */
2408 : /* None */
2409 : /* -------------------------------------------------------------------- */
2410 14 : if (nGType == 0)
2411 : {
2412 0 : if (pnBytesConsumed)
2413 0 : *pnBytesConsumed = 4;
2414 : }
2415 :
2416 : /* -------------------------------------------------------------------- */
2417 : /* Point */
2418 : /* -------------------------------------------------------------------- */
2419 14 : else if (nGType == 1)
2420 : {
2421 3 : if (nBytes < nTupleSize * 8 + 8)
2422 0 : return OGRERR_NOT_ENOUGH_DATA;
2423 :
2424 3 : double adfTuple[4] = {0.0, 0.0, 0.0, 0.0};
2425 3 : memcpy(adfTuple, pabyData + 8, nTupleSize * 8);
2426 : #ifdef CPL_MSB
2427 : for (int iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++)
2428 : CPL_SWAP64PTR(adfTuple + iOrdinal);
2429 : #endif
2430 3 : if (nTupleSize > 2)
2431 1 : poGeom = new OGRPoint(adfTuple[0], adfTuple[1], adfTuple[2]);
2432 : else
2433 2 : poGeom = new OGRPoint(adfTuple[0], adfTuple[1]);
2434 :
2435 3 : if (pnBytesConsumed)
2436 1 : *pnBytesConsumed = 8 + nTupleSize * 8;
2437 : }
2438 :
2439 : /* -------------------------------------------------------------------- */
2440 : /* LineString */
2441 : /* -------------------------------------------------------------------- */
2442 11 : else if (nGType == 2)
2443 : {
2444 2 : if (nBytes < 12)
2445 0 : return OGRERR_NOT_ENOUGH_DATA;
2446 :
2447 2 : GInt32 nPointCount = 0;
2448 2 : memcpy(&nPointCount, pabyData + 8, 4);
2449 2 : CPL_LSBPTR32(&nPointCount);
2450 :
2451 2 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
2452 0 : return OGRERR_CORRUPT_DATA;
2453 :
2454 2 : if (nBytes - 12 < nTupleSize * 8 * nPointCount)
2455 0 : return OGRERR_NOT_ENOUGH_DATA;
2456 :
2457 2 : OGRLineString *poLS = new OGRLineString();
2458 2 : poGeom = poLS;
2459 2 : poLS->setNumPoints(nPointCount);
2460 :
2461 4 : for (int iPoint = 0; iPoint < nPointCount; iPoint++)
2462 : {
2463 2 : double adfTuple[4] = {0.0, 0.0, 0.0, 0.0};
2464 2 : memcpy(adfTuple, pabyData + 12 + 8 * nTupleSize * iPoint,
2465 2 : nTupleSize * 8);
2466 : #ifdef CPL_MSB
2467 : for (int iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++)
2468 : CPL_SWAP64PTR(adfTuple + iOrdinal);
2469 : #endif
2470 2 : if (nTupleSize > 2)
2471 0 : poLS->setPoint(iPoint, adfTuple[0], adfTuple[1], adfTuple[2]);
2472 : else
2473 2 : poLS->setPoint(iPoint, adfTuple[0], adfTuple[1]);
2474 : }
2475 :
2476 2 : if (pnBytesConsumed)
2477 0 : *pnBytesConsumed = 12 + nTupleSize * 8 * nPointCount;
2478 : }
2479 :
2480 : /* -------------------------------------------------------------------- */
2481 : /* Polygon */
2482 : /* -------------------------------------------------------------------- */
2483 9 : else if (nGType == 3)
2484 : {
2485 4 : if (nBytes < 12)
2486 0 : return OGRERR_NOT_ENOUGH_DATA;
2487 :
2488 4 : GInt32 nRingCount = 0;
2489 4 : memcpy(&nRingCount, pabyData + 8, 4);
2490 4 : CPL_LSBPTR32(&nRingCount);
2491 :
2492 4 : if (nRingCount < 0 || nRingCount > INT_MAX / 4)
2493 0 : return OGRERR_CORRUPT_DATA;
2494 :
2495 : // Each ring takes at least 4 bytes.
2496 4 : if (nBytes - 12 < nRingCount * 4)
2497 0 : return OGRERR_NOT_ENOUGH_DATA;
2498 :
2499 4 : int nNextByte = 12;
2500 :
2501 4 : OGRPolygon *poPoly = new OGRPolygon();
2502 4 : poGeom = poPoly;
2503 :
2504 10 : for (int iRing = 0; iRing < nRingCount; iRing++)
2505 : {
2506 6 : if (nBytes - nNextByte < 4)
2507 : {
2508 0 : delete poGeom;
2509 0 : return OGRERR_NOT_ENOUGH_DATA;
2510 : }
2511 :
2512 6 : GInt32 nPointCount = 0;
2513 6 : memcpy(&nPointCount, pabyData + nNextByte, 4);
2514 6 : CPL_LSBPTR32(&nPointCount);
2515 :
2516 6 : if (nPointCount < 0 || nPointCount > INT_MAX / (nTupleSize * 8))
2517 : {
2518 0 : delete poGeom;
2519 0 : return OGRERR_CORRUPT_DATA;
2520 : }
2521 :
2522 6 : nNextByte += 4;
2523 :
2524 6 : if (nBytes - nNextByte < nTupleSize * 8 * nPointCount)
2525 : {
2526 0 : delete poGeom;
2527 0 : return OGRERR_NOT_ENOUGH_DATA;
2528 : }
2529 :
2530 6 : OGRLinearRing *poLR = new OGRLinearRing();
2531 6 : poLR->setNumPoints(nPointCount);
2532 :
2533 12 : for (int iPoint = 0; iPoint < nPointCount; iPoint++)
2534 : {
2535 6 : double adfTuple[4] = {0.0, 0.0, 0.0, 0.0};
2536 6 : memcpy(adfTuple, pabyData + nNextByte, nTupleSize * 8);
2537 6 : nNextByte += nTupleSize * 8;
2538 :
2539 : #ifdef CPL_MSB
2540 : for (int iOrdinal = 0; iOrdinal < nTupleSize; iOrdinal++)
2541 : CPL_SWAP64PTR(adfTuple + iOrdinal);
2542 : #endif
2543 6 : if (nTupleSize > 2)
2544 0 : poLR->setPoint(iPoint, adfTuple[0], adfTuple[1],
2545 : adfTuple[2]);
2546 : else
2547 6 : poLR->setPoint(iPoint, adfTuple[0], adfTuple[1]);
2548 : }
2549 :
2550 6 : poPoly->addRingDirectly(poLR);
2551 : }
2552 :
2553 4 : if (pnBytesConsumed)
2554 2 : *pnBytesConsumed = nNextByte;
2555 : }
2556 :
2557 : /* -------------------------------------------------------------------- */
2558 : /* GeometryCollections of various kinds. */
2559 : /* -------------------------------------------------------------------- */
2560 5 : else if (nGType == 4 // MultiPoint
2561 5 : || nGType == 5 // MultiLineString
2562 5 : || nGType == 6 // MultiPolygon
2563 3 : || nGType == 7) // MultiGeometry
2564 : {
2565 5 : if (nBytes < 8)
2566 3 : return OGRERR_NOT_ENOUGH_DATA;
2567 :
2568 5 : GInt32 nGeomCount = 0;
2569 5 : memcpy(&nGeomCount, pabyData + 4, 4);
2570 5 : CPL_LSBPTR32(&nGeomCount);
2571 :
2572 5 : if (nGeomCount < 0 || nGeomCount > INT_MAX / 4)
2573 0 : return OGRERR_CORRUPT_DATA;
2574 :
2575 : // Each geometry takes at least 4 bytes.
2576 5 : if (nBytes - 8 < 4 * nGeomCount)
2577 2 : return OGRERR_NOT_ENOUGH_DATA;
2578 :
2579 3 : OGRGeometryCollection *poGC = nullptr;
2580 3 : if (nGType == 4)
2581 0 : poGC = new OGRMultiPoint();
2582 3 : else if (nGType == 5)
2583 0 : poGC = new OGRMultiLineString();
2584 3 : else if (nGType == 6)
2585 1 : poGC = new OGRMultiPolygon();
2586 2 : else if (nGType == 7)
2587 2 : poGC = new OGRGeometryCollection();
2588 :
2589 3 : int nBytesUsed = 8;
2590 :
2591 5 : for (int iGeom = 0; iGeom < nGeomCount; iGeom++)
2592 : {
2593 3 : int nThisGeomSize = 0;
2594 3 : OGRGeometry *poThisGeom = nullptr;
2595 :
2596 6 : const OGRErr eErr = createFromFgfInternal(
2597 3 : pabyData + nBytesUsed, poSR, &poThisGeom, nBytes - nBytesUsed,
2598 : &nThisGeomSize, nRecLevel + 1);
2599 3 : if (eErr != OGRERR_NONE)
2600 : {
2601 0 : delete poGC;
2602 1 : return eErr;
2603 : }
2604 :
2605 3 : nBytesUsed += nThisGeomSize;
2606 3 : if (poThisGeom != nullptr)
2607 : {
2608 3 : const OGRErr eErr2 = poGC->addGeometryDirectly(poThisGeom);
2609 3 : if (eErr2 != OGRERR_NONE)
2610 : {
2611 1 : delete poGC;
2612 1 : delete poThisGeom;
2613 1 : return eErr2;
2614 : }
2615 : }
2616 : }
2617 :
2618 2 : poGeom = poGC;
2619 2 : if (pnBytesConsumed)
2620 2 : *pnBytesConsumed = nBytesUsed;
2621 : }
2622 :
2623 : /* -------------------------------------------------------------------- */
2624 : /* Currently unsupported geometry. */
2625 : /* */
2626 : /* We need to add 10/11/12/13 curve types in some fashion. */
2627 : /* -------------------------------------------------------------------- */
2628 : else
2629 : {
2630 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
2631 : }
2632 :
2633 : /* -------------------------------------------------------------------- */
2634 : /* Assign spatial reference system. */
2635 : /* -------------------------------------------------------------------- */
2636 11 : if (poGeom != nullptr && poSR)
2637 0 : poGeom->assignSpatialReference(poSR);
2638 11 : *ppoReturn = poGeom;
2639 :
2640 11 : return OGRERR_NONE;
2641 : }
2642 :
2643 : /************************************************************************/
2644 : /* OGR_G_CreateFromFgf() */
2645 : /************************************************************************/
2646 :
2647 : /**
2648 : * \brief Create a geometry object of the appropriate type from its FGF
2649 : * (FDO Geometry Format) binary representation.
2650 : *
2651 : * See OGRGeometryFactory::createFromFgf() */
2652 0 : OGRErr CPL_DLL OGR_G_CreateFromFgf(const void *pabyData,
2653 : OGRSpatialReferenceH hSRS,
2654 : OGRGeometryH *phGeometry, int nBytes,
2655 : int *pnBytesConsumed)
2656 :
2657 : {
2658 0 : return OGRGeometryFactory::createFromFgf(
2659 : pabyData, OGRSpatialReference::FromHandle(hSRS),
2660 0 : reinterpret_cast<OGRGeometry **>(phGeometry), nBytes, pnBytesConsumed);
2661 : }
2662 :
2663 : /************************************************************************/
2664 : /* SplitLineStringAtDateline() */
2665 : /************************************************************************/
2666 :
2667 8 : static void SplitLineStringAtDateline(OGRGeometryCollection *poMulti,
2668 : const OGRLineString *poLS,
2669 : double dfDateLineOffset, double dfXOffset)
2670 : {
2671 8 : const double dfLeftBorderX = 180 - dfDateLineOffset;
2672 8 : const double dfRightBorderX = -180 + dfDateLineOffset;
2673 8 : const double dfDiffSpace = 360 - dfDateLineOffset;
2674 :
2675 8 : const bool bIs3D = poLS->getCoordinateDimension() == 3;
2676 8 : OGRLineString *poNewLS = new OGRLineString();
2677 8 : poMulti->addGeometryDirectly(poNewLS);
2678 35 : for (int i = 0; i < poLS->getNumPoints(); i++)
2679 : {
2680 27 : const double dfX = poLS->getX(i) + dfXOffset;
2681 27 : if (i > 0 && fabs(dfX - (poLS->getX(i - 1) + dfXOffset)) > dfDiffSpace)
2682 : {
2683 9 : double dfX1 = poLS->getX(i - 1) + dfXOffset;
2684 9 : double dfY1 = poLS->getY(i - 1);
2685 9 : double dfZ1 = poLS->getY(i - 1);
2686 9 : double dfX2 = poLS->getX(i) + dfXOffset;
2687 9 : double dfY2 = poLS->getY(i);
2688 9 : double dfZ2 = poLS->getY(i);
2689 :
2690 8 : if (dfX1 > -180 && dfX1 < dfRightBorderX && dfX2 == 180 &&
2691 0 : i + 1 < poLS->getNumPoints() &&
2692 17 : poLS->getX(i + 1) + dfXOffset > -180 &&
2693 0 : poLS->getX(i + 1) + dfXOffset < dfRightBorderX)
2694 : {
2695 0 : if (bIs3D)
2696 0 : poNewLS->addPoint(-180, poLS->getY(i), poLS->getZ(i));
2697 : else
2698 0 : poNewLS->addPoint(-180, poLS->getY(i));
2699 :
2700 0 : i++;
2701 :
2702 0 : if (bIs3D)
2703 0 : poNewLS->addPoint(poLS->getX(i) + dfXOffset, poLS->getY(i),
2704 : poLS->getZ(i));
2705 : else
2706 0 : poNewLS->addPoint(poLS->getX(i) + dfXOffset, poLS->getY(i));
2707 0 : continue;
2708 : }
2709 4 : else if (dfX1 > dfLeftBorderX && dfX1 < 180 && dfX2 == -180 &&
2710 0 : i + 1 < poLS->getNumPoints() &&
2711 13 : poLS->getX(i + 1) + dfXOffset > dfLeftBorderX &&
2712 0 : poLS->getX(i + 1) + dfXOffset < 180)
2713 : {
2714 0 : if (bIs3D)
2715 0 : poNewLS->addPoint(180, poLS->getY(i), poLS->getZ(i));
2716 : else
2717 0 : poNewLS->addPoint(180, poLS->getY(i));
2718 :
2719 0 : i++;
2720 :
2721 0 : if (bIs3D)
2722 0 : poNewLS->addPoint(poLS->getX(i) + dfXOffset, poLS->getY(i),
2723 : poLS->getZ(i));
2724 : else
2725 0 : poNewLS->addPoint(poLS->getX(i) + dfXOffset, poLS->getY(i));
2726 0 : continue;
2727 : }
2728 :
2729 9 : if (dfX1 < dfRightBorderX && dfX2 > dfLeftBorderX)
2730 : {
2731 5 : std::swap(dfX1, dfX2);
2732 5 : std::swap(dfY1, dfY2);
2733 5 : std::swap(dfZ1, dfZ2);
2734 : }
2735 9 : if (dfX1 > dfLeftBorderX && dfX2 < dfRightBorderX)
2736 9 : dfX2 += 360;
2737 :
2738 9 : if (dfX1 <= 180 && dfX2 >= 180 && dfX1 < dfX2)
2739 : {
2740 9 : const double dfRatio = (180 - dfX1) / (dfX2 - dfX1);
2741 9 : const double dfY = dfRatio * dfY2 + (1 - dfRatio) * dfY1;
2742 9 : const double dfZ = dfRatio * dfZ2 + (1 - dfRatio) * dfZ1;
2743 : double dfNewX =
2744 9 : poLS->getX(i - 1) + dfXOffset > dfLeftBorderX ? 180 : -180;
2745 18 : if (poNewLS->getNumPoints() == 0 ||
2746 18 : poNewLS->getX(poNewLS->getNumPoints() - 1) != dfNewX ||
2747 2 : poNewLS->getY(poNewLS->getNumPoints() - 1) != dfY)
2748 : {
2749 7 : if (bIs3D)
2750 0 : poNewLS->addPoint(dfNewX, dfY, dfZ);
2751 : else
2752 7 : poNewLS->addPoint(dfNewX, dfY);
2753 : }
2754 9 : poNewLS = new OGRLineString();
2755 9 : if (bIs3D)
2756 0 : poNewLS->addPoint(
2757 0 : poLS->getX(i - 1) + dfXOffset > dfLeftBorderX ? -180
2758 : : 180,
2759 : dfY, dfZ);
2760 : else
2761 9 : poNewLS->addPoint(
2762 9 : poLS->getX(i - 1) + dfXOffset > dfLeftBorderX ? -180
2763 : : 180,
2764 : dfY);
2765 9 : poMulti->addGeometryDirectly(poNewLS);
2766 : }
2767 : else
2768 : {
2769 0 : poNewLS = new OGRLineString();
2770 0 : poMulti->addGeometryDirectly(poNewLS);
2771 : }
2772 : }
2773 27 : if (bIs3D)
2774 0 : poNewLS->addPoint(dfX, poLS->getY(i), poLS->getZ(i));
2775 : else
2776 27 : poNewLS->addPoint(dfX, poLS->getY(i));
2777 : }
2778 8 : }
2779 :
2780 : /************************************************************************/
2781 : /* FixPolygonCoordinatesAtDateLine() */
2782 : /************************************************************************/
2783 :
2784 : #ifdef HAVE_GEOS
2785 5 : static void FixPolygonCoordinatesAtDateLine(OGRPolygon *poPoly,
2786 : double dfDateLineOffset)
2787 : {
2788 5 : const double dfLeftBorderX = 180 - dfDateLineOffset;
2789 5 : const double dfRightBorderX = -180 + dfDateLineOffset;
2790 5 : const double dfDiffSpace = 360 - dfDateLineOffset;
2791 :
2792 10 : for (int iPart = 0; iPart < 1 + poPoly->getNumInteriorRings(); iPart++)
2793 : {
2794 5 : OGRLineString *poLS = (iPart == 0) ? poPoly->getExteriorRing()
2795 5 : : poPoly->getInteriorRing(iPart - 1);
2796 5 : bool bGoEast = false;
2797 5 : const bool bIs3D = poLS->getCoordinateDimension() == 3;
2798 41 : for (int i = 1; i < poLS->getNumPoints(); i++)
2799 : {
2800 36 : double dfX = poLS->getX(i);
2801 36 : const double dfPrevX = poLS->getX(i - 1);
2802 36 : const double dfDiffLong = fabs(dfX - dfPrevX);
2803 36 : if (dfDiffLong > dfDiffSpace)
2804 : {
2805 21 : if ((dfPrevX > dfLeftBorderX && dfX < dfRightBorderX) ||
2806 6 : (dfX < 0 && bGoEast))
2807 : {
2808 18 : dfX += 360;
2809 18 : bGoEast = true;
2810 18 : if (bIs3D)
2811 0 : poLS->setPoint(i, dfX, poLS->getY(i), poLS->getZ(i));
2812 : else
2813 18 : poLS->setPoint(i, dfX, poLS->getY(i));
2814 : }
2815 3 : else if (dfPrevX < dfRightBorderX && dfX > dfLeftBorderX)
2816 : {
2817 10 : for (int j = i - 1; j >= 0; j--)
2818 : {
2819 7 : dfX = poLS->getX(j);
2820 7 : if (dfX < 0)
2821 : {
2822 7 : if (bIs3D)
2823 0 : poLS->setPoint(j, dfX + 360, poLS->getY(j),
2824 : poLS->getZ(j));
2825 : else
2826 7 : poLS->setPoint(j, dfX + 360, poLS->getY(j));
2827 : }
2828 : }
2829 3 : bGoEast = false;
2830 : }
2831 : else
2832 : {
2833 0 : bGoEast = false;
2834 : }
2835 : }
2836 : }
2837 : }
2838 5 : }
2839 : #endif
2840 :
2841 : /************************************************************************/
2842 : /* AddOffsetToLon() */
2843 : /************************************************************************/
2844 :
2845 17 : static void AddOffsetToLon(OGRGeometry *poGeom, double dfOffset)
2846 : {
2847 17 : switch (wkbFlatten(poGeom->getGeometryType()))
2848 : {
2849 7 : case wkbPolygon:
2850 : {
2851 14 : for (auto poSubGeom : *(poGeom->toPolygon()))
2852 : {
2853 7 : AddOffsetToLon(poSubGeom, dfOffset);
2854 : }
2855 :
2856 7 : break;
2857 : }
2858 :
2859 0 : case wkbMultiLineString:
2860 : case wkbMultiPolygon:
2861 : case wkbGeometryCollection:
2862 : {
2863 0 : for (auto poSubGeom : *(poGeom->toGeometryCollection()))
2864 : {
2865 0 : AddOffsetToLon(poSubGeom, dfOffset);
2866 : }
2867 :
2868 0 : break;
2869 : }
2870 :
2871 10 : case wkbLineString:
2872 : {
2873 10 : OGRLineString *poLineString = poGeom->toLineString();
2874 10 : const int nPointCount = poLineString->getNumPoints();
2875 10 : const int nCoordDim = poLineString->getCoordinateDimension();
2876 63 : for (int iPoint = 0; iPoint < nPointCount; iPoint++)
2877 : {
2878 53 : if (nCoordDim == 2)
2879 106 : poLineString->setPoint(
2880 53 : iPoint, poLineString->getX(iPoint) + dfOffset,
2881 : poLineString->getY(iPoint));
2882 : else
2883 0 : poLineString->setPoint(
2884 0 : iPoint, poLineString->getX(iPoint) + dfOffset,
2885 : poLineString->getY(iPoint), poLineString->getZ(iPoint));
2886 : }
2887 10 : break;
2888 : }
2889 :
2890 0 : default:
2891 0 : break;
2892 : }
2893 17 : }
2894 :
2895 : /************************************************************************/
2896 : /* AddSimpleGeomToMulti() */
2897 : /************************************************************************/
2898 :
2899 : #ifdef HAVE_GEOS
2900 12 : static void AddSimpleGeomToMulti(OGRGeometryCollection *poMulti,
2901 : const OGRGeometry *poGeom)
2902 : {
2903 12 : switch (wkbFlatten(poGeom->getGeometryType()))
2904 : {
2905 12 : case wkbPolygon:
2906 : case wkbLineString:
2907 12 : poMulti->addGeometry(poGeom);
2908 12 : break;
2909 :
2910 0 : case wkbMultiLineString:
2911 : case wkbMultiPolygon:
2912 : case wkbGeometryCollection:
2913 : {
2914 0 : for (const auto poSubGeom : *(poGeom->toGeometryCollection()))
2915 : {
2916 0 : AddSimpleGeomToMulti(poMulti, poSubGeom);
2917 : }
2918 0 : break;
2919 : }
2920 :
2921 0 : default:
2922 0 : break;
2923 : }
2924 12 : }
2925 : #endif // #ifdef HAVE_GEOS
2926 :
2927 : /************************************************************************/
2928 : /* WrapPointDateLine() */
2929 : /************************************************************************/
2930 :
2931 14 : static void WrapPointDateLine(OGRPoint *poPoint)
2932 : {
2933 14 : if (poPoint->getX() > 180)
2934 : {
2935 2 : poPoint->setX(fmod(poPoint->getX() + 180, 360) - 180);
2936 : }
2937 12 : else if (poPoint->getX() < -180)
2938 : {
2939 3 : poPoint->setX(-(fmod(-poPoint->getX() + 180, 360) - 180));
2940 : }
2941 14 : }
2942 :
2943 : /************************************************************************/
2944 : /* CutGeometryOnDateLineAndAddToMulti() */
2945 : /************************************************************************/
2946 :
2947 69 : static void CutGeometryOnDateLineAndAddToMulti(OGRGeometryCollection *poMulti,
2948 : const OGRGeometry *poGeom,
2949 : double dfDateLineOffset)
2950 : {
2951 69 : const OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
2952 69 : switch (eGeomType)
2953 : {
2954 1 : case wkbPoint:
2955 : {
2956 1 : auto poPoint = poGeom->toPoint()->clone();
2957 1 : WrapPointDateLine(poPoint);
2958 1 : poMulti->addGeometryDirectly(poPoint);
2959 1 : break;
2960 : }
2961 :
2962 54 : case wkbPolygon:
2963 : case wkbLineString:
2964 : {
2965 54 : bool bSplitLineStringAtDateline = false;
2966 54 : OGREnvelope oEnvelope;
2967 :
2968 54 : poGeom->getEnvelope(&oEnvelope);
2969 54 : const bool bAroundMinus180 = (oEnvelope.MinX < -180.0);
2970 :
2971 : // Naive heuristics... Place to improve.
2972 : #ifdef HAVE_GEOS
2973 54 : std::unique_ptr<OGRGeometry> poDupGeom;
2974 54 : bool bWrapDateline = false;
2975 : #endif
2976 :
2977 54 : const double dfLeftBorderX = 180 - dfDateLineOffset;
2978 54 : const double dfRightBorderX = -180 + dfDateLineOffset;
2979 54 : const double dfDiffSpace = 360 - dfDateLineOffset;
2980 :
2981 54 : const double dfXOffset = (bAroundMinus180) ? 360.0 : 0.0;
2982 54 : if (oEnvelope.MinX < -180 || oEnvelope.MaxX > 180 ||
2983 52 : (oEnvelope.MinX + dfXOffset > dfLeftBorderX &&
2984 12 : oEnvelope.MaxX + dfXOffset > 180))
2985 : {
2986 : #ifndef HAVE_GEOS
2987 : CPLError(CE_Failure, CPLE_NotSupported,
2988 : "GEOS support not enabled.");
2989 : #else
2990 2 : bWrapDateline = true;
2991 : #endif
2992 : }
2993 : else
2994 : {
2995 : auto poLS = eGeomType == wkbPolygon
2996 52 : ? poGeom->toPolygon()->getExteriorRing()
2997 14 : : poGeom->toLineString();
2998 52 : if (poLS)
2999 : {
3000 52 : double dfMaxSmallDiffLong = 0;
3001 52 : bool bHasBigDiff = false;
3002 : // Detect big gaps in longitude.
3003 294 : for (int i = 1; i < poLS->getNumPoints(); i++)
3004 : {
3005 242 : const double dfPrevX = poLS->getX(i - 1) + dfXOffset;
3006 242 : const double dfX = poLS->getX(i) + dfXOffset;
3007 242 : const double dfDiffLong = fabs(dfX - dfPrevX);
3008 :
3009 242 : if (dfDiffLong > dfDiffSpace &&
3010 11 : ((dfX > dfLeftBorderX &&
3011 10 : dfPrevX < dfRightBorderX) ||
3012 10 : (dfPrevX > dfLeftBorderX && dfX < dfRightBorderX)))
3013 21 : bHasBigDiff = true;
3014 221 : else if (dfDiffLong > dfMaxSmallDiffLong)
3015 55 : dfMaxSmallDiffLong = dfDiffLong;
3016 : }
3017 52 : if (bHasBigDiff && dfMaxSmallDiffLong < dfDateLineOffset)
3018 : {
3019 13 : if (eGeomType == wkbLineString)
3020 8 : bSplitLineStringAtDateline = true;
3021 : else
3022 : {
3023 : #ifndef HAVE_GEOS
3024 : CPLError(CE_Failure, CPLE_NotSupported,
3025 : "GEOS support not enabled.");
3026 : #else
3027 5 : poDupGeom.reset(poGeom->clone());
3028 5 : FixPolygonCoordinatesAtDateLine(
3029 : poDupGeom->toPolygon(), dfDateLineOffset);
3030 :
3031 5 : OGREnvelope sEnvelope;
3032 5 : poDupGeom->getEnvelope(&sEnvelope);
3033 5 : bWrapDateline = sEnvelope.MinX != sEnvelope.MaxX;
3034 : #endif
3035 : }
3036 : }
3037 : }
3038 : }
3039 :
3040 54 : if (bSplitLineStringAtDateline)
3041 : {
3042 8 : SplitLineStringAtDateline(poMulti, poGeom->toLineString(),
3043 : dfDateLineOffset,
3044 : (bAroundMinus180) ? 360.0 : 0.0);
3045 : }
3046 : #ifdef HAVE_GEOS
3047 46 : else if (bWrapDateline)
3048 : {
3049 : const OGRGeometry *poWorkGeom =
3050 6 : poDupGeom ? poDupGeom.get() : poGeom;
3051 6 : OGRGeometry *poRectangle1 = nullptr;
3052 6 : OGRGeometry *poRectangle2 = nullptr;
3053 6 : const char *pszWKT1 =
3054 : !bAroundMinus180
3055 6 : ? "POLYGON((-180 90,180 90,180 -90,-180 -90,-180 90))"
3056 : : "POLYGON((180 90,-180 90,-180 -90,180 -90,180 90))";
3057 6 : const char *pszWKT2 =
3058 : !bAroundMinus180
3059 6 : ? "POLYGON((180 90,360 90,360 -90,180 -90,180 90))"
3060 : : "POLYGON((-180 90,-360 90,-360 -90,-180 -90,-180 "
3061 : "90))";
3062 6 : OGRGeometryFactory::createFromWkt(pszWKT1, nullptr,
3063 : &poRectangle1);
3064 6 : OGRGeometryFactory::createFromWkt(pszWKT2, nullptr,
3065 : &poRectangle2);
3066 : auto poGeom1 = std::unique_ptr<OGRGeometry>(
3067 12 : poWorkGeom->Intersection(poRectangle1));
3068 : auto poGeom2 = std::unique_ptr<OGRGeometry>(
3069 12 : poWorkGeom->Intersection(poRectangle2));
3070 6 : delete poRectangle1;
3071 6 : delete poRectangle2;
3072 :
3073 6 : if (poGeom1 != nullptr && poGeom2 != nullptr)
3074 : {
3075 6 : AddSimpleGeomToMulti(poMulti, poGeom1.get());
3076 6 : AddOffsetToLon(poGeom2.get(),
3077 : !bAroundMinus180 ? -360.0 : 360.0);
3078 6 : AddSimpleGeomToMulti(poMulti, poGeom2.get());
3079 : }
3080 : else
3081 : {
3082 0 : AddSimpleGeomToMulti(poMulti, poGeom);
3083 : }
3084 : }
3085 : #endif
3086 : else
3087 : {
3088 40 : poMulti->addGeometry(poGeom);
3089 : }
3090 54 : break;
3091 : }
3092 :
3093 14 : case wkbMultiLineString:
3094 : case wkbMultiPolygon:
3095 : case wkbGeometryCollection:
3096 : {
3097 45 : for (const auto poSubGeom : *(poGeom->toGeometryCollection()))
3098 : {
3099 31 : CutGeometryOnDateLineAndAddToMulti(poMulti, poSubGeom,
3100 : dfDateLineOffset);
3101 : }
3102 14 : break;
3103 : }
3104 :
3105 0 : default:
3106 0 : break;
3107 : }
3108 69 : }
3109 :
3110 : #ifdef HAVE_GEOS
3111 :
3112 : /************************************************************************/
3113 : /* RemovePoint() */
3114 : /************************************************************************/
3115 :
3116 9 : static void RemovePoint(OGRGeometry *poGeom, OGRPoint *poPoint)
3117 : {
3118 9 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3119 9 : switch (eType)
3120 : {
3121 4 : case wkbLineString:
3122 : {
3123 4 : OGRLineString *poLS = poGeom->toLineString();
3124 4 : const bool bIs3D = (poLS->getCoordinateDimension() == 3);
3125 4 : int j = 0;
3126 32 : for (int i = 0; i < poLS->getNumPoints(); i++)
3127 : {
3128 30 : if (poLS->getX(i) != poPoint->getX() ||
3129 2 : poLS->getY(i) != poPoint->getY())
3130 : {
3131 26 : if (i > j)
3132 : {
3133 4 : if (bIs3D)
3134 : {
3135 0 : poLS->setPoint(j, poLS->getX(i), poLS->getY(i),
3136 : poLS->getZ(i));
3137 : }
3138 : else
3139 : {
3140 4 : poLS->setPoint(j, poLS->getX(i), poLS->getY(i));
3141 : }
3142 : }
3143 26 : j++;
3144 : }
3145 : }
3146 4 : poLS->setNumPoints(j);
3147 4 : break;
3148 : }
3149 :
3150 4 : case wkbPolygon:
3151 : {
3152 4 : OGRPolygon *poPoly = poGeom->toPolygon();
3153 4 : if (poPoly->getExteriorRing() != nullptr)
3154 : {
3155 4 : RemovePoint(poPoly->getExteriorRing(), poPoint);
3156 4 : for (int i = 0; i < poPoly->getNumInteriorRings(); ++i)
3157 : {
3158 0 : RemovePoint(poPoly->getInteriorRing(i), poPoint);
3159 : }
3160 : }
3161 4 : break;
3162 : }
3163 :
3164 1 : case wkbMultiLineString:
3165 : case wkbMultiPolygon:
3166 : case wkbGeometryCollection:
3167 : {
3168 1 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
3169 3 : for (int i = 0; i < poGC->getNumGeometries(); ++i)
3170 : {
3171 2 : RemovePoint(poGC->getGeometryRef(i), poPoint);
3172 : }
3173 1 : break;
3174 : }
3175 :
3176 0 : default:
3177 0 : break;
3178 : }
3179 9 : }
3180 :
3181 : /************************************************************************/
3182 : /* GetDist() */
3183 : /************************************************************************/
3184 :
3185 30172 : static double GetDist(double dfDeltaX, double dfDeltaY)
3186 : {
3187 30172 : return sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
3188 : }
3189 :
3190 : /************************************************************************/
3191 : /* AlterPole() */
3192 : /* */
3193 : /* Replace and point at the pole by points really close to the pole, */
3194 : /* but on the previous and later segments. */
3195 : /************************************************************************/
3196 :
3197 5 : static void AlterPole(OGRGeometry *poGeom, OGRPoint *poPole,
3198 : bool bIsRing = false)
3199 : {
3200 5 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3201 5 : switch (eType)
3202 : {
3203 2 : case wkbLineString:
3204 : {
3205 2 : if (!bIsRing)
3206 0 : return;
3207 2 : OGRLineString *poLS = poGeom->toLineString();
3208 2 : const int nNumPoints = poLS->getNumPoints();
3209 2 : if (nNumPoints >= 4)
3210 : {
3211 2 : const bool bIs3D = (poLS->getCoordinateDimension() == 3);
3212 4 : std::vector<OGRRawPoint> aoPoints;
3213 4 : std::vector<double> adfZ;
3214 2 : bool bMustClose = false;
3215 10 : for (int i = 0; i < nNumPoints; i++)
3216 : {
3217 8 : const double dfX = poLS->getX(i);
3218 8 : const double dfY = poLS->getY(i);
3219 8 : if (dfX == poPole->getX() && dfY == poPole->getY())
3220 : {
3221 : // Replace the pole by points really close to it
3222 2 : if (i == 0)
3223 0 : bMustClose = true;
3224 2 : if (i == nNumPoints - 1)
3225 0 : continue;
3226 2 : const int iBefore = i > 0 ? i - 1 : nNumPoints - 2;
3227 2 : double dfXBefore = poLS->getX(iBefore);
3228 2 : double dfYBefore = poLS->getY(iBefore);
3229 : double dfNorm =
3230 2 : GetDist(dfXBefore - dfX, dfYBefore - dfY);
3231 2 : double dfXInterp =
3232 2 : dfX + (dfXBefore - dfX) / dfNorm * 1.0e-7;
3233 2 : double dfYInterp =
3234 2 : dfY + (dfYBefore - dfY) / dfNorm * 1.0e-7;
3235 2 : OGRRawPoint oPoint;
3236 2 : oPoint.x = dfXInterp;
3237 2 : oPoint.y = dfYInterp;
3238 2 : aoPoints.push_back(oPoint);
3239 2 : adfZ.push_back(poLS->getZ(i));
3240 :
3241 2 : const int iAfter = i + 1;
3242 2 : double dfXAfter = poLS->getX(iAfter);
3243 2 : double dfYAfter = poLS->getY(iAfter);
3244 2 : dfNorm = GetDist(dfXAfter - dfX, dfYAfter - dfY);
3245 2 : dfXInterp = dfX + (dfXAfter - dfX) / dfNorm * 1e-7;
3246 2 : dfYInterp = dfY + (dfYAfter - dfY) / dfNorm * 1e-7;
3247 2 : oPoint.x = dfXInterp;
3248 2 : oPoint.y = dfYInterp;
3249 2 : aoPoints.push_back(oPoint);
3250 2 : adfZ.push_back(poLS->getZ(i));
3251 : }
3252 : else
3253 : {
3254 6 : OGRRawPoint oPoint;
3255 6 : oPoint.x = dfX;
3256 6 : oPoint.y = dfY;
3257 6 : aoPoints.push_back(oPoint);
3258 6 : adfZ.push_back(poLS->getZ(i));
3259 : }
3260 : }
3261 2 : if (bMustClose)
3262 : {
3263 0 : aoPoints.push_back(aoPoints[0]);
3264 0 : adfZ.push_back(adfZ[0]);
3265 : }
3266 :
3267 4 : poLS->setPoints(static_cast<int>(aoPoints.size()),
3268 2 : &(aoPoints[0]), bIs3D ? &adfZ[0] : nullptr);
3269 : }
3270 2 : break;
3271 : }
3272 :
3273 2 : case wkbPolygon:
3274 : {
3275 2 : OGRPolygon *poPoly = poGeom->toPolygon();
3276 2 : if (poPoly->getExteriorRing() != nullptr)
3277 : {
3278 2 : AlterPole(poPoly->getExteriorRing(), poPole, true);
3279 2 : for (int i = 0; i < poPoly->getNumInteriorRings(); ++i)
3280 : {
3281 0 : AlterPole(poPoly->getInteriorRing(i), poPole, true);
3282 : }
3283 : }
3284 2 : break;
3285 : }
3286 :
3287 1 : case wkbMultiLineString:
3288 : case wkbMultiPolygon:
3289 : case wkbGeometryCollection:
3290 : {
3291 1 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
3292 2 : for (int i = 0; i < poGC->getNumGeometries(); ++i)
3293 : {
3294 1 : AlterPole(poGC->getGeometryRef(i), poPole);
3295 : }
3296 1 : break;
3297 : }
3298 :
3299 0 : default:
3300 0 : break;
3301 : }
3302 : }
3303 :
3304 : /************************************************************************/
3305 : /* IsPolarToWGS84() */
3306 : /* */
3307 : /* Returns true if poCT transforms from a projection that includes one */
3308 : /* of the pole in a continuous way. */
3309 : /************************************************************************/
3310 :
3311 30 : static bool IsPolarToWGS84(OGRCoordinateTransformation *poCT,
3312 : OGRCoordinateTransformation *poRevCT,
3313 : bool &bIsNorthPolarOut)
3314 : {
3315 30 : bool bIsNorthPolar = false;
3316 30 : bool bIsSouthPolar = false;
3317 30 : double x = 0.0;
3318 30 : double y = 90.0;
3319 :
3320 30 : const bool bBackupEmitErrors = poCT->GetEmitErrors();
3321 30 : poRevCT->SetEmitErrors(false);
3322 30 : poCT->SetEmitErrors(false);
3323 :
3324 30 : if (poRevCT->Transform(1, &x, &y) &&
3325 : // Surprisingly, pole south projects correctly back &
3326 : // forth for antarctic polar stereographic. Therefore, check that
3327 : // the projected value is not too big.
3328 30 : fabs(x) < 1e10 && fabs(y) < 1e10)
3329 : {
3330 27 : double x_tab[] = {x, x - 1e5, x + 1e5};
3331 27 : double y_tab[] = {y, y - 1e5, y + 1e5};
3332 27 : if (poCT->Transform(3, x_tab, y_tab) &&
3333 27 : fabs(y_tab[0] - (90.0)) < 1e-10 &&
3334 80 : fabs(x_tab[2] - x_tab[1]) > 170 &&
3335 26 : fabs(y_tab[2] - y_tab[1]) < 1e-10)
3336 : {
3337 26 : bIsNorthPolar = true;
3338 : }
3339 : }
3340 :
3341 30 : x = 0.0;
3342 30 : y = -90.0;
3343 30 : if (poRevCT->Transform(1, &x, &y) && fabs(x) < 1e10 && fabs(y) < 1e10)
3344 : {
3345 24 : double x_tab[] = {x, x - 1e5, x + 1e5};
3346 24 : double y_tab[] = {y, y - 1e5, y + 1e5};
3347 24 : if (poCT->Transform(3, x_tab, y_tab) &&
3348 24 : fabs(y_tab[0] - (-90.0)) < 1e-10 &&
3349 71 : fabs(x_tab[2] - x_tab[1]) > 170 &&
3350 23 : fabs(y_tab[2] - y_tab[1]) < 1e-10)
3351 : {
3352 23 : bIsSouthPolar = true;
3353 : }
3354 : }
3355 :
3356 30 : poCT->SetEmitErrors(bBackupEmitErrors);
3357 :
3358 30 : if (bIsNorthPolar && bIsSouthPolar)
3359 : {
3360 22 : bIsNorthPolar = false;
3361 22 : bIsSouthPolar = false;
3362 : }
3363 :
3364 30 : bIsNorthPolarOut = bIsNorthPolar;
3365 30 : return bIsNorthPolar || bIsSouthPolar;
3366 : }
3367 :
3368 : /************************************************************************/
3369 : /* TransformBeforePolarToWGS84() */
3370 : /* */
3371 : /* Transform the geometry (by intersection), so as to cut each geometry */
3372 : /* that crosses the pole, in 2 parts. Do also tricks for geometries */
3373 : /* that just touch the pole. */
3374 : /************************************************************************/
3375 :
3376 : static OGRGeometry *
3377 6 : TransformBeforePolarToWGS84(OGRCoordinateTransformation *poRevCT,
3378 : bool bIsNorthPolar, OGRGeometry *poDstGeom,
3379 : bool &bNeedPostCorrectionOut)
3380 : {
3381 6 : const int nSign = (bIsNorthPolar) ? 1 : -1;
3382 :
3383 : // Does the geometry fully contains the pole ? */
3384 6 : double dfXPole = 0.0;
3385 6 : double dfYPole = nSign * 90.0;
3386 6 : poRevCT->Transform(1, &dfXPole, &dfYPole);
3387 12 : OGRPoint oPole(dfXPole, dfYPole);
3388 6 : const bool bContainsPole = CPL_TO_BOOL(poDstGeom->Contains(&oPole));
3389 :
3390 6 : const double EPS = 1e-9;
3391 :
3392 : // Does the geometry touches the pole and intersects the antimeridian ?
3393 6 : double dfNearPoleAntiMeridianX = 180.0;
3394 6 : double dfNearPoleAntiMeridianY = nSign * (90.0 - EPS);
3395 6 : poRevCT->Transform(1, &dfNearPoleAntiMeridianX, &dfNearPoleAntiMeridianY);
3396 : OGRPoint oNearPoleAntimeridian(dfNearPoleAntiMeridianX,
3397 12 : dfNearPoleAntiMeridianY);
3398 : const bool bContainsNearPoleAntimeridian =
3399 6 : CPL_TO_BOOL(poDstGeom->Contains(&oNearPoleAntimeridian));
3400 :
3401 : // Does the geometry touches the pole (but not intersect the antimeridian) ?
3402 10 : const bool bRegularTouchesPole = !bContainsPole &&
3403 9 : !bContainsNearPoleAntimeridian &&
3404 3 : CPL_TO_BOOL(poDstGeom->Touches(&oPole));
3405 :
3406 : // Create a polygon of nearly a full hemisphere, but excluding the anti
3407 : // meridian and the pole.
3408 6 : OGRPolygon oCutter;
3409 6 : OGRLinearRing *poRing = new OGRLinearRing();
3410 6 : poRing->addPoint(180.0 - EPS, 0);
3411 6 : poRing->addPoint(180.0 - EPS, nSign * (90.0 - EPS));
3412 : // If the geometry doesn't contain the pole, then we add it to the cutter
3413 : // geometry, but will later remove it completely (geometry touching the
3414 : // pole but intersecting the antimeridian), or will replace it by 2
3415 : // close points (geometry touching the pole without intersecting the
3416 : // antimeridian)
3417 6 : if (!bContainsPole)
3418 4 : poRing->addPoint(180.0, nSign * 90);
3419 6 : poRing->addPoint(-180.0 + EPS, nSign * (90.0 - EPS));
3420 6 : poRing->addPoint(-180.0 + EPS, 0);
3421 6 : poRing->addPoint(180.0 - EPS, 0);
3422 6 : oCutter.addRingDirectly(poRing);
3423 :
3424 6 : if (oCutter.transform(poRevCT) == OGRERR_NONE &&
3425 : // Check that longitudes +/- 180 are continuous
3426 : // in the polar projection
3427 10 : fabs(poRing->getX(0) - poRing->getX(poRing->getNumPoints() - 2)) < 1 &&
3428 4 : (bContainsPole || bContainsNearPoleAntimeridian || bRegularTouchesPole))
3429 : {
3430 5 : if (bContainsPole || bContainsNearPoleAntimeridian)
3431 : {
3432 3 : OGRGeometry *poNewGeom = poDstGeom->Difference(&oCutter);
3433 3 : if (poNewGeom)
3434 : {
3435 3 : if (bContainsNearPoleAntimeridian)
3436 3 : RemovePoint(poNewGeom, &oPole);
3437 3 : delete poDstGeom;
3438 3 : poDstGeom = poNewGeom;
3439 : }
3440 : }
3441 :
3442 5 : if (bRegularTouchesPole)
3443 : {
3444 2 : AlterPole(poDstGeom, &oPole);
3445 : }
3446 :
3447 5 : bNeedPostCorrectionOut = true;
3448 : }
3449 12 : return poDstGeom;
3450 : }
3451 :
3452 : /************************************************************************/
3453 : /* IsAntimeridianProjToWGS84() */
3454 : /* */
3455 : /* Returns true if poCT transforms from a projection that includes the */
3456 : /* antimeridian in a continuous way. */
3457 : /************************************************************************/
3458 :
3459 10058 : static bool IsAntimeridianProjToWGS84(OGRCoordinateTransformation *poCT,
3460 : OGRCoordinateTransformation *poRevCT,
3461 : OGRGeometry *poDstGeometry)
3462 : {
3463 10058 : const bool bBackupEmitErrors = poCT->GetEmitErrors();
3464 10058 : poRevCT->SetEmitErrors(false);
3465 10058 : poCT->SetEmitErrors(false);
3466 :
3467 : // Find a reasonable latitude for the geometry
3468 10058 : OGREnvelope sEnvelope;
3469 10058 : poDstGeometry->getEnvelope(&sEnvelope);
3470 20116 : OGRPoint pMean(sEnvelope.MinX, (sEnvelope.MinY + sEnvelope.MaxY) / 2);
3471 10058 : if (pMean.transform(poCT) != OGRERR_NONE)
3472 : {
3473 0 : poCT->SetEmitErrors(bBackupEmitErrors);
3474 0 : return false;
3475 : }
3476 10058 : const double dfMeanLat = pMean.getY();
3477 :
3478 : // Check that close points on each side of the antimeridian in (long, lat)
3479 : // project to close points in the source projection, and check that they
3480 : // roundtrip correctly.
3481 10058 : const double EPS = 1.0e-8;
3482 10058 : double x1 = 180 - EPS;
3483 10058 : double y1 = dfMeanLat;
3484 10058 : double x2 = -180 + EPS;
3485 10058 : double y2 = dfMeanLat;
3486 30174 : if (!poRevCT->Transform(1, &x1, &y1) || !poRevCT->Transform(1, &x2, &y2) ||
3487 20113 : GetDist(x2 - x1, y2 - y1) > 1 || !poCT->Transform(1, &x1, &y1) ||
3488 20110 : !poCT->Transform(1, &x2, &y2) ||
3489 30171 : GetDist(x1 - (180 - EPS), y1 - dfMeanLat) > 2 * EPS ||
3490 10055 : GetDist(x2 - (-180 + EPS), y2 - dfMeanLat) > 2 * EPS)
3491 : {
3492 3 : poCT->SetEmitErrors(bBackupEmitErrors);
3493 3 : return false;
3494 : }
3495 :
3496 10055 : poCT->SetEmitErrors(bBackupEmitErrors);
3497 :
3498 10055 : return true;
3499 : }
3500 :
3501 : /************************************************************************/
3502 : /* CollectPointsOnAntimeridian() */
3503 : /* */
3504 : /* Collect points that are the intersection of the lines of the geometry*/
3505 : /* with the antimeridian. */
3506 : /************************************************************************/
3507 :
3508 19 : static void CollectPointsOnAntimeridian(OGRGeometry *poGeom,
3509 : OGRCoordinateTransformation *poCT,
3510 : OGRCoordinateTransformation *poRevCT,
3511 : std::vector<OGRRawPoint> &aoPoints)
3512 : {
3513 19 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3514 19 : switch (eType)
3515 : {
3516 9 : case wkbLineString:
3517 : {
3518 9 : OGRLineString *poLS = poGeom->toLineString();
3519 9 : const int nNumPoints = poLS->getNumPoints();
3520 40 : for (int i = 0; i < nNumPoints - 1; i++)
3521 : {
3522 31 : const double dfX = poLS->getX(i);
3523 31 : const double dfY = poLS->getY(i);
3524 31 : const double dfX2 = poLS->getX(i + 1);
3525 31 : const double dfY2 = poLS->getY(i + 1);
3526 31 : double dfXTrans = dfX;
3527 31 : double dfYTrans = dfY;
3528 31 : double dfX2Trans = dfX2;
3529 31 : double dfY2Trans = dfY2;
3530 31 : poCT->Transform(1, &dfXTrans, &dfYTrans);
3531 31 : poCT->Transform(1, &dfX2Trans, &dfY2Trans);
3532 : // Are we crossing the antimeridian ? (detecting by inversion of
3533 : // sign of X)
3534 31 : if ((dfX2 - dfX) * (dfX2Trans - dfXTrans) < 0 ||
3535 14 : (dfX == dfX2 && dfX2Trans * dfXTrans < 0 &&
3536 1 : fabs(fabs(dfXTrans) - 180) < 10 &&
3537 1 : fabs(fabs(dfX2Trans) - 180) < 10))
3538 : {
3539 15 : double dfXStart = dfX;
3540 15 : double dfYStart = dfY;
3541 15 : double dfXEnd = dfX2;
3542 15 : double dfYEnd = dfY2;
3543 15 : double dfXStartTrans = dfXTrans;
3544 15 : double dfXEndTrans = dfX2Trans;
3545 15 : int iIter = 0;
3546 15 : const double EPS = 1e-8;
3547 : // Find point of the segment intersecting the antimeridian
3548 : // by dichotomy
3549 397 : for (;
3550 412 : iIter < 50 && (fabs(fabs(dfXStartTrans) - 180) > EPS ||
3551 23 : fabs(fabs(dfXEndTrans) - 180) > EPS);
3552 : ++iIter)
3553 : {
3554 397 : double dfXMid = (dfXStart + dfXEnd) / 2;
3555 397 : double dfYMid = (dfYStart + dfYEnd) / 2;
3556 397 : double dfXMidTrans = dfXMid;
3557 397 : double dfYMidTrans = dfYMid;
3558 397 : poCT->Transform(1, &dfXMidTrans, &dfYMidTrans);
3559 397 : if ((dfXMid - dfXStart) *
3560 397 : (dfXMidTrans - dfXStartTrans) <
3561 203 : 0 ||
3562 22 : (dfXMid == dfXStart &&
3563 22 : dfXMidTrans * dfXStartTrans < 0))
3564 : {
3565 202 : dfXEnd = dfXMid;
3566 202 : dfYEnd = dfYMid;
3567 202 : dfXEndTrans = dfXMidTrans;
3568 : }
3569 : else
3570 : {
3571 195 : dfXStart = dfXMid;
3572 195 : dfYStart = dfYMid;
3573 195 : dfXStartTrans = dfXMidTrans;
3574 : }
3575 : }
3576 15 : if (iIter < 50)
3577 : {
3578 15 : OGRRawPoint oPoint;
3579 15 : oPoint.x = (dfXStart + dfXEnd) / 2;
3580 15 : oPoint.y = (dfYStart + dfYEnd) / 2;
3581 15 : poCT->Transform(1, &(oPoint.x), &(oPoint.y));
3582 15 : oPoint.x = 180.0;
3583 15 : aoPoints.push_back(oPoint);
3584 : }
3585 : }
3586 : }
3587 9 : break;
3588 : }
3589 :
3590 6 : case wkbPolygon:
3591 : {
3592 6 : OGRPolygon *poPoly = poGeom->toPolygon();
3593 6 : if (poPoly->getExteriorRing() != nullptr)
3594 : {
3595 6 : CollectPointsOnAntimeridian(poPoly->getExteriorRing(), poCT,
3596 : poRevCT, aoPoints);
3597 6 : for (int i = 0; i < poPoly->getNumInteriorRings(); ++i)
3598 : {
3599 0 : CollectPointsOnAntimeridian(poPoly->getInteriorRing(i),
3600 : poCT, poRevCT, aoPoints);
3601 : }
3602 : }
3603 6 : break;
3604 : }
3605 :
3606 4 : case wkbMultiLineString:
3607 : case wkbMultiPolygon:
3608 : case wkbGeometryCollection:
3609 : {
3610 4 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
3611 8 : for (int i = 0; i < poGC->getNumGeometries(); ++i)
3612 : {
3613 4 : CollectPointsOnAntimeridian(poGC->getGeometryRef(i), poCT,
3614 : poRevCT, aoPoints);
3615 : }
3616 4 : break;
3617 : }
3618 :
3619 0 : default:
3620 0 : break;
3621 : }
3622 19 : }
3623 :
3624 : /************************************************************************/
3625 : /* SortPointsByAscendingY() */
3626 : /************************************************************************/
3627 :
3628 : struct SortPointsByAscendingY
3629 : {
3630 8 : bool operator()(const OGRRawPoint &a, const OGRRawPoint &b)
3631 : {
3632 8 : return a.y < b.y;
3633 : }
3634 : };
3635 :
3636 : /************************************************************************/
3637 : /* TransformBeforeAntimeridianToWGS84() */
3638 : /* */
3639 : /* Transform the geometry (by intersection), so as to cut each geometry */
3640 : /* that crosses the antimeridian, in 2 parts. */
3641 : /************************************************************************/
3642 :
3643 10055 : static OGRGeometry *TransformBeforeAntimeridianToWGS84(
3644 : OGRCoordinateTransformation *poCT, OGRCoordinateTransformation *poRevCT,
3645 : OGRGeometry *poDstGeom, bool &bNeedPostCorrectionOut)
3646 : {
3647 10055 : OGREnvelope sEnvelope;
3648 10055 : poDstGeom->getEnvelope(&sEnvelope);
3649 20110 : OGRPoint pMean(sEnvelope.MinX, (sEnvelope.MinY + sEnvelope.MaxY) / 2);
3650 10055 : pMean.transform(poCT);
3651 10055 : const double dfMeanLat = pMean.getY();
3652 10055 : pMean.setX(180.0);
3653 10055 : pMean.setY(dfMeanLat);
3654 10055 : pMean.transform(poRevCT);
3655 : // Check if the antimeridian crosses the bbox of our geometry
3656 20078 : if (!(pMean.getX() >= sEnvelope.MinX && pMean.getY() >= sEnvelope.MinY &&
3657 10023 : pMean.getX() <= sEnvelope.MaxX && pMean.getY() <= sEnvelope.MaxY))
3658 : {
3659 10046 : return poDstGeom;
3660 : }
3661 :
3662 : // Collect points that are the intersection of the lines of the geometry
3663 : // with the antimeridian
3664 18 : std::vector<OGRRawPoint> aoPoints;
3665 9 : CollectPointsOnAntimeridian(poDstGeom, poCT, poRevCT, aoPoints);
3666 9 : if (aoPoints.empty())
3667 0 : return poDstGeom;
3668 :
3669 : SortPointsByAscendingY sortFunc;
3670 9 : std::sort(aoPoints.begin(), aoPoints.end(), sortFunc);
3671 :
3672 9 : const double EPS = 1e-9;
3673 :
3674 : // Build a very thin polygon cutting the antimeridian at our points
3675 9 : OGRLinearRing *poLR = new OGRLinearRing;
3676 : {
3677 9 : double x = 180.0 - EPS;
3678 9 : double y = aoPoints[0].y - EPS;
3679 9 : poRevCT->Transform(1, &x, &y);
3680 9 : poLR->addPoint(x, y);
3681 : }
3682 24 : for (const auto &oPoint : aoPoints)
3683 : {
3684 15 : double x = 180.0 - EPS;
3685 15 : double y = oPoint.y;
3686 15 : poRevCT->Transform(1, &x, &y);
3687 15 : poLR->addPoint(x, y);
3688 : }
3689 : {
3690 9 : double x = 180.0 - EPS;
3691 9 : double y = aoPoints.back().y + EPS;
3692 9 : poRevCT->Transform(1, &x, &y);
3693 9 : poLR->addPoint(x, y);
3694 : }
3695 : {
3696 9 : double x = 180.0 + EPS;
3697 9 : double y = aoPoints.back().y + EPS;
3698 9 : poRevCT->Transform(1, &x, &y);
3699 9 : poLR->addPoint(x, y);
3700 : }
3701 24 : for (size_t i = aoPoints.size(); i > 0;)
3702 : {
3703 15 : --i;
3704 15 : const OGRRawPoint &oPoint = aoPoints[i];
3705 15 : double x = 180.0 + EPS;
3706 15 : double y = oPoint.y;
3707 15 : poRevCT->Transform(1, &x, &y);
3708 15 : poLR->addPoint(x, y);
3709 : }
3710 : {
3711 9 : double x = 180.0 + EPS;
3712 9 : double y = aoPoints[0].y - EPS;
3713 9 : poRevCT->Transform(1, &x, &y);
3714 9 : poLR->addPoint(x, y);
3715 : }
3716 9 : poLR->closeRings();
3717 :
3718 9 : OGRPolygon oPolyToCut;
3719 9 : oPolyToCut.addRingDirectly(poLR);
3720 :
3721 : #if DEBUG_VERBOSE
3722 : char *pszWKT = NULL;
3723 : oPolyToCut.exportToWkt(&pszWKT);
3724 : CPLDebug("OGR", "Geometry to cut: %s", pszWKT);
3725 : CPLFree(pszWKT);
3726 : #endif
3727 :
3728 : // Get the geometry without the antimeridian
3729 9 : OGRGeometry *poInter = poDstGeom->Difference(&oPolyToCut);
3730 9 : if (poInter != nullptr)
3731 : {
3732 9 : delete poDstGeom;
3733 9 : poDstGeom = poInter;
3734 9 : bNeedPostCorrectionOut = true;
3735 : }
3736 :
3737 9 : return poDstGeom;
3738 : }
3739 :
3740 : /************************************************************************/
3741 : /* SnapCoordsCloseToLatLongBounds() */
3742 : /* */
3743 : /* This function snaps points really close to the antimerdian or poles */
3744 : /* to their exact longitudes/latitudes. */
3745 : /************************************************************************/
3746 :
3747 53 : static void SnapCoordsCloseToLatLongBounds(OGRGeometry *poGeom)
3748 : {
3749 53 : const OGRwkbGeometryType eType = wkbFlatten(poGeom->getGeometryType());
3750 53 : switch (eType)
3751 : {
3752 24 : case wkbLineString:
3753 : {
3754 24 : OGRLineString *poLS = poGeom->toLineString();
3755 24 : const double EPS = 1e-8;
3756 153 : for (int i = 0; i < poLS->getNumPoints(); i++)
3757 : {
3758 258 : OGRPoint p;
3759 129 : poLS->getPoint(i, &p);
3760 129 : if (fabs(p.getX() - 180.0) < EPS)
3761 : {
3762 34 : p.setX(180.0);
3763 34 : poLS->setPoint(i, &p);
3764 : }
3765 95 : else if (fabs(p.getX() - -180.0) < EPS)
3766 : {
3767 29 : p.setX(-180.0);
3768 29 : poLS->setPoint(i, &p);
3769 : }
3770 :
3771 129 : if (fabs(p.getY() - 90.0) < EPS)
3772 : {
3773 8 : p.setY(90.0);
3774 8 : poLS->setPoint(i, &p);
3775 : }
3776 121 : else if (fabs(p.getY() - -90.0) < EPS)
3777 : {
3778 2 : p.setY(-90.0);
3779 2 : poLS->setPoint(i, &p);
3780 : }
3781 : }
3782 24 : break;
3783 : }
3784 :
3785 18 : case wkbPolygon:
3786 : {
3787 18 : OGRPolygon *poPoly = poGeom->toPolygon();
3788 18 : if (poPoly->getExteriorRing() != nullptr)
3789 : {
3790 18 : SnapCoordsCloseToLatLongBounds(poPoly->getExteriorRing());
3791 18 : for (int i = 0; i < poPoly->getNumInteriorRings(); ++i)
3792 : {
3793 0 : SnapCoordsCloseToLatLongBounds(poPoly->getInteriorRing(i));
3794 : }
3795 : }
3796 18 : break;
3797 : }
3798 :
3799 11 : case wkbMultiLineString:
3800 : case wkbMultiPolygon:
3801 : case wkbGeometryCollection:
3802 : {
3803 11 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
3804 32 : for (int i = 0; i < poGC->getNumGeometries(); ++i)
3805 : {
3806 21 : SnapCoordsCloseToLatLongBounds(poGC->getGeometryRef(i));
3807 : }
3808 11 : break;
3809 : }
3810 :
3811 0 : default:
3812 0 : break;
3813 : }
3814 53 : }
3815 :
3816 : #endif
3817 :
3818 : /************************************************************************/
3819 : /* TransformWithOptionsCache::Private */
3820 : /************************************************************************/
3821 :
3822 : struct OGRGeometryFactory::TransformWithOptionsCache::Private
3823 : {
3824 : std::unique_ptr<OGRCoordinateTransformation> poRevCT{};
3825 : bool bIsPolar = false;
3826 : bool bIsNorthPolar = false;
3827 : };
3828 :
3829 : /************************************************************************/
3830 : /* TransformWithOptionsCache() */
3831 : /************************************************************************/
3832 :
3833 916 : OGRGeometryFactory::TransformWithOptionsCache::TransformWithOptionsCache()
3834 916 : : d(new Private())
3835 : {
3836 916 : }
3837 :
3838 : /************************************************************************/
3839 : /* ~TransformWithOptionsCache() */
3840 : /************************************************************************/
3841 :
3842 916 : OGRGeometryFactory::TransformWithOptionsCache::~TransformWithOptionsCache()
3843 : {
3844 916 : }
3845 :
3846 : /************************************************************************/
3847 : /* transformWithOptions() */
3848 : /************************************************************************/
3849 :
3850 : /** Transform a geometry.
3851 : * @param poSrcGeom source geometry
3852 : * @param poCT coordinate transformation object, or NULL.
3853 : * @param papszOptions options. Including WRAPDATELINE=YES and DATELINEOFFSET=.
3854 : * @param cache Cache. May increase performance if persisted between invocations
3855 : * @return (new) transformed geometry.
3856 : */
3857 10167 : OGRGeometry *OGRGeometryFactory::transformWithOptions(
3858 : const OGRGeometry *poSrcGeom, OGRCoordinateTransformation *poCT,
3859 : char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache)
3860 : {
3861 10167 : OGRGeometry *poDstGeom = poSrcGeom->clone();
3862 10167 : if (poCT != nullptr)
3863 : {
3864 : #ifdef HAVE_GEOS
3865 10132 : bool bNeedPostCorrection = false;
3866 :
3867 10132 : auto poSourceCRS = poCT->GetSourceCS();
3868 10132 : auto poTargetCRS = poCT->GetTargetCS();
3869 10126 : if (poSourceCRS != nullptr && poTargetCRS != nullptr &&
3870 20258 : poSourceCRS->IsProjected() && poTargetCRS->IsGeographic())
3871 : {
3872 20132 : OGRSpatialReference oSRSWGS84;
3873 10066 : oSRSWGS84.SetWellKnownGeogCS("WGS84");
3874 10066 : oSRSWGS84.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3875 10066 : if (poTargetCRS->IsSame(&oSRSWGS84))
3876 : {
3877 20099 : if (cache.d->poRevCT == nullptr ||
3878 10035 : !cache.d->poRevCT->GetTargetCS()->IsSame(poSourceCRS))
3879 : {
3880 30 : cache.d->poRevCT.reset(OGRCreateCoordinateTransformation(
3881 : &oSRSWGS84, poSourceCRS));
3882 30 : cache.d->bIsNorthPolar = false;
3883 30 : cache.d->bIsPolar = false;
3884 90 : if (cache.d->poRevCT &&
3885 30 : IsPolarToWGS84(poCT, cache.d->poRevCT.get(),
3886 60 : cache.d->bIsNorthPolar))
3887 : {
3888 5 : cache.d->bIsPolar = true;
3889 : }
3890 : }
3891 10064 : auto poRevCT = cache.d->poRevCT.get();
3892 10064 : if (poRevCT != nullptr)
3893 : {
3894 10064 : if (cache.d->bIsPolar)
3895 : {
3896 6 : poDstGeom = TransformBeforePolarToWGS84(
3897 6 : poRevCT, cache.d->bIsNorthPolar, poDstGeom,
3898 : bNeedPostCorrection);
3899 : }
3900 10058 : else if (IsAntimeridianProjToWGS84(poCT, poRevCT,
3901 : poDstGeom))
3902 : {
3903 10055 : poDstGeom = TransformBeforeAntimeridianToWGS84(
3904 : poCT, poRevCT, poDstGeom, bNeedPostCorrection);
3905 : }
3906 : }
3907 : }
3908 : }
3909 : #endif
3910 10132 : OGRErr eErr = poDstGeom->transform(poCT);
3911 10132 : if (eErr != OGRERR_NONE)
3912 : {
3913 0 : delete poDstGeom;
3914 0 : return nullptr;
3915 : }
3916 : #ifdef HAVE_GEOS
3917 10132 : if (bNeedPostCorrection)
3918 : {
3919 14 : SnapCoordsCloseToLatLongBounds(poDstGeom);
3920 : }
3921 : #endif
3922 : }
3923 :
3924 10167 : if (CPLTestBool(CSLFetchNameValueDef(papszOptions, "WRAPDATELINE", "NO")))
3925 : {
3926 101 : if (poDstGeom->getSpatialReference() &&
3927 49 : !poDstGeom->getSpatialReference()->IsGeographic())
3928 : {
3929 : static bool bHasWarned = false;
3930 0 : if (!bHasWarned)
3931 : {
3932 0 : CPLError(
3933 : CE_Warning, CPLE_AppDefined,
3934 : "WRAPDATELINE is without effect when reprojecting to a "
3935 : "non-geographic CRS");
3936 0 : bHasWarned = true;
3937 : }
3938 0 : return poDstGeom;
3939 : }
3940 : // TODO and we should probably also test that the axis order + data axis
3941 : // mapping is long-lat...
3942 : const OGRwkbGeometryType eType =
3943 52 : wkbFlatten(poDstGeom->getGeometryType());
3944 52 : if (eType == wkbPoint)
3945 : {
3946 9 : OGRPoint *poDstPoint = poDstGeom->toPoint();
3947 9 : WrapPointDateLine(poDstPoint);
3948 : }
3949 43 : else if (eType == wkbMultiPoint)
3950 : {
3951 5 : for (auto *poDstPoint : *(poDstGeom->toMultiPoint()))
3952 : {
3953 4 : WrapPointDateLine(poDstPoint);
3954 : }
3955 : }
3956 : else
3957 : {
3958 42 : OGREnvelope sEnvelope;
3959 42 : poDstGeom->getEnvelope(&sEnvelope);
3960 42 : if (sEnvelope.MinX >= -360.0 && sEnvelope.MaxX <= -180.0)
3961 2 : AddOffsetToLon(poDstGeom, 360.0);
3962 40 : else if (sEnvelope.MinX >= 180.0 && sEnvelope.MaxX <= 360.0)
3963 2 : AddOffsetToLon(poDstGeom, -360.0);
3964 : else
3965 : {
3966 : OGRwkbGeometryType eNewType;
3967 38 : if (eType == wkbPolygon || eType == wkbMultiPolygon)
3968 27 : eNewType = wkbMultiPolygon;
3969 11 : else if (eType == wkbLineString || eType == wkbMultiLineString)
3970 10 : eNewType = wkbMultiLineString;
3971 : else
3972 1 : eNewType = wkbGeometryCollection;
3973 :
3974 38 : OGRGeometry *poMultiGeom = createGeometry(eNewType);
3975 : OGRGeometryCollection *poMulti =
3976 38 : poMultiGeom->toGeometryCollection();
3977 :
3978 38 : double dfDateLineOffset = CPLAtofM(
3979 : CSLFetchNameValueDef(papszOptions, "DATELINEOFFSET", "10"));
3980 38 : if (dfDateLineOffset <= 0.0 || dfDateLineOffset >= 360.0)
3981 0 : dfDateLineOffset = 10.0;
3982 :
3983 38 : CutGeometryOnDateLineAndAddToMulti(poMulti, poDstGeom,
3984 : dfDateLineOffset);
3985 :
3986 38 : if (poMulti->getNumGeometries() == 0)
3987 : {
3988 0 : delete poMultiGeom;
3989 : }
3990 38 : else if (poMulti->getNumGeometries() == 1)
3991 : {
3992 13 : delete poDstGeom;
3993 13 : poDstGeom = poMulti->getGeometryRef(0)->clone();
3994 13 : delete poMultiGeom;
3995 : }
3996 : else
3997 : {
3998 25 : delete poDstGeom;
3999 25 : poDstGeom = poMultiGeom;
4000 : }
4001 : }
4002 : }
4003 : }
4004 :
4005 10167 : return poDstGeom;
4006 : }
4007 :
4008 : /************************************************************************/
4009 : /* OGRGeomTransformer() */
4010 : /************************************************************************/
4011 :
4012 : struct OGRGeomTransformer
4013 : {
4014 : std::unique_ptr<OGRCoordinateTransformation> poCT{};
4015 : OGRGeometryFactory::TransformWithOptionsCache cache{};
4016 : CPLStringList aosOptions{};
4017 :
4018 6 : OGRGeomTransformer() = default;
4019 : OGRGeomTransformer(const OGRGeomTransformer &) = delete;
4020 : OGRGeomTransformer &operator=(const OGRGeomTransformer &) = delete;
4021 : };
4022 :
4023 : /************************************************************************/
4024 : /* OGR_GeomTransformer_Create() */
4025 : /************************************************************************/
4026 :
4027 : /** Create a geometry transformer.
4028 : *
4029 : * This is a enhanced version of OGR_G_Transform().
4030 : *
4031 : * When reprojecting geometries from a Polar Stereographic projection or a
4032 : * projection naturally crossing the antimeridian (like UTM Zone 60) to a
4033 : * geographic CRS, it will cut geometries along the antimeridian. So a
4034 : * LineString might be returned as a MultiLineString.
4035 : *
4036 : * The WRAPDATELINE=YES option might be specified for circumstances to correct
4037 : * geometries that incorrectly go from a longitude on a side of the antimeridian
4038 : * to the other side, like a LINESTRING(-179 0,179 0) will be transformed to
4039 : * a MULTILINESTRING ((-179 0,-180 0),(180 0,179 0)). For that use case, hCT
4040 : * might be NULL.
4041 : *
4042 : * @param hCT Coordinate transformation object (will be cloned) or NULL.
4043 : * @param papszOptions NULL terminated list of options, or NULL.
4044 : * Supported options are:
4045 : * <ul>
4046 : * <li>WRAPDATELINE=YES</li>
4047 : * <li>DATELINEOFFSET=longitude_gap_in_degree. Defaults
4048 : * to 10.</li>
4049 : * </ul>
4050 : * @return transformer object to free with OGR_GeomTransformer_Destroy()
4051 : * @since GDAL 3.1
4052 : */
4053 6 : OGRGeomTransformerH OGR_GeomTransformer_Create(OGRCoordinateTransformationH hCT,
4054 : CSLConstList papszOptions)
4055 : {
4056 6 : OGRGeomTransformer *transformer = new OGRGeomTransformer;
4057 6 : if (hCT)
4058 : {
4059 3 : transformer->poCT.reset(
4060 3 : OGRCoordinateTransformation::FromHandle(hCT)->Clone());
4061 : }
4062 6 : transformer->aosOptions.Assign(CSLDuplicate(papszOptions));
4063 6 : return transformer;
4064 : }
4065 :
4066 : /************************************************************************/
4067 : /* OGR_GeomTransformer_Transform() */
4068 : /************************************************************************/
4069 :
4070 : /** Transforms a geometry.
4071 : *
4072 : * @param hTransformer transformer object.
4073 : * @param hGeom Source geometry.
4074 : * @return a new geometry (or NULL) to destroy with OGR_G_DestroyGeometry()
4075 : * @since GDAL 3.1
4076 : */
4077 6 : OGRGeometryH OGR_GeomTransformer_Transform(OGRGeomTransformerH hTransformer,
4078 : OGRGeometryH hGeom)
4079 : {
4080 6 : VALIDATE_POINTER1(hTransformer, "OGR_GeomTransformer_Transform", nullptr);
4081 6 : VALIDATE_POINTER1(hGeom, "OGR_GeomTransformer_Transform", nullptr);
4082 :
4083 12 : return OGRGeometry::ToHandle(OGRGeometryFactory::transformWithOptions(
4084 6 : OGRGeometry::FromHandle(hGeom), hTransformer->poCT.get(),
4085 12 : hTransformer->aosOptions.List(), hTransformer->cache));
4086 : }
4087 :
4088 : /************************************************************************/
4089 : /* OGR_GeomTransformer_Destroy() */
4090 : /************************************************************************/
4091 :
4092 : /** Destroy a geometry transformer allocated with OGR_GeomTransformer_Create()
4093 : *
4094 : * @param hTransformer transformer object.
4095 : * @since GDAL 3.1
4096 : */
4097 6 : void OGR_GeomTransformer_Destroy(OGRGeomTransformerH hTransformer)
4098 : {
4099 6 : delete hTransformer;
4100 6 : }
4101 :
4102 : /************************************************************************/
4103 : /* OGRGF_GetDefaultStepSize() */
4104 : /************************************************************************/
4105 :
4106 4300 : static double OGRGF_GetDefaultStepSize()
4107 : {
4108 : // coverity[tainted_data]
4109 4300 : return CPLAtofM(CPLGetConfigOption("OGR_ARC_STEPSIZE", "4"));
4110 : }
4111 :
4112 : /************************************************************************/
4113 : /* DISTANCE() */
4114 : /************************************************************************/
4115 :
4116 311641 : static inline double DISTANCE(double x1, double y1, double x2, double y2)
4117 : {
4118 311641 : return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
4119 : }
4120 :
4121 : /************************************************************************/
4122 : /* approximateArcAngles() */
4123 : /************************************************************************/
4124 :
4125 : /**
4126 : * Stroke arc to linestring.
4127 : *
4128 : * Stroke an arc of a circle to a linestring based on a center
4129 : * point, radius, start angle and end angle, all angles in degrees.
4130 : *
4131 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
4132 : * used. This is currently 4 degrees unless the user has overridden the
4133 : * value with the OGR_ARC_STEPSIZE configuration variable.
4134 : *
4135 : * If the OGR_ARC_MAX_GAP configuration variable is set, the straight-line
4136 : * distance between adjacent pairs of interpolated points will be limited to
4137 : * the specified distance. If the distance between a pair of points exceeds
4138 : * this maximum, additional points are interpolated between the two points.
4139 : *
4140 : * @see CPLSetConfigOption()
4141 : *
4142 : * @param dfCenterX center X
4143 : * @param dfCenterY center Y
4144 : * @param dfZ center Z
4145 : * @param dfPrimaryRadius X radius of ellipse.
4146 : * @param dfSecondaryRadius Y radius of ellipse.
4147 : * @param dfRotation rotation of the ellipse clockwise.
4148 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
4149 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
4150 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
4151 : * arc, zero to use the default setting.
4152 : * @param bUseMaxGap Optional: whether to honor OGR_ARC_MAX_GAP.
4153 : *
4154 : * @return OGRLineString geometry representing an approximation of the arc.
4155 : *
4156 : * @since OGR 1.8.0
4157 : */
4158 :
4159 111 : OGRGeometry *OGRGeometryFactory::approximateArcAngles(
4160 : double dfCenterX, double dfCenterY, double dfZ, double dfPrimaryRadius,
4161 : double dfSecondaryRadius, double dfRotation, double dfStartAngle,
4162 : double dfEndAngle, double dfMaxAngleStepSizeDegrees,
4163 : const bool bUseMaxGap /* = false */)
4164 :
4165 : {
4166 111 : OGRLineString *poLine = new OGRLineString();
4167 111 : const double dfRotationRadians = dfRotation * M_PI / 180.0;
4168 :
4169 : // Support default arc step setting.
4170 111 : if (dfMaxAngleStepSizeDegrees < 1e-6)
4171 : {
4172 110 : dfMaxAngleStepSizeDegrees = OGRGF_GetDefaultStepSize();
4173 : }
4174 :
4175 : // Determine maximum interpolation gap. This is the largest straight-line
4176 : // distance allowed between pairs of interpolated points. Default zero,
4177 : // meaning no gap.
4178 : // coverity[tainted_data]
4179 : const double dfMaxInterpolationGap =
4180 111 : bUseMaxGap ? CPLAtofM(CPLGetConfigOption("OGR_ARC_MAX_GAP", "0")) : 0.0;
4181 :
4182 : // Is this a full circle?
4183 111 : const bool bIsFullCircle = fabs(dfEndAngle - dfStartAngle) == 360.0;
4184 :
4185 : // Switch direction.
4186 111 : dfStartAngle *= -1;
4187 111 : dfEndAngle *= -1;
4188 :
4189 : // Figure out the number of slices to make this into.
4190 : int nVertexCount =
4191 222 : std::max(2, static_cast<int>(ceil(fabs(dfEndAngle - dfStartAngle) /
4192 111 : dfMaxAngleStepSizeDegrees) +
4193 111 : 1));
4194 111 : const double dfSlice = (dfEndAngle - dfStartAngle) / (nVertexCount - 1);
4195 :
4196 : // If it is a full circle we will work out the last point separately.
4197 111 : if (bIsFullCircle)
4198 : {
4199 51 : nVertexCount--;
4200 : }
4201 :
4202 : /* -------------------------------------------------------------------- */
4203 : /* Compute the interpolated points. */
4204 : /* -------------------------------------------------------------------- */
4205 111 : double dfLastX = 0.0;
4206 111 : double dfLastY = 0.0;
4207 111 : int nTotalAddPoints = 0;
4208 6876 : for (int iPoint = 0; iPoint < nVertexCount; iPoint++)
4209 : {
4210 6765 : const double dfAngleOnEllipse =
4211 6765 : (dfStartAngle + iPoint * dfSlice) * M_PI / 180.0;
4212 :
4213 : // Compute position on the unrotated ellipse.
4214 6765 : const double dfEllipseX = cos(dfAngleOnEllipse) * dfPrimaryRadius;
4215 6765 : const double dfEllipseY = sin(dfAngleOnEllipse) * dfSecondaryRadius;
4216 :
4217 : // Is this point too far from the previous point?
4218 6765 : if (iPoint && dfMaxInterpolationGap != 0.0)
4219 : {
4220 : const double dfDistFromLast =
4221 1 : DISTANCE(dfLastX, dfLastY, dfEllipseX, dfEllipseY);
4222 :
4223 1 : if (dfDistFromLast > dfMaxInterpolationGap)
4224 : {
4225 1 : const int nAddPoints =
4226 1 : static_cast<int>(dfDistFromLast / dfMaxInterpolationGap);
4227 1 : const double dfAddSlice = dfSlice / (nAddPoints + 1);
4228 :
4229 : // Interpolate additional points
4230 3 : for (int iAddPoint = 0; iAddPoint < nAddPoints; iAddPoint++)
4231 : {
4232 2 : const double dfAddAngleOnEllipse =
4233 2 : (dfStartAngle + (iPoint - 1) * dfSlice +
4234 2 : (iAddPoint + 1) * dfAddSlice) *
4235 : (M_PI / 180.0);
4236 :
4237 2 : poLine->setPoint(
4238 2 : iPoint + nTotalAddPoints + iAddPoint,
4239 2 : cos(dfAddAngleOnEllipse) * dfPrimaryRadius,
4240 2 : sin(dfAddAngleOnEllipse) * dfSecondaryRadius, dfZ);
4241 : }
4242 :
4243 1 : nTotalAddPoints += nAddPoints;
4244 : }
4245 : }
4246 :
4247 6765 : poLine->setPoint(iPoint + nTotalAddPoints, dfEllipseX, dfEllipseY, dfZ);
4248 6765 : dfLastX = dfEllipseX;
4249 6765 : dfLastY = dfEllipseY;
4250 : }
4251 :
4252 : /* -------------------------------------------------------------------- */
4253 : /* Rotate and translate the ellipse. */
4254 : /* -------------------------------------------------------------------- */
4255 111 : nVertexCount = poLine->getNumPoints();
4256 6878 : for (int iPoint = 0; iPoint < nVertexCount; iPoint++)
4257 : {
4258 6767 : const double dfEllipseX = poLine->getX(iPoint);
4259 6767 : const double dfEllipseY = poLine->getY(iPoint);
4260 :
4261 : // Rotate this position around the center of the ellipse.
4262 6767 : const double dfArcX = dfCenterX + dfEllipseX * cos(dfRotationRadians) +
4263 6767 : dfEllipseY * sin(dfRotationRadians);
4264 6767 : const double dfArcY = dfCenterY - dfEllipseX * sin(dfRotationRadians) +
4265 6767 : dfEllipseY * cos(dfRotationRadians);
4266 :
4267 6767 : poLine->setPoint(iPoint, dfArcX, dfArcY, dfZ);
4268 : }
4269 :
4270 : /* -------------------------------------------------------------------- */
4271 : /* If we're asked to make a full circle, ensure the start and */
4272 : /* end points coincide exactly, in spite of any rounding error. */
4273 : /* -------------------------------------------------------------------- */
4274 111 : if (bIsFullCircle)
4275 : {
4276 102 : OGRPoint oPoint;
4277 51 : poLine->getPoint(0, &oPoint);
4278 51 : poLine->setPoint(nVertexCount, &oPoint);
4279 : }
4280 :
4281 111 : return poLine;
4282 : }
4283 :
4284 : /************************************************************************/
4285 : /* OGR_G_ApproximateArcAngles() */
4286 : /************************************************************************/
4287 :
4288 : /**
4289 : * Stroke arc to linestring.
4290 : *
4291 : * Stroke an arc of a circle to a linestring based on a center
4292 : * point, radius, start angle and end angle, all angles in degrees.
4293 : *
4294 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
4295 : * used. This is currently 4 degrees unless the user has overridden the
4296 : * value with the OGR_ARC_STEPSIZE configuration variable.
4297 : *
4298 : * @see CPLSetConfigOption()
4299 : *
4300 : * @param dfCenterX center X
4301 : * @param dfCenterY center Y
4302 : * @param dfZ center Z
4303 : * @param dfPrimaryRadius X radius of ellipse.
4304 : * @param dfSecondaryRadius Y radius of ellipse.
4305 : * @param dfRotation rotation of the ellipse clockwise.
4306 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
4307 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
4308 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
4309 : * arc, zero to use the default setting.
4310 : *
4311 : * @return OGRLineString geometry representing an approximation of the arc.
4312 : *
4313 : * @since OGR 1.8.0
4314 : */
4315 :
4316 1 : OGRGeometryH CPL_DLL OGR_G_ApproximateArcAngles(
4317 : double dfCenterX, double dfCenterY, double dfZ, double dfPrimaryRadius,
4318 : double dfSecondaryRadius, double dfRotation, double dfStartAngle,
4319 : double dfEndAngle, double dfMaxAngleStepSizeDegrees)
4320 :
4321 : {
4322 1 : return OGRGeometry::ToHandle(OGRGeometryFactory::approximateArcAngles(
4323 : dfCenterX, dfCenterY, dfZ, dfPrimaryRadius, dfSecondaryRadius,
4324 1 : dfRotation, dfStartAngle, dfEndAngle, dfMaxAngleStepSizeDegrees));
4325 : }
4326 :
4327 : /************************************************************************/
4328 : /* forceToLineString() */
4329 : /************************************************************************/
4330 :
4331 : /**
4332 : * \brief Convert to line string.
4333 : *
4334 : * Tries to force the provided geometry to be a line string. This nominally
4335 : * effects a change on multilinestrings.
4336 : * In GDAL 2.0, for polygons or curvepolygons that have a single exterior ring,
4337 : * it will return the ring. For circular strings or compound curves, it will
4338 : * return an approximated line string.
4339 : *
4340 : * The passed in geometry is
4341 : * consumed and a new one returned (or potentially the same one).
4342 : *
4343 : * @param poGeom the input geometry - ownership is passed to the method.
4344 : * @param bOnlyInOrder flag that, if set to FALSE, indicate that the order of
4345 : * points in a linestring might be reversed if it enables
4346 : * to match the extremity of another linestring. If set
4347 : * to TRUE, the start of a linestring must match the end
4348 : * of another linestring.
4349 : * @return new geometry.
4350 : */
4351 :
4352 171 : OGRGeometry *OGRGeometryFactory::forceToLineString(OGRGeometry *poGeom,
4353 : bool bOnlyInOrder)
4354 :
4355 : {
4356 171 : if (poGeom == nullptr)
4357 2 : return nullptr;
4358 :
4359 169 : const OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
4360 :
4361 : /* -------------------------------------------------------------------- */
4362 : /* If this is already a LineString, nothing to do */
4363 : /* -------------------------------------------------------------------- */
4364 169 : if (eGeomType == wkbLineString)
4365 : {
4366 : // Except if it is a linearring.
4367 25 : poGeom = OGRCurve::CastToLineString(poGeom->toCurve());
4368 :
4369 25 : return poGeom;
4370 : }
4371 :
4372 : /* -------------------------------------------------------------------- */
4373 : /* If it is a polygon with a single ring, return it */
4374 : /* -------------------------------------------------------------------- */
4375 144 : if (eGeomType == wkbPolygon || eGeomType == wkbCurvePolygon)
4376 : {
4377 30 : OGRCurvePolygon *poCP = poGeom->toCurvePolygon();
4378 30 : if (poCP->getNumInteriorRings() == 0)
4379 : {
4380 28 : OGRCurve *poRing = poCP->stealExteriorRingCurve();
4381 28 : delete poCP;
4382 28 : return forceToLineString(poRing);
4383 : }
4384 2 : return poGeom;
4385 : }
4386 :
4387 : /* -------------------------------------------------------------------- */
4388 : /* If it is a curve line, call CurveToLine() */
4389 : /* -------------------------------------------------------------------- */
4390 114 : if (eGeomType == wkbCircularString || eGeomType == wkbCompoundCurve)
4391 : {
4392 66 : OGRGeometry *poNewGeom = poGeom->toCurve()->CurveToLine();
4393 66 : delete poGeom;
4394 66 : return poNewGeom;
4395 : }
4396 :
4397 48 : if (eGeomType != wkbGeometryCollection && eGeomType != wkbMultiLineString &&
4398 : eGeomType != wkbMultiCurve)
4399 17 : return poGeom;
4400 :
4401 : // Build an aggregated linestring from all the linestrings in the container.
4402 31 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4403 31 : if (poGeom->hasCurveGeometry())
4404 : {
4405 : OGRGeometryCollection *poNewGC =
4406 7 : poGC->getLinearGeometry()->toGeometryCollection();
4407 7 : delete poGC;
4408 7 : poGC = poNewGC;
4409 : }
4410 :
4411 31 : if (poGC->getNumGeometries() == 0)
4412 : {
4413 3 : poGeom = new OGRLineString();
4414 3 : poGeom->assignSpatialReference(poGC->getSpatialReference());
4415 3 : delete poGC;
4416 3 : return poGeom;
4417 : }
4418 :
4419 28 : int iGeom0 = 0;
4420 69 : while (iGeom0 < poGC->getNumGeometries())
4421 : {
4422 41 : if (wkbFlatten(poGC->getGeometryRef(iGeom0)->getGeometryType()) !=
4423 : wkbLineString)
4424 : {
4425 12 : iGeom0++;
4426 26 : continue;
4427 : }
4428 :
4429 : OGRLineString *poLineString0 =
4430 29 : poGC->getGeometryRef(iGeom0)->toLineString();
4431 29 : if (poLineString0->getNumPoints() < 2)
4432 : {
4433 14 : iGeom0++;
4434 14 : continue;
4435 : }
4436 :
4437 30 : OGRPoint pointStart0;
4438 15 : poLineString0->StartPoint(&pointStart0);
4439 30 : OGRPoint pointEnd0;
4440 15 : poLineString0->EndPoint(&pointEnd0);
4441 :
4442 15 : int iGeom1 = iGeom0 + 1; // Used after for.
4443 17 : for (; iGeom1 < poGC->getNumGeometries(); iGeom1++)
4444 : {
4445 6 : if (wkbFlatten(poGC->getGeometryRef(iGeom1)->getGeometryType()) !=
4446 : wkbLineString)
4447 1 : continue;
4448 :
4449 : OGRLineString *poLineString1 =
4450 6 : poGC->getGeometryRef(iGeom1)->toLineString();
4451 6 : if (poLineString1->getNumPoints() < 2)
4452 1 : continue;
4453 :
4454 5 : OGRPoint pointStart1;
4455 5 : poLineString1->StartPoint(&pointStart1);
4456 5 : OGRPoint pointEnd1;
4457 5 : poLineString1->EndPoint(&pointEnd1);
4458 :
4459 5 : if (!bOnlyInOrder && (pointEnd0.Equals(&pointEnd1) ||
4460 0 : pointStart0.Equals(&pointStart1)))
4461 : {
4462 0 : poLineString1->reversePoints();
4463 0 : poLineString1->StartPoint(&pointStart1);
4464 0 : poLineString1->EndPoint(&pointEnd1);
4465 : }
4466 :
4467 5 : if (pointEnd0.Equals(&pointStart1))
4468 : {
4469 4 : poLineString0->addSubLineString(poLineString1, 1);
4470 4 : poGC->removeGeometry(iGeom1);
4471 4 : break;
4472 : }
4473 :
4474 1 : if (pointEnd1.Equals(&pointStart0))
4475 : {
4476 0 : poLineString1->addSubLineString(poLineString0, 1);
4477 0 : poGC->removeGeometry(iGeom0);
4478 0 : break;
4479 : }
4480 : }
4481 :
4482 15 : if (iGeom1 == poGC->getNumGeometries())
4483 : {
4484 14 : iGeom0++;
4485 : }
4486 : }
4487 :
4488 28 : if (poGC->getNumGeometries() == 1)
4489 : {
4490 20 : OGRGeometry *poSingleGeom = poGC->getGeometryRef(0);
4491 20 : poGC->removeGeometry(0, FALSE);
4492 20 : delete poGC;
4493 :
4494 20 : return poSingleGeom;
4495 : }
4496 :
4497 8 : return poGC;
4498 : }
4499 :
4500 : /************************************************************************/
4501 : /* OGR_G_ForceToLineString() */
4502 : /************************************************************************/
4503 :
4504 : /**
4505 : * \brief Convert to line string.
4506 : *
4507 : * This function is the same as the C++ method
4508 : * OGRGeometryFactory::forceToLineString().
4509 : *
4510 : * @param hGeom handle to the geometry to convert (ownership surrendered).
4511 : * @return the converted geometry (ownership to caller).
4512 : *
4513 : * @since GDAL/OGR 1.10.0
4514 : */
4515 :
4516 58 : OGRGeometryH OGR_G_ForceToLineString(OGRGeometryH hGeom)
4517 :
4518 : {
4519 58 : return OGRGeometry::ToHandle(
4520 58 : OGRGeometryFactory::forceToLineString(OGRGeometry::FromHandle(hGeom)));
4521 : }
4522 :
4523 : /************************************************************************/
4524 : /* forceTo() */
4525 : /************************************************************************/
4526 :
4527 : /**
4528 : * \brief Convert to another geometry type
4529 : *
4530 : * Tries to force the provided geometry to the specified geometry type.
4531 : *
4532 : * It can promote 'single' geometry type to their corresponding collection type
4533 : * (see OGR_GT_GetCollection()) or the reverse. non-linear geometry type to
4534 : * their corresponding linear geometry type (see OGR_GT_GetLinear()), by
4535 : * possibly approximating circular arcs they may contain. Regarding conversion
4536 : * from linear geometry types to curve geometry types, only "wrapping" will be
4537 : * done. No attempt to retrieve potential circular arcs by de-approximating
4538 : * stroking will be done. For that, OGRGeometry::getCurveGeometry() can be used.
4539 : *
4540 : * The passed in geometry is consumed and a new one returned (or potentially the
4541 : * same one).
4542 : *
4543 : * Starting with GDAL 3.9, this method honours the dimensionality of eTargetType.
4544 : *
4545 : * @param poGeom the input geometry - ownership is passed to the method.
4546 : * @param eTargetType target output geometry type.
4547 : * @param papszOptions options as a null-terminated list of strings or NULL.
4548 : * @return new geometry, or nullptr in case of error.
4549 : *
4550 : * @since GDAL 2.0
4551 : */
4552 :
4553 2242 : OGRGeometry *OGRGeometryFactory::forceTo(OGRGeometry *poGeom,
4554 : OGRwkbGeometryType eTargetType,
4555 : const char *const *papszOptions)
4556 : {
4557 2242 : if (poGeom == nullptr)
4558 0 : return poGeom;
4559 :
4560 2242 : const OGRwkbGeometryType eTargetTypeFlat = wkbFlatten(eTargetType);
4561 2242 : if (eTargetTypeFlat == wkbUnknown)
4562 252 : return poGeom;
4563 :
4564 1990 : if (poGeom->IsEmpty())
4565 : {
4566 277 : OGRGeometry *poRet = createGeometry(eTargetType);
4567 277 : if (poRet)
4568 : {
4569 277 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4570 277 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4571 277 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4572 : }
4573 277 : delete poGeom;
4574 277 : return poRet;
4575 : }
4576 :
4577 1713 : OGRwkbGeometryType eType = poGeom->getGeometryType();
4578 1713 : OGRwkbGeometryType eTypeFlat = wkbFlatten(eType);
4579 :
4580 1713 : if (eTargetTypeFlat != eTargetType && (eType == eTypeFlat))
4581 : {
4582 57 : auto poGeomNew = forceTo(poGeom, eTargetTypeFlat, papszOptions);
4583 57 : if (poGeomNew)
4584 : {
4585 57 : poGeomNew->set3D(OGR_GT_HasZ(eTargetType));
4586 57 : poGeomNew->setMeasured(OGR_GT_HasM(eTargetType));
4587 : }
4588 57 : return poGeomNew;
4589 : }
4590 :
4591 1656 : if (eTypeFlat == eTargetTypeFlat)
4592 : {
4593 502 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4594 502 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4595 502 : return poGeom;
4596 : }
4597 :
4598 1154 : eType = eTypeFlat;
4599 :
4600 1279 : if (OGR_GT_IsSubClassOf(eType, wkbPolyhedralSurface) &&
4601 125 : (eTargetTypeFlat == wkbMultiSurface ||
4602 : eTargetTypeFlat == wkbGeometryCollection))
4603 : {
4604 55 : OGRwkbGeometryType eTempGeomType = wkbMultiPolygon;
4605 55 : if (OGR_GT_HasZ(eTargetType))
4606 51 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4607 55 : if (OGR_GT_HasM(eTargetType))
4608 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4609 55 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4610 55 : eTargetType, papszOptions);
4611 : }
4612 :
4613 1099 : if (OGR_GT_IsSubClassOf(eType, wkbGeometryCollection) &&
4614 : eTargetTypeFlat == wkbGeometryCollection)
4615 : {
4616 111 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4617 111 : auto poRet = OGRGeometryCollection::CastToGeometryCollection(poGC);
4618 111 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4619 111 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4620 111 : return poRet;
4621 : }
4622 :
4623 988 : if (eType == wkbTriangle && eTargetTypeFlat == wkbPolyhedralSurface)
4624 : {
4625 1 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4626 1 : poPS->assignSpatialReference(poGeom->getSpatialReference());
4627 1 : poPS->addGeometryDirectly(OGRTriangle::CastToPolygon(poGeom));
4628 1 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4629 1 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4630 1 : return poPS;
4631 : }
4632 987 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbPolyhedralSurface)
4633 : {
4634 3 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4635 3 : poPS->assignSpatialReference(poGeom->getSpatialReference());
4636 3 : poPS->addGeometryDirectly(poGeom);
4637 3 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4638 3 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4639 3 : return poPS;
4640 : }
4641 984 : else if (eType == wkbMultiPolygon &&
4642 : eTargetTypeFlat == wkbPolyhedralSurface)
4643 : {
4644 2 : OGRMultiPolygon *poMP = poGeom->toMultiPolygon();
4645 2 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4646 4 : for (int i = 0; i < poMP->getNumGeometries(); ++i)
4647 : {
4648 2 : poPS->addGeometry(poMP->getGeometryRef(i));
4649 : }
4650 2 : delete poGeom;
4651 2 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4652 2 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4653 2 : return poPS;
4654 : }
4655 982 : else if (eType == wkbTIN && eTargetTypeFlat == wkbPolyhedralSurface)
4656 : {
4657 1 : poGeom = OGRTriangulatedSurface::CastToPolyhedralSurface(
4658 : poGeom->toTriangulatedSurface());
4659 : }
4660 981 : else if (eType == wkbCurvePolygon &&
4661 : eTargetTypeFlat == wkbPolyhedralSurface)
4662 : {
4663 1 : OGRwkbGeometryType eTempGeomType = wkbPolygon;
4664 1 : if (OGR_GT_HasZ(eTargetType))
4665 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4666 1 : if (OGR_GT_HasM(eTargetType))
4667 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4668 1 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4669 1 : eTargetType, papszOptions);
4670 : }
4671 980 : else if (eType == wkbMultiSurface &&
4672 : eTargetTypeFlat == wkbPolyhedralSurface)
4673 : {
4674 1 : OGRwkbGeometryType eTempGeomType = wkbMultiPolygon;
4675 1 : if (OGR_GT_HasZ(eTargetType))
4676 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4677 1 : if (OGR_GT_HasM(eTargetType))
4678 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4679 1 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4680 1 : eTargetType, papszOptions);
4681 : }
4682 :
4683 979 : else if (eType == wkbTriangle && eTargetTypeFlat == wkbTIN)
4684 : {
4685 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4686 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4687 1 : poTS->addGeometryDirectly(poGeom);
4688 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4689 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4690 1 : return poTS;
4691 : }
4692 978 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbTIN)
4693 : {
4694 4 : OGRPolygon *poPoly = poGeom->toPolygon();
4695 4 : OGRLinearRing *poLR = poPoly->getExteriorRing();
4696 7 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4697 3 : poPoly->getNumInteriorRings() == 0))
4698 : {
4699 1 : return poGeom;
4700 : }
4701 3 : OGRErr eErr = OGRERR_NONE;
4702 3 : OGRTriangle *poTriangle = new OGRTriangle(*poPoly, eErr);
4703 3 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4704 3 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4705 3 : poTS->addGeometryDirectly(poTriangle);
4706 3 : delete poGeom;
4707 3 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4708 3 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4709 3 : return poTS;
4710 : }
4711 974 : else if (eType == wkbMultiPolygon && eTargetTypeFlat == wkbTIN)
4712 : {
4713 1 : OGRMultiPolygon *poMP = poGeom->toMultiPolygon();
4714 2 : for (const auto poPoly : *poMP)
4715 : {
4716 1 : const OGRLinearRing *poLR = poPoly->getExteriorRing();
4717 2 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4718 1 : poPoly->getNumInteriorRings() == 0))
4719 : {
4720 0 : return poGeom;
4721 : }
4722 : }
4723 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4724 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4725 2 : for (const auto poPoly : *poMP)
4726 : {
4727 1 : OGRErr eErr = OGRERR_NONE;
4728 1 : poTS->addGeometryDirectly(new OGRTriangle(*poPoly, eErr));
4729 : }
4730 1 : delete poGeom;
4731 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4732 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4733 1 : return poTS;
4734 : }
4735 973 : else if (eType == wkbPolyhedralSurface && eTargetTypeFlat == wkbTIN)
4736 : {
4737 2 : OGRPolyhedralSurface *poPS = poGeom->toPolyhedralSurface();
4738 3 : for (const auto poPoly : *poPS)
4739 : {
4740 2 : const OGRLinearRing *poLR = poPoly->getExteriorRing();
4741 3 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4742 1 : poPoly->getNumInteriorRings() == 0))
4743 : {
4744 1 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4745 1 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4746 1 : return poGeom;
4747 : }
4748 : }
4749 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4750 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4751 2 : for (const auto poPoly : *poPS)
4752 : {
4753 1 : OGRErr eErr = OGRERR_NONE;
4754 1 : poTS->addGeometryDirectly(new OGRTriangle(*poPoly, eErr));
4755 : }
4756 1 : delete poGeom;
4757 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4758 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4759 1 : return poTS;
4760 : }
4761 :
4762 971 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbTriangle)
4763 : {
4764 7 : OGRPolygon *poPoly = poGeom->toPolygon();
4765 7 : OGRLinearRing *poLR = poPoly->getExteriorRing();
4766 13 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4767 6 : poPoly->getNumInteriorRings() == 0))
4768 : {
4769 1 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4770 1 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4771 1 : return poGeom;
4772 : }
4773 6 : OGRErr eErr = OGRERR_NONE;
4774 6 : OGRTriangle *poTriangle = new OGRTriangle(*poPoly, eErr);
4775 6 : delete poGeom;
4776 6 : poTriangle->set3D(OGR_GT_HasZ(eTargetType));
4777 6 : poTriangle->setMeasured(OGR_GT_HasM(eTargetType));
4778 6 : return poTriangle;
4779 : }
4780 :
4781 965 : if (eTargetTypeFlat == wkbTriangle || eTargetTypeFlat == wkbTIN ||
4782 : eTargetTypeFlat == wkbPolyhedralSurface)
4783 : {
4784 9 : OGRwkbGeometryType eTempGeomType = wkbPolygon;
4785 9 : if (OGR_GT_HasZ(eTargetType))
4786 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4787 9 : if (OGR_GT_HasM(eTargetType))
4788 1 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4789 9 : OGRGeometry *poPoly = forceTo(poGeom, eTempGeomType, papszOptions);
4790 9 : if (poPoly == poGeom)
4791 0 : return poGeom;
4792 9 : return forceTo(poPoly, eTargetType, papszOptions);
4793 : }
4794 :
4795 956 : if (eType == wkbTriangle && eTargetTypeFlat == wkbGeometryCollection)
4796 : {
4797 1 : OGRGeometryCollection *poGC = new OGRGeometryCollection();
4798 1 : poGC->assignSpatialReference(poGeom->getSpatialReference());
4799 1 : poGC->addGeometryDirectly(poGeom);
4800 1 : poGC->set3D(OGR_GT_HasZ(eTargetType));
4801 1 : poGC->setMeasured(OGR_GT_HasM(eTargetType));
4802 1 : return poGC;
4803 : }
4804 :
4805 : // Promote single to multi.
4806 1645 : if (!OGR_GT_IsSubClassOf(eType, wkbGeometryCollection) &&
4807 690 : OGR_GT_IsSubClassOf(OGR_GT_GetCollection(eType), eTargetType))
4808 : {
4809 197 : OGRGeometry *poRet = createGeometry(eTargetType);
4810 197 : if (poRet == nullptr)
4811 : {
4812 0 : delete poGeom;
4813 0 : return nullptr;
4814 : }
4815 197 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4816 197 : if (eType == wkbLineString)
4817 42 : poGeom = OGRCurve::CastToLineString(poGeom->toCurve());
4818 197 : poRet->toGeometryCollection()->addGeometryDirectly(poGeom);
4819 197 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4820 197 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4821 197 : return poRet;
4822 : }
4823 :
4824 758 : const bool bIsCurve = CPL_TO_BOOL(OGR_GT_IsCurve(eType));
4825 758 : if (bIsCurve && eTargetTypeFlat == wkbCompoundCurve)
4826 : {
4827 30 : auto poRet = OGRCurve::CastToCompoundCurve(poGeom->toCurve());
4828 30 : if (poRet)
4829 : {
4830 28 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4831 28 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4832 : }
4833 30 : return poRet;
4834 : }
4835 728 : else if (bIsCurve && eTargetTypeFlat == wkbCurvePolygon)
4836 : {
4837 26 : OGRCurve *poCurve = poGeom->toCurve();
4838 26 : if (poCurve->getNumPoints() >= 3 && poCurve->get_IsClosed())
4839 : {
4840 18 : OGRCurvePolygon *poCP = new OGRCurvePolygon();
4841 18 : if (poCP->addRingDirectly(poCurve) == OGRERR_NONE)
4842 : {
4843 18 : poCP->assignSpatialReference(poGeom->getSpatialReference());
4844 18 : poCP->set3D(OGR_GT_HasZ(eTargetType));
4845 18 : poCP->setMeasured(OGR_GT_HasM(eTargetType));
4846 18 : return poCP;
4847 : }
4848 : else
4849 : {
4850 0 : delete poCP;
4851 : }
4852 8 : }
4853 : }
4854 778 : else if (eType == wkbLineString &&
4855 76 : OGR_GT_IsSubClassOf(eTargetType, wkbMultiSurface))
4856 : {
4857 23 : OGRGeometry *poTmp = forceTo(poGeom, wkbPolygon, papszOptions);
4858 23 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4859 15 : return forceTo(poTmp, eTargetType, papszOptions);
4860 : }
4861 679 : else if (bIsCurve && eTargetTypeFlat == wkbMultiSurface)
4862 : {
4863 10 : OGRGeometry *poTmp = forceTo(poGeom, wkbCurvePolygon, papszOptions);
4864 10 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4865 10 : return forceTo(poTmp, eTargetType, papszOptions);
4866 : }
4867 669 : else if (bIsCurve && eTargetTypeFlat == wkbMultiPolygon)
4868 : {
4869 13 : OGRGeometry *poTmp = forceTo(poGeom, wkbPolygon, papszOptions);
4870 13 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4871 13 : return forceTo(poTmp, eTargetType, papszOptions);
4872 : }
4873 656 : else if (eType == wkbTriangle && eTargetTypeFlat == wkbCurvePolygon)
4874 : {
4875 1 : auto poRet = OGRSurface::CastToCurvePolygon(
4876 : OGRTriangle::CastToPolygon(poGeom)->toSurface());
4877 1 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4878 1 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4879 1 : return poRet;
4880 : }
4881 655 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbCurvePolygon)
4882 : {
4883 18 : auto poRet = OGRSurface::CastToCurvePolygon(poGeom->toPolygon());
4884 18 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4885 18 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4886 18 : return poRet;
4887 : }
4888 637 : else if (OGR_GT_IsSubClassOf(eType, wkbCurvePolygon) &&
4889 : eTargetTypeFlat == wkbCompoundCurve)
4890 : {
4891 15 : OGRCurvePolygon *poPoly = poGeom->toCurvePolygon();
4892 15 : if (poPoly->getNumInteriorRings() == 0)
4893 : {
4894 14 : OGRCurve *poRet = poPoly->stealExteriorRingCurve();
4895 14 : if (poRet)
4896 14 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4897 14 : delete poPoly;
4898 14 : return forceTo(poRet, eTargetType, papszOptions);
4899 : }
4900 : }
4901 622 : else if (eType == wkbMultiPolygon && eTargetTypeFlat == wkbMultiSurface)
4902 : {
4903 : auto poRet =
4904 10 : OGRMultiPolygon::CastToMultiSurface(poGeom->toMultiPolygon());
4905 10 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4906 10 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4907 10 : return poRet;
4908 : }
4909 612 : else if (eType == wkbMultiLineString && eTargetTypeFlat == wkbMultiCurve)
4910 : {
4911 : auto poRet =
4912 8 : OGRMultiLineString::CastToMultiCurve(poGeom->toMultiLineString());
4913 8 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4914 8 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4915 8 : return poRet;
4916 : }
4917 604 : else if (OGR_GT_IsSubClassOf(eType, wkbGeometryCollection))
4918 : {
4919 247 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4920 247 : if (poGC->getNumGeometries() == 1)
4921 : {
4922 156 : OGRGeometry *poSubGeom = poGC->getGeometryRef(0);
4923 156 : if (poSubGeom)
4924 : {
4925 156 : poSubGeom->assignSpatialReference(
4926 156 : poGeom->getSpatialReference());
4927 156 : poGC->removeGeometry(0, FALSE);
4928 : OGRGeometry *poRet =
4929 156 : forceTo(poSubGeom->clone(), eTargetType, papszOptions);
4930 156 : if (OGR_GT_IsSubClassOf(wkbFlatten(poRet->getGeometryType()),
4931 156 : eTargetType))
4932 : {
4933 121 : delete poGC;
4934 121 : delete poSubGeom;
4935 121 : return poRet;
4936 : }
4937 35 : poGC->addGeometryDirectly(poSubGeom);
4938 35 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4939 35 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4940 35 : delete poRet;
4941 : }
4942 : }
4943 : }
4944 471 : else if (OGR_GT_IsSubClassOf(eType, wkbCurvePolygon) &&
4945 114 : (OGR_GT_IsSubClassOf(eTargetType, wkbMultiSurface) ||
4946 102 : OGR_GT_IsSubClassOf(eTargetType, wkbMultiCurve)))
4947 : {
4948 43 : OGRCurvePolygon *poCP = poGeom->toCurvePolygon();
4949 43 : if (poCP->getNumInteriorRings() == 0)
4950 : {
4951 41 : OGRCurve *poRing = poCP->getExteriorRingCurve();
4952 41 : poRing->assignSpatialReference(poGeom->getSpatialReference());
4953 41 : OGRwkbGeometryType eRingType = poRing->getGeometryType();
4954 41 : OGRGeometry *poRingDup = poRing->clone();
4955 41 : OGRGeometry *poRet = forceTo(poRingDup, eTargetType, papszOptions);
4956 57 : if (poRet->getGeometryType() != eRingType &&
4957 16 : !(eTypeFlat == wkbPolygon &&
4958 : eTargetTypeFlat == wkbMultiLineString))
4959 : {
4960 29 : delete poCP;
4961 29 : return poRet;
4962 : }
4963 : else
4964 : {
4965 12 : delete poRet;
4966 : }
4967 : }
4968 : }
4969 :
4970 475 : if (eTargetTypeFlat == wkbLineString)
4971 : {
4972 85 : poGeom = forceToLineString(poGeom);
4973 85 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4974 85 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4975 : }
4976 390 : else if (eTargetTypeFlat == wkbPolygon)
4977 : {
4978 99 : poGeom = forceToPolygon(poGeom);
4979 99 : if (poGeom)
4980 : {
4981 99 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4982 99 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4983 : }
4984 : }
4985 291 : else if (eTargetTypeFlat == wkbMultiPolygon)
4986 : {
4987 112 : poGeom = forceToMultiPolygon(poGeom);
4988 112 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4989 112 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4990 : }
4991 179 : else if (eTargetTypeFlat == wkbMultiLineString)
4992 : {
4993 37 : poGeom = forceToMultiLineString(poGeom);
4994 37 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4995 37 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4996 : }
4997 142 : else if (eTargetTypeFlat == wkbMultiPoint)
4998 : {
4999 22 : poGeom = forceToMultiPoint(poGeom);
5000 22 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5001 22 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5002 : }
5003 :
5004 475 : return poGeom;
5005 : }
5006 :
5007 : /************************************************************************/
5008 : /* OGR_G_ForceTo() */
5009 : /************************************************************************/
5010 :
5011 : /**
5012 : * \brief Convert to another geometry type
5013 : *
5014 : * This function is the same as the C++ method OGRGeometryFactory::forceTo().
5015 : *
5016 : * @param hGeom the input geometry - ownership is passed to the method.
5017 : * @param eTargetType target output geometry type.
5018 : * @param papszOptions options as a null-terminated list of strings or NULL.
5019 : * @return new geometry.
5020 : *
5021 : * @since GDAL 2.0
5022 : */
5023 :
5024 848 : OGRGeometryH OGR_G_ForceTo(OGRGeometryH hGeom, OGRwkbGeometryType eTargetType,
5025 : char **papszOptions)
5026 :
5027 : {
5028 848 : return OGRGeometry::ToHandle(OGRGeometryFactory::forceTo(
5029 848 : OGRGeometry::FromHandle(hGeom), eTargetType, papszOptions));
5030 : }
5031 :
5032 : /************************************************************************/
5033 : /* GetCurveParameters() */
5034 : /************************************************************************/
5035 :
5036 : /**
5037 : * \brief Returns the parameter of an arc circle.
5038 : *
5039 : * Angles are return in radians, with trigonometic convention (counter clock
5040 : * wise)
5041 : *
5042 : * @param x0 x of first point
5043 : * @param y0 y of first point
5044 : * @param x1 x of intermediate point
5045 : * @param y1 y of intermediate point
5046 : * @param x2 x of final point
5047 : * @param y2 y of final point
5048 : * @param R radius (output)
5049 : * @param cx x of arc center (output)
5050 : * @param cy y of arc center (output)
5051 : * @param alpha0 angle between center and first point, in radians (output)
5052 : * @param alpha1 angle between center and intermediate point, in radians
5053 : * (output)
5054 : * @param alpha2 angle between center and final point, in radians (output)
5055 : * @return TRUE if the points are not aligned and define an arc circle.
5056 : *
5057 : * @since GDAL 2.0
5058 : */
5059 :
5060 187552 : int OGRGeometryFactory::GetCurveParameters(double x0, double y0, double x1,
5061 : double y1, double x2, double y2,
5062 : double &R, double &cx, double &cy,
5063 : double &alpha0, double &alpha1,
5064 : double &alpha2)
5065 : {
5066 187552 : if (CPLIsNan(x0) || CPLIsNan(y0) || CPLIsNan(x1) || CPLIsNan(y1) ||
5067 187552 : CPLIsNan(x2) || CPLIsNan(y2))
5068 : {
5069 0 : return FALSE;
5070 : }
5071 :
5072 : // Circle.
5073 187552 : if (x0 == x2 && y0 == y2)
5074 : {
5075 153 : if (x0 != x1 || y0 != y1)
5076 : {
5077 152 : cx = (x0 + x1) / 2;
5078 152 : cy = (y0 + y1) / 2;
5079 152 : R = DISTANCE(cx, cy, x0, y0);
5080 : // Arbitrarily pick counter-clock-wise order (like PostGIS does).
5081 152 : alpha0 = atan2(y0 - cy, x0 - cx);
5082 152 : alpha1 = alpha0 + M_PI;
5083 152 : alpha2 = alpha0 + 2 * M_PI;
5084 152 : return TRUE;
5085 : }
5086 : else
5087 : {
5088 1 : return FALSE;
5089 : }
5090 : }
5091 :
5092 187399 : double dx01 = x1 - x0;
5093 187399 : double dy01 = y1 - y0;
5094 187399 : double dx12 = x2 - x1;
5095 187399 : double dy12 = y2 - y1;
5096 :
5097 : // Normalize above values so as to make sure we don't end up with
5098 : // computing a difference of too big values.
5099 187399 : double dfScale = fabs(dx01);
5100 187399 : if (fabs(dy01) > dfScale)
5101 93210 : dfScale = fabs(dy01);
5102 187399 : if (fabs(dx12) > dfScale)
5103 47004 : dfScale = fabs(dx12);
5104 187399 : if (fabs(dy12) > dfScale)
5105 46301 : dfScale = fabs(dy12);
5106 187399 : const double dfInvScale = 1.0 / dfScale;
5107 187399 : dx01 *= dfInvScale;
5108 187399 : dy01 *= dfInvScale;
5109 187399 : dx12 *= dfInvScale;
5110 187399 : dy12 *= dfInvScale;
5111 :
5112 187399 : const double det = dx01 * dy12 - dx12 * dy01;
5113 187399 : if (fabs(det) < 1.0e-8 || CPLIsNan(det))
5114 : {
5115 143 : return FALSE;
5116 : }
5117 187256 : const double x01_mid = (x0 + x1) * dfInvScale;
5118 187256 : const double x12_mid = (x1 + x2) * dfInvScale;
5119 187256 : const double y01_mid = (y0 + y1) * dfInvScale;
5120 187256 : const double y12_mid = (y1 + y2) * dfInvScale;
5121 187256 : const double c01 = dx01 * x01_mid + dy01 * y01_mid;
5122 187256 : const double c12 = dx12 * x12_mid + dy12 * y12_mid;
5123 187256 : cx = 0.5 * dfScale * (c01 * dy12 - c12 * dy01) / det;
5124 187256 : cy = 0.5 * dfScale * (-c01 * dx12 + c12 * dx01) / det;
5125 :
5126 187256 : alpha0 = atan2((y0 - cy) * dfInvScale, (x0 - cx) * dfInvScale);
5127 187256 : alpha1 = atan2((y1 - cy) * dfInvScale, (x1 - cx) * dfInvScale);
5128 187256 : alpha2 = atan2((y2 - cy) * dfInvScale, (x2 - cx) * dfInvScale);
5129 187256 : R = DISTANCE(cx, cy, x0, y0);
5130 :
5131 : // If det is negative, the orientation if clockwise.
5132 187256 : if (det < 0)
5133 : {
5134 95252 : if (alpha1 > alpha0)
5135 1331 : alpha1 -= 2 * M_PI;
5136 95252 : if (alpha2 > alpha1)
5137 3424 : alpha2 -= 2 * M_PI;
5138 : }
5139 : else
5140 : {
5141 92004 : if (alpha1 < alpha0)
5142 1308 : alpha1 += 2 * M_PI;
5143 92004 : if (alpha2 < alpha1)
5144 3136 : alpha2 += 2 * M_PI;
5145 : }
5146 :
5147 187256 : CPLAssert((alpha0 <= alpha1 && alpha1 <= alpha2) ||
5148 : (alpha0 >= alpha1 && alpha1 >= alpha2));
5149 :
5150 187256 : return TRUE;
5151 : }
5152 :
5153 : /************************************************************************/
5154 : /* OGRGeometryFactoryStrokeArc() */
5155 : /************************************************************************/
5156 :
5157 4307 : static void OGRGeometryFactoryStrokeArc(OGRLineString *poLine, double cx,
5158 : double cy, double R, double z0,
5159 : double z1, int bHasZ, double alpha0,
5160 : double alpha1, double dfStep,
5161 : int bStealthConstraints)
5162 : {
5163 4307 : const int nSign = dfStep > 0 ? 1 : -1;
5164 :
5165 : // Constant angle between all points, so as to not depend on winding order.
5166 4307 : const double dfNumSteps = fabs((alpha1 - alpha0) / dfStep) + 0.5;
5167 4307 : if (dfNumSteps >= std::numeric_limits<int>::max() ||
5168 4307 : dfNumSteps <= std::numeric_limits<int>::min() || CPLIsNan(dfNumSteps))
5169 : {
5170 0 : CPLError(CE_Warning, CPLE_AppDefined,
5171 : "OGRGeometryFactoryStrokeArc: bogus steps: "
5172 : "%lf %lf %lf %lf",
5173 : alpha0, alpha1, dfStep, dfNumSteps);
5174 0 : return;
5175 : }
5176 :
5177 4307 : int nSteps = static_cast<int>(dfNumSteps);
5178 4307 : if (bStealthConstraints)
5179 : {
5180 : // We need at least 6 intermediate vertex, and if more additional
5181 : // multiples of 2.
5182 4109 : if (nSteps < 1 + 6)
5183 104 : nSteps = 1 + 6;
5184 : else
5185 4005 : nSteps = 1 + 6 + 2 * ((nSteps - (1 + 6) + (2 - 1)) / 2);
5186 : }
5187 198 : else if (nSteps < 4)
5188 : {
5189 194 : nSteps = 4;
5190 : }
5191 4307 : dfStep = nSign * fabs((alpha1 - alpha0) / nSteps);
5192 4307 : double alpha = alpha0 + dfStep;
5193 :
5194 230429 : for (; (alpha - alpha1) * nSign < -1e-8; alpha += dfStep)
5195 : {
5196 226122 : const double dfX = cx + R * cos(alpha);
5197 226122 : const double dfY = cy + R * sin(alpha);
5198 226122 : if (bHasZ)
5199 : {
5200 9104 : const double z =
5201 9104 : z0 + (z1 - z0) * (alpha - alpha0) / (alpha1 - alpha0);
5202 9104 : poLine->addPoint(dfX, dfY, z);
5203 : }
5204 : else
5205 : {
5206 217018 : poLine->addPoint(dfX, dfY);
5207 : }
5208 : }
5209 : }
5210 :
5211 : /************************************************************************/
5212 : /* OGRGF_SetHiddenValue() */
5213 : /************************************************************************/
5214 :
5215 : // TODO(schwehr): Cleanup these static constants.
5216 : constexpr int HIDDEN_ALPHA_WIDTH = 32;
5217 : constexpr GUInt32 HIDDEN_ALPHA_SCALE =
5218 : static_cast<GUInt32>((static_cast<GUIntBig>(1) << HIDDEN_ALPHA_WIDTH) - 2);
5219 : constexpr int HIDDEN_ALPHA_HALF_WIDTH = (HIDDEN_ALPHA_WIDTH / 2);
5220 : constexpr int HIDDEN_ALPHA_HALF_MASK = (1 << HIDDEN_ALPHA_HALF_WIDTH) - 1;
5221 :
5222 : // Encode 16-bit nValue in the 8-lsb of dfX and dfY.
5223 :
5224 : #ifdef CPL_LSB
5225 : constexpr int DOUBLE_LSB_OFFSET = 0;
5226 : #else
5227 : constexpr int DOUBLE_LSB_OFFSET = 7;
5228 : #endif
5229 :
5230 225988 : static void OGRGF_SetHiddenValue(GUInt16 nValue, double &dfX, double &dfY)
5231 : {
5232 225988 : GByte abyData[8] = {};
5233 :
5234 225988 : memcpy(abyData, &dfX, sizeof(double));
5235 225988 : abyData[DOUBLE_LSB_OFFSET] = static_cast<GByte>(nValue & 0xFF);
5236 225988 : memcpy(&dfX, abyData, sizeof(double));
5237 :
5238 225988 : memcpy(abyData, &dfY, sizeof(double));
5239 225988 : abyData[DOUBLE_LSB_OFFSET] = static_cast<GByte>(nValue >> 8);
5240 225988 : memcpy(&dfY, abyData, sizeof(double));
5241 225988 : }
5242 :
5243 : /************************************************************************/
5244 : /* OGRGF_GetHiddenValue() */
5245 : /************************************************************************/
5246 :
5247 : // Decode 16-bit nValue from the 8-lsb of dfX and dfY.
5248 182380 : static GUInt16 OGRGF_GetHiddenValue(double dfX, double dfY)
5249 : {
5250 182380 : GByte abyData[8] = {};
5251 182380 : memcpy(abyData, &dfX, sizeof(double));
5252 182380 : GUInt16 nValue = abyData[DOUBLE_LSB_OFFSET];
5253 182380 : memcpy(abyData, &dfY, sizeof(double));
5254 182380 : nValue |= (abyData[DOUBLE_LSB_OFFSET] << 8);
5255 :
5256 182380 : return nValue;
5257 : }
5258 :
5259 : /************************************************************************/
5260 : /* OGRGF_NeedSwithArcOrder() */
5261 : /************************************************************************/
5262 :
5263 : // We need to define a full ordering between starting point and ending point
5264 : // whatever it is.
5265 9432 : static bool OGRGF_NeedSwithArcOrder(double x0, double y0, double x2, double y2)
5266 : {
5267 9432 : return x0 < x2 || (x0 == x2 && y0 < y2);
5268 : }
5269 :
5270 : /************************************************************************/
5271 : /* curveToLineString() */
5272 : /************************************************************************/
5273 :
5274 : /* clang-format off */
5275 : /**
5276 : * \brief Converts an arc circle into an approximate line string
5277 : *
5278 : * The arc circle is defined by a first point, an intermediate point and a
5279 : * final point.
5280 : *
5281 : * The provided dfMaxAngleStepSizeDegrees is a hint. The discretization
5282 : * algorithm may pick a slightly different value.
5283 : *
5284 : * So as to avoid gaps when rendering curve polygons that share common arcs,
5285 : * this method is guaranteed to return a line with reversed vertex if called
5286 : * with inverted first and final point, and identical intermediate point.
5287 : *
5288 : * @param x0 x of first point
5289 : * @param y0 y of first point
5290 : * @param z0 z of first point
5291 : * @param x1 x of intermediate point
5292 : * @param y1 y of intermediate point
5293 : * @param z1 z of intermediate point
5294 : * @param x2 x of final point
5295 : * @param y2 y of final point
5296 : * @param z2 z of final point
5297 : * @param bHasZ TRUE if z must be taken into account
5298 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
5299 : * arc, zero to use the default setting.
5300 : * @param papszOptions options as a null-terminated list of strings or NULL.
5301 : * Recognized options:
5302 : * <ul>
5303 : * <li>ADD_INTERMEDIATE_POINT=STEALTH/YES/NO (Default to STEALTH).
5304 : * Determine if and how the intermediate point must be output in the
5305 : * linestring. If set to STEALTH, no explicit intermediate point is
5306 : * added but its properties are encoded in low significant bits of
5307 : * intermediate points and OGRGeometryFactory::curveFromLineString() can
5308 : * decode them. This is the best compromise for round-tripping in OGR
5309 : * and better results with PostGIS
5310 : * <a href="http://postgis.org/docs/ST_LineToCurve.html">ST_LineToCurve()</a>.
5311 : * If set to YES, the intermediate point is explicitly added to the
5312 : * linestring. If set to NO, the intermediate point is not explicitly
5313 : * added.
5314 : * </li>
5315 : * </ul>
5316 : *
5317 : * @return the converted geometry (ownership to caller).
5318 : *
5319 : * @since GDAL 2.0
5320 : */
5321 : /* clang-format on */
5322 :
5323 6365 : OGRLineString *OGRGeometryFactory::curveToLineString(
5324 : double x0, double y0, double z0, double x1, double y1, double z1, double x2,
5325 : double y2, double z2, int bHasZ, double dfMaxAngleStepSizeDegrees,
5326 : const char *const *papszOptions)
5327 : {
5328 : // So as to make sure the same curve followed in both direction results
5329 : // in perfectly(=binary identical) symmetrical points.
5330 6365 : if (OGRGF_NeedSwithArcOrder(x0, y0, x2, y2))
5331 : {
5332 : OGRLineString *poLS =
5333 2156 : curveToLineString(x2, y2, z2, x1, y1, z1, x0, y0, z0, bHasZ,
5334 : dfMaxAngleStepSizeDegrees, papszOptions);
5335 2156 : poLS->reversePoints();
5336 2156 : return poLS;
5337 : }
5338 :
5339 4209 : double R = 0.0;
5340 4209 : double cx = 0.0;
5341 4209 : double cy = 0.0;
5342 4209 : double alpha0 = 0.0;
5343 4209 : double alpha1 = 0.0;
5344 4209 : double alpha2 = 0.0;
5345 :
5346 4209 : OGRLineString *poLine = new OGRLineString();
5347 4209 : bool bIsArc = true;
5348 4209 : if (!GetCurveParameters(x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1,
5349 : alpha2))
5350 : {
5351 97 : bIsArc = false;
5352 97 : cx = 0.0;
5353 97 : cy = 0.0;
5354 97 : R = 0.0;
5355 97 : alpha0 = 0.0;
5356 97 : alpha1 = 0.0;
5357 97 : alpha2 = 0.0;
5358 : }
5359 :
5360 4209 : const int nSign = alpha1 >= alpha0 ? 1 : -1;
5361 :
5362 : // support default arc step setting.
5363 4209 : if (dfMaxAngleStepSizeDegrees < 1e-6)
5364 : {
5365 4190 : dfMaxAngleStepSizeDegrees = OGRGF_GetDefaultStepSize();
5366 : }
5367 :
5368 4209 : double dfStep = dfMaxAngleStepSizeDegrees / 180 * M_PI;
5369 4209 : if (dfStep <= 0.01 / 180 * M_PI)
5370 : {
5371 0 : CPLDebug("OGR", "Too small arc step size: limiting to 0.01 degree.");
5372 0 : dfStep = 0.01 / 180 * M_PI;
5373 : }
5374 :
5375 4209 : dfStep *= nSign;
5376 :
5377 4209 : if (bHasZ)
5378 252 : poLine->addPoint(x0, y0, z0);
5379 : else
5380 3957 : poLine->addPoint(x0, y0);
5381 :
5382 4209 : bool bAddIntermediatePoint = false;
5383 4209 : bool bStealth = true;
5384 4215 : for (const char *const *papszIter = papszOptions; papszIter && *papszIter;
5385 : papszIter++)
5386 : {
5387 6 : char *pszKey = nullptr;
5388 6 : const char *pszValue = CPLParseNameValue(*papszIter, &pszKey);
5389 6 : if (pszKey != nullptr && EQUAL(pszKey, "ADD_INTERMEDIATE_POINT"))
5390 : {
5391 4 : if (EQUAL(pszValue, "YES") || EQUAL(pszValue, "TRUE") ||
5392 3 : EQUAL(pszValue, "ON"))
5393 : {
5394 1 : bAddIntermediatePoint = true;
5395 1 : bStealth = false;
5396 : }
5397 3 : else if (EQUAL(pszValue, "NO") || EQUAL(pszValue, "FALSE") ||
5398 1 : EQUAL(pszValue, "OFF"))
5399 : {
5400 2 : bAddIntermediatePoint = false;
5401 2 : bStealth = false;
5402 : }
5403 : else if (EQUAL(pszValue, "STEALTH"))
5404 : {
5405 : // default.
5406 : }
5407 : }
5408 : else
5409 : {
5410 2 : CPLError(CE_Warning, CPLE_NotSupported, "Unsupported option: %s",
5411 : *papszIter);
5412 : }
5413 6 : CPLFree(pszKey);
5414 : }
5415 :
5416 4209 : if (!bIsArc || bAddIntermediatePoint)
5417 : {
5418 98 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z0, z1, bHasZ, alpha0,
5419 : alpha1, dfStep, FALSE);
5420 :
5421 98 : if (bHasZ)
5422 25 : poLine->addPoint(x1, y1, z1);
5423 : else
5424 73 : poLine->addPoint(x1, y1);
5425 :
5426 98 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z1, z2, bHasZ, alpha1,
5427 : alpha2, dfStep, FALSE);
5428 : }
5429 : else
5430 : {
5431 4111 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z0, z2, bHasZ, alpha0,
5432 : alpha2, dfStep, bStealth);
5433 :
5434 4111 : if (bStealth && poLine->getNumPoints() > 6)
5435 : {
5436 : // 'Hide' the angle of the intermediate point in the 8
5437 : // low-significant bits of the x, y of the first 2 computed points
5438 : // (so 32 bits), then put 0xFF, and on the last couple points put
5439 : // again the angle but in reverse order, so that overall the
5440 : // low-significant bits of all the points are symmetrical w.r.t the
5441 : // mid-point.
5442 4109 : const double dfRatio = (alpha1 - alpha0) / (alpha2 - alpha0);
5443 4109 : double dfAlphaRatio = 0.5 + HIDDEN_ALPHA_SCALE * dfRatio;
5444 4109 : if (dfAlphaRatio < 0.0)
5445 : {
5446 0 : CPLError(CE_Warning, CPLE_AppDefined, "AlphaRation < 0: %lf",
5447 : dfAlphaRatio);
5448 0 : dfAlphaRatio *= -1;
5449 : }
5450 8218 : else if (dfAlphaRatio >= std::numeric_limits<GUInt32>::max() ||
5451 4109 : CPLIsNan(dfAlphaRatio))
5452 : {
5453 0 : CPLError(CE_Warning, CPLE_AppDefined,
5454 : "AlphaRatio too large: %lf", dfAlphaRatio);
5455 0 : dfAlphaRatio = std::numeric_limits<GUInt32>::max();
5456 : }
5457 4109 : const GUInt32 nAlphaRatio = static_cast<GUInt32>(dfAlphaRatio);
5458 4109 : const GUInt16 nAlphaRatioLow = nAlphaRatio & HIDDEN_ALPHA_HALF_MASK;
5459 4109 : const GUInt16 nAlphaRatioHigh =
5460 4109 : nAlphaRatio >> HIDDEN_ALPHA_HALF_WIDTH;
5461 : // printf("alpha0=%f, alpha1=%f, alpha2=%f, dfRatio=%f, "/*ok*/
5462 : // "nAlphaRatio = %u\n",
5463 : // alpha0, alpha1, alpha2, dfRatio, nAlphaRatio);
5464 :
5465 4109 : CPLAssert(((poLine->getNumPoints() - 1 - 6) % 2) == 0);
5466 :
5467 117103 : for (int i = 1; i + 1 < poLine->getNumPoints(); i += 2)
5468 : {
5469 112994 : GUInt16 nVal = 0xFFFF;
5470 :
5471 112994 : double dfX = poLine->getX(i);
5472 112994 : double dfY = poLine->getY(i);
5473 112994 : if (i == 1)
5474 4109 : nVal = nAlphaRatioLow;
5475 108885 : else if (i == poLine->getNumPoints() - 2)
5476 4109 : nVal = nAlphaRatioHigh;
5477 112994 : OGRGF_SetHiddenValue(nVal, dfX, dfY);
5478 112994 : poLine->setPoint(i, dfX, dfY);
5479 :
5480 112994 : dfX = poLine->getX(i + 1);
5481 112994 : dfY = poLine->getY(i + 1);
5482 112994 : if (i == 1)
5483 4109 : nVal = nAlphaRatioHigh;
5484 108885 : else if (i == poLine->getNumPoints() - 2)
5485 4109 : nVal = nAlphaRatioLow;
5486 112994 : OGRGF_SetHiddenValue(nVal, dfX, dfY);
5487 112994 : poLine->setPoint(i + 1, dfX, dfY);
5488 : }
5489 : }
5490 : }
5491 :
5492 4209 : if (bHasZ)
5493 252 : poLine->addPoint(x2, y2, z2);
5494 : else
5495 3957 : poLine->addPoint(x2, y2);
5496 :
5497 4209 : return poLine;
5498 : }
5499 :
5500 : /************************************************************************/
5501 : /* OGRGF_FixAngle() */
5502 : /************************************************************************/
5503 :
5504 : // Fix dfAngle by offsets of 2 PI so that it lies between dfAngleStart and
5505 : // dfAngleStop, whatever their respective order.
5506 181205 : static double OGRGF_FixAngle(double dfAngleStart, double dfAngleStop,
5507 : double dfAngle)
5508 : {
5509 181205 : if (dfAngleStart < dfAngleStop)
5510 : {
5511 124138 : while (dfAngle <= dfAngleStart + 1e-8)
5512 35045 : dfAngle += 2 * M_PI;
5513 : }
5514 : else
5515 : {
5516 126980 : while (dfAngle >= dfAngleStart - 1e-8)
5517 34868 : dfAngle -= 2 * M_PI;
5518 : }
5519 181205 : return dfAngle;
5520 : }
5521 :
5522 : /************************************************************************/
5523 : /* OGRGF_DetectArc() */
5524 : /************************************************************************/
5525 :
5526 : // #define VERBOSE_DEBUG_CURVEFROMLINESTRING
5527 :
5528 12211 : static inline bool IS_ALMOST_INTEGER(double x)
5529 : {
5530 12211 : const double val = fabs(x - floor(x + 0.5));
5531 12211 : return val < 1.0e-8;
5532 : }
5533 :
5534 3468 : static int OGRGF_DetectArc(const OGRLineString *poLS, int i,
5535 : OGRCompoundCurve *&poCC, OGRCircularString *&poCS,
5536 : OGRLineString *&poLSNew)
5537 : {
5538 3468 : if (i + 3 >= poLS->getNumPoints())
5539 305 : return -1;
5540 :
5541 6326 : OGRPoint p0;
5542 6326 : OGRPoint p1;
5543 6326 : OGRPoint p2;
5544 3163 : poLS->getPoint(i, &p0);
5545 3163 : poLS->getPoint(i + 1, &p1);
5546 3163 : poLS->getPoint(i + 2, &p2);
5547 3163 : double R_1 = 0.0;
5548 3163 : double cx_1 = 0.0;
5549 3163 : double cy_1 = 0.0;
5550 3163 : double alpha0_1 = 0.0;
5551 3163 : double alpha1_1 = 0.0;
5552 3163 : double alpha2_1 = 0.0;
5553 6319 : if (!(OGRGeometryFactory::GetCurveParameters(
5554 : p0.getX(), p0.getY(), p1.getX(), p1.getY(), p2.getX(), p2.getY(),
5555 : R_1, cx_1, cy_1, alpha0_1, alpha1_1, alpha2_1) &&
5556 3156 : fabs(alpha2_1 - alpha0_1) < 2.0 * 20.0 / 180.0 * M_PI))
5557 : {
5558 24 : return -1;
5559 : }
5560 :
5561 3139 : const double dfDeltaAlpha10 = alpha1_1 - alpha0_1;
5562 3139 : const double dfDeltaAlpha21 = alpha2_1 - alpha1_1;
5563 : const double dfMaxDeltaAlpha =
5564 3139 : std::max(fabs(dfDeltaAlpha10), fabs(dfDeltaAlpha21));
5565 : GUInt32 nAlphaRatioRef =
5566 3139 : OGRGF_GetHiddenValue(p1.getX(), p1.getY()) |
5567 3139 : (OGRGF_GetHiddenValue(p2.getX(), p2.getY()) << HIDDEN_ALPHA_HALF_WIDTH);
5568 3139 : bool bFoundFFFFFFFFPattern = false;
5569 3139 : bool bFoundReversedAlphaRatioRef = false;
5570 3139 : bool bValidAlphaRatio = nAlphaRatioRef > 0 && nAlphaRatioRef < 0xFFFFFFFF;
5571 3139 : int nCountValidAlphaRatio = 1;
5572 :
5573 3139 : double dfScale = std::max(1.0, R_1);
5574 3139 : dfScale = std::max(dfScale, fabs(cx_1));
5575 3139 : dfScale = std::max(dfScale, fabs(cy_1));
5576 3139 : dfScale = pow(10.0, ceil(log10(dfScale)));
5577 3139 : const double dfInvScale = 1.0 / dfScale;
5578 :
5579 3139 : const int bInitialConstantStep =
5580 3139 : (fabs(dfDeltaAlpha10 - dfDeltaAlpha21) / dfMaxDeltaAlpha) < 1.0e-4;
5581 3139 : const double dfDeltaEpsilon =
5582 3139 : bInitialConstantStep ? dfMaxDeltaAlpha * 1e-4 : dfMaxDeltaAlpha / 10;
5583 :
5584 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5585 : printf("----------------------------\n"); /*ok*/
5586 : printf("Curve beginning at offset i = %d\n", i); /*ok*/
5587 : printf("Initial alpha ratio = %u\n", nAlphaRatioRef); /*ok*/
5588 : /*ok*/ printf("Initial R = %.16g, cx = %.16g, cy = %.16g\n", R_1, cx_1,
5589 : cy_1);
5590 : printf("dfScale = %f\n", dfScale); /*ok*/
5591 : printf("bInitialConstantStep = %d, " /*ok*/
5592 : "fabs(dfDeltaAlpha10 - dfDeltaAlpha21)=%.8g, "
5593 : "dfMaxDeltaAlpha = %.8f, "
5594 : "dfDeltaEpsilon = %.8f (%.8f)\n",
5595 : bInitialConstantStep, fabs(dfDeltaAlpha10 - dfDeltaAlpha21),
5596 : dfMaxDeltaAlpha, dfDeltaEpsilon, 1.0 / 180.0 * M_PI);
5597 : #endif
5598 3139 : int iMidPoint = -1;
5599 3139 : double dfLastValidAlpha = alpha2_1;
5600 :
5601 3139 : double dfLastLogRelDiff = 0;
5602 :
5603 6278 : OGRPoint p3;
5604 3139 : int j = i + 1; // Used after for.
5605 182753 : for (; j + 2 < poLS->getNumPoints(); j++)
5606 : {
5607 179713 : poLS->getPoint(j, &p1);
5608 179713 : poLS->getPoint(j + 1, &p2);
5609 179713 : poLS->getPoint(j + 2, &p3);
5610 179713 : double R_2 = 0.0;
5611 179713 : double cx_2 = 0.0;
5612 179713 : double cy_2 = 0.0;
5613 179713 : double alpha0_2 = 0.0;
5614 179713 : double alpha1_2 = 0.0;
5615 179713 : double alpha2_2 = 0.0;
5616 : // Check that the new candidate arc shares the same
5617 : // radius, center and winding order.
5618 179713 : if (!(OGRGeometryFactory::GetCurveParameters(
5619 : p1.getX(), p1.getY(), p2.getX(), p2.getY(), p3.getX(),
5620 : p3.getY(), R_2, cx_2, cy_2, alpha0_2, alpha1_2, alpha2_2)))
5621 : {
5622 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5623 : printf("End of curve at j=%d\n : straight line", j); /*ok*/
5624 : #endif
5625 99 : break;
5626 : }
5627 :
5628 179705 : const double dfRelDiffR = fabs(R_1 - R_2) * dfInvScale;
5629 179705 : const double dfRelDiffCx = fabs(cx_1 - cx_2) * dfInvScale;
5630 179705 : const double dfRelDiffCy = fabs(cy_1 - cy_2) * dfInvScale;
5631 :
5632 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5633 : printf("j=%d: R = %.16g, cx = %.16g, cy = %.16g, " /*ok*/
5634 : "rel_diff_R=%.8g rel_diff_cx=%.8g rel_diff_cy=%.8g\n",
5635 : j, R_2, cx_2, cy_2, dfRelDiffR, dfRelDiffCx, dfRelDiffCy);
5636 : #endif
5637 :
5638 179705 : if (dfRelDiffR > 1.0e-7 || dfRelDiffCx > 1.0e-7 ||
5639 179636 : dfRelDiffCy > 1.0e-7 ||
5640 179636 : dfDeltaAlpha10 * (alpha1_2 - alpha0_2) < 0.0)
5641 : {
5642 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5643 : printf("End of curve at j=%d\n", j); /*ok*/
5644 : #endif
5645 : break;
5646 : }
5647 :
5648 179636 : if (dfRelDiffR > 0.0 && dfRelDiffCx > 0.0 && dfRelDiffCy > 0.0)
5649 : {
5650 : const double dfLogRelDiff = std::min(
5651 359242 : std::min(fabs(log10(dfRelDiffR)), fabs(log10(dfRelDiffCx))),
5652 179621 : fabs(log10(dfRelDiffCy)));
5653 : // printf("dfLogRelDiff = %f, dfLastLogRelDiff=%f, "/*ok*/
5654 : // "dfLogRelDiff - dfLastLogRelDiff=%f\n",
5655 : // dfLogRelDiff, dfLastLogRelDiff,
5656 : // dfLogRelDiff - dfLastLogRelDiff);
5657 179621 : if (dfLogRelDiff > 0.0 && dfLastLogRelDiff >= 8.0 &&
5658 2 : dfLogRelDiff <= 8.0 && dfLogRelDiff < dfLastLogRelDiff - 2.0)
5659 : {
5660 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5661 : printf("End of curve at j=%d. Significant different in " /*ok*/
5662 : "relative error w.r.t previous points\n",
5663 : j);
5664 : #endif
5665 2 : break;
5666 : }
5667 179619 : dfLastLogRelDiff = dfLogRelDiff;
5668 : }
5669 :
5670 179634 : const double dfStep10 = fabs(alpha1_2 - alpha0_2);
5671 179634 : const double dfStep21 = fabs(alpha2_2 - alpha1_2);
5672 : // Check that the angle step is consistent with the original step.
5673 179634 : if (!(dfStep10 < 2.0 * dfMaxDeltaAlpha &&
5674 179634 : dfStep21 < 2.0 * dfMaxDeltaAlpha))
5675 : {
5676 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5677 : printf("End of curve at j=%d: dfStep10=%f, dfStep21=%f, " /*ok*/
5678 : "2*dfMaxDeltaAlpha=%f\n",
5679 : j, dfStep10, dfStep21, 2 * dfMaxDeltaAlpha);
5680 : #endif
5681 : break;
5682 : }
5683 :
5684 179633 : if (bValidAlphaRatio && j > i + 1 && (i % 2) != (j % 2))
5685 : {
5686 : const GUInt32 nAlphaRatioReversed =
5687 88051 : (OGRGF_GetHiddenValue(p1.getX(), p1.getY())
5688 176102 : << HIDDEN_ALPHA_HALF_WIDTH) |
5689 88051 : (OGRGF_GetHiddenValue(p2.getX(), p2.getY()));
5690 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5691 : printf("j=%d, nAlphaRatioReversed = %u\n", /*ok*/
5692 : j, nAlphaRatioReversed);
5693 : #endif
5694 88051 : if (!bFoundFFFFFFFFPattern && nAlphaRatioReversed == 0xFFFFFFFF)
5695 : {
5696 3067 : bFoundFFFFFFFFPattern = true;
5697 3067 : nCountValidAlphaRatio++;
5698 : }
5699 84984 : else if (bFoundFFFFFFFFPattern && !bFoundReversedAlphaRatioRef &&
5700 : nAlphaRatioReversed == 0xFFFFFFFF)
5701 : {
5702 81890 : nCountValidAlphaRatio++;
5703 : }
5704 3094 : else if (bFoundFFFFFFFFPattern && !bFoundReversedAlphaRatioRef &&
5705 : nAlphaRatioReversed == nAlphaRatioRef)
5706 : {
5707 3067 : bFoundReversedAlphaRatioRef = true;
5708 3067 : nCountValidAlphaRatio++;
5709 : }
5710 : else
5711 : {
5712 27 : if (bInitialConstantStep &&
5713 26 : fabs(dfLastValidAlpha - alpha0_1) >= M_PI &&
5714 : nCountValidAlphaRatio > 10)
5715 : {
5716 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5717 : printf("End of curve at j=%d: " /*ok*/
5718 : "fabs(dfLastValidAlpha - alpha0_1)=%f, "
5719 : "nCountValidAlphaRatio=%d\n",
5720 : j, fabs(dfLastValidAlpha - alpha0_1),
5721 : nCountValidAlphaRatio);
5722 : #endif
5723 19 : if (dfLastValidAlpha - alpha0_1 > 0)
5724 : {
5725 21 : while (dfLastValidAlpha - alpha0_1 - dfMaxDeltaAlpha -
5726 14 : M_PI >
5727 14 : -dfMaxDeltaAlpha / 10)
5728 : {
5729 7 : dfLastValidAlpha -= dfMaxDeltaAlpha;
5730 7 : j--;
5731 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5732 : printf(/*ok*/
5733 : "--> corrected as fabs(dfLastValidAlpha - "
5734 : "alpha0_1)=%f, j=%d\n",
5735 : fabs(dfLastValidAlpha - alpha0_1), j);
5736 : #endif
5737 : }
5738 : }
5739 : else
5740 : {
5741 36 : while (dfLastValidAlpha - alpha0_1 + dfMaxDeltaAlpha +
5742 24 : M_PI <
5743 24 : dfMaxDeltaAlpha / 10)
5744 : {
5745 12 : dfLastValidAlpha += dfMaxDeltaAlpha;
5746 12 : j--;
5747 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5748 : printf(/*ok*/
5749 : "--> corrected as fabs(dfLastValidAlpha - "
5750 : "alpha0_1)=%f, j=%d\n",
5751 : fabs(dfLastValidAlpha - alpha0_1), j);
5752 : #endif
5753 : }
5754 : }
5755 19 : poLS->getPoint(j + 1, &p2);
5756 19 : break;
5757 : }
5758 :
5759 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5760 : printf("j=%d, nAlphaRatioReversed = %u --> inconsistent " /*ok*/
5761 : "values across arc. Don't use it\n",
5762 : j, nAlphaRatioReversed);
5763 : #endif
5764 8 : bValidAlphaRatio = false;
5765 : }
5766 : }
5767 :
5768 : // Correct current end angle, consistently with start angle.
5769 179614 : dfLastValidAlpha = OGRGF_FixAngle(alpha0_1, alpha1_1, alpha2_2);
5770 :
5771 : // Try to detect the precise intermediate point of the
5772 : // arc circle by detecting irregular angle step
5773 : // This is OK if we don't detect the right point or fail
5774 : // to detect it.
5775 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5776 : printf("j=%d A(0,1)-maxDelta=%.8f A(1,2)-maxDelta=%.8f " /*ok*/
5777 : "x1=%.8f y1=%.8f x2=%.8f y2=%.8f x3=%.8f y3=%.8f\n",
5778 : j, fabs(dfStep10 - dfMaxDeltaAlpha),
5779 : fabs(dfStep21 - dfMaxDeltaAlpha), p1.getX(), p1.getY(),
5780 : p2.getX(), p2.getY(), p3.getX(), p3.getY());
5781 : #endif
5782 179614 : if (j > i + 1 && iMidPoint < 0 && dfDeltaEpsilon < 1.0 / 180.0 * M_PI)
5783 : {
5784 176143 : if (fabs(dfStep10 - dfMaxDeltaAlpha) > dfDeltaEpsilon)
5785 8 : iMidPoint = j + ((bInitialConstantStep) ? 0 : 1);
5786 176135 : else if (fabs(dfStep21 - dfMaxDeltaAlpha) > dfDeltaEpsilon)
5787 4 : iMidPoint = j + ((bInitialConstantStep) ? 1 : 2);
5788 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5789 : if (iMidPoint >= 0)
5790 : {
5791 : OGRPoint pMid;
5792 : poLS->getPoint(iMidPoint, &pMid);
5793 : printf("Midpoint detected at j = %d, iMidPoint = %d, " /*ok*/
5794 : "x=%.8f y=%.8f\n",
5795 : j, iMidPoint, pMid.getX(), pMid.getY());
5796 : }
5797 : #endif
5798 : }
5799 : }
5800 :
5801 : // Take a minimum threshold of consecutive points
5802 : // on the arc to avoid false positives.
5803 3139 : if (j < i + 3)
5804 61 : return -1;
5805 :
5806 3078 : bValidAlphaRatio &= bFoundFFFFFFFFPattern && bFoundReversedAlphaRatioRef;
5807 :
5808 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5809 : printf("bValidAlphaRatio=%d bFoundFFFFFFFFPattern=%d, " /*ok*/
5810 : "bFoundReversedAlphaRatioRef=%d\n",
5811 : static_cast<int>(bValidAlphaRatio),
5812 : static_cast<int>(bFoundFFFFFFFFPattern),
5813 : static_cast<int>(bFoundReversedAlphaRatioRef));
5814 : printf("alpha0_1=%f dfLastValidAlpha=%f\n", /*ok*/
5815 : alpha0_1, dfLastValidAlpha);
5816 : #endif
5817 :
5818 3078 : if (poLSNew != nullptr)
5819 : {
5820 11 : double dfScale2 = std::max(1.0, fabs(p0.getX()));
5821 11 : dfScale2 = std::max(dfScale2, fabs(p0.getY()));
5822 : // Not strictly necessary, but helps having 'clean' lines without
5823 : // duplicated points.
5824 11 : constexpr double dfToleranceEps =
5825 : OGRCompoundCurve::DEFAULT_TOLERANCE_EPSILON;
5826 11 : if (fabs(poLSNew->getX(poLSNew->getNumPoints() - 1) - p0.getX()) >
5827 12 : dfToleranceEps * dfScale2 ||
5828 1 : fabs(poLSNew->getY(poLSNew->getNumPoints() - 1) - p0.getY()) >
5829 1 : dfToleranceEps * dfScale2)
5830 10 : poLSNew->addPoint(&p0);
5831 11 : if (poLSNew->getNumPoints() >= 2)
5832 : {
5833 10 : if (poCC == nullptr)
5834 3 : poCC = new OGRCompoundCurve();
5835 10 : poCC->addCurveDirectly(poLSNew);
5836 : }
5837 : else
5838 1 : delete poLSNew;
5839 11 : poLSNew = nullptr;
5840 : }
5841 :
5842 3078 : if (poCS == nullptr)
5843 : {
5844 3054 : poCS = new OGRCircularString();
5845 3054 : poCS->addPoint(&p0);
5846 : }
5847 :
5848 3078 : OGRPoint *poFinalPoint = (j + 2 >= poLS->getNumPoints()) ? &p3 : &p2;
5849 :
5850 3078 : double dfXMid = 0.0;
5851 3078 : double dfYMid = 0.0;
5852 3078 : double dfZMid = 0.0;
5853 3078 : if (bValidAlphaRatio)
5854 : {
5855 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5856 : printf("Using alpha ratio...\n"); /*ok*/
5857 : #endif
5858 3067 : double dfAlphaMid = 0.0;
5859 3067 : if (OGRGF_NeedSwithArcOrder(p0.getX(), p0.getY(), poFinalPoint->getX(),
5860 : poFinalPoint->getY()))
5861 : {
5862 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5863 : printf("Switching angles\n"); /*ok*/
5864 : #endif
5865 1543 : dfAlphaMid = dfLastValidAlpha + nAlphaRatioRef *
5866 1543 : (alpha0_1 - dfLastValidAlpha) /
5867 : HIDDEN_ALPHA_SCALE;
5868 1543 : dfAlphaMid = OGRGF_FixAngle(alpha0_1, dfLastValidAlpha, dfAlphaMid);
5869 : }
5870 : else
5871 : {
5872 1524 : dfAlphaMid = alpha0_1 + nAlphaRatioRef *
5873 1524 : (dfLastValidAlpha - alpha0_1) /
5874 : HIDDEN_ALPHA_SCALE;
5875 : }
5876 :
5877 3067 : dfXMid = cx_1 + R_1 * cos(dfAlphaMid);
5878 3067 : dfYMid = cy_1 + R_1 * sin(dfAlphaMid);
5879 :
5880 3067 : if (poLS->getCoordinateDimension() == 3)
5881 : {
5882 2 : double dfLastAlpha = 0.0;
5883 2 : double dfLastZ = 0.0;
5884 2 : int k = i; // Used after for.
5885 48 : for (; k < j + 2; k++)
5886 : {
5887 48 : OGRPoint p;
5888 48 : poLS->getPoint(k, &p);
5889 48 : double dfAlpha = atan2(p.getY() - cy_1, p.getX() - cx_1);
5890 48 : dfAlpha = OGRGF_FixAngle(alpha0_1, dfLastValidAlpha, dfAlpha);
5891 48 : if (k > i &&
5892 46 : ((dfAlpha < dfLastValidAlpha && dfAlphaMid < dfAlpha) ||
5893 23 : (dfAlpha > dfLastValidAlpha && dfAlphaMid > dfAlpha)))
5894 : {
5895 2 : const double dfRatio =
5896 2 : (dfAlphaMid - dfLastAlpha) / (dfAlpha - dfLastAlpha);
5897 2 : dfZMid = (1 - dfRatio) * dfLastZ + dfRatio * p.getZ();
5898 2 : break;
5899 : }
5900 46 : dfLastAlpha = dfAlpha;
5901 46 : dfLastZ = p.getZ();
5902 : }
5903 2 : if (k == j + 2)
5904 0 : dfZMid = dfLastZ;
5905 2 : if (IS_ALMOST_INTEGER(dfZMid))
5906 2 : dfZMid = static_cast<int>(floor(dfZMid + 0.5));
5907 : }
5908 :
5909 : // A few rounding strategies in case the mid point was at "exact"
5910 : // coordinates.
5911 3067 : if (R_1 > 1e-5)
5912 : {
5913 : const bool bStartEndInteger =
5914 9161 : IS_ALMOST_INTEGER(p0.getX()) && IS_ALMOST_INTEGER(p0.getY()) &&
5915 9161 : IS_ALMOST_INTEGER(poFinalPoint->getX()) &&
5916 3048 : IS_ALMOST_INTEGER(poFinalPoint->getY());
5917 3061 : if (bStartEndInteger &&
5918 3048 : fabs(dfXMid - floor(dfXMid + 0.5)) / dfScale < 1e-4 &&
5919 3029 : fabs(dfYMid - floor(dfYMid + 0.5)) / dfScale < 1e-4)
5920 : {
5921 3029 : dfXMid = static_cast<int>(floor(dfXMid + 0.5));
5922 3029 : dfYMid = static_cast<int>(floor(dfYMid + 0.5));
5923 : // Sometimes rounding to closest is not best approach
5924 : // Try neighbouring integers to look for the one that
5925 : // minimize the error w.r.t to the arc center
5926 : // But only do that if the radius is greater than
5927 : // the magnitude of the delta that we will try!
5928 : double dfBestRError =
5929 3029 : fabs(R_1 - DISTANCE(dfXMid, dfYMid, cx_1, cy_1));
5930 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5931 : printf("initial_error=%f\n", dfBestRError); /*ok*/
5932 : #endif
5933 3029 : int iBestX = 0;
5934 3029 : int iBestY = 0;
5935 3029 : if (dfBestRError > 0.001 && R_1 > 2)
5936 : {
5937 3 : int nSearchRadius = 1;
5938 : // Extend the search radius if the arc circle radius
5939 : // is much higher than the coordinate values.
5940 : double dfMaxCoords =
5941 3 : std::max(fabs(p0.getX()), fabs(p0.getY()));
5942 3 : dfMaxCoords = std::max(dfMaxCoords, poFinalPoint->getX());
5943 3 : dfMaxCoords = std::max(dfMaxCoords, poFinalPoint->getY());
5944 3 : dfMaxCoords = std::max(dfMaxCoords, dfXMid);
5945 3 : dfMaxCoords = std::max(dfMaxCoords, dfYMid);
5946 3 : if (R_1 > dfMaxCoords * 1000)
5947 3 : nSearchRadius = 100;
5948 0 : else if (R_1 > dfMaxCoords * 10)
5949 0 : nSearchRadius = 10;
5950 606 : for (int iY = -nSearchRadius; iY <= nSearchRadius; iY++)
5951 : {
5952 121806 : for (int iX = -nSearchRadius; iX <= nSearchRadius; iX++)
5953 : {
5954 121203 : const double dfCandidateX = dfXMid + iX;
5955 121203 : const double dfCandidateY = dfYMid + iY;
5956 121203 : if (fabs(dfCandidateX - p0.getX()) < 1e-8 &&
5957 0 : fabs(dfCandidateY - p0.getY()) < 1e-8)
5958 0 : continue;
5959 121203 : if (fabs(dfCandidateX - poFinalPoint->getX()) <
5960 121203 : 1e-8 &&
5961 0 : fabs(dfCandidateY - poFinalPoint->getY()) <
5962 : 1e-8)
5963 0 : continue;
5964 : const double dfRError =
5965 121203 : fabs(R_1 - DISTANCE(dfCandidateX, dfCandidateY,
5966 121203 : cx_1, cy_1));
5967 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5968 : printf("x=%d y=%d error=%f besterror=%f\n", /*ok*/
5969 : static_cast<int>(dfXMid + iX),
5970 : static_cast<int>(dfYMid + iY), dfRError,
5971 : dfBestRError);
5972 : #endif
5973 121203 : if (dfRError < dfBestRError)
5974 : {
5975 20 : iBestX = iX;
5976 20 : iBestY = iY;
5977 20 : dfBestRError = dfRError;
5978 : }
5979 : }
5980 : }
5981 : }
5982 3029 : dfXMid += iBestX;
5983 3029 : dfYMid += iBestY;
5984 : }
5985 : else
5986 : {
5987 : // Limit the number of significant figures in decimal
5988 : // representation.
5989 32 : if (fabs(dfXMid) < 100000000.0)
5990 : {
5991 32 : dfXMid =
5992 32 : static_cast<GIntBig>(floor(dfXMid * 100000000 + 0.5)) /
5993 : 100000000.0;
5994 : }
5995 32 : if (fabs(dfYMid) < 100000000.0)
5996 : {
5997 32 : dfYMid =
5998 32 : static_cast<GIntBig>(floor(dfYMid * 100000000 + 0.5)) /
5999 : 100000000.0;
6000 : }
6001 : }
6002 : }
6003 :
6004 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6005 : printf("dfAlphaMid=%f, x_mid = %f, y_mid = %f\n", /*ok*/
6006 : dfLastValidAlpha, dfXMid, dfYMid);
6007 : #endif
6008 : }
6009 :
6010 : // If this is a full circle of a non-polygonal zone, we must
6011 : // use a 5-point representation to keep the winding order.
6012 3089 : if (p0.Equals(poFinalPoint) &&
6013 11 : !EQUAL(poLS->getGeometryName(), "LINEARRING"))
6014 : {
6015 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6016 : printf("Full circle of a non-polygonal zone\n"); /*ok*/
6017 : #endif
6018 1 : poLS->getPoint((i + j + 2) / 4, &p1);
6019 1 : poCS->addPoint(&p1);
6020 1 : if (bValidAlphaRatio)
6021 : {
6022 1 : p1.setX(dfXMid);
6023 1 : p1.setY(dfYMid);
6024 1 : if (poLS->getCoordinateDimension() == 3)
6025 0 : p1.setZ(dfZMid);
6026 : }
6027 : else
6028 : {
6029 0 : poLS->getPoint((i + j + 1) / 2, &p1);
6030 : }
6031 1 : poCS->addPoint(&p1);
6032 1 : poLS->getPoint(3 * (i + j + 2) / 4, &p1);
6033 1 : poCS->addPoint(&p1);
6034 : }
6035 :
6036 3077 : else if (bValidAlphaRatio)
6037 : {
6038 3066 : p1.setX(dfXMid);
6039 3066 : p1.setY(dfYMid);
6040 3066 : if (poLS->getCoordinateDimension() == 3)
6041 2 : p1.setZ(dfZMid);
6042 3066 : poCS->addPoint(&p1);
6043 : }
6044 :
6045 : // If we have found a candidate for a precise intermediate
6046 : // point, use it.
6047 11 : else if (iMidPoint >= 1 && iMidPoint < j)
6048 : {
6049 3 : poLS->getPoint(iMidPoint, &p1);
6050 3 : poCS->addPoint(&p1);
6051 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6052 : printf("Using detected midpoint...\n"); /*ok*/
6053 : printf("x_mid = %f, y_mid = %f\n", p1.getX(), p1.getY()); /*ok*/
6054 : #endif
6055 : }
6056 : // Otherwise pick up the mid point between both extremities.
6057 : else
6058 : {
6059 8 : poLS->getPoint((i + j + 1) / 2, &p1);
6060 8 : poCS->addPoint(&p1);
6061 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6062 : printf("Pickup 'random' midpoint at index=%d...\n", /*ok*/
6063 : (i + j + 1) / 2);
6064 : printf("x_mid = %f, y_mid = %f\n", p1.getX(), p1.getY()); /*ok*/
6065 : #endif
6066 : }
6067 3078 : poCS->addPoint(poFinalPoint);
6068 :
6069 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6070 : printf("----------------------------\n"); /*ok*/
6071 : #endif
6072 :
6073 3078 : if (j + 2 >= poLS->getNumPoints())
6074 3040 : return -2;
6075 38 : return j + 1;
6076 : }
6077 :
6078 : /************************************************************************/
6079 : /* curveFromLineString() */
6080 : /************************************************************************/
6081 :
6082 : /**
6083 : * \brief Try to convert a linestring approximating curves into a curve.
6084 : *
6085 : * This method can return a COMPOUNDCURVE, a CIRCULARSTRING or a LINESTRING.
6086 : *
6087 : * This method is the reverse of curveFromLineString().
6088 : *
6089 : * @param poLS handle to the geometry to convert.
6090 : * @param papszOptions options as a null-terminated list of strings.
6091 : * Unused for now. Must be set to NULL.
6092 : *
6093 : * @return the converted geometry (ownership to caller).
6094 : *
6095 : * @since GDAL 2.0
6096 : */
6097 :
6098 3190 : OGRCurve *OGRGeometryFactory::curveFromLineString(
6099 : const OGRLineString *poLS, CPL_UNUSED const char *const *papszOptions)
6100 : {
6101 3190 : OGRCompoundCurve *poCC = nullptr;
6102 3190 : OGRCircularString *poCS = nullptr;
6103 3190 : OGRLineString *poLSNew = nullptr;
6104 3190 : const int nLSNumPoints = poLS->getNumPoints();
6105 3190 : const bool bIsClosed = nLSNumPoints >= 4 && poLS->get_IsClosed();
6106 3618 : for (int i = 0; i < nLSNumPoints; /* nothing */)
6107 : {
6108 3468 : const int iNewI = OGRGF_DetectArc(poLS, i, poCC, poCS, poLSNew);
6109 3468 : if (iNewI == -2)
6110 3040 : break;
6111 428 : if (iNewI >= 0)
6112 : {
6113 38 : i = iNewI;
6114 38 : continue;
6115 : }
6116 :
6117 390 : if (poCS != nullptr)
6118 : {
6119 14 : if (poCC == nullptr)
6120 5 : poCC = new OGRCompoundCurve();
6121 14 : poCC->addCurveDirectly(poCS);
6122 14 : poCS = nullptr;
6123 : }
6124 :
6125 390 : OGRPoint p;
6126 390 : poLS->getPoint(i, &p);
6127 390 : if (poLSNew == nullptr)
6128 : {
6129 160 : poLSNew = new OGRLineString();
6130 160 : poLSNew->addPoint(&p);
6131 : }
6132 : // Not strictly necessary, but helps having 'clean' lines without
6133 : // duplicated points.
6134 : else
6135 : {
6136 230 : double dfScale = std::max(1.0, fabs(p.getX()));
6137 230 : dfScale = std::max(dfScale, fabs(p.getY()));
6138 230 : if (bIsClosed && i == nLSNumPoints - 1)
6139 7 : dfScale = 0;
6140 230 : constexpr double dfToleranceEps =
6141 : OGRCompoundCurve::DEFAULT_TOLERANCE_EPSILON;
6142 230 : if (fabs(poLSNew->getX(poLSNew->getNumPoints() - 1) - p.getX()) >
6143 239 : dfToleranceEps * dfScale ||
6144 9 : fabs(poLSNew->getY(poLSNew->getNumPoints() - 1) - p.getY()) >
6145 9 : dfToleranceEps * dfScale)
6146 : {
6147 229 : poLSNew->addPoint(&p);
6148 : }
6149 : }
6150 :
6151 390 : i++;
6152 : }
6153 :
6154 3190 : OGRCurve *poRet = nullptr;
6155 :
6156 3190 : if (poLSNew != nullptr && poLSNew->getNumPoints() < 2)
6157 : {
6158 1 : delete poLSNew;
6159 1 : poLSNew = nullptr;
6160 1 : if (poCC != nullptr)
6161 : {
6162 1 : if (poCC->getNumCurves() == 1)
6163 : {
6164 1 : poRet = poCC->stealCurve(0);
6165 1 : delete poCC;
6166 1 : poCC = nullptr;
6167 : }
6168 : else
6169 0 : poRet = poCC;
6170 : }
6171 : else
6172 0 : poRet = poLS->clone();
6173 : }
6174 3189 : else if (poCC != nullptr)
6175 : {
6176 7 : if (poLSNew)
6177 6 : poCC->addCurveDirectly(poLSNew);
6178 : else
6179 1 : poCC->addCurveDirectly(poCS);
6180 7 : poRet = poCC;
6181 : }
6182 3182 : else if (poLSNew != nullptr)
6183 142 : poRet = poLSNew;
6184 3040 : else if (poCS != nullptr)
6185 3039 : poRet = poCS;
6186 : else
6187 1 : poRet = poLS->clone();
6188 :
6189 3190 : poRet->assignSpatialReference(poLS->getSpatialReference());
6190 :
6191 3190 : return poRet;
6192 : }
6193 :
6194 : /************************************************************************/
6195 : /* createFromGeoJson( const char* ) */
6196 : /************************************************************************/
6197 :
6198 : /**
6199 : * @brief Create geometry from GeoJson fragment.
6200 : * @param pszJsonString The GeoJSON fragment for the geometry.
6201 : * @param nSize (new in GDAL 3.4) Optional length of the string
6202 : * if it is not null-terminated
6203 : * @return a geometry on success, or NULL on error.
6204 : * @since GDAL 2.3
6205 : */
6206 0 : OGRGeometry *OGRGeometryFactory::createFromGeoJson(const char *pszJsonString,
6207 : int nSize)
6208 : {
6209 0 : CPLJSONDocument oDocument;
6210 0 : if (!oDocument.LoadMemory(reinterpret_cast<const GByte *>(pszJsonString),
6211 : nSize))
6212 : {
6213 0 : return nullptr;
6214 : }
6215 :
6216 0 : return createFromGeoJson(oDocument.GetRoot());
6217 : }
6218 :
6219 : /************************************************************************/
6220 : /* createFromGeoJson( const CPLJSONObject& ) */
6221 : /************************************************************************/
6222 :
6223 : /**
6224 : * @brief Create geometry from GeoJson fragment.
6225 : * @param oJsonObject The JSONObject class describes the GeoJSON geometry.
6226 : * @return a geometry on success, or NULL on error.
6227 : * @since GDAL 2.3
6228 : */
6229 : OGRGeometry *
6230 0 : OGRGeometryFactory::createFromGeoJson(const CPLJSONObject &oJsonObject)
6231 : {
6232 0 : if (!oJsonObject.IsValid())
6233 : {
6234 0 : return nullptr;
6235 : }
6236 :
6237 : // TODO: Move from GeoJSON driver functions create geometry here, and
6238 : // replace json-c specific json_object to CPLJSONObject
6239 0 : return OGRGeoJSONReadGeometry(
6240 0 : static_cast<json_object *>(oJsonObject.GetInternalHandle()));
6241 : }
|
