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