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