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 57316 : OGRErr OGRGeometryFactory::createFromWkb(const void *pabyData,
90 : const OGRSpatialReference *poSR,
91 : OGRGeometry **ppoReturn, size_t nBytes,
92 : OGRwkbVariant eWkbVariant)
93 :
94 : {
95 57316 : size_t nBytesConsumedOutIgnored = 0;
96 57316 : return createFromWkb(pabyData, poSR, ppoReturn, nBytes, eWkbVariant,
97 114628 : 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 94974 : 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 94974 : const GByte *l_pabyData = static_cast<const GByte *>(pabyData);
142 94974 : nBytesConsumedOut = 0;
143 94974 : *ppoReturn = nullptr;
144 :
145 94974 : if (nBytes < 9 && nBytes != static_cast<size_t>(-1))
146 1387 : 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 93587 : const int nByteOrder = DB2_V72_FIX_BYTE_ORDER(*l_pabyData);
153 93587 : if (nByteOrder != wkbXDR && nByteOrder != wkbNDR)
154 : {
155 354 : CPLDebug("OGR",
156 : "OGRGeometryFactory::createFromWkb() - got corrupt data.\n"
157 : "%02X%02X%02X%02X%02X%02X%02X%02X%02X",
158 354 : l_pabyData[0], l_pabyData[1], l_pabyData[2], l_pabyData[3],
159 354 : l_pabyData[4], l_pabyData[5], l_pabyData[6], l_pabyData[7],
160 354 : l_pabyData[8]);
161 354 : 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 93233 : OGRwkbGeometryType eGeometryType = wkbUnknown;
171 : const OGRErr err =
172 93233 : OGRReadWKBGeometryType(l_pabyData, eWkbVariant, &eGeometryType);
173 :
174 93231 : 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 92668 : OGRGeometry *poGeom = createGeometry(eGeometryType);
182 :
183 92669 : if (poGeom == nullptr)
184 0 : return OGRERR_UNSUPPORTED_GEOMETRY_TYPE;
185 :
186 : /* -------------------------------------------------------------------- */
187 : /* Import from binary. */
188 : /* -------------------------------------------------------------------- */
189 185339 : const OGRErr eErr = poGeom->importFromWkb(l_pabyData, nBytes, eWkbVariant,
190 92669 : nBytesConsumedOut);
191 92670 : if (eErr != OGRERR_NONE)
192 : {
193 7315 : delete poGeom;
194 7315 : return eErr;
195 : }
196 :
197 : /* -------------------------------------------------------------------- */
198 : /* Assign spatial reference system. */
199 : /* -------------------------------------------------------------------- */
200 :
201 88916 : 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 85353 : poGeom->assignSpatialReference(poSR);
209 85353 : *ppoReturn = poGeom;
210 :
211 85353 : 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 124169 : OGRErr OGRGeometryFactory::createFromWkt(const char **ppszData,
340 : const OGRSpatialReference *poSR,
341 : OGRGeometry **ppoReturn)
342 :
343 : {
344 124169 : const char *pszInput = *ppszData;
345 124169 : *ppoReturn = nullptr;
346 :
347 : /* -------------------------------------------------------------------- */
348 : /* Get the first token, which should be the geometry type. */
349 : /* -------------------------------------------------------------------- */
350 124169 : char szToken[OGR_WKT_TOKEN_MAX] = {};
351 124169 : if (OGRWktReadToken(pszInput, szToken) == nullptr)
352 0 : return OGRERR_CORRUPT_DATA;
353 :
354 : /* -------------------------------------------------------------------- */
355 : /* Instantiate a geometry of the appropriate type. */
356 : /* -------------------------------------------------------------------- */
357 124169 : OGRGeometry *poGeom = nullptr;
358 124169 : if (STARTS_WITH_CI(szToken, "POINT"))
359 : {
360 97737 : poGeom = new OGRPoint();
361 : }
362 26432 : else if (STARTS_WITH_CI(szToken, "LINESTRING"))
363 : {
364 1673 : poGeom = new OGRLineString();
365 : }
366 24759 : else if (STARTS_WITH_CI(szToken, "POLYGON"))
367 : {
368 16122 : poGeom = new OGRPolygon();
369 : }
370 8637 : else if (STARTS_WITH_CI(szToken, "TRIANGLE"))
371 : {
372 62 : poGeom = new OGRTriangle();
373 : }
374 8575 : else if (STARTS_WITH_CI(szToken, "GEOMETRYCOLLECTION"))
375 : {
376 543 : poGeom = new OGRGeometryCollection();
377 : }
378 8032 : else if (STARTS_WITH_CI(szToken, "MULTIPOLYGON"))
379 : {
380 980 : poGeom = new OGRMultiPolygon();
381 : }
382 7052 : else if (STARTS_WITH_CI(szToken, "MULTIPOINT"))
383 : {
384 616 : poGeom = new OGRMultiPoint();
385 : }
386 6436 : else if (STARTS_WITH_CI(szToken, "MULTILINESTRING"))
387 : {
388 667 : poGeom = new OGRMultiLineString();
389 : }
390 5769 : else if (STARTS_WITH_CI(szToken, "CIRCULARSTRING"))
391 : {
392 3510 : 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 123093 : const OGRErr eErr = poGeom->importFromWkt(&pszInput);
430 :
431 : /* -------------------------------------------------------------------- */
432 : /* Assign spatial reference system. */
433 : /* -------------------------------------------------------------------- */
434 123093 : if (eErr == OGRERR_NONE)
435 : {
436 127249 : if (poGeom->hasCurveGeometry() &&
437 4395 : CPLTestBool(CPLGetConfigOption("OGR_STROKE_CURVE", "FALSE")))
438 : {
439 9 : OGRGeometry *poNewGeom = poGeom->getLinearGeometry();
440 9 : delete poGeom;
441 9 : poGeom = poNewGeom;
442 : }
443 122854 : poGeom->assignSpatialReference(poSR);
444 122854 : *ppoReturn = poGeom;
445 122854 : *ppszData = pszInput;
446 : }
447 : else
448 : {
449 239 : delete poGeom;
450 : }
451 :
452 123093 : 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 5747 : OGRErr OGRGeometryFactory::createFromWkt(const char *pszData,
477 : const OGRSpatialReference *poSR,
478 : OGRGeometry **ppoReturn)
479 :
480 : {
481 5747 : 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 116653 : OGRErr CPL_DLL OGR_G_CreateFromWkt(char **ppszData, OGRSpatialReferenceH hSRS,
541 : OGRGeometryH *phGeometry)
542 :
543 : {
544 116653 : return OGRGeometryFactory::createFromWkt(
545 : const_cast<const char **>(ppszData),
546 116653 : OGRSpatialReference::FromHandle(hSRS),
547 116653 : 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 258895 : OGRGeometryFactory::createGeometry(OGRwkbGeometryType eGeometryType)
571 :
572 : {
573 258895 : OGRGeometry *poGeom = nullptr;
574 258895 : switch (wkbFlatten(eGeometryType))
575 : {
576 182994 : case wkbPoint:
577 365988 : poGeom = new (std::nothrow) OGRPoint();
578 182994 : break;
579 :
580 7626 : case wkbLineString:
581 15252 : poGeom = new (std::nothrow) OGRLineString();
582 7626 : break;
583 :
584 27600 : case wkbPolygon:
585 55199 : poGeom = new (std::nothrow) OGRPolygon();
586 27599 : break;
587 :
588 1977 : case wkbGeometryCollection:
589 3954 : poGeom = new (std::nothrow) OGRGeometryCollection();
590 1977 : break;
591 :
592 2818 : case wkbMultiPolygon:
593 5636 : poGeom = new (std::nothrow) OGRMultiPolygon();
594 2818 : 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 115 : case wkbLinearRing:
605 230 : poGeom = new (std::nothrow) OGRLinearRing();
606 115 : 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 258895 : if (poGeom)
648 : {
649 258895 : if (OGR_GT_HasZ(eGeometryType))
650 64736 : poGeom->set3D(true);
651 258894 : if (OGR_GT_HasM(eGeometryType))
652 59821 : poGeom->setMeasured(true);
653 : }
654 258894 : 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 163121 : OGRGeometryH OGR_G_CreateGeometry(OGRwkbGeometryType eGeometryType)
677 :
678 : {
679 163121 : return OGRGeometry::ToHandle(
680 163121 : 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 291511 : void OGR_G_DestroyGeometry(OGRGeometryH hGeom)
720 :
721 : {
722 291511 : delete OGRGeometry::FromHandle(hGeom);
723 291511 : }
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 60178 : sPolyExtended() = default;
1498 111890 : 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 48212 : OGRGeometry *OGRGeometryFactory::organizePolygons(OGRGeometry **papoPolygons,
1592 : int nPolygonCount,
1593 : int *pbIsValidGeometry,
1594 : const char **papszOptions)
1595 : {
1596 48212 : if (nPolygonCount == 0)
1597 : {
1598 4 : if (pbIsValidGeometry)
1599 0 : *pbIsValidGeometry = TRUE;
1600 :
1601 4 : return new OGRPolygon();
1602 : }
1603 :
1604 48208 : OGRGeometry *geom = nullptr;
1605 48208 : OrganizePolygonMethod method = METHOD_NORMAL;
1606 48208 : bool bHasCurves = false;
1607 :
1608 : /* -------------------------------------------------------------------- */
1609 : /* Trivial case of a single polygon. */
1610 : /* -------------------------------------------------------------------- */
1611 48208 : 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 14770 : bool bUseFastVersion = true;
1640 14770 : 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 29540 : std::vector<sPolyExtended> asPolyEx;
1665 14770 : asPolyEx.reserve(nPolygonCount);
1666 :
1667 14770 : bool bValidTopology = true;
1668 14770 : bool bMixedUpGeometries = false;
1669 14770 : bool bFoundCCW = false;
1670 :
1671 14770 : const char *pszMethodValue = CSLFetchNameValue(papszOptions, "METHOD");
1672 : const char *pszMethodValueOption =
1673 14770 : CPLGetConfigOption("OGR_ORGANIZE_POLYGONS", nullptr);
1674 14770 : if (pszMethodValueOption != nullptr && pszMethodValueOption[0] != '\0')
1675 13944 : pszMethodValue = pszMethodValueOption;
1676 :
1677 14770 : if (pszMethodValue != nullptr)
1678 : {
1679 14286 : if (EQUAL(pszMethodValue, "SKIP"))
1680 : {
1681 13948 : method = METHOD_SKIP;
1682 13948 : bMixedUpGeometries = true;
1683 : }
1684 338 : else if (EQUAL(pszMethodValue, "ONLY_CCW"))
1685 : {
1686 282 : 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 14770 : int nCountCWPolygon = 0;
1701 14770 : int indexOfCWPolygon = -1;
1702 :
1703 74951 : for (int i = 0; i < nPolygonCount; i++)
1704 : {
1705 : OGRwkbGeometryType eType =
1706 60181 : wkbFlatten(papoPolygons[i]->getGeometryType());
1707 :
1708 60181 : 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 120356 : sPolyExtended sPolyEx;
1718 :
1719 60178 : sPolyEx.nInitialIndex = i;
1720 60178 : sPolyEx.poGeometry = papoPolygons[i];
1721 60178 : sPolyEx.poPolygon = papoPolygons[i]->toCurvePolygon();
1722 :
1723 60178 : papoPolygons[i]->getEnvelope(&sPolyEx.sEnvelope);
1724 :
1725 60178 : if (eType == wkbCurvePolygon)
1726 33 : bHasCurves = true;
1727 60178 : if (!sPolyEx.poPolygon->IsEmpty() &&
1728 120356 : sPolyEx.poPolygon->getNumInteriorRings() == 0 &&
1729 60178 : sPolyEx.poPolygon->getExteriorRingCurve()->getNumPoints() >= 4)
1730 : {
1731 60176 : if (method != METHOD_CCW_INNER_JUST_AFTER_CW_OUTER)
1732 60176 : sPolyEx.dfArea = sPolyEx.poPolygon->get_Area();
1733 60176 : sPolyEx.poExteriorRing = sPolyEx.poPolygon->getExteriorRingCurve();
1734 60176 : sPolyEx.poExteriorRing->StartPoint(&sPolyEx.poAPoint);
1735 60176 : if (eType == wkbPolygon)
1736 : {
1737 60143 : sPolyEx.bIsCW = CPL_TO_BOOL(
1738 60143 : sPolyEx.poExteriorRing->toLinearRing()->isClockwise());
1739 60143 : 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 60176 : if (sPolyEx.bIsCW)
1751 : {
1752 17133 : indexOfCWPolygon = i;
1753 17133 : nCountCWPolygon++;
1754 : }
1755 60176 : if (!bFoundCCW)
1756 29580 : 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 60178 : 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 14770 : if ((method == METHOD_ONLY_CCW ||
1778 282 : method == METHOD_CCW_INNER_JUST_AFTER_CW_OUTER) &&
1779 101 : nCountCWPolygon == 1 && bUseFastVersion)
1780 : {
1781 101 : OGRCurvePolygon *poCP = asPolyEx[indexOfCWPolygon].poPolygon;
1782 353 : for (int i = 0; i < static_cast<int>(asPolyEx.size()); i++)
1783 : {
1784 252 : if (i != indexOfCWPolygon)
1785 : {
1786 151 : poCP->addRingDirectly(
1787 151 : asPolyEx[i].poPolygon->stealExteriorRingCurve());
1788 151 : delete asPolyEx[i].poPolygon;
1789 : }
1790 : }
1791 :
1792 101 : if (pbIsValidGeometry)
1793 99 : *pbIsValidGeometry = TRUE;
1794 101 : 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 1935 : 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 1935 : size_t nSize = 0;
2260 1935 : unsigned char *pabyBuf = nullptr;
2261 1935 : OGRGeometry *poGeometry = nullptr;
2262 :
2263 : // Special case as POINT EMPTY cannot be translated to WKB.
2264 1992 : if (GEOSGeomTypeId_r(hGEOSCtxt, geosGeom) == GEOS_POINT &&
2265 57 : GEOSisEmpty_r(hGEOSCtxt, geosGeom))
2266 4 : return new OGRPoint();
2267 :
2268 : const int nCoordDim =
2269 1929 : GEOSGeom_getCoordinateDimension_r(hGEOSCtxt, geosGeom);
2270 1930 : GEOSWKBWriter *wkbwriter = GEOSWKBWriter_create_r(hGEOSCtxt);
2271 1931 : GEOSWKBWriter_setOutputDimension_r(hGEOSCtxt, wkbwriter, nCoordDim);
2272 1931 : pabyBuf = GEOSWKBWriter_write_r(hGEOSCtxt, wkbwriter, geosGeom, &nSize);
2273 1929 : GEOSWKBWriter_destroy_r(hGEOSCtxt, wkbwriter);
2274 :
2275 1930 : if (pabyBuf == nullptr || nSize == 0)
2276 : {
2277 1 : return nullptr;
2278 : }
2279 :
2280 1929 : if (OGRGeometryFactory::createFromWkb(pabyBuf, nullptr, &poGeometry,
2281 1931 : static_cast<int>(nSize)) !=
2282 : OGRERR_NONE)
2283 : {
2284 0 : poGeometry = nullptr;
2285 : }
2286 :
2287 1931 : GEOSFree_r(hGEOSCtxt, pabyBuf);
2288 :
2289 1928 : 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 31105 : bool OGRGeometryFactory::haveGEOS()
2308 :
2309 : {
2310 : #ifndef HAVE_GEOS
2311 : return false;
2312 : #else
2313 31105 : 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 352 : OGRErr OGRGeometryFactory::createFromFgf(const void *pabyData,
2347 : OGRSpatialReference *poSR,
2348 : OGRGeometry **ppoReturn, int nBytes,
2349 : int *pnBytesConsumed)
2350 :
2351 : {
2352 352 : return createFromFgfInternal(static_cast<const GByte *>(pabyData), poSR,
2353 352 : ppoReturn, nBytes, pnBytesConsumed, 0);
2354 : }
2355 :
2356 : /************************************************************************/
2357 : /* createFromFgfInternal() */
2358 : /************************************************************************/
2359 :
2360 355 : 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 355 : 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 355 : *ppoReturn = nullptr;
2374 :
2375 355 : if (nBytes < 4)
2376 109 : return OGRERR_NOT_ENOUGH_DATA;
2377 :
2378 : /* -------------------------------------------------------------------- */
2379 : /* Decode the geometry type. */
2380 : /* -------------------------------------------------------------------- */
2381 246 : GInt32 nGType = 0;
2382 246 : memcpy(&nGType, pabyData + 0, 4);
2383 246 : CPL_LSBPTR32(&nGType);
2384 :
2385 246 : if (nGType < 0 || nGType > 13)
2386 232 : 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 964 : OGRGeometryFactory::TransformWithOptionsCache::TransformWithOptionsCache()
3862 964 : : d(new Private())
3863 : {
3864 964 : }
3865 :
3866 : /************************************************************************/
3867 : /* ~TransformWithOptionsCache() */
3868 : /************************************************************************/
3869 :
3870 964 : OGRGeometryFactory::TransformWithOptionsCache::~TransformWithOptionsCache()
3871 : {
3872 964 : }
3873 :
3874 : /************************************************************************/
3875 : /* isTransformWithOptionsRegularTransform() */
3876 : /************************************************************************/
3877 :
3878 : //! @cond Doxygen_Suppress
3879 : /*static */
3880 58 : bool OGRGeometryFactory::isTransformWithOptionsRegularTransform(
3881 : [[maybe_unused]] const OGRSpatialReference *poSourceCRS,
3882 : [[maybe_unused]] const OGRSpatialReference *poTargetCRS,
3883 : CSLConstList papszOptions)
3884 : {
3885 58 : if (papszOptions)
3886 11 : return false;
3887 :
3888 : #ifdef HAVE_GEOS
3889 84 : if (poSourceCRS->IsProjected() && poTargetCRS->IsGeographic() &&
3890 115 : poTargetCRS->GetAxisMappingStrategy() == OAMS_TRADITIONAL_GIS_ORDER &&
3891 : // check that angular units is degree
3892 31 : std::fabs(poTargetCRS->GetAngularUnits(nullptr) -
3893 31 : CPLAtof(SRS_UA_DEGREE_CONV)) <=
3894 31 : 1e-8 * CPLAtof(SRS_UA_DEGREE_CONV))
3895 : {
3896 31 : double dfWestLong = 0.0;
3897 31 : double dfSouthLat = 0.0;
3898 31 : double dfEastLong = 0.0;
3899 31 : double dfNorthLat = 0.0;
3900 31 : if (poSourceCRS->GetAreaOfUse(&dfWestLong, &dfSouthLat, &dfEastLong,
3901 59 : &dfNorthLat, nullptr) &&
3902 28 : !(dfSouthLat == -90.0 || dfNorthLat == 90.0 ||
3903 28 : dfWestLong == -180.0 || dfEastLong == 180.0 ||
3904 21 : dfWestLong > dfEastLong))
3905 : {
3906 : // Not a global geographic CRS
3907 21 : return true;
3908 : }
3909 10 : return false;
3910 : }
3911 : #endif
3912 :
3913 16 : return true;
3914 : }
3915 :
3916 : //! @endcond
3917 :
3918 : /************************************************************************/
3919 : /* transformWithOptions() */
3920 : /************************************************************************/
3921 :
3922 : /** Transform a geometry.
3923 : * @param poSrcGeom source geometry
3924 : * @param poCT coordinate transformation object, or NULL.
3925 : * @param papszOptions options. Including WRAPDATELINE=YES and DATELINEOFFSET=.
3926 : * @param cache Cache. May increase performance if persisted between invocations
3927 : * @return (new) transformed geometry.
3928 : */
3929 181 : OGRGeometry *OGRGeometryFactory::transformWithOptions(
3930 : const OGRGeometry *poSrcGeom, OGRCoordinateTransformation *poCT,
3931 : char **papszOptions, CPL_UNUSED const TransformWithOptionsCache &cache)
3932 : {
3933 362 : auto poDstGeom = std::unique_ptr<OGRGeometry>(poSrcGeom->clone());
3934 181 : if (poCT)
3935 : {
3936 : #ifdef HAVE_GEOS
3937 146 : bool bNeedPostCorrection = false;
3938 146 : const auto poSourceCRS = poCT->GetSourceCS();
3939 146 : const auto poTargetCRS = poCT->GetTargetCS();
3940 146 : const auto eSrcGeomType = wkbFlatten(poSrcGeom->getGeometryType());
3941 : // Check if we are transforming from projected coordinates to
3942 : // geographic coordinates, with a chance that there might be polar or
3943 : // anti-meridian discontinuities. If so, create the inverse transform.
3944 272 : if (eSrcGeomType != wkbPoint && eSrcGeomType != wkbMultiPoint &&
3945 126 : (poSourceCRS != cache.d->poSourceCRS ||
3946 80 : poTargetCRS != cache.d->poTargetCRS || poCT != cache.d->poCT))
3947 : {
3948 46 : cache.d->clear();
3949 46 : cache.d->poSourceCRS = poSourceCRS;
3950 46 : cache.d->poTargetCRS = poTargetCRS;
3951 46 : cache.d->poCT = poCT;
3952 92 : if (poSourceCRS && poTargetCRS &&
3953 46 : !isTransformWithOptionsRegularTransform(
3954 : poSourceCRS, poTargetCRS, papszOptions))
3955 : {
3956 20 : cache.d->poRevCT.reset(OGRCreateCoordinateTransformation(
3957 : poTargetCRS, poSourceCRS));
3958 20 : cache.d->bIsNorthPolar = false;
3959 20 : cache.d->bIsPolar = false;
3960 20 : cache.d->poRevCT.reset(poCT->GetInverse());
3961 60 : if (cache.d->poRevCT &&
3962 20 : IsPolarToGeographic(poCT, cache.d->poRevCT.get(),
3963 40 : cache.d->bIsNorthPolar))
3964 : {
3965 5 : cache.d->bIsPolar = true;
3966 : }
3967 : }
3968 : }
3969 :
3970 146 : if (auto poRevCT = cache.d->poRevCT.get())
3971 : {
3972 23 : if (cache.d->bIsPolar)
3973 : {
3974 12 : poDstGeom = TransformBeforePolarToGeographic(
3975 12 : poRevCT, cache.d->bIsNorthPolar, std::move(poDstGeom),
3976 6 : bNeedPostCorrection);
3977 : }
3978 17 : else if (IsAntimeridianProjToGeographic(poCT, poRevCT,
3979 : poDstGeom.get()))
3980 : {
3981 30 : poDstGeom = TransformBeforeAntimeridianToGeographic(
3982 30 : poCT, poRevCT, std::move(poDstGeom), bNeedPostCorrection);
3983 : }
3984 : }
3985 : #endif
3986 146 : OGRErr eErr = poDstGeom->transform(poCT);
3987 146 : if (eErr != OGRERR_NONE)
3988 : {
3989 0 : return nullptr;
3990 : }
3991 : #ifdef HAVE_GEOS
3992 146 : if (bNeedPostCorrection)
3993 : {
3994 16 : SnapCoordsCloseToLatLongBounds(poDstGeom.get());
3995 : }
3996 : #endif
3997 : }
3998 :
3999 181 : if (CPLTestBool(CSLFetchNameValueDef(papszOptions, "WRAPDATELINE", "NO")))
4000 : {
4001 101 : if (poDstGeom->getSpatialReference() &&
4002 49 : !poDstGeom->getSpatialReference()->IsGeographic())
4003 : {
4004 : static bool bHasWarned = false;
4005 0 : if (!bHasWarned)
4006 : {
4007 0 : CPLError(
4008 : CE_Warning, CPLE_AppDefined,
4009 : "WRAPDATELINE is without effect when reprojecting to a "
4010 : "non-geographic CRS");
4011 0 : bHasWarned = true;
4012 : }
4013 0 : return poDstGeom.release();
4014 : }
4015 : // TODO and we should probably also test that the axis order + data axis
4016 : // mapping is long-lat...
4017 : const OGRwkbGeometryType eType =
4018 52 : wkbFlatten(poDstGeom->getGeometryType());
4019 52 : if (eType == wkbPoint)
4020 : {
4021 9 : OGRPoint *poDstPoint = poDstGeom->toPoint();
4022 9 : WrapPointDateLine(poDstPoint);
4023 : }
4024 43 : else if (eType == wkbMultiPoint)
4025 : {
4026 5 : for (auto *poDstPoint : *(poDstGeom->toMultiPoint()))
4027 : {
4028 4 : WrapPointDateLine(poDstPoint);
4029 : }
4030 : }
4031 : else
4032 : {
4033 42 : OGREnvelope sEnvelope;
4034 42 : poDstGeom->getEnvelope(&sEnvelope);
4035 42 : if (sEnvelope.MinX >= -360.0 && sEnvelope.MaxX <= -180.0)
4036 2 : AddOffsetToLon(poDstGeom.get(), 360.0);
4037 40 : else if (sEnvelope.MinX >= 180.0 && sEnvelope.MaxX <= 360.0)
4038 2 : AddOffsetToLon(poDstGeom.get(), -360.0);
4039 : else
4040 : {
4041 : OGRwkbGeometryType eNewType;
4042 38 : if (eType == wkbPolygon || eType == wkbMultiPolygon)
4043 27 : eNewType = wkbMultiPolygon;
4044 11 : else if (eType == wkbLineString || eType == wkbMultiLineString)
4045 10 : eNewType = wkbMultiLineString;
4046 : else
4047 1 : eNewType = wkbGeometryCollection;
4048 :
4049 : auto poMulti = std::unique_ptr<OGRGeometryCollection>(
4050 76 : createGeometry(eNewType)->toGeometryCollection());
4051 :
4052 38 : double dfDateLineOffset = CPLAtofM(
4053 : CSLFetchNameValueDef(papszOptions, "DATELINEOFFSET", "10"));
4054 38 : if (dfDateLineOffset <= 0.0 || dfDateLineOffset >= 360.0)
4055 0 : dfDateLineOffset = 10.0;
4056 :
4057 38 : CutGeometryOnDateLineAndAddToMulti(
4058 38 : poMulti.get(), poDstGeom.get(), dfDateLineOffset);
4059 :
4060 38 : if (poMulti->getNumGeometries() == 0)
4061 : {
4062 : // do nothing
4063 : }
4064 39 : else if (poMulti->getNumGeometries() == 1 &&
4065 1 : (eType == wkbPolygon || eType == wkbLineString))
4066 : {
4067 12 : poDstGeom = poMulti->stealGeometry(0);
4068 : }
4069 : else
4070 : {
4071 26 : poDstGeom = std::move(poMulti);
4072 : }
4073 : }
4074 : }
4075 : }
4076 :
4077 181 : return poDstGeom.release();
4078 : }
4079 :
4080 : /************************************************************************/
4081 : /* OGRGeomTransformer() */
4082 : /************************************************************************/
4083 :
4084 : struct OGRGeomTransformer
4085 : {
4086 : std::unique_ptr<OGRCoordinateTransformation> poCT{};
4087 : OGRGeometryFactory::TransformWithOptionsCache cache{};
4088 : CPLStringList aosOptions{};
4089 :
4090 6 : OGRGeomTransformer() = default;
4091 : OGRGeomTransformer(const OGRGeomTransformer &) = delete;
4092 : OGRGeomTransformer &operator=(const OGRGeomTransformer &) = delete;
4093 : };
4094 :
4095 : /************************************************************************/
4096 : /* OGR_GeomTransformer_Create() */
4097 : /************************************************************************/
4098 :
4099 : /** Create a geometry transformer.
4100 : *
4101 : * This is a enhanced version of OGR_G_Transform().
4102 : *
4103 : * When reprojecting geometries from a Polar Stereographic projection or a
4104 : * projection naturally crossing the antimeridian (like UTM Zone 60) to a
4105 : * geographic CRS, it will cut geometries along the antimeridian. So a
4106 : * LineString might be returned as a MultiLineString.
4107 : *
4108 : * The WRAPDATELINE=YES option might be specified for circumstances to correct
4109 : * geometries that incorrectly go from a longitude on a side of the antimeridian
4110 : * to the other side, like a LINESTRING(-179 0,179 0) will be transformed to
4111 : * a MULTILINESTRING ((-179 0,-180 0),(180 0,179 0)). For that use case, hCT
4112 : * might be NULL.
4113 : *
4114 : * @param hCT Coordinate transformation object (will be cloned) or NULL.
4115 : * @param papszOptions NULL terminated list of options, or NULL.
4116 : * Supported options are:
4117 : * <ul>
4118 : * <li>WRAPDATELINE=YES</li>
4119 : * <li>DATELINEOFFSET=longitude_gap_in_degree. Defaults
4120 : * to 10.</li>
4121 : * </ul>
4122 : * @return transformer object to free with OGR_GeomTransformer_Destroy()
4123 : * @since GDAL 3.1
4124 : */
4125 6 : OGRGeomTransformerH OGR_GeomTransformer_Create(OGRCoordinateTransformationH hCT,
4126 : CSLConstList papszOptions)
4127 : {
4128 6 : OGRGeomTransformer *transformer = new OGRGeomTransformer;
4129 6 : if (hCT)
4130 : {
4131 3 : transformer->poCT.reset(
4132 3 : OGRCoordinateTransformation::FromHandle(hCT)->Clone());
4133 : }
4134 6 : transformer->aosOptions.Assign(CSLDuplicate(papszOptions));
4135 6 : return transformer;
4136 : }
4137 :
4138 : /************************************************************************/
4139 : /* OGR_GeomTransformer_Transform() */
4140 : /************************************************************************/
4141 :
4142 : /** Transforms a geometry.
4143 : *
4144 : * @param hTransformer transformer object.
4145 : * @param hGeom Source geometry.
4146 : * @return a new geometry (or NULL) to destroy with OGR_G_DestroyGeometry()
4147 : * @since GDAL 3.1
4148 : */
4149 6 : OGRGeometryH OGR_GeomTransformer_Transform(OGRGeomTransformerH hTransformer,
4150 : OGRGeometryH hGeom)
4151 : {
4152 6 : VALIDATE_POINTER1(hTransformer, "OGR_GeomTransformer_Transform", nullptr);
4153 6 : VALIDATE_POINTER1(hGeom, "OGR_GeomTransformer_Transform", nullptr);
4154 :
4155 12 : return OGRGeometry::ToHandle(OGRGeometryFactory::transformWithOptions(
4156 6 : OGRGeometry::FromHandle(hGeom), hTransformer->poCT.get(),
4157 12 : hTransformer->aosOptions.List(), hTransformer->cache));
4158 : }
4159 :
4160 : /************************************************************************/
4161 : /* OGR_GeomTransformer_Destroy() */
4162 : /************************************************************************/
4163 :
4164 : /** Destroy a geometry transformer allocated with OGR_GeomTransformer_Create()
4165 : *
4166 : * @param hTransformer transformer object.
4167 : * @since GDAL 3.1
4168 : */
4169 6 : void OGR_GeomTransformer_Destroy(OGRGeomTransformerH hTransformer)
4170 : {
4171 6 : delete hTransformer;
4172 6 : }
4173 :
4174 : /************************************************************************/
4175 : /* OGRGF_GetDefaultStepSize() */
4176 : /************************************************************************/
4177 :
4178 4314 : static double OGRGF_GetDefaultStepSize()
4179 : {
4180 : // coverity[tainted_data]
4181 4314 : return CPLAtofM(CPLGetConfigOption("OGR_ARC_STEPSIZE", "4"));
4182 : }
4183 :
4184 : /************************************************************************/
4185 : /* DISTANCE() */
4186 : /************************************************************************/
4187 :
4188 309996 : static inline double DISTANCE(double x1, double y1, double x2, double y2)
4189 : {
4190 309996 : return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
4191 : }
4192 :
4193 : /************************************************************************/
4194 : /* approximateArcAngles() */
4195 : /************************************************************************/
4196 :
4197 : /**
4198 : * Stroke arc to linestring.
4199 : *
4200 : * Stroke an arc of a circle to a linestring based on a center
4201 : * point, radius, start angle and end angle, all angles in degrees.
4202 : *
4203 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
4204 : * used. This is currently 4 degrees unless the user has overridden the
4205 : * value with the OGR_ARC_STEPSIZE configuration variable.
4206 : *
4207 : * If the OGR_ARC_MAX_GAP configuration variable is set, the straight-line
4208 : * distance between adjacent pairs of interpolated points will be limited to
4209 : * the specified distance. If the distance between a pair of points exceeds
4210 : * this maximum, additional points are interpolated between the two points.
4211 : *
4212 : * @see CPLSetConfigOption()
4213 : *
4214 : * @param dfCenterX center X
4215 : * @param dfCenterY center Y
4216 : * @param dfZ center Z
4217 : * @param dfPrimaryRadius X radius of ellipse.
4218 : * @param dfSecondaryRadius Y radius of ellipse.
4219 : * @param dfRotation rotation of the ellipse clockwise.
4220 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
4221 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
4222 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
4223 : * arc, zero to use the default setting.
4224 : * @param bUseMaxGap Optional: whether to honor OGR_ARC_MAX_GAP.
4225 : *
4226 : * @return OGRLineString geometry representing an approximation of the arc.
4227 : *
4228 : * @since OGR 1.8.0
4229 : */
4230 :
4231 117 : OGRGeometry *OGRGeometryFactory::approximateArcAngles(
4232 : double dfCenterX, double dfCenterY, double dfZ, double dfPrimaryRadius,
4233 : double dfSecondaryRadius, double dfRotation, double dfStartAngle,
4234 : double dfEndAngle, double dfMaxAngleStepSizeDegrees,
4235 : const bool bUseMaxGap /* = false */)
4236 :
4237 : {
4238 117 : OGRLineString *poLine = new OGRLineString();
4239 117 : const double dfRotationRadians = dfRotation * M_PI / 180.0;
4240 :
4241 : // Support default arc step setting.
4242 117 : if (dfMaxAngleStepSizeDegrees < 1e-6)
4243 : {
4244 116 : dfMaxAngleStepSizeDegrees = OGRGF_GetDefaultStepSize();
4245 : }
4246 :
4247 : // Determine maximum interpolation gap. This is the largest straight-line
4248 : // distance allowed between pairs of interpolated points. Default zero,
4249 : // meaning no gap.
4250 : // coverity[tainted_data]
4251 : const double dfMaxInterpolationGap =
4252 117 : bUseMaxGap ? CPLAtofM(CPLGetConfigOption("OGR_ARC_MAX_GAP", "0")) : 0.0;
4253 :
4254 : // Is this a full circle?
4255 117 : const bool bIsFullCircle = fabs(dfEndAngle - dfStartAngle) == 360.0;
4256 :
4257 : // Switch direction.
4258 117 : dfStartAngle *= -1;
4259 117 : dfEndAngle *= -1;
4260 :
4261 : // Figure out the number of slices to make this into.
4262 : int nVertexCount =
4263 234 : std::max(2, static_cast<int>(ceil(fabs(dfEndAngle - dfStartAngle) /
4264 117 : dfMaxAngleStepSizeDegrees) +
4265 117 : 1));
4266 117 : const double dfSlice = (dfEndAngle - dfStartAngle) / (nVertexCount - 1);
4267 :
4268 : // If it is a full circle we will work out the last point separately.
4269 117 : if (bIsFullCircle)
4270 : {
4271 51 : nVertexCount--;
4272 : }
4273 :
4274 : /* -------------------------------------------------------------------- */
4275 : /* Compute the interpolated points. */
4276 : /* -------------------------------------------------------------------- */
4277 117 : double dfLastX = 0.0;
4278 117 : double dfLastY = 0.0;
4279 117 : int nTotalAddPoints = 0;
4280 6980 : for (int iPoint = 0; iPoint < nVertexCount; iPoint++)
4281 : {
4282 6863 : const double dfAngleOnEllipse =
4283 6863 : (dfStartAngle + iPoint * dfSlice) * M_PI / 180.0;
4284 :
4285 : // Compute position on the unrotated ellipse.
4286 6863 : const double dfEllipseX = cos(dfAngleOnEllipse) * dfPrimaryRadius;
4287 6863 : const double dfEllipseY = sin(dfAngleOnEllipse) * dfSecondaryRadius;
4288 :
4289 : // Is this point too far from the previous point?
4290 6863 : if (iPoint && dfMaxInterpolationGap != 0.0)
4291 : {
4292 : const double dfDistFromLast =
4293 1 : DISTANCE(dfLastX, dfLastY, dfEllipseX, dfEllipseY);
4294 :
4295 1 : if (dfDistFromLast > dfMaxInterpolationGap)
4296 : {
4297 1 : const int nAddPoints =
4298 1 : static_cast<int>(dfDistFromLast / dfMaxInterpolationGap);
4299 1 : const double dfAddSlice = dfSlice / (nAddPoints + 1);
4300 :
4301 : // Interpolate additional points
4302 3 : for (int iAddPoint = 0; iAddPoint < nAddPoints; iAddPoint++)
4303 : {
4304 2 : const double dfAddAngleOnEllipse =
4305 2 : (dfStartAngle + (iPoint - 1) * dfSlice +
4306 2 : (iAddPoint + 1) * dfAddSlice) *
4307 : (M_PI / 180.0);
4308 :
4309 2 : poLine->setPoint(
4310 2 : iPoint + nTotalAddPoints + iAddPoint,
4311 2 : cos(dfAddAngleOnEllipse) * dfPrimaryRadius,
4312 2 : sin(dfAddAngleOnEllipse) * dfSecondaryRadius, dfZ);
4313 : }
4314 :
4315 1 : nTotalAddPoints += nAddPoints;
4316 : }
4317 : }
4318 :
4319 6863 : poLine->setPoint(iPoint + nTotalAddPoints, dfEllipseX, dfEllipseY, dfZ);
4320 6863 : dfLastX = dfEllipseX;
4321 6863 : dfLastY = dfEllipseY;
4322 : }
4323 :
4324 : /* -------------------------------------------------------------------- */
4325 : /* Rotate and translate the ellipse. */
4326 : /* -------------------------------------------------------------------- */
4327 117 : nVertexCount = poLine->getNumPoints();
4328 6982 : for (int iPoint = 0; iPoint < nVertexCount; iPoint++)
4329 : {
4330 6865 : const double dfEllipseX = poLine->getX(iPoint);
4331 6865 : const double dfEllipseY = poLine->getY(iPoint);
4332 :
4333 : // Rotate this position around the center of the ellipse.
4334 6865 : const double dfArcX = dfCenterX + dfEllipseX * cos(dfRotationRadians) +
4335 6865 : dfEllipseY * sin(dfRotationRadians);
4336 6865 : const double dfArcY = dfCenterY - dfEllipseX * sin(dfRotationRadians) +
4337 6865 : dfEllipseY * cos(dfRotationRadians);
4338 :
4339 6865 : poLine->setPoint(iPoint, dfArcX, dfArcY, dfZ);
4340 : }
4341 :
4342 : /* -------------------------------------------------------------------- */
4343 : /* If we're asked to make a full circle, ensure the start and */
4344 : /* end points coincide exactly, in spite of any rounding error. */
4345 : /* -------------------------------------------------------------------- */
4346 117 : if (bIsFullCircle)
4347 : {
4348 102 : OGRPoint oPoint;
4349 51 : poLine->getPoint(0, &oPoint);
4350 51 : poLine->setPoint(nVertexCount, &oPoint);
4351 : }
4352 :
4353 117 : return poLine;
4354 : }
4355 :
4356 : /************************************************************************/
4357 : /* OGR_G_ApproximateArcAngles() */
4358 : /************************************************************************/
4359 :
4360 : /**
4361 : * Stroke arc to linestring.
4362 : *
4363 : * Stroke an arc of a circle to a linestring based on a center
4364 : * point, radius, start angle and end angle, all angles in degrees.
4365 : *
4366 : * If the dfMaxAngleStepSizeDegrees is zero, then a default value will be
4367 : * used. This is currently 4 degrees unless the user has overridden the
4368 : * value with the OGR_ARC_STEPSIZE configuration variable.
4369 : *
4370 : * @see CPLSetConfigOption()
4371 : *
4372 : * @param dfCenterX center X
4373 : * @param dfCenterY center Y
4374 : * @param dfZ center Z
4375 : * @param dfPrimaryRadius X radius of ellipse.
4376 : * @param dfSecondaryRadius Y radius of ellipse.
4377 : * @param dfRotation rotation of the ellipse clockwise.
4378 : * @param dfStartAngle angle to first point on arc (clockwise of X-positive)
4379 : * @param dfEndAngle angle to last point on arc (clockwise of X-positive)
4380 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
4381 : * arc, zero to use the default setting.
4382 : *
4383 : * @return OGRLineString geometry representing an approximation of the arc.
4384 : *
4385 : * @since OGR 1.8.0
4386 : */
4387 :
4388 1 : OGRGeometryH CPL_DLL OGR_G_ApproximateArcAngles(
4389 : double dfCenterX, double dfCenterY, double dfZ, double dfPrimaryRadius,
4390 : double dfSecondaryRadius, double dfRotation, double dfStartAngle,
4391 : double dfEndAngle, double dfMaxAngleStepSizeDegrees)
4392 :
4393 : {
4394 1 : return OGRGeometry::ToHandle(OGRGeometryFactory::approximateArcAngles(
4395 : dfCenterX, dfCenterY, dfZ, dfPrimaryRadius, dfSecondaryRadius,
4396 1 : dfRotation, dfStartAngle, dfEndAngle, dfMaxAngleStepSizeDegrees));
4397 : }
4398 :
4399 : /************************************************************************/
4400 : /* forceToLineString() */
4401 : /************************************************************************/
4402 :
4403 : /**
4404 : * \brief Convert to line string.
4405 : *
4406 : * Tries to force the provided geometry to be a line string. This nominally
4407 : * effects a change on multilinestrings.
4408 : * In GDAL 2.0, for polygons or curvepolygons that have a single exterior ring,
4409 : * it will return the ring. For circular strings or compound curves, it will
4410 : * return an approximated line string.
4411 : *
4412 : * The passed in geometry is
4413 : * consumed and a new one returned (or potentially the same one).
4414 : *
4415 : * @param poGeom the input geometry - ownership is passed to the method.
4416 : * @param bOnlyInOrder flag that, if set to FALSE, indicate that the order of
4417 : * points in a linestring might be reversed if it enables
4418 : * to match the extremity of another linestring. If set
4419 : * to TRUE, the start of a linestring must match the end
4420 : * of another linestring.
4421 : * @return new geometry.
4422 : */
4423 :
4424 173 : OGRGeometry *OGRGeometryFactory::forceToLineString(OGRGeometry *poGeom,
4425 : bool bOnlyInOrder)
4426 :
4427 : {
4428 173 : if (poGeom == nullptr)
4429 2 : return nullptr;
4430 :
4431 171 : const OGRwkbGeometryType eGeomType = wkbFlatten(poGeom->getGeometryType());
4432 :
4433 : /* -------------------------------------------------------------------- */
4434 : /* If this is already a LineString, nothing to do */
4435 : /* -------------------------------------------------------------------- */
4436 171 : if (eGeomType == wkbLineString)
4437 : {
4438 : // Except if it is a linearring.
4439 25 : poGeom = OGRCurve::CastToLineString(poGeom->toCurve());
4440 :
4441 25 : return poGeom;
4442 : }
4443 :
4444 : /* -------------------------------------------------------------------- */
4445 : /* If it is a polygon with a single ring, return it */
4446 : /* -------------------------------------------------------------------- */
4447 146 : if (eGeomType == wkbPolygon || eGeomType == wkbCurvePolygon)
4448 : {
4449 30 : OGRCurvePolygon *poCP = poGeom->toCurvePolygon();
4450 30 : if (poCP->getNumInteriorRings() == 0)
4451 : {
4452 28 : OGRCurve *poRing = poCP->stealExteriorRingCurve();
4453 28 : delete poCP;
4454 28 : return forceToLineString(poRing);
4455 : }
4456 2 : return poGeom;
4457 : }
4458 :
4459 : /* -------------------------------------------------------------------- */
4460 : /* If it is a curve line, call CurveToLine() */
4461 : /* -------------------------------------------------------------------- */
4462 116 : if (eGeomType == wkbCircularString || eGeomType == wkbCompoundCurve)
4463 : {
4464 68 : OGRGeometry *poNewGeom = poGeom->toCurve()->CurveToLine();
4465 68 : delete poGeom;
4466 68 : return poNewGeom;
4467 : }
4468 :
4469 48 : if (eGeomType != wkbGeometryCollection && eGeomType != wkbMultiLineString &&
4470 : eGeomType != wkbMultiCurve)
4471 17 : return poGeom;
4472 :
4473 : // Build an aggregated linestring from all the linestrings in the container.
4474 31 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4475 31 : if (poGeom->hasCurveGeometry())
4476 : {
4477 : OGRGeometryCollection *poNewGC =
4478 7 : poGC->getLinearGeometry()->toGeometryCollection();
4479 7 : delete poGC;
4480 7 : poGC = poNewGC;
4481 : }
4482 :
4483 31 : if (poGC->getNumGeometries() == 0)
4484 : {
4485 3 : poGeom = new OGRLineString();
4486 3 : poGeom->assignSpatialReference(poGC->getSpatialReference());
4487 3 : delete poGC;
4488 3 : return poGeom;
4489 : }
4490 :
4491 28 : int iGeom0 = 0;
4492 69 : while (iGeom0 < poGC->getNumGeometries())
4493 : {
4494 41 : if (wkbFlatten(poGC->getGeometryRef(iGeom0)->getGeometryType()) !=
4495 : wkbLineString)
4496 : {
4497 12 : iGeom0++;
4498 26 : continue;
4499 : }
4500 :
4501 : OGRLineString *poLineString0 =
4502 29 : poGC->getGeometryRef(iGeom0)->toLineString();
4503 29 : if (poLineString0->getNumPoints() < 2)
4504 : {
4505 14 : iGeom0++;
4506 14 : continue;
4507 : }
4508 :
4509 30 : OGRPoint pointStart0;
4510 15 : poLineString0->StartPoint(&pointStart0);
4511 30 : OGRPoint pointEnd0;
4512 15 : poLineString0->EndPoint(&pointEnd0);
4513 :
4514 15 : int iGeom1 = iGeom0 + 1; // Used after for.
4515 17 : for (; iGeom1 < poGC->getNumGeometries(); iGeom1++)
4516 : {
4517 6 : if (wkbFlatten(poGC->getGeometryRef(iGeom1)->getGeometryType()) !=
4518 : wkbLineString)
4519 1 : continue;
4520 :
4521 : OGRLineString *poLineString1 =
4522 6 : poGC->getGeometryRef(iGeom1)->toLineString();
4523 6 : if (poLineString1->getNumPoints() < 2)
4524 1 : continue;
4525 :
4526 5 : OGRPoint pointStart1;
4527 5 : poLineString1->StartPoint(&pointStart1);
4528 5 : OGRPoint pointEnd1;
4529 5 : poLineString1->EndPoint(&pointEnd1);
4530 :
4531 5 : if (!bOnlyInOrder && (pointEnd0.Equals(&pointEnd1) ||
4532 0 : pointStart0.Equals(&pointStart1)))
4533 : {
4534 0 : poLineString1->reversePoints();
4535 0 : poLineString1->StartPoint(&pointStart1);
4536 0 : poLineString1->EndPoint(&pointEnd1);
4537 : }
4538 :
4539 5 : if (pointEnd0.Equals(&pointStart1))
4540 : {
4541 4 : poLineString0->addSubLineString(poLineString1, 1);
4542 4 : poGC->removeGeometry(iGeom1);
4543 4 : break;
4544 : }
4545 :
4546 1 : if (pointEnd1.Equals(&pointStart0))
4547 : {
4548 0 : poLineString1->addSubLineString(poLineString0, 1);
4549 0 : poGC->removeGeometry(iGeom0);
4550 0 : break;
4551 : }
4552 : }
4553 :
4554 15 : if (iGeom1 == poGC->getNumGeometries())
4555 : {
4556 14 : iGeom0++;
4557 : }
4558 : }
4559 :
4560 28 : if (poGC->getNumGeometries() == 1)
4561 : {
4562 20 : OGRGeometry *poSingleGeom = poGC->getGeometryRef(0);
4563 20 : poGC->removeGeometry(0, FALSE);
4564 20 : delete poGC;
4565 :
4566 20 : return poSingleGeom;
4567 : }
4568 :
4569 8 : return poGC;
4570 : }
4571 :
4572 : /************************************************************************/
4573 : /* OGR_G_ForceToLineString() */
4574 : /************************************************************************/
4575 :
4576 : /**
4577 : * \brief Convert to line string.
4578 : *
4579 : * This function is the same as the C++ method
4580 : * OGRGeometryFactory::forceToLineString().
4581 : *
4582 : * @param hGeom handle to the geometry to convert (ownership surrendered).
4583 : * @return the converted geometry (ownership to caller).
4584 : *
4585 : * @since GDAL/OGR 1.10.0
4586 : */
4587 :
4588 60 : OGRGeometryH OGR_G_ForceToLineString(OGRGeometryH hGeom)
4589 :
4590 : {
4591 60 : return OGRGeometry::ToHandle(
4592 60 : OGRGeometryFactory::forceToLineString(OGRGeometry::FromHandle(hGeom)));
4593 : }
4594 :
4595 : /************************************************************************/
4596 : /* forceTo() */
4597 : /************************************************************************/
4598 :
4599 : /**
4600 : * \brief Convert to another geometry type
4601 : *
4602 : * Tries to force the provided geometry to the specified geometry type.
4603 : *
4604 : * It can promote 'single' geometry type to their corresponding collection type
4605 : * (see OGR_GT_GetCollection()) or the reverse. non-linear geometry type to
4606 : * their corresponding linear geometry type (see OGR_GT_GetLinear()), by
4607 : * possibly approximating circular arcs they may contain. Regarding conversion
4608 : * from linear geometry types to curve geometry types, only "wrapping" will be
4609 : * done. No attempt to retrieve potential circular arcs by de-approximating
4610 : * stroking will be done. For that, OGRGeometry::getCurveGeometry() can be used.
4611 : *
4612 : * The passed in geometry is consumed and a new one returned (or potentially the
4613 : * same one).
4614 : *
4615 : * Starting with GDAL 3.9, this method honours the dimensionality of eTargetType.
4616 : *
4617 : * @param poGeom the input geometry - ownership is passed to the method.
4618 : * @param eTargetType target output geometry type.
4619 : * @param papszOptions options as a null-terminated list of strings or NULL.
4620 : * @return new geometry, or nullptr in case of error.
4621 : *
4622 : * @since GDAL 2.0
4623 : */
4624 :
4625 4646 : OGRGeometry *OGRGeometryFactory::forceTo(OGRGeometry *poGeom,
4626 : OGRwkbGeometryType eTargetType,
4627 : const char *const *papszOptions)
4628 : {
4629 4646 : if (poGeom == nullptr)
4630 0 : return poGeom;
4631 :
4632 4646 : const OGRwkbGeometryType eTargetTypeFlat = wkbFlatten(eTargetType);
4633 4646 : if (eTargetTypeFlat == wkbUnknown)
4634 252 : return poGeom;
4635 :
4636 4394 : if (poGeom->IsEmpty())
4637 : {
4638 277 : OGRGeometry *poRet = createGeometry(eTargetType);
4639 277 : if (poRet)
4640 : {
4641 277 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4642 277 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4643 277 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4644 : }
4645 277 : delete poGeom;
4646 277 : return poRet;
4647 : }
4648 :
4649 4117 : OGRwkbGeometryType eType = poGeom->getGeometryType();
4650 4117 : OGRwkbGeometryType eTypeFlat = wkbFlatten(eType);
4651 :
4652 4117 : if (eTargetTypeFlat != eTargetType && (eType == eTypeFlat))
4653 : {
4654 57 : auto poGeomNew = forceTo(poGeom, eTargetTypeFlat, papszOptions);
4655 57 : if (poGeomNew)
4656 : {
4657 57 : poGeomNew->set3D(OGR_GT_HasZ(eTargetType));
4658 57 : poGeomNew->setMeasured(OGR_GT_HasM(eTargetType));
4659 : }
4660 57 : return poGeomNew;
4661 : }
4662 :
4663 4060 : if (eTypeFlat == eTargetTypeFlat)
4664 : {
4665 505 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4666 505 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4667 505 : return poGeom;
4668 : }
4669 :
4670 3555 : eType = eTypeFlat;
4671 :
4672 5272 : if (OGR_GT_IsSubClassOf(eType, wkbPolyhedralSurface) &&
4673 1717 : (eTargetTypeFlat == wkbMultiSurface ||
4674 : eTargetTypeFlat == wkbGeometryCollection))
4675 : {
4676 851 : OGRwkbGeometryType eTempGeomType = wkbMultiPolygon;
4677 851 : if (OGR_GT_HasZ(eTargetType))
4678 847 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4679 851 : if (OGR_GT_HasM(eTargetType))
4680 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4681 851 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4682 851 : eTargetType, papszOptions);
4683 : }
4684 :
4685 2704 : if (OGR_GT_IsSubClassOf(eType, wkbGeometryCollection) &&
4686 : eTargetTypeFlat == wkbGeometryCollection)
4687 : {
4688 909 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4689 909 : auto poRet = OGRGeometryCollection::CastToGeometryCollection(poGC);
4690 909 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4691 909 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4692 909 : return poRet;
4693 : }
4694 :
4695 1795 : if (eType == wkbTriangle && eTargetTypeFlat == wkbPolyhedralSurface)
4696 : {
4697 1 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4698 1 : poPS->assignSpatialReference(poGeom->getSpatialReference());
4699 1 : poPS->addGeometryDirectly(OGRTriangle::CastToPolygon(poGeom));
4700 1 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4701 1 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4702 1 : return poPS;
4703 : }
4704 1794 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbPolyhedralSurface)
4705 : {
4706 3 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4707 3 : poPS->assignSpatialReference(poGeom->getSpatialReference());
4708 3 : poPS->addGeometryDirectly(poGeom);
4709 3 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4710 3 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4711 3 : return poPS;
4712 : }
4713 1791 : else if (eType == wkbMultiPolygon &&
4714 : eTargetTypeFlat == wkbPolyhedralSurface)
4715 : {
4716 2 : OGRMultiPolygon *poMP = poGeom->toMultiPolygon();
4717 2 : OGRPolyhedralSurface *poPS = new OGRPolyhedralSurface();
4718 4 : for (int i = 0; i < poMP->getNumGeometries(); ++i)
4719 : {
4720 2 : poPS->addGeometry(poMP->getGeometryRef(i));
4721 : }
4722 2 : delete poGeom;
4723 2 : poPS->set3D(OGR_GT_HasZ(eTargetType));
4724 2 : poPS->setMeasured(OGR_GT_HasM(eTargetType));
4725 2 : return poPS;
4726 : }
4727 1789 : else if (eType == wkbTIN && eTargetTypeFlat == wkbPolyhedralSurface)
4728 : {
4729 1 : poGeom = OGRTriangulatedSurface::CastToPolyhedralSurface(
4730 : poGeom->toTriangulatedSurface());
4731 : }
4732 1788 : else if (eType == wkbCurvePolygon &&
4733 : eTargetTypeFlat == wkbPolyhedralSurface)
4734 : {
4735 1 : OGRwkbGeometryType eTempGeomType = wkbPolygon;
4736 1 : if (OGR_GT_HasZ(eTargetType))
4737 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4738 1 : if (OGR_GT_HasM(eTargetType))
4739 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4740 1 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4741 1 : eTargetType, papszOptions);
4742 : }
4743 1787 : else if (eType == wkbMultiSurface &&
4744 : eTargetTypeFlat == wkbPolyhedralSurface)
4745 : {
4746 1 : OGRwkbGeometryType eTempGeomType = wkbMultiPolygon;
4747 1 : if (OGR_GT_HasZ(eTargetType))
4748 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4749 1 : if (OGR_GT_HasM(eTargetType))
4750 0 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4751 1 : return forceTo(forceTo(poGeom, eTempGeomType, papszOptions),
4752 1 : eTargetType, papszOptions);
4753 : }
4754 :
4755 1786 : else if (eType == wkbTriangle && eTargetTypeFlat == wkbTIN)
4756 : {
4757 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4758 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4759 1 : poTS->addGeometryDirectly(poGeom);
4760 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4761 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4762 1 : return poTS;
4763 : }
4764 1785 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbTIN)
4765 : {
4766 4 : OGRPolygon *poPoly = poGeom->toPolygon();
4767 4 : OGRLinearRing *poLR = poPoly->getExteriorRing();
4768 7 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4769 3 : poPoly->getNumInteriorRings() == 0))
4770 : {
4771 1 : return poGeom;
4772 : }
4773 3 : OGRErr eErr = OGRERR_NONE;
4774 3 : OGRTriangle *poTriangle = new OGRTriangle(*poPoly, eErr);
4775 3 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4776 3 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4777 3 : poTS->addGeometryDirectly(poTriangle);
4778 3 : delete poGeom;
4779 3 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4780 3 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4781 3 : return poTS;
4782 : }
4783 1781 : else if (eType == wkbMultiPolygon && eTargetTypeFlat == wkbTIN)
4784 : {
4785 1 : OGRMultiPolygon *poMP = poGeom->toMultiPolygon();
4786 2 : for (const auto poPoly : *poMP)
4787 : {
4788 1 : const OGRLinearRing *poLR = poPoly->getExteriorRing();
4789 2 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4790 1 : poPoly->getNumInteriorRings() == 0))
4791 : {
4792 0 : return poGeom;
4793 : }
4794 : }
4795 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4796 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4797 2 : for (const auto poPoly : *poMP)
4798 : {
4799 1 : OGRErr eErr = OGRERR_NONE;
4800 1 : poTS->addGeometryDirectly(new OGRTriangle(*poPoly, eErr));
4801 : }
4802 1 : delete poGeom;
4803 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4804 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4805 1 : return poTS;
4806 : }
4807 1780 : else if (eType == wkbPolyhedralSurface && eTargetTypeFlat == wkbTIN)
4808 : {
4809 2 : OGRPolyhedralSurface *poPS = poGeom->toPolyhedralSurface();
4810 3 : for (const auto poPoly : *poPS)
4811 : {
4812 2 : const OGRLinearRing *poLR = poPoly->getExteriorRing();
4813 3 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4814 1 : poPoly->getNumInteriorRings() == 0))
4815 : {
4816 1 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4817 1 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4818 1 : return poGeom;
4819 : }
4820 : }
4821 1 : OGRTriangulatedSurface *poTS = new OGRTriangulatedSurface();
4822 1 : poTS->assignSpatialReference(poGeom->getSpatialReference());
4823 2 : for (const auto poPoly : *poPS)
4824 : {
4825 1 : OGRErr eErr = OGRERR_NONE;
4826 1 : poTS->addGeometryDirectly(new OGRTriangle(*poPoly, eErr));
4827 : }
4828 1 : delete poGeom;
4829 1 : poTS->set3D(OGR_GT_HasZ(eTargetType));
4830 1 : poTS->setMeasured(OGR_GT_HasM(eTargetType));
4831 1 : return poTS;
4832 : }
4833 :
4834 1778 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbTriangle)
4835 : {
4836 7 : OGRPolygon *poPoly = poGeom->toPolygon();
4837 7 : OGRLinearRing *poLR = poPoly->getExteriorRing();
4838 13 : if (!(poLR != nullptr && poLR->getNumPoints() == 4 &&
4839 6 : poPoly->getNumInteriorRings() == 0))
4840 : {
4841 1 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
4842 1 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
4843 1 : return poGeom;
4844 : }
4845 6 : OGRErr eErr = OGRERR_NONE;
4846 6 : OGRTriangle *poTriangle = new OGRTriangle(*poPoly, eErr);
4847 6 : delete poGeom;
4848 6 : poTriangle->set3D(OGR_GT_HasZ(eTargetType));
4849 6 : poTriangle->setMeasured(OGR_GT_HasM(eTargetType));
4850 6 : return poTriangle;
4851 : }
4852 :
4853 1772 : if (eTargetTypeFlat == wkbTriangle || eTargetTypeFlat == wkbTIN ||
4854 : eTargetTypeFlat == wkbPolyhedralSurface)
4855 : {
4856 9 : OGRwkbGeometryType eTempGeomType = wkbPolygon;
4857 9 : if (OGR_GT_HasZ(eTargetType))
4858 0 : eTempGeomType = OGR_GT_SetZ(eTempGeomType);
4859 9 : if (OGR_GT_HasM(eTargetType))
4860 1 : eTempGeomType = OGR_GT_SetM(eTempGeomType);
4861 9 : OGRGeometry *poPoly = forceTo(poGeom, eTempGeomType, papszOptions);
4862 9 : if (poPoly == poGeom)
4863 0 : return poGeom;
4864 9 : return forceTo(poPoly, eTargetType, papszOptions);
4865 : }
4866 :
4867 1763 : if (eType == wkbTriangle && eTargetTypeFlat == wkbGeometryCollection)
4868 : {
4869 1 : OGRGeometryCollection *poGC = new OGRGeometryCollection();
4870 1 : poGC->assignSpatialReference(poGeom->getSpatialReference());
4871 1 : poGC->addGeometryDirectly(poGeom);
4872 1 : poGC->set3D(OGR_GT_HasZ(eTargetType));
4873 1 : poGC->setMeasured(OGR_GT_HasM(eTargetType));
4874 1 : return poGC;
4875 : }
4876 :
4877 : // Promote single to multi.
4878 3252 : if (!OGR_GT_IsSubClassOf(eType, wkbGeometryCollection) &&
4879 1490 : OGR_GT_IsSubClassOf(OGR_GT_GetCollection(eType), eTargetType))
4880 : {
4881 201 : OGRGeometry *poRet = createGeometry(eTargetType);
4882 201 : if (poRet == nullptr)
4883 : {
4884 0 : delete poGeom;
4885 0 : return nullptr;
4886 : }
4887 201 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4888 201 : if (eType == wkbLineString)
4889 42 : poGeom = OGRCurve::CastToLineString(poGeom->toCurve());
4890 201 : poRet->toGeometryCollection()->addGeometryDirectly(poGeom);
4891 201 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4892 201 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4893 201 : return poRet;
4894 : }
4895 :
4896 1561 : const bool bIsCurve = CPL_TO_BOOL(OGR_GT_IsCurve(eType));
4897 1561 : if (bIsCurve && eTargetTypeFlat == wkbCompoundCurve)
4898 : {
4899 30 : auto poRet = OGRCurve::CastToCompoundCurve(poGeom->toCurve());
4900 30 : if (poRet)
4901 : {
4902 28 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4903 28 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4904 : }
4905 30 : return poRet;
4906 : }
4907 1531 : else if (bIsCurve && eTargetTypeFlat == wkbCurvePolygon)
4908 : {
4909 26 : OGRCurve *poCurve = poGeom->toCurve();
4910 26 : if (poCurve->getNumPoints() >= 3 && poCurve->get_IsClosed())
4911 : {
4912 18 : OGRCurvePolygon *poCP = new OGRCurvePolygon();
4913 18 : if (poCP->addRingDirectly(poCurve) == OGRERR_NONE)
4914 : {
4915 18 : poCP->assignSpatialReference(poGeom->getSpatialReference());
4916 18 : poCP->set3D(OGR_GT_HasZ(eTargetType));
4917 18 : poCP->setMeasured(OGR_GT_HasM(eTargetType));
4918 18 : return poCP;
4919 : }
4920 : else
4921 : {
4922 0 : delete poCP;
4923 : }
4924 8 : }
4925 : }
4926 1581 : else if (eType == wkbLineString &&
4927 76 : OGR_GT_IsSubClassOf(eTargetType, wkbMultiSurface))
4928 : {
4929 23 : OGRGeometry *poTmp = forceTo(poGeom, wkbPolygon, papszOptions);
4930 23 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4931 15 : return forceTo(poTmp, eTargetType, papszOptions);
4932 : }
4933 1482 : else if (bIsCurve && eTargetTypeFlat == wkbMultiSurface)
4934 : {
4935 10 : OGRGeometry *poTmp = forceTo(poGeom, wkbCurvePolygon, papszOptions);
4936 10 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4937 10 : return forceTo(poTmp, eTargetType, papszOptions);
4938 : }
4939 1472 : else if (bIsCurve && eTargetTypeFlat == wkbMultiPolygon)
4940 : {
4941 13 : OGRGeometry *poTmp = forceTo(poGeom, wkbPolygon, papszOptions);
4942 13 : if (wkbFlatten(poTmp->getGeometryType()) != eType)
4943 13 : return forceTo(poTmp, eTargetType, papszOptions);
4944 : }
4945 1459 : else if (eType == wkbTriangle && eTargetTypeFlat == wkbCurvePolygon)
4946 : {
4947 1 : auto poRet = OGRSurface::CastToCurvePolygon(
4948 : OGRTriangle::CastToPolygon(poGeom)->toSurface());
4949 1 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4950 1 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4951 1 : return poRet;
4952 : }
4953 1458 : else if (eType == wkbPolygon && eTargetTypeFlat == wkbCurvePolygon)
4954 : {
4955 18 : auto poRet = OGRSurface::CastToCurvePolygon(poGeom->toPolygon());
4956 18 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4957 18 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4958 18 : return poRet;
4959 : }
4960 1440 : else if (OGR_GT_IsSubClassOf(eType, wkbCurvePolygon) &&
4961 : eTargetTypeFlat == wkbCompoundCurve)
4962 : {
4963 15 : OGRCurvePolygon *poPoly = poGeom->toCurvePolygon();
4964 15 : if (poPoly->getNumInteriorRings() == 0)
4965 : {
4966 14 : OGRCurve *poRet = poPoly->stealExteriorRingCurve();
4967 14 : if (poRet)
4968 14 : poRet->assignSpatialReference(poGeom->getSpatialReference());
4969 14 : delete poPoly;
4970 14 : return forceTo(poRet, eTargetType, papszOptions);
4971 : }
4972 : }
4973 1425 : else if (eType == wkbMultiPolygon && eTargetTypeFlat == wkbMultiSurface)
4974 : {
4975 : auto poRet =
4976 14 : OGRMultiPolygon::CastToMultiSurface(poGeom->toMultiPolygon());
4977 14 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4978 14 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4979 14 : return poRet;
4980 : }
4981 1411 : else if (eType == wkbMultiLineString && eTargetTypeFlat == wkbMultiCurve)
4982 : {
4983 : auto poRet =
4984 9 : OGRMultiLineString::CastToMultiCurve(poGeom->toMultiLineString());
4985 9 : poRet->set3D(OGR_GT_HasZ(eTargetType));
4986 9 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
4987 9 : return poRet;
4988 : }
4989 1402 : else if (OGR_GT_IsSubClassOf(eType, wkbGeometryCollection))
4990 : {
4991 249 : OGRGeometryCollection *poGC = poGeom->toGeometryCollection();
4992 249 : if (poGC->getNumGeometries() == 1)
4993 : {
4994 158 : OGRGeometry *poSubGeom = poGC->getGeometryRef(0);
4995 158 : if (poSubGeom)
4996 : {
4997 158 : poSubGeom->assignSpatialReference(
4998 158 : poGeom->getSpatialReference());
4999 158 : poGC->removeGeometry(0, FALSE);
5000 : OGRGeometry *poRet =
5001 158 : forceTo(poSubGeom->clone(), eTargetType, papszOptions);
5002 158 : if (OGR_GT_IsSubClassOf(wkbFlatten(poRet->getGeometryType()),
5003 158 : eTargetType))
5004 : {
5005 123 : delete poGC;
5006 123 : delete poSubGeom;
5007 123 : return poRet;
5008 : }
5009 35 : poGC->addGeometryDirectly(poSubGeom);
5010 35 : poRet->set3D(OGR_GT_HasZ(eTargetType));
5011 35 : poRet->setMeasured(OGR_GT_HasM(eTargetType));
5012 35 : delete poRet;
5013 : }
5014 : }
5015 : }
5016 1267 : else if (OGR_GT_IsSubClassOf(eType, wkbCurvePolygon) &&
5017 114 : (OGR_GT_IsSubClassOf(eTargetType, wkbMultiSurface) ||
5018 102 : OGR_GT_IsSubClassOf(eTargetType, wkbMultiCurve)))
5019 : {
5020 43 : OGRCurvePolygon *poCP = poGeom->toCurvePolygon();
5021 43 : if (poCP->getNumInteriorRings() == 0)
5022 : {
5023 41 : OGRCurve *poRing = poCP->getExteriorRingCurve();
5024 41 : poRing->assignSpatialReference(poGeom->getSpatialReference());
5025 41 : OGRwkbGeometryType eRingType = poRing->getGeometryType();
5026 41 : OGRGeometry *poRingDup = poRing->clone();
5027 41 : OGRGeometry *poRet = forceTo(poRingDup, eTargetType, papszOptions);
5028 57 : if (poRet->getGeometryType() != eRingType &&
5029 16 : !(eTypeFlat == wkbPolygon &&
5030 : eTargetTypeFlat == wkbMultiLineString))
5031 : {
5032 29 : delete poCP;
5033 29 : return poRet;
5034 : }
5035 : else
5036 : {
5037 12 : delete poRet;
5038 : }
5039 : }
5040 : }
5041 :
5042 1271 : if (eTargetTypeFlat == wkbLineString)
5043 : {
5044 85 : poGeom = forceToLineString(poGeom);
5045 85 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5046 85 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5047 : }
5048 1186 : else if (eTargetTypeFlat == wkbPolygon)
5049 : {
5050 99 : poGeom = forceToPolygon(poGeom);
5051 99 : if (poGeom)
5052 : {
5053 99 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5054 99 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5055 : }
5056 : }
5057 1087 : else if (eTargetTypeFlat == wkbMultiPolygon)
5058 : {
5059 908 : poGeom = forceToMultiPolygon(poGeom);
5060 908 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5061 908 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5062 : }
5063 179 : else if (eTargetTypeFlat == wkbMultiLineString)
5064 : {
5065 37 : poGeom = forceToMultiLineString(poGeom);
5066 37 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5067 37 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5068 : }
5069 142 : else if (eTargetTypeFlat == wkbMultiPoint)
5070 : {
5071 22 : poGeom = forceToMultiPoint(poGeom);
5072 22 : poGeom->set3D(OGR_GT_HasZ(eTargetType));
5073 22 : poGeom->setMeasured(OGR_GT_HasM(eTargetType));
5074 : }
5075 :
5076 1271 : return poGeom;
5077 : }
5078 :
5079 : /************************************************************************/
5080 : /* OGR_G_ForceTo() */
5081 : /************************************************************************/
5082 :
5083 : /**
5084 : * \brief Convert to another geometry type
5085 : *
5086 : * This function is the same as the C++ method OGRGeometryFactory::forceTo().
5087 : *
5088 : * @param hGeom the input geometry - ownership is passed to the method.
5089 : * @param eTargetType target output geometry type.
5090 : * @param papszOptions options as a null-terminated list of strings or NULL.
5091 : * @return new geometry.
5092 : *
5093 : * @since GDAL 2.0
5094 : */
5095 :
5096 848 : OGRGeometryH OGR_G_ForceTo(OGRGeometryH hGeom, OGRwkbGeometryType eTargetType,
5097 : char **papszOptions)
5098 :
5099 : {
5100 848 : return OGRGeometry::ToHandle(OGRGeometryFactory::forceTo(
5101 848 : OGRGeometry::FromHandle(hGeom), eTargetType, papszOptions));
5102 : }
5103 :
5104 : /************************************************************************/
5105 : /* GetCurveParameters() */
5106 : /************************************************************************/
5107 :
5108 : /**
5109 : * \brief Returns the parameter of an arc circle.
5110 : *
5111 : * Angles are return in radians, with trigonometic convention (counter clock
5112 : * wise)
5113 : *
5114 : * @param x0 x of first point
5115 : * @param y0 y of first point
5116 : * @param x1 x of intermediate point
5117 : * @param y1 y of intermediate point
5118 : * @param x2 x of final point
5119 : * @param y2 y of final point
5120 : * @param R radius (output)
5121 : * @param cx x of arc center (output)
5122 : * @param cy y of arc center (output)
5123 : * @param alpha0 angle between center and first point, in radians (output)
5124 : * @param alpha1 angle between center and intermediate point, in radians
5125 : * (output)
5126 : * @param alpha2 angle between center and final point, in radians (output)
5127 : * @return TRUE if the points are not aligned and define an arc circle.
5128 : *
5129 : * @since GDAL 2.0
5130 : */
5131 :
5132 185879 : int OGRGeometryFactory::GetCurveParameters(double x0, double y0, double x1,
5133 : double y1, double x2, double y2,
5134 : double &R, double &cx, double &cy,
5135 : double &alpha0, double &alpha1,
5136 : double &alpha2)
5137 : {
5138 557637 : if (std::isnan(x0) || std::isnan(y0) || std::isnan(x1) || std::isnan(y1) ||
5139 557637 : std::isnan(x2) || std::isnan(y2))
5140 : {
5141 0 : return FALSE;
5142 : }
5143 :
5144 : // Circle.
5145 185879 : if (x0 == x2 && y0 == y2)
5146 : {
5147 153 : if (x0 != x1 || y0 != y1)
5148 : {
5149 152 : cx = (x0 + x1) / 2;
5150 152 : cy = (y0 + y1) / 2;
5151 152 : R = DISTANCE(cx, cy, x0, y0);
5152 : // Arbitrarily pick counter-clock-wise order (like PostGIS does).
5153 152 : alpha0 = atan2(y0 - cy, x0 - cx);
5154 152 : alpha1 = alpha0 + M_PI;
5155 152 : alpha2 = alpha0 + 2 * M_PI;
5156 152 : return TRUE;
5157 : }
5158 : else
5159 : {
5160 1 : return FALSE;
5161 : }
5162 : }
5163 :
5164 185726 : double dx01 = x1 - x0;
5165 185726 : double dy01 = y1 - y0;
5166 185726 : double dx12 = x2 - x1;
5167 185726 : double dy12 = y2 - y1;
5168 :
5169 : // Normalize above values so as to make sure we don't end up with
5170 : // computing a difference of too big values.
5171 185726 : double dfScale = fabs(dx01);
5172 185726 : if (fabs(dy01) > dfScale)
5173 93070 : dfScale = fabs(dy01);
5174 185726 : if (fabs(dx12) > dfScale)
5175 46253 : dfScale = fabs(dx12);
5176 185726 : if (fabs(dy12) > dfScale)
5177 46514 : dfScale = fabs(dy12);
5178 185726 : const double dfInvScale = 1.0 / dfScale;
5179 185726 : dx01 *= dfInvScale;
5180 185726 : dy01 *= dfInvScale;
5181 185726 : dx12 *= dfInvScale;
5182 185726 : dy12 *= dfInvScale;
5183 :
5184 185726 : const double det = dx01 * dy12 - dx12 * dy01;
5185 185726 : if (fabs(det) < 1.0e-8 || std::isnan(det))
5186 : {
5187 123 : return FALSE;
5188 : }
5189 185603 : const double x01_mid = (x0 + x1) * dfInvScale;
5190 185603 : const double x12_mid = (x1 + x2) * dfInvScale;
5191 185603 : const double y01_mid = (y0 + y1) * dfInvScale;
5192 185603 : const double y12_mid = (y1 + y2) * dfInvScale;
5193 185603 : const double c01 = dx01 * x01_mid + dy01 * y01_mid;
5194 185603 : const double c12 = dx12 * x12_mid + dy12 * y12_mid;
5195 185603 : cx = 0.5 * dfScale * (c01 * dy12 - c12 * dy01) / det;
5196 185603 : cy = 0.5 * dfScale * (-c01 * dx12 + c12 * dx01) / det;
5197 :
5198 185603 : alpha0 = atan2((y0 - cy) * dfInvScale, (x0 - cx) * dfInvScale);
5199 185603 : alpha1 = atan2((y1 - cy) * dfInvScale, (x1 - cx) * dfInvScale);
5200 185603 : alpha2 = atan2((y2 - cy) * dfInvScale, (x2 - cx) * dfInvScale);
5201 185603 : R = DISTANCE(cx, cy, x0, y0);
5202 :
5203 : // If det is negative, the orientation if clockwise.
5204 185603 : if (det < 0)
5205 : {
5206 94471 : if (alpha1 > alpha0)
5207 1295 : alpha1 -= 2 * M_PI;
5208 94471 : if (alpha2 > alpha1)
5209 3349 : alpha2 -= 2 * M_PI;
5210 : }
5211 : else
5212 : {
5213 91132 : if (alpha1 < alpha0)
5214 1275 : alpha1 += 2 * M_PI;
5215 91132 : if (alpha2 < alpha1)
5216 3145 : alpha2 += 2 * M_PI;
5217 : }
5218 :
5219 185603 : CPLAssert((alpha0 <= alpha1 && alpha1 <= alpha2) ||
5220 : (alpha0 >= alpha1 && alpha1 >= alpha2));
5221 :
5222 185603 : return TRUE;
5223 : }
5224 :
5225 : /************************************************************************/
5226 : /* OGRGeometryFactoryStrokeArc() */
5227 : /************************************************************************/
5228 :
5229 4307 : static void OGRGeometryFactoryStrokeArc(OGRLineString *poLine, double cx,
5230 : double cy, double R, double z0,
5231 : double z1, int bHasZ, double alpha0,
5232 : double alpha1, double dfStep,
5233 : int bStealthConstraints)
5234 : {
5235 4307 : const int nSign = dfStep > 0 ? 1 : -1;
5236 :
5237 : // Constant angle between all points, so as to not depend on winding order.
5238 4307 : const double dfNumSteps = fabs((alpha1 - alpha0) / dfStep) + 0.5;
5239 4307 : if (dfNumSteps >= std::numeric_limits<int>::max() ||
5240 4307 : dfNumSteps <= std::numeric_limits<int>::min() || std::isnan(dfNumSteps))
5241 : {
5242 0 : CPLError(CE_Warning, CPLE_AppDefined,
5243 : "OGRGeometryFactoryStrokeArc: bogus steps: "
5244 : "%lf %lf %lf %lf",
5245 : alpha0, alpha1, dfStep, dfNumSteps);
5246 0 : return;
5247 : }
5248 :
5249 4307 : int nSteps = static_cast<int>(dfNumSteps);
5250 4307 : if (bStealthConstraints)
5251 : {
5252 : // We need at least 6 intermediate vertex, and if more additional
5253 : // multiples of 2.
5254 4125 : if (nSteps < 1 + 6)
5255 120 : nSteps = 1 + 6;
5256 : else
5257 4005 : nSteps = 1 + 6 + 2 * ((nSteps - (1 + 6) + (2 - 1)) / 2);
5258 : }
5259 182 : else if (nSteps < 4)
5260 : {
5261 178 : nSteps = 4;
5262 : }
5263 4307 : dfStep = nSign * fabs((alpha1 - alpha0) / nSteps);
5264 4307 : double alpha = alpha0 + dfStep;
5265 :
5266 229149 : for (; (alpha - alpha1) * nSign < -1e-8; alpha += dfStep)
5267 : {
5268 224842 : const double dfX = cx + R * cos(alpha);
5269 224842 : const double dfY = cy + R * sin(alpha);
5270 224842 : if (bHasZ)
5271 : {
5272 9104 : const double z =
5273 9104 : z0 + (z1 - z0) * (alpha - alpha0) / (alpha1 - alpha0);
5274 9104 : poLine->addPoint(dfX, dfY, z);
5275 : }
5276 : else
5277 : {
5278 215738 : poLine->addPoint(dfX, dfY);
5279 : }
5280 : }
5281 : }
5282 :
5283 : /************************************************************************/
5284 : /* OGRGF_SetHiddenValue() */
5285 : /************************************************************************/
5286 :
5287 : // TODO(schwehr): Cleanup these static constants.
5288 : constexpr int HIDDEN_ALPHA_WIDTH = 32;
5289 : constexpr GUInt32 HIDDEN_ALPHA_SCALE =
5290 : static_cast<GUInt32>((static_cast<GUIntBig>(1) << HIDDEN_ALPHA_WIDTH) - 2);
5291 : constexpr int HIDDEN_ALPHA_HALF_WIDTH = (HIDDEN_ALPHA_WIDTH / 2);
5292 : constexpr int HIDDEN_ALPHA_HALF_MASK = (1 << HIDDEN_ALPHA_HALF_WIDTH) - 1;
5293 :
5294 : // Encode 16-bit nValue in the 8-lsb of dfX and dfY.
5295 :
5296 : #ifdef CPL_LSB
5297 : constexpr int DOUBLE_LSB_OFFSET = 0;
5298 : #else
5299 : constexpr int DOUBLE_LSB_OFFSET = 7;
5300 : #endif
5301 :
5302 224708 : static void OGRGF_SetHiddenValue(GUInt16 nValue, double &dfX, double &dfY)
5303 : {
5304 224708 : GByte abyData[8] = {};
5305 :
5306 224708 : memcpy(abyData, &dfX, sizeof(double));
5307 224708 : abyData[DOUBLE_LSB_OFFSET] = static_cast<GByte>(nValue & 0xFF);
5308 224708 : memcpy(&dfX, abyData, sizeof(double));
5309 :
5310 224708 : memcpy(abyData, &dfY, sizeof(double));
5311 224708 : abyData[DOUBLE_LSB_OFFSET] = static_cast<GByte>(nValue >> 8);
5312 224708 : memcpy(&dfY, abyData, sizeof(double));
5313 224708 : }
5314 :
5315 : /************************************************************************/
5316 : /* OGRGF_GetHiddenValue() */
5317 : /************************************************************************/
5318 :
5319 : // Decode 16-bit nValue from the 8-lsb of dfX and dfY.
5320 180748 : static GUInt16 OGRGF_GetHiddenValue(double dfX, double dfY)
5321 : {
5322 180748 : GByte abyData[8] = {};
5323 180748 : memcpy(abyData, &dfX, sizeof(double));
5324 180748 : GUInt16 nValue = abyData[DOUBLE_LSB_OFFSET];
5325 180748 : memcpy(abyData, &dfY, sizeof(double));
5326 180748 : nValue |= (abyData[DOUBLE_LSB_OFFSET] << 8);
5327 :
5328 180748 : return nValue;
5329 : }
5330 :
5331 : /************************************************************************/
5332 : /* OGRGF_NeedSwithArcOrder() */
5333 : /************************************************************************/
5334 :
5335 : // We need to define a full ordering between starting point and ending point
5336 : // whatever it is.
5337 9433 : static bool OGRGF_NeedSwithArcOrder(double x0, double y0, double x2, double y2)
5338 : {
5339 9433 : return x0 < x2 || (x0 == x2 && y0 < y2);
5340 : }
5341 :
5342 : /************************************************************************/
5343 : /* curveToLineString() */
5344 : /************************************************************************/
5345 :
5346 : /* clang-format off */
5347 : /**
5348 : * \brief Converts an arc circle into an approximate line string
5349 : *
5350 : * The arc circle is defined by a first point, an intermediate point and a
5351 : * final point.
5352 : *
5353 : * The provided dfMaxAngleStepSizeDegrees is a hint. The discretization
5354 : * algorithm may pick a slightly different value.
5355 : *
5356 : * So as to avoid gaps when rendering curve polygons that share common arcs,
5357 : * this method is guaranteed to return a line with reversed vertex if called
5358 : * with inverted first and final point, and identical intermediate point.
5359 : *
5360 : * @param x0 x of first point
5361 : * @param y0 y of first point
5362 : * @param z0 z of first point
5363 : * @param x1 x of intermediate point
5364 : * @param y1 y of intermediate point
5365 : * @param z1 z of intermediate point
5366 : * @param x2 x of final point
5367 : * @param y2 y of final point
5368 : * @param z2 z of final point
5369 : * @param bHasZ TRUE if z must be taken into account
5370 : * @param dfMaxAngleStepSizeDegrees the largest step in degrees along the
5371 : * arc, zero to use the default setting.
5372 : * @param papszOptions options as a null-terminated list of strings or NULL.
5373 : * Recognized options:
5374 : * <ul>
5375 : * <li>ADD_INTERMEDIATE_POINT=STEALTH/YES/NO (Default to STEALTH).
5376 : * Determine if and how the intermediate point must be output in the
5377 : * linestring. If set to STEALTH, no explicit intermediate point is
5378 : * added but its properties are encoded in low significant bits of
5379 : * intermediate points and OGRGeometryFactory::curveFromLineString() can
5380 : * decode them. This is the best compromise for round-tripping in OGR
5381 : * and better results with PostGIS
5382 : * <a href="http://postgis.org/docs/ST_LineToCurve.html">ST_LineToCurve()</a>.
5383 : * If set to YES, the intermediate point is explicitly added to the
5384 : * linestring. If set to NO, the intermediate point is not explicitly
5385 : * added.
5386 : * </li>
5387 : * </ul>
5388 : *
5389 : * @return the converted geometry (ownership to caller).
5390 : *
5391 : * @since GDAL 2.0
5392 : */
5393 : /* clang-format on */
5394 :
5395 6358 : OGRLineString *OGRGeometryFactory::curveToLineString(
5396 : double x0, double y0, double z0, double x1, double y1, double z1, double x2,
5397 : double y2, double z2, int bHasZ, double dfMaxAngleStepSizeDegrees,
5398 : const char *const *papszOptions)
5399 : {
5400 : // So as to make sure the same curve followed in both direction results
5401 : // in perfectly(=binary identical) symmetrical points.
5402 6358 : if (OGRGF_NeedSwithArcOrder(x0, y0, x2, y2))
5403 : {
5404 : OGRLineString *poLS =
5405 2141 : curveToLineString(x2, y2, z2, x1, y1, z1, x0, y0, z0, bHasZ,
5406 : dfMaxAngleStepSizeDegrees, papszOptions);
5407 2141 : poLS->reversePoints();
5408 2141 : return poLS;
5409 : }
5410 :
5411 4217 : double R = 0.0;
5412 4217 : double cx = 0.0;
5413 4217 : double cy = 0.0;
5414 4217 : double alpha0 = 0.0;
5415 4217 : double alpha1 = 0.0;
5416 4217 : double alpha2 = 0.0;
5417 :
5418 4217 : OGRLineString *poLine = new OGRLineString();
5419 4217 : bool bIsArc = true;
5420 4217 : if (!GetCurveParameters(x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1,
5421 : alpha2))
5422 : {
5423 89 : bIsArc = false;
5424 89 : cx = 0.0;
5425 89 : cy = 0.0;
5426 89 : R = 0.0;
5427 89 : alpha0 = 0.0;
5428 89 : alpha1 = 0.0;
5429 89 : alpha2 = 0.0;
5430 : }
5431 :
5432 4217 : const int nSign = alpha1 >= alpha0 ? 1 : -1;
5433 :
5434 : // support default arc step setting.
5435 4217 : if (dfMaxAngleStepSizeDegrees < 1e-6)
5436 : {
5437 4198 : dfMaxAngleStepSizeDegrees = OGRGF_GetDefaultStepSize();
5438 : }
5439 :
5440 4217 : double dfStep = dfMaxAngleStepSizeDegrees / 180 * M_PI;
5441 4217 : if (dfStep <= 0.01 / 180 * M_PI)
5442 : {
5443 0 : CPLDebug("OGR", "Too small arc step size: limiting to 0.01 degree.");
5444 0 : dfStep = 0.01 / 180 * M_PI;
5445 : }
5446 :
5447 4217 : dfStep *= nSign;
5448 :
5449 4217 : if (bHasZ)
5450 252 : poLine->addPoint(x0, y0, z0);
5451 : else
5452 3965 : poLine->addPoint(x0, y0);
5453 :
5454 4217 : bool bAddIntermediatePoint = false;
5455 4217 : bool bStealth = true;
5456 4223 : for (const char *const *papszIter = papszOptions; papszIter && *papszIter;
5457 : papszIter++)
5458 : {
5459 6 : char *pszKey = nullptr;
5460 6 : const char *pszValue = CPLParseNameValue(*papszIter, &pszKey);
5461 6 : if (pszKey != nullptr && EQUAL(pszKey, "ADD_INTERMEDIATE_POINT"))
5462 : {
5463 4 : if (EQUAL(pszValue, "YES") || EQUAL(pszValue, "TRUE") ||
5464 3 : EQUAL(pszValue, "ON"))
5465 : {
5466 1 : bAddIntermediatePoint = true;
5467 1 : bStealth = false;
5468 : }
5469 3 : else if (EQUAL(pszValue, "NO") || EQUAL(pszValue, "FALSE") ||
5470 1 : EQUAL(pszValue, "OFF"))
5471 : {
5472 2 : bAddIntermediatePoint = false;
5473 2 : bStealth = false;
5474 : }
5475 : else if (EQUAL(pszValue, "STEALTH"))
5476 : {
5477 : // default.
5478 : }
5479 : }
5480 : else
5481 : {
5482 2 : CPLError(CE_Warning, CPLE_NotSupported, "Unsupported option: %s",
5483 : *papszIter);
5484 : }
5485 6 : CPLFree(pszKey);
5486 : }
5487 :
5488 4217 : if (!bIsArc || bAddIntermediatePoint)
5489 : {
5490 90 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z0, z1, bHasZ, alpha0,
5491 : alpha1, dfStep, FALSE);
5492 :
5493 90 : if (bHasZ)
5494 25 : poLine->addPoint(x1, y1, z1);
5495 : else
5496 65 : poLine->addPoint(x1, y1);
5497 :
5498 90 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z1, z2, bHasZ, alpha1,
5499 : alpha2, dfStep, FALSE);
5500 : }
5501 : else
5502 : {
5503 4127 : OGRGeometryFactoryStrokeArc(poLine, cx, cy, R, z0, z2, bHasZ, alpha0,
5504 : alpha2, dfStep, bStealth);
5505 :
5506 4127 : if (bStealth && poLine->getNumPoints() > 6)
5507 : {
5508 : // 'Hide' the angle of the intermediate point in the 8
5509 : // low-significant bits of the x, y of the first 2 computed points
5510 : // (so 32 bits), then put 0xFF, and on the last couple points put
5511 : // again the angle but in reverse order, so that overall the
5512 : // low-significant bits of all the points are symmetrical w.r.t the
5513 : // mid-point.
5514 4125 : const double dfRatio = (alpha1 - alpha0) / (alpha2 - alpha0);
5515 4125 : double dfAlphaRatio = 0.5 + HIDDEN_ALPHA_SCALE * dfRatio;
5516 4125 : if (dfAlphaRatio < 0.0)
5517 : {
5518 0 : CPLError(CE_Warning, CPLE_AppDefined, "AlphaRation < 0: %lf",
5519 : dfAlphaRatio);
5520 0 : dfAlphaRatio *= -1;
5521 : }
5522 8250 : else if (dfAlphaRatio >= std::numeric_limits<GUInt32>::max() ||
5523 4125 : std::isnan(dfAlphaRatio))
5524 : {
5525 0 : CPLError(CE_Warning, CPLE_AppDefined,
5526 : "AlphaRatio too large: %lf", dfAlphaRatio);
5527 0 : dfAlphaRatio = std::numeric_limits<GUInt32>::max();
5528 : }
5529 4125 : const GUInt32 nAlphaRatio = static_cast<GUInt32>(dfAlphaRatio);
5530 4125 : const GUInt16 nAlphaRatioLow = nAlphaRatio & HIDDEN_ALPHA_HALF_MASK;
5531 4125 : const GUInt16 nAlphaRatioHigh =
5532 4125 : nAlphaRatio >> HIDDEN_ALPHA_HALF_WIDTH;
5533 : // printf("alpha0=%f, alpha1=%f, alpha2=%f, dfRatio=%f, "/*ok*/
5534 : // "nAlphaRatio = %u\n",
5535 : // alpha0, alpha1, alpha2, dfRatio, nAlphaRatio);
5536 :
5537 4125 : CPLAssert(((poLine->getNumPoints() - 1 - 6) % 2) == 0);
5538 :
5539 116479 : for (int i = 1; i + 1 < poLine->getNumPoints(); i += 2)
5540 : {
5541 112354 : GUInt16 nVal = 0xFFFF;
5542 :
5543 112354 : double dfX = poLine->getX(i);
5544 112354 : double dfY = poLine->getY(i);
5545 112354 : if (i == 1)
5546 4125 : nVal = nAlphaRatioLow;
5547 108229 : else if (i == poLine->getNumPoints() - 2)
5548 4125 : nVal = nAlphaRatioHigh;
5549 112354 : OGRGF_SetHiddenValue(nVal, dfX, dfY);
5550 112354 : poLine->setPoint(i, dfX, dfY);
5551 :
5552 112354 : dfX = poLine->getX(i + 1);
5553 112354 : dfY = poLine->getY(i + 1);
5554 112354 : if (i == 1)
5555 4125 : nVal = nAlphaRatioHigh;
5556 108229 : else if (i == poLine->getNumPoints() - 2)
5557 4125 : nVal = nAlphaRatioLow;
5558 112354 : OGRGF_SetHiddenValue(nVal, dfX, dfY);
5559 112354 : poLine->setPoint(i + 1, dfX, dfY);
5560 : }
5561 : }
5562 : }
5563 :
5564 4217 : if (bHasZ)
5565 252 : poLine->addPoint(x2, y2, z2);
5566 : else
5567 3965 : poLine->addPoint(x2, y2);
5568 :
5569 4217 : return poLine;
5570 : }
5571 :
5572 : /************************************************************************/
5573 : /* OGRGF_FixAngle() */
5574 : /************************************************************************/
5575 :
5576 : // Fix dfAngle by offsets of 2 PI so that it lies between dfAngleStart and
5577 : // dfAngleStop, whatever their respective order.
5578 179548 : static double OGRGF_FixAngle(double dfAngleStart, double dfAngleStop,
5579 : double dfAngle)
5580 : {
5581 179548 : if (dfAngleStart < dfAngleStop)
5582 : {
5583 122937 : while (dfAngle <= dfAngleStart + 1e-8)
5584 34732 : dfAngle += 2 * M_PI;
5585 : }
5586 : else
5587 : {
5588 125233 : while (dfAngle >= dfAngleStart - 1e-8)
5589 33890 : dfAngle -= 2 * M_PI;
5590 : }
5591 179548 : return dfAngle;
5592 : }
5593 :
5594 : /************************************************************************/
5595 : /* OGRGF_DetectArc() */
5596 : /************************************************************************/
5597 :
5598 : // #define VERBOSE_DEBUG_CURVEFROMLINESTRING
5599 :
5600 12243 : static inline bool IS_ALMOST_INTEGER(double x)
5601 : {
5602 12243 : const double val = fabs(x - floor(x + 0.5));
5603 12243 : return val < 1.0e-8;
5604 : }
5605 :
5606 3476 : static int OGRGF_DetectArc(const OGRLineString *poLS, int i,
5607 : OGRCompoundCurve *&poCC, OGRCircularString *&poCS,
5608 : OGRLineString *&poLSNew)
5609 : {
5610 3476 : if (i + 3 >= poLS->getNumPoints())
5611 305 : return -1;
5612 :
5613 6342 : OGRPoint p0;
5614 6342 : OGRPoint p1;
5615 6342 : OGRPoint p2;
5616 3171 : poLS->getPoint(i, &p0);
5617 3171 : poLS->getPoint(i + 1, &p1);
5618 3171 : poLS->getPoint(i + 2, &p2);
5619 3171 : double R_1 = 0.0;
5620 3171 : double cx_1 = 0.0;
5621 3171 : double cy_1 = 0.0;
5622 3171 : double alpha0_1 = 0.0;
5623 3171 : double alpha1_1 = 0.0;
5624 3171 : double alpha2_1 = 0.0;
5625 6335 : if (!(OGRGeometryFactory::GetCurveParameters(
5626 : p0.getX(), p0.getY(), p1.getX(), p1.getY(), p2.getX(), p2.getY(),
5627 : R_1, cx_1, cy_1, alpha0_1, alpha1_1, alpha2_1) &&
5628 3164 : fabs(alpha2_1 - alpha0_1) < 2.0 * 20.0 / 180.0 * M_PI))
5629 : {
5630 24 : return -1;
5631 : }
5632 :
5633 3147 : const double dfDeltaAlpha10 = alpha1_1 - alpha0_1;
5634 3147 : const double dfDeltaAlpha21 = alpha2_1 - alpha1_1;
5635 : const double dfMaxDeltaAlpha =
5636 3147 : std::max(fabs(dfDeltaAlpha10), fabs(dfDeltaAlpha21));
5637 : GUInt32 nAlphaRatioRef =
5638 3147 : OGRGF_GetHiddenValue(p1.getX(), p1.getY()) |
5639 3147 : (OGRGF_GetHiddenValue(p2.getX(), p2.getY()) << HIDDEN_ALPHA_HALF_WIDTH);
5640 3147 : bool bFoundFFFFFFFFPattern = false;
5641 3147 : bool bFoundReversedAlphaRatioRef = false;
5642 3147 : bool bValidAlphaRatio = nAlphaRatioRef > 0 && nAlphaRatioRef < 0xFFFFFFFF;
5643 3147 : int nCountValidAlphaRatio = 1;
5644 :
5645 3147 : double dfScale = std::max(1.0, R_1);
5646 3147 : dfScale = std::max(dfScale, fabs(cx_1));
5647 3147 : dfScale = std::max(dfScale, fabs(cy_1));
5648 3147 : dfScale = pow(10.0, ceil(log10(dfScale)));
5649 3147 : const double dfInvScale = 1.0 / dfScale;
5650 :
5651 3147 : const int bInitialConstantStep =
5652 3147 : (fabs(dfDeltaAlpha10 - dfDeltaAlpha21) / dfMaxDeltaAlpha) < 1.0e-4;
5653 3147 : const double dfDeltaEpsilon =
5654 3147 : bInitialConstantStep ? dfMaxDeltaAlpha * 1e-4 : dfMaxDeltaAlpha / 10;
5655 :
5656 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5657 : printf("----------------------------\n"); /*ok*/
5658 : printf("Curve beginning at offset i = %d\n", i); /*ok*/
5659 : printf("Initial alpha ratio = %u\n", nAlphaRatioRef); /*ok*/
5660 : /*ok*/ printf("Initial R = %.16g, cx = %.16g, cy = %.16g\n", R_1, cx_1,
5661 : cy_1);
5662 : printf("dfScale = %f\n", dfScale); /*ok*/
5663 : printf("bInitialConstantStep = %d, " /*ok*/
5664 : "fabs(dfDeltaAlpha10 - dfDeltaAlpha21)=%.8g, "
5665 : "dfMaxDeltaAlpha = %.8f, "
5666 : "dfDeltaEpsilon = %.8f (%.8f)\n",
5667 : bInitialConstantStep, fabs(dfDeltaAlpha10 - dfDeltaAlpha21),
5668 : dfMaxDeltaAlpha, dfDeltaEpsilon, 1.0 / 180.0 * M_PI);
5669 : #endif
5670 3147 : int iMidPoint = -1;
5671 3147 : double dfLastValidAlpha = alpha2_1;
5672 :
5673 3147 : double dfLastLogRelDiff = 0;
5674 :
5675 6294 : OGRPoint p3;
5676 3147 : int j = i + 1; // Used after for.
5677 181121 : for (; j + 2 < poLS->getNumPoints(); j++)
5678 : {
5679 178073 : poLS->getPoint(j, &p1);
5680 178073 : poLS->getPoint(j + 1, &p2);
5681 178073 : poLS->getPoint(j + 2, &p3);
5682 178073 : double R_2 = 0.0;
5683 178073 : double cx_2 = 0.0;
5684 178073 : double cy_2 = 0.0;
5685 178073 : double alpha0_2 = 0.0;
5686 178073 : double alpha1_2 = 0.0;
5687 178073 : double alpha2_2 = 0.0;
5688 : // Check that the new candidate arc shares the same
5689 : // radius, center and winding order.
5690 178073 : if (!(OGRGeometryFactory::GetCurveParameters(
5691 : p1.getX(), p1.getY(), p2.getX(), p2.getY(), p3.getX(),
5692 : p3.getY(), R_2, cx_2, cy_2, alpha0_2, alpha1_2, alpha2_2)))
5693 : {
5694 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5695 : printf("End of curve at j=%d\n : straight line", j); /*ok*/
5696 : #endif
5697 99 : break;
5698 : }
5699 :
5700 178065 : const double dfRelDiffR = fabs(R_1 - R_2) * dfInvScale;
5701 178065 : const double dfRelDiffCx = fabs(cx_1 - cx_2) * dfInvScale;
5702 178065 : const double dfRelDiffCy = fabs(cy_1 - cy_2) * dfInvScale;
5703 :
5704 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5705 : printf("j=%d: R = %.16g, cx = %.16g, cy = %.16g, " /*ok*/
5706 : "rel_diff_R=%.8g rel_diff_cx=%.8g rel_diff_cy=%.8g\n",
5707 : j, R_2, cx_2, cy_2, dfRelDiffR, dfRelDiffCx, dfRelDiffCy);
5708 : #endif
5709 :
5710 178065 : if (dfRelDiffR > 1.0e-7 || dfRelDiffCx > 1.0e-7 ||
5711 177996 : dfRelDiffCy > 1.0e-7 ||
5712 177996 : dfDeltaAlpha10 * (alpha1_2 - alpha0_2) < 0.0)
5713 : {
5714 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5715 : printf("End of curve at j=%d\n", j); /*ok*/
5716 : #endif
5717 : break;
5718 : }
5719 :
5720 177996 : if (dfRelDiffR > 0.0 && dfRelDiffCx > 0.0 && dfRelDiffCy > 0.0)
5721 : {
5722 : const double dfLogRelDiff = std::min(
5723 355964 : std::min(fabs(log10(dfRelDiffR)), fabs(log10(dfRelDiffCx))),
5724 177982 : fabs(log10(dfRelDiffCy)));
5725 : // printf("dfLogRelDiff = %f, dfLastLogRelDiff=%f, "/*ok*/
5726 : // "dfLogRelDiff - dfLastLogRelDiff=%f\n",
5727 : // dfLogRelDiff, dfLastLogRelDiff,
5728 : // dfLogRelDiff - dfLastLogRelDiff);
5729 177982 : if (dfLogRelDiff > 0.0 && dfLastLogRelDiff >= 8.0 &&
5730 2 : dfLogRelDiff <= 8.0 && dfLogRelDiff < dfLastLogRelDiff - 2.0)
5731 : {
5732 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5733 : printf("End of curve at j=%d. Significant different in " /*ok*/
5734 : "relative error w.r.t previous points\n",
5735 : j);
5736 : #endif
5737 2 : break;
5738 : }
5739 177980 : dfLastLogRelDiff = dfLogRelDiff;
5740 : }
5741 :
5742 177994 : const double dfStep10 = fabs(alpha1_2 - alpha0_2);
5743 177994 : const double dfStep21 = fabs(alpha2_2 - alpha1_2);
5744 : // Check that the angle step is consistent with the original step.
5745 177994 : if (!(dfStep10 < 2.0 * dfMaxDeltaAlpha &&
5746 177994 : dfStep21 < 2.0 * dfMaxDeltaAlpha))
5747 : {
5748 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5749 : printf("End of curve at j=%d: dfStep10=%f, dfStep21=%f, " /*ok*/
5750 : "2*dfMaxDeltaAlpha=%f\n",
5751 : j, dfStep10, dfStep21, 2 * dfMaxDeltaAlpha);
5752 : #endif
5753 : break;
5754 : }
5755 :
5756 177993 : if (bValidAlphaRatio && j > i + 1 && (i % 2) != (j % 2))
5757 : {
5758 : const GUInt32 nAlphaRatioReversed =
5759 87227 : (OGRGF_GetHiddenValue(p1.getX(), p1.getY())
5760 174454 : << HIDDEN_ALPHA_HALF_WIDTH) |
5761 87227 : (OGRGF_GetHiddenValue(p2.getX(), p2.getY()));
5762 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5763 : printf("j=%d, nAlphaRatioReversed = %u\n", /*ok*/
5764 : j, nAlphaRatioReversed);
5765 : #endif
5766 87227 : if (!bFoundFFFFFFFFPattern && nAlphaRatioReversed == 0xFFFFFFFF)
5767 : {
5768 3075 : bFoundFFFFFFFFPattern = true;
5769 3075 : nCountValidAlphaRatio++;
5770 : }
5771 84152 : else if (bFoundFFFFFFFFPattern && !bFoundReversedAlphaRatioRef &&
5772 : nAlphaRatioReversed == 0xFFFFFFFF)
5773 : {
5774 81050 : nCountValidAlphaRatio++;
5775 : }
5776 3102 : else if (bFoundFFFFFFFFPattern && !bFoundReversedAlphaRatioRef &&
5777 : nAlphaRatioReversed == nAlphaRatioRef)
5778 : {
5779 3075 : bFoundReversedAlphaRatioRef = true;
5780 3075 : nCountValidAlphaRatio++;
5781 : }
5782 : else
5783 : {
5784 27 : if (bInitialConstantStep &&
5785 26 : fabs(dfLastValidAlpha - alpha0_1) >= M_PI &&
5786 : nCountValidAlphaRatio > 10)
5787 : {
5788 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5789 : printf("End of curve at j=%d: " /*ok*/
5790 : "fabs(dfLastValidAlpha - alpha0_1)=%f, "
5791 : "nCountValidAlphaRatio=%d\n",
5792 : j, fabs(dfLastValidAlpha - alpha0_1),
5793 : nCountValidAlphaRatio);
5794 : #endif
5795 19 : if (dfLastValidAlpha - alpha0_1 > 0)
5796 : {
5797 21 : while (dfLastValidAlpha - alpha0_1 - dfMaxDeltaAlpha -
5798 14 : M_PI >
5799 14 : -dfMaxDeltaAlpha / 10)
5800 : {
5801 7 : dfLastValidAlpha -= dfMaxDeltaAlpha;
5802 7 : j--;
5803 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5804 : printf(/*ok*/
5805 : "--> corrected as fabs(dfLastValidAlpha - "
5806 : "alpha0_1)=%f, j=%d\n",
5807 : fabs(dfLastValidAlpha - alpha0_1), j);
5808 : #endif
5809 : }
5810 : }
5811 : else
5812 : {
5813 36 : while (dfLastValidAlpha - alpha0_1 + dfMaxDeltaAlpha +
5814 24 : M_PI <
5815 24 : dfMaxDeltaAlpha / 10)
5816 : {
5817 12 : dfLastValidAlpha += dfMaxDeltaAlpha;
5818 12 : j--;
5819 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5820 : printf(/*ok*/
5821 : "--> corrected as fabs(dfLastValidAlpha - "
5822 : "alpha0_1)=%f, j=%d\n",
5823 : fabs(dfLastValidAlpha - alpha0_1), j);
5824 : #endif
5825 : }
5826 : }
5827 19 : poLS->getPoint(j + 1, &p2);
5828 19 : break;
5829 : }
5830 :
5831 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5832 : printf("j=%d, nAlphaRatioReversed = %u --> inconsistent " /*ok*/
5833 : "values across arc. Don't use it\n",
5834 : j, nAlphaRatioReversed);
5835 : #endif
5836 8 : bValidAlphaRatio = false;
5837 : }
5838 : }
5839 :
5840 : // Correct current end angle, consistently with start angle.
5841 177974 : dfLastValidAlpha = OGRGF_FixAngle(alpha0_1, alpha1_1, alpha2_2);
5842 :
5843 : // Try to detect the precise intermediate point of the
5844 : // arc circle by detecting irregular angle step
5845 : // This is OK if we don't detect the right point or fail
5846 : // to detect it.
5847 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5848 : printf("j=%d A(0,1)-maxDelta=%.8f A(1,2)-maxDelta=%.8f " /*ok*/
5849 : "x1=%.8f y1=%.8f x2=%.8f y2=%.8f x3=%.8f y3=%.8f\n",
5850 : j, fabs(dfStep10 - dfMaxDeltaAlpha),
5851 : fabs(dfStep21 - dfMaxDeltaAlpha), p1.getX(), p1.getY(),
5852 : p2.getX(), p2.getY(), p3.getX(), p3.getY());
5853 : #endif
5854 177974 : if (j > i + 1 && iMidPoint < 0 && dfDeltaEpsilon < 1.0 / 180.0 * M_PI)
5855 : {
5856 174495 : if (fabs(dfStep10 - dfMaxDeltaAlpha) > dfDeltaEpsilon)
5857 8 : iMidPoint = j + ((bInitialConstantStep) ? 0 : 1);
5858 174487 : else if (fabs(dfStep21 - dfMaxDeltaAlpha) > dfDeltaEpsilon)
5859 4 : iMidPoint = j + ((bInitialConstantStep) ? 1 : 2);
5860 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5861 : if (iMidPoint >= 0)
5862 : {
5863 : OGRPoint pMid;
5864 : poLS->getPoint(iMidPoint, &pMid);
5865 : printf("Midpoint detected at j = %d, iMidPoint = %d, " /*ok*/
5866 : "x=%.8f y=%.8f\n",
5867 : j, iMidPoint, pMid.getX(), pMid.getY());
5868 : }
5869 : #endif
5870 : }
5871 : }
5872 :
5873 : // Take a minimum threshold of consecutive points
5874 : // on the arc to avoid false positives.
5875 3147 : if (j < i + 3)
5876 61 : return -1;
5877 :
5878 3086 : bValidAlphaRatio &= bFoundFFFFFFFFPattern && bFoundReversedAlphaRatioRef;
5879 :
5880 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5881 : printf("bValidAlphaRatio=%d bFoundFFFFFFFFPattern=%d, " /*ok*/
5882 : "bFoundReversedAlphaRatioRef=%d\n",
5883 : static_cast<int>(bValidAlphaRatio),
5884 : static_cast<int>(bFoundFFFFFFFFPattern),
5885 : static_cast<int>(bFoundReversedAlphaRatioRef));
5886 : printf("alpha0_1=%f dfLastValidAlpha=%f\n", /*ok*/
5887 : alpha0_1, dfLastValidAlpha);
5888 : #endif
5889 :
5890 3086 : if (poLSNew != nullptr)
5891 : {
5892 11 : double dfScale2 = std::max(1.0, fabs(p0.getX()));
5893 11 : dfScale2 = std::max(dfScale2, fabs(p0.getY()));
5894 : // Not strictly necessary, but helps having 'clean' lines without
5895 : // duplicated points.
5896 11 : constexpr double dfToleranceEps =
5897 : OGRCompoundCurve::DEFAULT_TOLERANCE_EPSILON;
5898 11 : if (fabs(poLSNew->getX(poLSNew->getNumPoints() - 1) - p0.getX()) >
5899 12 : dfToleranceEps * dfScale2 ||
5900 1 : fabs(poLSNew->getY(poLSNew->getNumPoints() - 1) - p0.getY()) >
5901 1 : dfToleranceEps * dfScale2)
5902 10 : poLSNew->addPoint(&p0);
5903 11 : if (poLSNew->getNumPoints() >= 2)
5904 : {
5905 10 : if (poCC == nullptr)
5906 3 : poCC = new OGRCompoundCurve();
5907 10 : poCC->addCurveDirectly(poLSNew);
5908 : }
5909 : else
5910 1 : delete poLSNew;
5911 11 : poLSNew = nullptr;
5912 : }
5913 :
5914 3086 : if (poCS == nullptr)
5915 : {
5916 3062 : poCS = new OGRCircularString();
5917 3062 : poCS->addPoint(&p0);
5918 : }
5919 :
5920 3086 : OGRPoint *poFinalPoint = (j + 2 >= poLS->getNumPoints()) ? &p3 : &p2;
5921 :
5922 3086 : double dfXMid = 0.0;
5923 3086 : double dfYMid = 0.0;
5924 3086 : double dfZMid = 0.0;
5925 3086 : if (bValidAlphaRatio)
5926 : {
5927 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5928 : printf("Using alpha ratio...\n"); /*ok*/
5929 : #endif
5930 3075 : double dfAlphaMid = 0.0;
5931 3075 : if (OGRGF_NeedSwithArcOrder(p0.getX(), p0.getY(), poFinalPoint->getX(),
5932 : poFinalPoint->getY()))
5933 : {
5934 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
5935 : printf("Switching angles\n"); /*ok*/
5936 : #endif
5937 1526 : dfAlphaMid = dfLastValidAlpha + nAlphaRatioRef *
5938 1526 : (alpha0_1 - dfLastValidAlpha) /
5939 : HIDDEN_ALPHA_SCALE;
5940 1526 : dfAlphaMid = OGRGF_FixAngle(alpha0_1, dfLastValidAlpha, dfAlphaMid);
5941 : }
5942 : else
5943 : {
5944 1549 : dfAlphaMid = alpha0_1 + nAlphaRatioRef *
5945 1549 : (dfLastValidAlpha - alpha0_1) /
5946 : HIDDEN_ALPHA_SCALE;
5947 : }
5948 :
5949 3075 : dfXMid = cx_1 + R_1 * cos(dfAlphaMid);
5950 3075 : dfYMid = cy_1 + R_1 * sin(dfAlphaMid);
5951 :
5952 3075 : if (poLS->getCoordinateDimension() == 3)
5953 : {
5954 2 : double dfLastAlpha = 0.0;
5955 2 : double dfLastZ = 0.0;
5956 2 : int k = i; // Used after for.
5957 48 : for (; k < j + 2; k++)
5958 : {
5959 48 : OGRPoint p;
5960 48 : poLS->getPoint(k, &p);
5961 48 : double dfAlpha = atan2(p.getY() - cy_1, p.getX() - cx_1);
5962 48 : dfAlpha = OGRGF_FixAngle(alpha0_1, dfLastValidAlpha, dfAlpha);
5963 48 : if (k > i &&
5964 46 : ((dfAlpha < dfLastValidAlpha && dfAlphaMid < dfAlpha) ||
5965 23 : (dfAlpha > dfLastValidAlpha && dfAlphaMid > dfAlpha)))
5966 : {
5967 2 : const double dfRatio =
5968 2 : (dfAlphaMid - dfLastAlpha) / (dfAlpha - dfLastAlpha);
5969 2 : dfZMid = (1 - dfRatio) * dfLastZ + dfRatio * p.getZ();
5970 2 : break;
5971 : }
5972 46 : dfLastAlpha = dfAlpha;
5973 46 : dfLastZ = p.getZ();
5974 : }
5975 2 : if (k == j + 2)
5976 0 : dfZMid = dfLastZ;
5977 2 : if (IS_ALMOST_INTEGER(dfZMid))
5978 2 : dfZMid = static_cast<int>(floor(dfZMid + 0.5));
5979 : }
5980 :
5981 : // A few rounding strategies in case the mid point was at "exact"
5982 : // coordinates.
5983 3075 : if (R_1 > 1e-5)
5984 : {
5985 : const bool bStartEndInteger =
5986 9185 : IS_ALMOST_INTEGER(p0.getX()) && IS_ALMOST_INTEGER(p0.getY()) &&
5987 9185 : IS_ALMOST_INTEGER(poFinalPoint->getX()) &&
5988 3056 : IS_ALMOST_INTEGER(poFinalPoint->getY());
5989 3069 : if (bStartEndInteger &&
5990 3056 : fabs(dfXMid - floor(dfXMid + 0.5)) / dfScale < 1e-4 &&
5991 3037 : fabs(dfYMid - floor(dfYMid + 0.5)) / dfScale < 1e-4)
5992 : {
5993 3037 : dfXMid = static_cast<int>(floor(dfXMid + 0.5));
5994 3037 : dfYMid = static_cast<int>(floor(dfYMid + 0.5));
5995 : // Sometimes rounding to closest is not best approach
5996 : // Try neighbouring integers to look for the one that
5997 : // minimize the error w.r.t to the arc center
5998 : // But only do that if the radius is greater than
5999 : // the magnitude of the delta that we will try!
6000 : double dfBestRError =
6001 3037 : fabs(R_1 - DISTANCE(dfXMid, dfYMid, cx_1, cy_1));
6002 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6003 : printf("initial_error=%f\n", dfBestRError); /*ok*/
6004 : #endif
6005 3037 : int iBestX = 0;
6006 3037 : int iBestY = 0;
6007 3037 : if (dfBestRError > 0.001 && R_1 > 2)
6008 : {
6009 3 : int nSearchRadius = 1;
6010 : // Extend the search radius if the arc circle radius
6011 : // is much higher than the coordinate values.
6012 : double dfMaxCoords =
6013 3 : std::max(fabs(p0.getX()), fabs(p0.getY()));
6014 3 : dfMaxCoords = std::max(dfMaxCoords, poFinalPoint->getX());
6015 3 : dfMaxCoords = std::max(dfMaxCoords, poFinalPoint->getY());
6016 3 : dfMaxCoords = std::max(dfMaxCoords, dfXMid);
6017 3 : dfMaxCoords = std::max(dfMaxCoords, dfYMid);
6018 3 : if (R_1 > dfMaxCoords * 1000)
6019 3 : nSearchRadius = 100;
6020 0 : else if (R_1 > dfMaxCoords * 10)
6021 0 : nSearchRadius = 10;
6022 606 : for (int iY = -nSearchRadius; iY <= nSearchRadius; iY++)
6023 : {
6024 121806 : for (int iX = -nSearchRadius; iX <= nSearchRadius; iX++)
6025 : {
6026 121203 : const double dfCandidateX = dfXMid + iX;
6027 121203 : const double dfCandidateY = dfYMid + iY;
6028 121203 : if (fabs(dfCandidateX - p0.getX()) < 1e-8 &&
6029 0 : fabs(dfCandidateY - p0.getY()) < 1e-8)
6030 0 : continue;
6031 121203 : if (fabs(dfCandidateX - poFinalPoint->getX()) <
6032 121203 : 1e-8 &&
6033 0 : fabs(dfCandidateY - poFinalPoint->getY()) <
6034 : 1e-8)
6035 0 : continue;
6036 : const double dfRError =
6037 121203 : fabs(R_1 - DISTANCE(dfCandidateX, dfCandidateY,
6038 121203 : cx_1, cy_1));
6039 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6040 : printf("x=%d y=%d error=%f besterror=%f\n", /*ok*/
6041 : static_cast<int>(dfXMid + iX),
6042 : static_cast<int>(dfYMid + iY), dfRError,
6043 : dfBestRError);
6044 : #endif
6045 121203 : if (dfRError < dfBestRError)
6046 : {
6047 20 : iBestX = iX;
6048 20 : iBestY = iY;
6049 20 : dfBestRError = dfRError;
6050 : }
6051 : }
6052 : }
6053 : }
6054 3037 : dfXMid += iBestX;
6055 3037 : dfYMid += iBestY;
6056 : }
6057 : else
6058 : {
6059 : // Limit the number of significant figures in decimal
6060 : // representation.
6061 32 : if (fabs(dfXMid) < 100000000.0)
6062 : {
6063 32 : dfXMid =
6064 32 : static_cast<GIntBig>(floor(dfXMid * 100000000 + 0.5)) /
6065 : 100000000.0;
6066 : }
6067 32 : if (fabs(dfYMid) < 100000000.0)
6068 : {
6069 32 : dfYMid =
6070 32 : static_cast<GIntBig>(floor(dfYMid * 100000000 + 0.5)) /
6071 : 100000000.0;
6072 : }
6073 : }
6074 : }
6075 :
6076 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6077 : printf("dfAlphaMid=%f, x_mid = %f, y_mid = %f\n", /*ok*/
6078 : dfLastValidAlpha, dfXMid, dfYMid);
6079 : #endif
6080 : }
6081 :
6082 : // If this is a full circle of a non-polygonal zone, we must
6083 : // use a 5-point representation to keep the winding order.
6084 3097 : if (p0.Equals(poFinalPoint) &&
6085 11 : !EQUAL(poLS->getGeometryName(), "LINEARRING"))
6086 : {
6087 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6088 : printf("Full circle of a non-polygonal zone\n"); /*ok*/
6089 : #endif
6090 1 : poLS->getPoint((i + j + 2) / 4, &p1);
6091 1 : poCS->addPoint(&p1);
6092 1 : if (bValidAlphaRatio)
6093 : {
6094 1 : p1.setX(dfXMid);
6095 1 : p1.setY(dfYMid);
6096 1 : if (poLS->getCoordinateDimension() == 3)
6097 0 : p1.setZ(dfZMid);
6098 : }
6099 : else
6100 : {
6101 0 : poLS->getPoint((i + j + 1) / 2, &p1);
6102 : }
6103 1 : poCS->addPoint(&p1);
6104 1 : poLS->getPoint(3 * (i + j + 2) / 4, &p1);
6105 1 : poCS->addPoint(&p1);
6106 : }
6107 :
6108 3085 : else if (bValidAlphaRatio)
6109 : {
6110 3074 : p1.setX(dfXMid);
6111 3074 : p1.setY(dfYMid);
6112 3074 : if (poLS->getCoordinateDimension() == 3)
6113 2 : p1.setZ(dfZMid);
6114 3074 : poCS->addPoint(&p1);
6115 : }
6116 :
6117 : // If we have found a candidate for a precise intermediate
6118 : // point, use it.
6119 11 : else if (iMidPoint >= 1 && iMidPoint < j)
6120 : {
6121 3 : poLS->getPoint(iMidPoint, &p1);
6122 3 : poCS->addPoint(&p1);
6123 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6124 : printf("Using detected midpoint...\n"); /*ok*/
6125 : printf("x_mid = %f, y_mid = %f\n", p1.getX(), p1.getY()); /*ok*/
6126 : #endif
6127 : }
6128 : // Otherwise pick up the mid point between both extremities.
6129 : else
6130 : {
6131 8 : poLS->getPoint((i + j + 1) / 2, &p1);
6132 8 : poCS->addPoint(&p1);
6133 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6134 : printf("Pickup 'random' midpoint at index=%d...\n", /*ok*/
6135 : (i + j + 1) / 2);
6136 : printf("x_mid = %f, y_mid = %f\n", p1.getX(), p1.getY()); /*ok*/
6137 : #endif
6138 : }
6139 3086 : poCS->addPoint(poFinalPoint);
6140 :
6141 : #ifdef VERBOSE_DEBUG_CURVEFROMLINESTRING
6142 : printf("----------------------------\n"); /*ok*/
6143 : #endif
6144 :
6145 3086 : if (j + 2 >= poLS->getNumPoints())
6146 3048 : return -2;
6147 38 : return j + 1;
6148 : }
6149 :
6150 : /************************************************************************/
6151 : /* curveFromLineString() */
6152 : /************************************************************************/
6153 :
6154 : /**
6155 : * \brief Try to convert a linestring approximating curves into a curve.
6156 : *
6157 : * This method can return a COMPOUNDCURVE, a CIRCULARSTRING or a LINESTRING.
6158 : *
6159 : * This method is the reverse of curveFromLineString().
6160 : *
6161 : * @param poLS handle to the geometry to convert.
6162 : * @param papszOptions options as a null-terminated list of strings.
6163 : * Unused for now. Must be set to NULL.
6164 : *
6165 : * @return the converted geometry (ownership to caller).
6166 : *
6167 : * @since GDAL 2.0
6168 : */
6169 :
6170 3198 : OGRCurve *OGRGeometryFactory::curveFromLineString(
6171 : const OGRLineString *poLS, CPL_UNUSED const char *const *papszOptions)
6172 : {
6173 3198 : OGRCompoundCurve *poCC = nullptr;
6174 3198 : OGRCircularString *poCS = nullptr;
6175 3198 : OGRLineString *poLSNew = nullptr;
6176 3198 : const int nLSNumPoints = poLS->getNumPoints();
6177 3198 : const bool bIsClosed = nLSNumPoints >= 4 && poLS->get_IsClosed();
6178 3626 : for (int i = 0; i < nLSNumPoints; /* nothing */)
6179 : {
6180 3476 : const int iNewI = OGRGF_DetectArc(poLS, i, poCC, poCS, poLSNew);
6181 3476 : if (iNewI == -2)
6182 3048 : break;
6183 428 : if (iNewI >= 0)
6184 : {
6185 38 : i = iNewI;
6186 38 : continue;
6187 : }
6188 :
6189 390 : if (poCS != nullptr)
6190 : {
6191 14 : if (poCC == nullptr)
6192 5 : poCC = new OGRCompoundCurve();
6193 14 : poCC->addCurveDirectly(poCS);
6194 14 : poCS = nullptr;
6195 : }
6196 :
6197 390 : OGRPoint p;
6198 390 : poLS->getPoint(i, &p);
6199 390 : if (poLSNew == nullptr)
6200 : {
6201 160 : poLSNew = new OGRLineString();
6202 160 : poLSNew->addPoint(&p);
6203 : }
6204 : // Not strictly necessary, but helps having 'clean' lines without
6205 : // duplicated points.
6206 : else
6207 : {
6208 230 : double dfScale = std::max(1.0, fabs(p.getX()));
6209 230 : dfScale = std::max(dfScale, fabs(p.getY()));
6210 230 : if (bIsClosed && i == nLSNumPoints - 1)
6211 7 : dfScale = 0;
6212 230 : constexpr double dfToleranceEps =
6213 : OGRCompoundCurve::DEFAULT_TOLERANCE_EPSILON;
6214 230 : if (fabs(poLSNew->getX(poLSNew->getNumPoints() - 1) - p.getX()) >
6215 239 : dfToleranceEps * dfScale ||
6216 9 : fabs(poLSNew->getY(poLSNew->getNumPoints() - 1) - p.getY()) >
6217 9 : dfToleranceEps * dfScale)
6218 : {
6219 229 : poLSNew->addPoint(&p);
6220 : }
6221 : }
6222 :
6223 390 : i++;
6224 : }
6225 :
6226 3198 : OGRCurve *poRet = nullptr;
6227 :
6228 3198 : if (poLSNew != nullptr && poLSNew->getNumPoints() < 2)
6229 : {
6230 1 : delete poLSNew;
6231 1 : poLSNew = nullptr;
6232 1 : if (poCC != nullptr)
6233 : {
6234 1 : if (poCC->getNumCurves() == 1)
6235 : {
6236 1 : poRet = poCC->stealCurve(0);
6237 1 : delete poCC;
6238 1 : poCC = nullptr;
6239 : }
6240 : else
6241 0 : poRet = poCC;
6242 : }
6243 : else
6244 0 : poRet = poLS->clone();
6245 : }
6246 3197 : else if (poCC != nullptr)
6247 : {
6248 7 : if (poLSNew)
6249 6 : poCC->addCurveDirectly(poLSNew);
6250 : else
6251 1 : poCC->addCurveDirectly(poCS);
6252 7 : poRet = poCC;
6253 : }
6254 3190 : else if (poLSNew != nullptr)
6255 142 : poRet = poLSNew;
6256 3048 : else if (poCS != nullptr)
6257 3047 : poRet = poCS;
6258 : else
6259 1 : poRet = poLS->clone();
6260 :
6261 3198 : poRet->assignSpatialReference(poLS->getSpatialReference());
6262 :
6263 3198 : return poRet;
6264 : }
6265 :
6266 : /************************************************************************/
6267 : /* createFromGeoJson( const char* ) */
6268 : /************************************************************************/
6269 :
6270 : /**
6271 : * @brief Create geometry from GeoJson fragment.
6272 : * @param pszJsonString The GeoJSON fragment for the geometry.
6273 : * @param nSize (new in GDAL 3.4) Optional length of the string
6274 : * if it is not null-terminated
6275 : * @return a geometry on success, or NULL on error.
6276 : * @since GDAL 2.3
6277 : */
6278 0 : OGRGeometry *OGRGeometryFactory::createFromGeoJson(const char *pszJsonString,
6279 : int nSize)
6280 : {
6281 0 : CPLJSONDocument oDocument;
6282 0 : if (!oDocument.LoadMemory(reinterpret_cast<const GByte *>(pszJsonString),
6283 : nSize))
6284 : {
6285 0 : return nullptr;
6286 : }
6287 :
6288 0 : return createFromGeoJson(oDocument.GetRoot());
6289 : }
6290 :
6291 : /************************************************************************/
6292 : /* createFromGeoJson( const CPLJSONObject& ) */
6293 : /************************************************************************/
6294 :
6295 : /**
6296 : * @brief Create geometry from GeoJson fragment.
6297 : * @param oJsonObject The JSONObject class describes the GeoJSON geometry.
6298 : * @return a geometry on success, or NULL on error.
6299 : * @since GDAL 2.3
6300 : */
6301 : OGRGeometry *
6302 0 : OGRGeometryFactory::createFromGeoJson(const CPLJSONObject &oJsonObject)
6303 : {
6304 0 : if (!oJsonObject.IsValid())
6305 : {
6306 0 : return nullptr;
6307 : }
6308 :
6309 : // TODO: Move from GeoJSON driver functions create geometry here, and
6310 : // replace json-c specific json_object to CPLJSONObject
6311 0 : return OGRGeoJSONReadGeometry(
6312 0 : static_cast<json_object *>(oJsonObject.GetInternalHandle()));
6313 : }
|
