Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: OpenGIS Simple Features Reference Implementation
4 : * Purpose: The OGRSCoordinateTransformation class.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2000, Frank Warmerdam
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "ogr_spatialref.h"
16 :
17 : #include <algorithm>
18 : #include <cassert>
19 : #include <cmath>
20 : #include <cstring>
21 : #include <limits>
22 : #include <list>
23 : #include <mutex>
24 :
25 : #include "cpl_conv.h"
26 : #include "cpl_error.h"
27 : #include "cpl_mem_cache.h"
28 : #include "cpl_string.h"
29 : #include "ogr_core.h"
30 : #include "ogr_srs_api.h"
31 : #include "ogr_proj_p.h"
32 : #include "ogrct_priv.h"
33 :
34 : #include "proj.h"
35 : #include "proj_experimental.h"
36 :
37 : #ifdef DEBUG_PERF
38 : static double g_dfTotalTimeCRStoCRS = 0;
39 : static double g_dfTotalTimeReprojection = 0;
40 :
41 : /************************************************************************/
42 : /* CPLGettimeofday() */
43 : /************************************************************************/
44 :
45 : #if defined(_WIN32) && !defined(__CYGWIN__)
46 : #include <sys/timeb.h>
47 :
48 : namespace
49 : {
50 : struct CPLTimeVal
51 : {
52 : time_t tv_sec; /* seconds */
53 : long tv_usec; /* and microseconds */
54 : };
55 : } // namespace
56 :
57 : static void CPLGettimeofday(struct CPLTimeVal *tp, void * /* timezonep*/)
58 : {
59 : struct _timeb theTime;
60 :
61 : _ftime(&theTime);
62 : tp->tv_sec = static_cast<time_t>(theTime.time);
63 : tp->tv_usec = theTime.millitm * 1000;
64 : }
65 : #else
66 : #include <sys/time.h> /* for gettimeofday() */
67 : #define CPLTimeVal timeval
68 : #define CPLGettimeofday(t, u) gettimeofday(t, u)
69 : #endif
70 :
71 : #endif // DEBUG_PERF
72 :
73 : // Cache of OGRProjCT objects
74 : static std::mutex g_oCTCacheMutex;
75 : class OGRProjCT;
76 : typedef std::string CTCacheKey;
77 : typedef std::unique_ptr<OGRProjCT> CTCacheValue;
78 : static lru11::Cache<CTCacheKey, CTCacheValue> *g_poCTCache = nullptr;
79 :
80 : /************************************************************************/
81 : /* OGRCoordinateTransformationOptions::Private */
82 : /************************************************************************/
83 :
84 5661 : struct OGRCoordinateTransformationOptions::Private
85 : {
86 : bool bHasAreaOfInterest = false;
87 : double dfWestLongitudeDeg = 0.0;
88 : double dfSouthLatitudeDeg = 0.0;
89 : double dfEastLongitudeDeg = 0.0;
90 : double dfNorthLatitudeDeg = 0.0;
91 :
92 : CPLString osCoordOperation{};
93 : bool bReverseCO = false;
94 :
95 : bool bAllowBallpark = true;
96 : double dfAccuracy = -1; // no constraint
97 :
98 : bool bOnlyBest = false;
99 : bool bOnlyBestOptionSet = false;
100 :
101 : bool bHasSourceCenterLong = false;
102 : double dfSourceCenterLong = 0.0;
103 :
104 : bool bHasTargetCenterLong = false;
105 : double dfTargetCenterLong = 0.0;
106 :
107 : bool bCheckWithInvertProj = false;
108 :
109 : Private();
110 2287 : Private(const Private &) = default;
111 : Private(Private &&) = default;
112 : Private &operator=(const Private &) = default;
113 : Private &operator=(Private &&) = default;
114 :
115 : std::string GetKey() const;
116 : void RefreshCheckWithInvertProj();
117 : };
118 :
119 : /************************************************************************/
120 : /* Private() */
121 : /************************************************************************/
122 :
123 11359 : OGRCoordinateTransformationOptions::Private::Private()
124 : {
125 11359 : RefreshCheckWithInvertProj();
126 11359 : }
127 :
128 : /************************************************************************/
129 : /* GetKey() */
130 : /************************************************************************/
131 :
132 10177 : std::string OGRCoordinateTransformationOptions::Private::GetKey() const
133 : {
134 10177 : std::string ret;
135 10177 : ret += std::to_string(static_cast<int>(bHasAreaOfInterest));
136 10177 : ret += std::to_string(dfWestLongitudeDeg);
137 10177 : ret += std::to_string(dfSouthLatitudeDeg);
138 10177 : ret += std::to_string(dfEastLongitudeDeg);
139 10177 : ret += std::to_string(dfNorthLatitudeDeg);
140 10177 : ret += osCoordOperation;
141 10177 : ret += std::to_string(static_cast<int>(bReverseCO));
142 10177 : ret += std::to_string(static_cast<int>(bAllowBallpark));
143 10177 : ret += std::to_string(dfAccuracy);
144 10177 : ret += std::to_string(static_cast<int>(bOnlyBestOptionSet));
145 10177 : ret += std::to_string(static_cast<int>(bOnlyBest));
146 10177 : ret += std::to_string(static_cast<int>(bHasSourceCenterLong));
147 10177 : ret += std::to_string(dfSourceCenterLong);
148 10177 : ret += std::to_string(static_cast<int>(bHasTargetCenterLong));
149 10177 : ret += std::to_string(dfTargetCenterLong);
150 10177 : ret += std::to_string(static_cast<int>(bCheckWithInvertProj));
151 10177 : return ret;
152 : }
153 :
154 : /************************************************************************/
155 : /* RefreshCheckWithInvertProj() */
156 : /************************************************************************/
157 :
158 13612 : void OGRCoordinateTransformationOptions::Private::RefreshCheckWithInvertProj()
159 : {
160 13612 : bCheckWithInvertProj =
161 13612 : CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO"));
162 13612 : }
163 :
164 : /************************************************************************/
165 : /* GetAsAProjRecognizableString() */
166 : /************************************************************************/
167 :
168 2796 : static char *GetAsAProjRecognizableString(const OGRSpatialReference *poSRS)
169 : {
170 2796 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
171 : // If there's a PROJ4 EXTENSION node in WKT1, then use
172 : // it. For example when dealing with "+proj=longlat +lon_wrap=180"
173 2796 : char *pszText = nullptr;
174 2796 : if (poSRS->GetExtension(nullptr, "PROJ4", nullptr))
175 : {
176 17 : poSRS->exportToProj4(&pszText);
177 17 : if (strstr(pszText, " +type=crs") == nullptr)
178 : {
179 34 : auto tmpText = std::string(pszText) + " +type=crs";
180 17 : CPLFree(pszText);
181 17 : pszText = CPLStrdup(tmpText.c_str());
182 : }
183 : }
184 2779 : else if (poSRS->IsEmpty())
185 : {
186 3 : pszText = CPLStrdup("");
187 : }
188 : else
189 : {
190 : // We export to PROJJSON rather than WKT2:2019 because PROJJSON
191 : // is a bit more verbose, which helps in situations like
192 : // https://github.com/OSGeo/gdal/issues/9732 /
193 : // https://github.com/OSGeo/PROJ/pull/4124 where we want to export
194 : // a DerivedProjectedCRS whose base ProjectedCRS has non-metre axis.
195 2776 : poSRS->exportToPROJJSON(&pszText, nullptr);
196 : }
197 :
198 5592 : return pszText;
199 : }
200 :
201 : /************************************************************************/
202 : /* GetTextRepresentation() */
203 : /************************************************************************/
204 :
205 12083 : static char *GetTextRepresentation(const OGRSpatialReference *poSRS)
206 : {
207 9289 : const auto CanUseAuthorityDef = [](const OGRSpatialReference *poSRS1,
208 : OGRSpatialReference *poSRSFromAuth,
209 : const char *pszAuth)
210 : {
211 18502 : if (EQUAL(pszAuth, "EPSG") &&
212 9213 : CPLTestBool(
213 : CPLGetConfigOption("OSR_CT_USE_DEFAULT_EPSG_TOWGS84", "NO")))
214 : {
215 : // We don't want by default to honour 'default' TOWGS84 terms that
216 : // come with the EPSG code because there might be a better
217 : // transformation from that Typical case if EPSG:31468 "DHDN /
218 : // 3-degree Gauss-Kruger zone 4" where the DHDN->TOWGS84
219 : // transformation can use the BETA2007.gsb grid instead of
220 : // TOWGS84[598.1,73.7,418.2,0.202,0.045,-2.455,6.7] But if the user
221 : // really wants it, it can set the OSR_CT_USE_DEFAULT_EPSG_TOWGS84
222 : // configuration option to YES
223 : double adfTOWGS84_1[7];
224 : double adfTOWGS84_2[7];
225 :
226 4 : poSRSFromAuth->AddGuessedTOWGS84();
227 :
228 4 : if (poSRS1->GetTOWGS84(adfTOWGS84_1) == OGRERR_NONE &&
229 6 : poSRSFromAuth->GetTOWGS84(adfTOWGS84_2) == OGRERR_NONE &&
230 2 : memcmp(adfTOWGS84_1, adfTOWGS84_2, sizeof(adfTOWGS84_1)) == 0)
231 : {
232 2 : return false;
233 : }
234 : }
235 9287 : return true;
236 : };
237 :
238 12083 : char *pszText = nullptr;
239 : // If we have a AUTH:CODE attached, use it to retrieve the full
240 : // definition in case a trip to WKT1 has lost the area of use.
241 : // unless OGR_CT_PREFER_OFFICIAL_SRS_DEF=NO (see
242 : // https://github.com/OSGeo/PROJ/issues/2955)
243 12083 : const char *pszAuth = poSRS->GetAuthorityName();
244 12083 : const char *pszCode = poSRS->GetAuthorityCode();
245 21423 : if (pszAuth && pszCode &&
246 9340 : CPLTestBool(
247 : CPLGetConfigOption("OGR_CT_PREFER_OFFICIAL_SRS_DEF", "YES")))
248 : {
249 18680 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
250 :
251 18680 : CPLString osAuthCode(pszAuth);
252 9340 : osAuthCode += ':';
253 9340 : osAuthCode += pszCode;
254 18680 : OGRSpatialReference oTmpSRS;
255 9340 : if (oTmpSRS.SetFromUserInput(osAuthCode) == OGRERR_NONE)
256 : {
257 9327 : const char *pszAuthAfter = oTmpSRS.GetAuthorityName();
258 9327 : const char *pszCodeAfter = oTmpSRS.GetAuthorityCode();
259 9327 : if (pszAuthAfter && pszCodeAfter && EQUAL(pszAuthAfter, pszAuth) &&
260 9326 : EQUAL(pszCodeAfter, pszCode))
261 : {
262 9326 : oTmpSRS.SetDataAxisToSRSAxisMapping(
263 : poSRS->GetDataAxisToSRSAxisMapping());
264 9326 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
265 : nullptr};
266 9326 : if (oTmpSRS.IsSame(poSRS, apszOptionsIsSame))
267 : {
268 9289 : if (CanUseAuthorityDef(poSRS, &oTmpSRS, pszAuth))
269 : {
270 9287 : pszText = CPLStrdup(osAuthCode);
271 : }
272 : }
273 : }
274 : }
275 : }
276 12083 : if (pszText == nullptr)
277 : {
278 2796 : pszText = GetAsAProjRecognizableString(poSRS);
279 : }
280 12083 : return pszText;
281 : }
282 :
283 : /************************************************************************/
284 : /* OGRCoordinateTransformationOptions() */
285 : /************************************************************************/
286 :
287 : /** \brief Constructs a new OGRCoordinateTransformationOptions.
288 : *
289 : * @since GDAL 3.0
290 : */
291 11359 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions()
292 11359 : : d(new Private())
293 : {
294 11359 : }
295 :
296 : /************************************************************************/
297 : /* OGRCoordinateTransformationOptions() */
298 : /************************************************************************/
299 :
300 : /** \brief Copy constructor
301 : *
302 : * @since GDAL 3.1
303 : */
304 2287 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions(
305 2287 : const OGRCoordinateTransformationOptions &other)
306 2287 : : d(new Private(*(other.d)))
307 : {
308 2287 : }
309 :
310 : /************************************************************************/
311 : /* operator =() */
312 : /************************************************************************/
313 :
314 : /** \brief Assignment operator
315 : *
316 : * @since GDAL 3.1
317 : */
318 : OGRCoordinateTransformationOptions &
319 5661 : OGRCoordinateTransformationOptions::operator=(
320 : const OGRCoordinateTransformationOptions &other)
321 : {
322 5661 : if (this != &other)
323 : {
324 5661 : *d = *(other.d);
325 : }
326 5661 : return *this;
327 : }
328 :
329 : /************************************************************************/
330 : /* OGRCoordinateTransformationOptions() */
331 : /************************************************************************/
332 :
333 : /** \brief Destroys a OGRCoordinateTransformationOptions.
334 : *
335 : * @since GDAL 3.0
336 : */
337 13179 : OGRCoordinateTransformationOptions::~OGRCoordinateTransformationOptions()
338 : {
339 13179 : }
340 :
341 : /************************************************************************/
342 : /* OCTNewCoordinateTransformationOptions() */
343 : /************************************************************************/
344 :
345 : /** \brief Create coordinate transformation options.
346 : *
347 : * To be freed with OCTDestroyCoordinateTransformationOptions()
348 : *
349 : * @since GDAL 3.0
350 : */
351 13 : OGRCoordinateTransformationOptionsH OCTNewCoordinateTransformationOptions(void)
352 : {
353 13 : return new OGRCoordinateTransformationOptions();
354 : }
355 :
356 : /************************************************************************/
357 : /* OCTDestroyCoordinateTransformationOptions() */
358 : /************************************************************************/
359 :
360 : /** \brief Destroy coordinate transformation options.
361 : *
362 : * @since GDAL 3.0
363 : */
364 11 : void OCTDestroyCoordinateTransformationOptions(
365 : OGRCoordinateTransformationOptionsH hOptions)
366 : {
367 11 : delete hOptions;
368 11 : }
369 :
370 : /************************************************************************/
371 : /* SetAreaOfInterest() */
372 : /************************************************************************/
373 :
374 : /** \brief Sets an area of interest.
375 : *
376 : * The west longitude is generally lower than the east longitude, except for
377 : * areas of interest that go across the anti-meridian.
378 : *
379 : * @param dfWestLongitudeDeg West longitude (in degree). Must be in [-180,180]
380 : * @param dfSouthLatitudeDeg South latitude (in degree). Must be in [-90,90]
381 : * @param dfEastLongitudeDeg East longitude (in degree). Must be in [-180,180]
382 : * @param dfNorthLatitudeDeg North latitude (in degree). Must be in [-90,90]
383 : * @return true in case of success.
384 : *
385 : * @since GDAL 3.0
386 : */
387 1232 : bool OGRCoordinateTransformationOptions::SetAreaOfInterest(
388 : double dfWestLongitudeDeg, double dfSouthLatitudeDeg,
389 : double dfEastLongitudeDeg, double dfNorthLatitudeDeg)
390 : {
391 1232 : if (std::fabs(dfWestLongitudeDeg) > 180)
392 : {
393 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfWestLongitudeDeg");
394 1 : return false;
395 : }
396 1231 : if (std::fabs(dfSouthLatitudeDeg) > 90)
397 : {
398 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfSouthLatitudeDeg");
399 1 : return false;
400 : }
401 1230 : if (std::fabs(dfEastLongitudeDeg) > 180)
402 : {
403 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfEastLongitudeDeg");
404 1 : return false;
405 : }
406 1229 : if (std::fabs(dfNorthLatitudeDeg) > 90)
407 : {
408 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfNorthLatitudeDeg");
409 1 : return false;
410 : }
411 1228 : if (dfSouthLatitudeDeg > dfNorthLatitudeDeg)
412 : {
413 0 : CPLError(CE_Failure, CPLE_AppDefined,
414 : "dfSouthLatitudeDeg should be lower than dfNorthLatitudeDeg");
415 0 : return false;
416 : }
417 1228 : d->bHasAreaOfInterest = true;
418 1228 : d->dfWestLongitudeDeg = dfWestLongitudeDeg;
419 1228 : d->dfSouthLatitudeDeg = dfSouthLatitudeDeg;
420 1228 : d->dfEastLongitudeDeg = dfEastLongitudeDeg;
421 1228 : d->dfNorthLatitudeDeg = dfNorthLatitudeDeg;
422 1228 : return true;
423 : }
424 :
425 : /************************************************************************/
426 : /* OCTCoordinateTransformationOptionsSetAreaOfInterest() */
427 : /************************************************************************/
428 :
429 : /** \brief Sets an area of interest.
430 : *
431 : * See OGRCoordinateTransformationOptions::SetAreaOfInterest()
432 : *
433 : * @since GDAL 3.0
434 : */
435 5 : int OCTCoordinateTransformationOptionsSetAreaOfInterest(
436 : OGRCoordinateTransformationOptionsH hOptions, double dfWestLongitudeDeg,
437 : double dfSouthLatitudeDeg, double dfEastLongitudeDeg,
438 : double dfNorthLatitudeDeg)
439 : {
440 5 : return hOptions->SetAreaOfInterest(dfWestLongitudeDeg, dfSouthLatitudeDeg,
441 5 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
442 : }
443 :
444 : /************************************************************************/
445 : /* SetCoordinateOperation() */
446 : /************************************************************************/
447 :
448 : /** \brief Sets a coordinate operation.
449 : *
450 : * This is a user override to be used instead of the normally computed pipeline.
451 : *
452 : * The pipeline must take into account the axis order of the source and target
453 : * SRS.
454 : *
455 : * The pipeline may be provided as a PROJ string (single step operation or
456 : * multiple step string starting with +proj=pipeline), a WKT2 string describing
457 : * a CoordinateOperation, or a "urn:ogc:def:coordinateOperation:EPSG::XXXX" URN
458 : *
459 : * @param pszCO PROJ or WKT string describing a coordinate operation
460 : * @param bReverseCO Whether the PROJ or WKT string should be evaluated in the
461 : * reverse path
462 : * @return true in case of success.
463 : *
464 : * @since GDAL 3.0
465 : */
466 22 : bool OGRCoordinateTransformationOptions::SetCoordinateOperation(
467 : const char *pszCO, bool bReverseCO)
468 : {
469 22 : d->osCoordOperation = pszCO ? pszCO : "";
470 22 : d->bReverseCO = bReverseCO;
471 22 : return true;
472 : }
473 :
474 : /************************************************************************/
475 : /* SetSourceCenterLong() */
476 : /************************************************************************/
477 :
478 : /*! @cond Doxygen_Suppress */
479 884 : void OGRCoordinateTransformationOptions::SetSourceCenterLong(
480 : double dfCenterLong)
481 : {
482 884 : d->dfSourceCenterLong = dfCenterLong;
483 884 : d->bHasSourceCenterLong = true;
484 884 : }
485 :
486 : /*! @endcond */
487 :
488 : /************************************************************************/
489 : /* SetTargetCenterLong() */
490 : /************************************************************************/
491 :
492 : /*! @cond Doxygen_Suppress */
493 0 : void OGRCoordinateTransformationOptions::SetTargetCenterLong(
494 : double dfCenterLong)
495 : {
496 0 : d->dfTargetCenterLong = dfCenterLong;
497 0 : d->bHasTargetCenterLong = true;
498 0 : }
499 :
500 : /*! @endcond */
501 :
502 : /************************************************************************/
503 : /* OCTCoordinateTransformationOptionsSetOperation() */
504 : /************************************************************************/
505 :
506 : /** \brief Sets a coordinate operation.
507 : *
508 : * See OGRCoordinateTransformationOptions::SetCoordinateTransformation()
509 : *
510 : * @since GDAL 3.0
511 : */
512 8 : int OCTCoordinateTransformationOptionsSetOperation(
513 : OGRCoordinateTransformationOptionsH hOptions, const char *pszCO,
514 : int bReverseCO)
515 : {
516 : // cppcheck-suppress knownConditionTrueFalse
517 8 : return hOptions->SetCoordinateOperation(pszCO, CPL_TO_BOOL(bReverseCO));
518 : }
519 :
520 : /************************************************************************/
521 : /* SetDesiredAccuracy() */
522 : /************************************************************************/
523 :
524 : /** \brief Sets the desired accuracy for coordinate operations.
525 : *
526 : * Only coordinate operations that offer an accuracy of at least the one
527 : * specified will be considered.
528 : *
529 : * An accuracy of 0 is valid and means a coordinate operation made only of one
530 : * or several conversions (map projections, unit conversion, etc.) Operations
531 : * involving ballpark transformations have an unknown accuracy, and will be
532 : * filtered out by any dfAccuracy >= 0 value.
533 : *
534 : * If this option is specified with PROJ < 8, the OGR_CT_OP_SELECTION
535 : * configuration option will default to BEST_ACCURACY.
536 : *
537 : * @param dfAccuracy accuracy in meters (or a negative value to disable this
538 : * filter)
539 : *
540 : * @since GDAL 3.3
541 : */
542 2 : bool OGRCoordinateTransformationOptions::SetDesiredAccuracy(double dfAccuracy)
543 : {
544 2 : d->dfAccuracy = dfAccuracy;
545 2 : return true;
546 : }
547 :
548 : /************************************************************************/
549 : /* OCTCoordinateTransformationOptionsSetDesiredAccuracy() */
550 : /************************************************************************/
551 :
552 : /** \brief Sets the desired accuracy for coordinate operations.
553 : *
554 : * See OGRCoordinateTransformationOptions::SetDesiredAccuracy()
555 : *
556 : * @since GDAL 3.3
557 : */
558 2 : int OCTCoordinateTransformationOptionsSetDesiredAccuracy(
559 : OGRCoordinateTransformationOptionsH hOptions, double dfAccuracy)
560 : {
561 : // cppcheck-suppress knownConditionTrueFalse
562 2 : return hOptions->SetDesiredAccuracy(dfAccuracy);
563 : }
564 :
565 : /************************************************************************/
566 : /* SetBallparkAllowed() */
567 : /************************************************************************/
568 :
569 : /** \brief Sets whether ballpark transformations are allowed.
570 : *
571 : * By default, PROJ may generate "ballpark transformations" (see
572 : * https://proj.org/glossary.html) when precise datum transformations are
573 : * missing. For high accuracy use cases, such transformations might not be
574 : * allowed.
575 : *
576 : * If this option is specified with PROJ < 8, the OGR_CT_OP_SELECTION
577 : * configuration option will default to BEST_ACCURACY.
578 : *
579 : * @param bAllowBallpark false to disable the user of ballpark transformations
580 : *
581 : * @since GDAL 3.3
582 : */
583 1296 : bool OGRCoordinateTransformationOptions::SetBallparkAllowed(bool bAllowBallpark)
584 : {
585 1296 : d->bAllowBallpark = bAllowBallpark;
586 1296 : return true;
587 : }
588 :
589 : /************************************************************************/
590 : /* OCTCoordinateTransformationOptionsSetBallparkAllowed() */
591 : /************************************************************************/
592 :
593 : /** \brief Sets whether ballpark transformations are allowed.
594 : *
595 : * See OGRCoordinateTransformationOptions::SetDesiredAccuracy()
596 : *
597 : * @since GDAL 3.3 and PROJ 8
598 : */
599 2 : int OCTCoordinateTransformationOptionsSetBallparkAllowed(
600 : OGRCoordinateTransformationOptionsH hOptions, int bAllowBallpark)
601 : {
602 : // cppcheck-suppress knownConditionTrueFalse
603 2 : return hOptions->SetBallparkAllowed(CPL_TO_BOOL(bAllowBallpark));
604 : }
605 :
606 : /************************************************************************/
607 : /* SetOnlyBest() */
608 : /************************************************************************/
609 :
610 : /** \brief Sets whether only the "best" operation should be used.
611 : *
612 : * By default (at least in the PROJ 9.x series), PROJ may use coordinate
613 : * operations that are not the "best" if resources (typically grids) needed
614 : * to use them are missing. It will then fallback to other coordinate operations
615 : * that have a lesser accuracy, for example using Helmert transformations,
616 : * or in the absence of such operations, to ones with potential very rough
617 : * accuracy, using "ballpark" transformations
618 : * (see https://proj.org/glossary.html).
619 : *
620 : * When calling this method with bOnlyBest = true, PROJ will only consider the
621 : * "best" operation, and error out (at Transform() time) if they cannot be
622 : * used.
623 : * This method may be used together with SetBallparkAllowed(false) to
624 : * only allow best operations that have a known accuracy.
625 : *
626 : * Note that this method has no effect on PROJ versions before 9.2.
627 : *
628 : * The default value for this option can be also set with the
629 : * PROJ_ONLY_BEST_DEFAULT environment variable, or with the "only_best_default"
630 : * setting of proj.ini. Calling SetOnlyBest() overrides such default value.
631 : *
632 : * @param bOnlyBest set to true to ask PROJ to use only the best operation(s)
633 : *
634 : * @since GDAL 3.8 and PROJ 9.2
635 : */
636 1 : bool OGRCoordinateTransformationOptions::SetOnlyBest(bool bOnlyBest)
637 : {
638 1 : d->bOnlyBest = bOnlyBest;
639 1 : d->bOnlyBestOptionSet = true;
640 1 : return true;
641 : }
642 :
643 : /************************************************************************/
644 : /* OCTCoordinateTransformationOptionsSetOnlyBest() */
645 : /************************************************************************/
646 :
647 : /** \brief Sets whether only the "best" operation(s) should be used.
648 : *
649 : * See OGRCoordinateTransformationOptions::SetOnlyBest()
650 : *
651 : * @since GDAL 3.8 and PROJ 9.2
652 : */
653 0 : int OCTCoordinateTransformationOptionsSetOnlyBest(
654 : OGRCoordinateTransformationOptionsH hOptions, bool bOnlyBest)
655 : {
656 : // cppcheck-suppress knownConditionTrueFalse
657 0 : return hOptions->SetOnlyBest(bOnlyBest);
658 : }
659 :
660 : /************************************************************************/
661 : /* OGRProjCT */
662 : /************************************************************************/
663 :
664 : //! @cond Doxygen_Suppress
665 : class OGRProjCT final : public OGRCoordinateTransformation
666 : {
667 : friend void
668 : OGRProjCTDifferentOperationsStart(OGRCoordinateTransformation *poCT);
669 : friend void
670 : OGRProjCTDifferentOperationsStop(OGRCoordinateTransformation *poCT);
671 : friend bool
672 : OGRProjCTDifferentOperationsUsed(OGRCoordinateTransformation *poCT);
673 :
674 : class PjPtr
675 : {
676 : PJ *m_pj = nullptr;
677 :
678 10629 : void reset()
679 : {
680 10629 : if (m_pj)
681 : {
682 4981 : proj_assign_context(m_pj, OSRGetProjTLSContext());
683 4981 : proj_destroy(m_pj);
684 : }
685 10629 : }
686 :
687 : public:
688 5661 : PjPtr() : m_pj(nullptr)
689 : {
690 5661 : }
691 :
692 0 : explicit PjPtr(PJ *pjIn) : m_pj(pjIn)
693 : {
694 0 : }
695 :
696 5218 : ~PjPtr()
697 5218 : {
698 5218 : reset();
699 5218 : }
700 :
701 22 : PjPtr(const PjPtr &other)
702 22 : : m_pj((other.m_pj != nullptr)
703 22 : ? (proj_clone(OSRGetProjTLSContext(), other.m_pj))
704 22 : : (nullptr))
705 : {
706 22 : }
707 :
708 : PjPtr(PjPtr &&other) : m_pj(other.m_pj)
709 : {
710 : other.m_pj = nullptr;
711 : }
712 :
713 : PjPtr &operator=(const PjPtr &other)
714 : {
715 : if (this != &other)
716 : {
717 : reset();
718 : m_pj = (other.m_pj != nullptr)
719 : ? (proj_clone(OSRGetProjTLSContext(), other.m_pj))
720 : : (nullptr);
721 : }
722 : return *this;
723 : }
724 :
725 5439 : PjPtr &operator=(PJ *pjIn)
726 : {
727 5439 : if (m_pj != pjIn)
728 : {
729 5411 : reset();
730 5411 : m_pj = pjIn;
731 : }
732 5439 : return *this;
733 : }
734 :
735 11711500 : operator PJ *()
736 : {
737 11711500 : return m_pj;
738 : }
739 :
740 4219 : operator const PJ *() const
741 : {
742 4219 : return m_pj;
743 : }
744 : };
745 :
746 : OGRSpatialReference *poSRSSource = nullptr;
747 : OGRAxisOrientation m_eSourceFirstAxisOrient = OAO_Other;
748 : bool bSourceLatLong = false;
749 : bool bSourceWrap = false;
750 : double dfSourceWrapLong = 0.0;
751 : bool bSourceIsDynamicCRS = false;
752 : double dfSourceCoordinateEpoch = 0.0;
753 : std::string m_osSrcSRS{}; // WKT, PROJ4 or AUTH:CODE
754 :
755 : OGRSpatialReference *poSRSTarget = nullptr;
756 : OGRAxisOrientation m_eTargetFirstAxisOrient = OAO_Other;
757 : bool bTargetLatLong = false;
758 : bool bTargetWrap = false;
759 : double dfTargetWrapLong = 0.0;
760 : bool bTargetIsDynamicCRS = false;
761 : double dfTargetCoordinateEpoch = 0.0;
762 : std::string m_osTargetSRS{}; // WKT, PROJ4 or AUTH:CODE
763 :
764 : bool bWebMercatorToWGS84LongLat = false;
765 :
766 : size_t nErrorCount = 0;
767 :
768 : double dfThreshold = 0.0;
769 :
770 : PjPtr m_pj{};
771 : bool m_bReversePj = false;
772 :
773 : bool m_bEmitErrors = true;
774 :
775 : bool bNoTransform = false;
776 :
777 : enum class Strategy
778 : {
779 : PROJ,
780 : BEST_ACCURACY,
781 : FIRST_MATCHING
782 : };
783 : Strategy m_eStrategy = Strategy::PROJ;
784 :
785 : bool
786 : ListCoordinateOperations(const char *pszSrcSRS, const char *pszTargetSRS,
787 : const OGRCoordinateTransformationOptions &options);
788 :
789 : struct Transformation
790 : {
791 : double minx = 0.0;
792 : double miny = 0.0;
793 : double maxx = 0.0;
794 : double maxy = 0.0;
795 : PjPtr pj{};
796 : CPLString osName{};
797 : CPLString osProjString{};
798 : double accuracy = 0.0;
799 :
800 0 : Transformation(double minxIn, double minyIn, double maxxIn,
801 : double maxyIn, PJ *pjIn, const CPLString &osNameIn,
802 : const CPLString &osProjStringIn, double accuracyIn)
803 0 : : minx(minxIn), miny(minyIn), maxx(maxxIn), maxy(maxyIn), pj(pjIn),
804 : osName(osNameIn), osProjString(osProjStringIn),
805 0 : accuracy(accuracyIn)
806 : {
807 0 : }
808 : };
809 :
810 : std::vector<Transformation> m_oTransformations{};
811 : int m_iCurTransformation = -1;
812 : OGRCoordinateTransformationOptions m_options{};
813 :
814 : bool m_recordDifferentOperationsUsed = false;
815 : std::string m_lastPjUsedPROJString{};
816 : bool m_differentOperationsUsed = false;
817 :
818 : void ComputeThreshold();
819 : void DetectWebMercatorToWGS84();
820 :
821 : OGRProjCT &operator=(const OGRProjCT &) = delete;
822 :
823 : static CTCacheKey
824 : MakeCacheKey(const OGRSpatialReference *poSRS1, const char *pszSrcSRS,
825 : const OGRSpatialReference *poSRS2, const char *pszTargetSRS,
826 : const OGRCoordinateTransformationOptions &options);
827 : bool ContainsNorthPole(const double xmin, const double ymin,
828 : const double xmax, const double ymax,
829 : bool lon_lat_order);
830 : bool ContainsSouthPole(const double xmin, const double ymin,
831 : const double xmax, const double ymax,
832 : bool lon_lat_order);
833 :
834 : public:
835 : OGRProjCT();
836 : ~OGRProjCT() override;
837 : OGRProjCT(const OGRProjCT &other);
838 :
839 : int Initialize(const OGRSpatialReference *poSource, const char *pszSrcSRS,
840 : const OGRSpatialReference *poTarget,
841 : const char *pszTargetSRS,
842 : const OGRCoordinateTransformationOptions &options);
843 :
844 : const OGRSpatialReference *GetSourceCS() const override;
845 : const OGRSpatialReference *GetTargetCS() const override;
846 :
847 : int Transform(size_t nCount, double *x, double *y, double *z, double *t,
848 : int *pabSuccess) override;
849 :
850 : int TransformWithErrorCodes(size_t nCount, double *x, double *y, double *z,
851 : double *t, int *panErrorCodes) override;
852 :
853 : int TransformBounds(const double xmin, const double ymin, const double xmax,
854 : const double ymax, double *out_xmin, double *out_ymin,
855 : double *out_xmax, double *out_ymax,
856 : const int densify_pts) override;
857 :
858 52 : bool GetEmitErrors() const override
859 : {
860 52 : return m_bEmitErrors;
861 : }
862 :
863 4748 : void SetEmitErrors(bool bEmitErrors) override
864 : {
865 4748 : m_bEmitErrors = bEmitErrors;
866 4748 : }
867 :
868 : OGRCoordinateTransformation *Clone() const override;
869 :
870 : OGRCoordinateTransformation *GetInverse() const override;
871 :
872 : static void InsertIntoCache(OGRProjCT *poCT);
873 :
874 : static OGRProjCT *
875 : FindFromCache(const OGRSpatialReference *poSource, const char *pszSrcSRS,
876 : const OGRSpatialReference *poTarget, const char *pszTargetSRS,
877 : const OGRCoordinateTransformationOptions &options);
878 : };
879 :
880 : /************************************************************************/
881 : /* OGRProjCTDifferentOperationsStart() */
882 : /************************************************************************/
883 :
884 657 : void OGRProjCTDifferentOperationsStart(OGRCoordinateTransformation *poCT)
885 : {
886 657 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
887 657 : if (poOGRCT)
888 : {
889 657 : poOGRCT->m_recordDifferentOperationsUsed = true;
890 657 : poOGRCT->m_differentOperationsUsed = false;
891 657 : poOGRCT->m_lastPjUsedPROJString.clear();
892 : }
893 : else
894 : {
895 0 : CPLError(CE_Failure, CPLE_AppDefined,
896 : "OGRProjCTDifferentOperationsStart() called with a non "
897 : "OGRProjCT instance");
898 : }
899 657 : }
900 :
901 : /************************************************************************/
902 : /* OGRProjCTDifferentOperationsStop() */
903 : /************************************************************************/
904 :
905 657 : void OGRProjCTDifferentOperationsStop(OGRCoordinateTransformation *poCT)
906 : {
907 657 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
908 657 : if (poOGRCT)
909 : {
910 657 : poOGRCT->m_recordDifferentOperationsUsed = false;
911 : }
912 : else
913 : {
914 0 : CPLError(CE_Failure, CPLE_AppDefined,
915 : "OGRProjCTDifferentOperationsStop() called with a non "
916 : "OGRProjCT instance");
917 : }
918 657 : }
919 :
920 : /************************************************************************/
921 : /* OGRProjCTDifferentOperationsUsed() */
922 : /************************************************************************/
923 :
924 657 : bool OGRProjCTDifferentOperationsUsed(OGRCoordinateTransformation *poCT)
925 : {
926 657 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
927 657 : if (poOGRCT)
928 : {
929 657 : return poOGRCT->m_differentOperationsUsed;
930 : }
931 : else
932 : {
933 0 : CPLError(CE_Failure, CPLE_AppDefined,
934 : "OGRProjCTDifferentOperationsReset() called with a non "
935 : "OGRProjCT instance");
936 0 : return false;
937 : }
938 : }
939 :
940 : //! @endcond
941 :
942 : /************************************************************************/
943 : /* ~OGRCoordinateTransformation() */
944 : /************************************************************************/
945 :
946 : OGRCoordinateTransformation::~OGRCoordinateTransformation() = default;
947 :
948 : /************************************************************************/
949 : /* OCTDestroyCoordinateTransformation() */
950 : /************************************************************************/
951 :
952 : /**
953 : * \brief OGRCoordinateTransformation destructor.
954 : *
955 : * This function is the same as OGRCoordinateTransformation::DestroyCT()
956 : *
957 : * @param hCT the object to delete
958 : */
959 :
960 : void CPL_STDCALL
961 365 : OCTDestroyCoordinateTransformation(OGRCoordinateTransformationH hCT)
962 :
963 : {
964 365 : OGRCoordinateTransformation::DestroyCT(
965 : OGRCoordinateTransformation::FromHandle(hCT));
966 365 : }
967 :
968 : /************************************************************************/
969 : /* DestroyCT() */
970 : /************************************************************************/
971 :
972 : /**
973 : * \brief OGRCoordinateTransformation destructor.
974 : *
975 : * This function is the same as
976 : * OGRCoordinateTransformation::~OGRCoordinateTransformation()
977 : * and OCTDestroyCoordinateTransformation()
978 : *
979 : * This static method will destroy a OGRCoordinateTransformation. It is
980 : * equivalent to calling delete on the object, but it ensures that the
981 : * deallocation is properly executed within the OGR libraries heap on
982 : * platforms where this can matter (win32).
983 : *
984 : * @param poCT the object to delete
985 : *
986 : */
987 :
988 4531 : void OGRCoordinateTransformation::DestroyCT(OGRCoordinateTransformation *poCT)
989 : {
990 4531 : auto poProjCT = dynamic_cast<OGRProjCT *>(poCT);
991 4531 : if (poProjCT)
992 : {
993 4525 : OGRProjCT::InsertIntoCache(poProjCT);
994 : }
995 : else
996 : {
997 6 : delete poCT;
998 : }
999 4531 : }
1000 :
1001 : /************************************************************************/
1002 : /* OGRCreateCoordinateTransformation() */
1003 : /************************************************************************/
1004 :
1005 : /**
1006 : * Create transformation object.
1007 : *
1008 : * This is the same as the C function OCTNewCoordinateTransformation().
1009 : *
1010 : * Input spatial reference system objects are assigned
1011 : * by copy (calling clone() method) and no ownership transfer occurs.
1012 : *
1013 : * The delete operator, or OCTDestroyCoordinateTransformation() should
1014 : * be used to destroy transformation objects.
1015 : *
1016 : * This will honour the axis order advertized by the source and target SRS,
1017 : * as well as their "data axis to SRS axis mapping".
1018 : * To have a behavior similar to GDAL < 3.0, the
1019 : * OGR_CT_FORCE_TRADITIONAL_GIS_ORDER configuration option can be set to YES.
1020 : *
1021 : * @param poSource source spatial reference system.
1022 : * @param poTarget target spatial reference system.
1023 : * @return NULL on failure or a ready to use transformation object.
1024 : */
1025 :
1026 : OGRCoordinateTransformation *
1027 4274 : OGRCreateCoordinateTransformation(const OGRSpatialReference *poSource,
1028 : const OGRSpatialReference *poTarget)
1029 :
1030 : {
1031 4274 : return OGRCreateCoordinateTransformation(
1032 8548 : poSource, poTarget, OGRCoordinateTransformationOptions());
1033 : }
1034 :
1035 : /**
1036 : * Create transformation object.
1037 : *
1038 : * This is the same as the C function OCTNewCoordinateTransformationEx().
1039 : *
1040 : * Input spatial reference system objects are assigned
1041 : * by copy (calling clone() method) and no ownership transfer occurs.
1042 : *
1043 : * The delete operator, or OCTDestroyCoordinateTransformation() should
1044 : * be used to destroy transformation objects.
1045 : *
1046 : * This will honour the axis order advertized by the source and target SRS,
1047 : * as well as their "data axis to SRS axis mapping".
1048 : * To have a behavior similar to GDAL < 3.0, the
1049 : * OGR_CT_FORCE_TRADITIONAL_GIS_ORDER configuration option can be set to YES.
1050 : *
1051 : * The source SRS and target SRS should generally not be NULL. This is only
1052 : * allowed if a custom coordinate operation is set through the hOptions
1053 : * argument.
1054 : *
1055 : * Starting with GDAL 3.0.3, the OGR_CT_OP_SELECTION configuration option can be
1056 : * set to PROJ (default if PROJ >= 6.3), BEST_ACCURACY or FIRST_MATCHING to
1057 : * decide of the strategy to select the operation to use among candidates, whose
1058 : * area of use is compatible with the points to transform. It is only taken into
1059 : * account if no user defined coordinate transformation pipeline has been
1060 : * specified. <ul> <li>PROJ means the default behavior used by PROJ
1061 : * proj_create_crs_to_crs(). In particular the operation to use among several
1062 : * initial candidates is evaluated for each point to transform.</li>
1063 : * <li>BEST_ACCURACY means the operation whose accuracy is best. It should be
1064 : * close to PROJ behavior, except that the operation to select is decided
1065 : * for the average point of the coordinates passed in a single Transform()
1066 : * call. Note: if the OGRCoordinateTransformationOptions::SetDesiredAccuracy()
1067 : * or OGRCoordinateTransformationOptions::SetBallparkAllowed() methods are
1068 : * called with PROJ < 8, this strategy will be selected instead of PROJ.
1069 : * </li>
1070 : * <li>FIRST_MATCHING is the operation ordered first in the list of candidates:
1071 : * it will not necessarily have the best accuracy, but generally a larger
1072 : * area of use. It is evaluated for the average point of the coordinates passed
1073 : * in a single Transform() call. This was the default behavior for GDAL 3.0.0 to
1074 : * 3.0.2</li>
1075 : * </ul>
1076 : *
1077 : * By default, if the source or target SRS definition refers to an official
1078 : * CRS through a code, GDAL will use the official definition if the official
1079 : * definition and the source/target SRS definition are equivalent. Note that
1080 : * TOWGS84[] clauses are ignored when checking equivalence. Starting with
1081 : * GDAL 3.4.1, if you set the OGR_CT_PREFER_OFFICIAL_SRS_DEF configuration
1082 : * option to NO, the source or target SRS definition will be always used.
1083 : *
1084 : * If options contains a user defined coordinate transformation pipeline, it
1085 : * will be unconditionally used.
1086 : * If options has an area of interest defined, it will be used to research the
1087 : * best fitting coordinate transformation (which will be used for all coordinate
1088 : * transformations, even if they don't fall into the declared area of interest)
1089 : * If no options are set, then a list of candidate coordinate operations will be
1090 : * researched, and at each call to Transform(), the best of those candidate
1091 : * regarding the centroid of the coordinate set will be dynamically selected.
1092 : *
1093 : * @param poSource source spatial reference system.
1094 : * @param poTarget target spatial reference system.
1095 : * @param options Coordinate transformation options.
1096 : * @return NULL on failure or a ready to use transformation object.
1097 : * @since GDAL 3.0
1098 : */
1099 :
1100 6054 : OGRCoordinateTransformation *OGRCreateCoordinateTransformation(
1101 : const OGRSpatialReference *poSource, const OGRSpatialReference *poTarget,
1102 : const OGRCoordinateTransformationOptions &options)
1103 :
1104 : {
1105 6054 : char *pszSrcSRS = poSource ? GetTextRepresentation(poSource) : nullptr;
1106 6054 : char *pszTargetSRS = poTarget ? GetTextRepresentation(poTarget) : nullptr;
1107 : // Try to find if we have a match in the case
1108 6054 : OGRProjCT *poCT = OGRProjCT::FindFromCache(poSource, pszSrcSRS, poTarget,
1109 : pszTargetSRS, options);
1110 6054 : if (poCT == nullptr)
1111 : {
1112 3764 : poCT = new OGRProjCT();
1113 3764 : if (!poCT->Initialize(poSource, pszSrcSRS, poTarget, pszTargetSRS,
1114 : options))
1115 : {
1116 29 : delete poCT;
1117 29 : poCT = nullptr;
1118 : }
1119 : }
1120 6054 : CPLFree(pszSrcSRS);
1121 6054 : CPLFree(pszTargetSRS);
1122 :
1123 6054 : return poCT;
1124 : }
1125 :
1126 : /************************************************************************/
1127 : /* OCTNewCoordinateTransformation() */
1128 : /************************************************************************/
1129 :
1130 : /**
1131 : * Create transformation object.
1132 : *
1133 : * This is the same as the C++ function OGRCreateCoordinateTransformation(const
1134 : * OGRSpatialReference *, const OGRSpatialReference *)
1135 : *
1136 : * Input spatial reference system objects are assigned
1137 : * by copy (calling clone() method) and no ownership transfer occurs.
1138 : *
1139 : * OCTDestroyCoordinateTransformation() should
1140 : * be used to destroy transformation objects.
1141 : *
1142 : * This will honour the axis order advertized by the source and target SRS,
1143 : * as well as their "data axis to SRS axis mapping".
1144 : * To have a behavior similar to GDAL < 3.0, the
1145 : * OGR_CT_FORCE_TRADITIONAL_GIS_ORDER configuration option can be set to YES.
1146 : *
1147 : * @param hSourceSRS source spatial reference system.
1148 : * @param hTargetSRS target spatial reference system.
1149 : * @return NULL on failure or a ready to use transformation object.
1150 : */
1151 :
1152 340 : OGRCoordinateTransformationH CPL_STDCALL OCTNewCoordinateTransformation(
1153 : OGRSpatialReferenceH hSourceSRS, OGRSpatialReferenceH hTargetSRS)
1154 :
1155 : {
1156 : return reinterpret_cast<OGRCoordinateTransformationH>(
1157 340 : OGRCreateCoordinateTransformation(
1158 : reinterpret_cast<OGRSpatialReference *>(hSourceSRS),
1159 340 : reinterpret_cast<OGRSpatialReference *>(hTargetSRS)));
1160 : }
1161 :
1162 : /************************************************************************/
1163 : /* OCTNewCoordinateTransformationEx() */
1164 : /************************************************************************/
1165 :
1166 : /**
1167 : * Create transformation object.
1168 : *
1169 : * This is the same as the C++ function OGRCreateCoordinateTransformation(const
1170 : * OGRSpatialReference *, const OGRSpatialReference *, const
1171 : * OGRCoordinateTransformationOptions& )
1172 : *
1173 : * Input spatial reference system objects are assigned
1174 : * by copy (calling clone() method) and no ownership transfer occurs.
1175 : *
1176 : * OCTDestroyCoordinateTransformation() should
1177 : * be used to destroy transformation objects.
1178 : *
1179 : * The source SRS and target SRS should generally not be NULL. This is only
1180 : * allowed if a custom coordinate operation is set through the hOptions
1181 : * argument.
1182 : *
1183 : * This will honour the axis order advertized by the source and target SRS,
1184 : * as well as their "data axis to SRS axis mapping".
1185 : * To have a behavior similar to GDAL < 3.0, the
1186 : * OGR_CT_FORCE_TRADITIONAL_GIS_ORDER configuration option can be set to YES.
1187 : *
1188 : * If options contains a user defined coordinate transformation pipeline, it
1189 : * will be unconditionally used.
1190 : * If options has an area of interest defined, it will be used to research the
1191 : * best fitting coordinate transformation (which will be used for all coordinate
1192 : * transformations, even if they don't fall into the declared area of interest)
1193 : * If no options are set, then a list of candidate coordinate operations will be
1194 : * researched, and at each call to Transform(), the best of those candidate
1195 : * regarding the centroid of the coordinate set will be dynamically selected.
1196 : *
1197 : * @param hSourceSRS source spatial reference system.
1198 : * @param hTargetSRS target spatial reference system.
1199 : * @param hOptions Coordinate transformation options.
1200 : * @return NULL on failure or a ready to use transformation object.
1201 : * @since GDAL 3.0
1202 : */
1203 :
1204 : OGRCoordinateTransformationH
1205 12 : OCTNewCoordinateTransformationEx(OGRSpatialReferenceH hSourceSRS,
1206 : OGRSpatialReferenceH hTargetSRS,
1207 : OGRCoordinateTransformationOptionsH hOptions)
1208 :
1209 : {
1210 : return reinterpret_cast<OGRCoordinateTransformationH>(
1211 12 : OGRCreateCoordinateTransformation(
1212 : reinterpret_cast<OGRSpatialReference *>(hSourceSRS),
1213 : reinterpret_cast<OGRSpatialReference *>(hTargetSRS),
1214 24 : hOptions ? *hOptions : OGRCoordinateTransformationOptions()));
1215 : }
1216 :
1217 : /************************************************************************/
1218 : /* OCTClone() */
1219 : /************************************************************************/
1220 :
1221 : /**
1222 : * Clone transformation object.
1223 : *
1224 : * This is the same as the C++ function OGRCreateCoordinateTransformation::Clone
1225 : *
1226 : * @return handle to transformation's clone or NULL on error,
1227 : * must be freed with OCTDestroyCoordinateTransformation
1228 : *
1229 : * @since GDAL 3.4
1230 : */
1231 :
1232 1 : OGRCoordinateTransformationH OCTClone(OGRCoordinateTransformationH hTransform)
1233 :
1234 : {
1235 1 : VALIDATE_POINTER1(hTransform, "OCTClone", nullptr);
1236 1 : return OGRCoordinateTransformation::ToHandle(
1237 2 : OGRCoordinateTransformation::FromHandle(hTransform)->Clone());
1238 : }
1239 :
1240 : /************************************************************************/
1241 : /* OCTGetSourceCS() */
1242 : /************************************************************************/
1243 :
1244 : /**
1245 : * Transformation's source coordinate system reference.
1246 : *
1247 : * This is the same as the C++ function
1248 : * OGRCreateCoordinateTransformation::GetSourceCS
1249 : *
1250 : * @return handle to transformation's source coordinate system or NULL if not
1251 : * present.
1252 : *
1253 : * The ownership of the returned SRS belongs to the transformation object, and
1254 : * the returned SRS should not be modified.
1255 : *
1256 : * @since GDAL 3.4
1257 : */
1258 :
1259 2 : OGRSpatialReferenceH OCTGetSourceCS(OGRCoordinateTransformationH hTransform)
1260 :
1261 : {
1262 2 : VALIDATE_POINTER1(hTransform, "OCTGetSourceCS", nullptr);
1263 2 : return OGRSpatialReference::ToHandle(const_cast<OGRSpatialReference *>(
1264 2 : OGRCoordinateTransformation::FromHandle(hTransform)->GetSourceCS()));
1265 : }
1266 :
1267 : /************************************************************************/
1268 : /* OCTGetTargetCS() */
1269 : /************************************************************************/
1270 :
1271 : /**
1272 : * Transformation's target coordinate system reference.
1273 : *
1274 : * This is the same as the C++ function
1275 : * OGRCreateCoordinateTransformation::GetTargetCS
1276 : *
1277 : * @return handle to transformation's target coordinate system or NULL if not
1278 : * present.
1279 : *
1280 : * The ownership of the returned SRS belongs to the transformation object, and
1281 : * the returned SRS should not be modified.
1282 : *
1283 : * @since GDAL 3.4
1284 : */
1285 :
1286 2 : OGRSpatialReferenceH OCTGetTargetCS(OGRCoordinateTransformationH hTransform)
1287 :
1288 : {
1289 2 : VALIDATE_POINTER1(hTransform, "OCTGetTargetCS", nullptr);
1290 2 : return OGRSpatialReference::ToHandle(const_cast<OGRSpatialReference *>(
1291 2 : OGRCoordinateTransformation::FromHandle(hTransform)->GetTargetCS()));
1292 : }
1293 :
1294 : /************************************************************************/
1295 : /* OCTGetInverse() */
1296 : /************************************************************************/
1297 :
1298 : /**
1299 : * Inverse transformation object.
1300 : *
1301 : * This is the same as the C++ function
1302 : * OGRCreateCoordinateTransformation::GetInverse
1303 : *
1304 : * @return handle to inverse transformation or NULL on error,
1305 : * must be freed with OCTDestroyCoordinateTransformation
1306 : *
1307 : * @since GDAL 3.4
1308 : */
1309 :
1310 : OGRCoordinateTransformationH CPL_DLL
1311 3 : OCTGetInverse(OGRCoordinateTransformationH hTransform)
1312 :
1313 : {
1314 3 : VALIDATE_POINTER1(hTransform, "OCTGetInverse", nullptr);
1315 3 : return OGRCoordinateTransformation::ToHandle(
1316 6 : OGRCoordinateTransformation::FromHandle(hTransform)->GetInverse());
1317 : }
1318 :
1319 : /************************************************************************/
1320 : /* OGRProjCT() */
1321 : /************************************************************************/
1322 :
1323 : //! @cond Doxygen_Suppress
1324 5661 : OGRProjCT::OGRProjCT()
1325 : {
1326 5661 : }
1327 :
1328 : /************************************************************************/
1329 : /* OGRProjCT(const OGRProjCT& other) */
1330 : /************************************************************************/
1331 :
1332 22 : OGRProjCT::OGRProjCT(const OGRProjCT &other)
1333 22 : : poSRSSource((other.poSRSSource != nullptr) ? (other.poSRSSource->Clone())
1334 : : (nullptr)),
1335 22 : m_eSourceFirstAxisOrient(other.m_eSourceFirstAxisOrient),
1336 22 : bSourceLatLong(other.bSourceLatLong), bSourceWrap(other.bSourceWrap),
1337 22 : dfSourceWrapLong(other.dfSourceWrapLong),
1338 22 : bSourceIsDynamicCRS(other.bSourceIsDynamicCRS),
1339 22 : dfSourceCoordinateEpoch(other.dfSourceCoordinateEpoch),
1340 22 : m_osSrcSRS(other.m_osSrcSRS),
1341 22 : poSRSTarget((other.poSRSTarget != nullptr) ? (other.poSRSTarget->Clone())
1342 : : (nullptr)),
1343 22 : m_eTargetFirstAxisOrient(other.m_eTargetFirstAxisOrient),
1344 22 : bTargetLatLong(other.bTargetLatLong), bTargetWrap(other.bTargetWrap),
1345 22 : dfTargetWrapLong(other.dfTargetWrapLong),
1346 22 : bTargetIsDynamicCRS(other.bTargetIsDynamicCRS),
1347 22 : dfTargetCoordinateEpoch(other.dfTargetCoordinateEpoch),
1348 22 : m_osTargetSRS(other.m_osTargetSRS),
1349 22 : bWebMercatorToWGS84LongLat(other.bWebMercatorToWGS84LongLat),
1350 22 : nErrorCount(other.nErrorCount), dfThreshold(other.dfThreshold),
1351 44 : m_pj(other.m_pj), m_bReversePj(other.m_bReversePj),
1352 22 : m_bEmitErrors(other.m_bEmitErrors), bNoTransform(other.bNoTransform),
1353 22 : m_eStrategy(other.m_eStrategy),
1354 22 : m_oTransformations(other.m_oTransformations),
1355 22 : m_iCurTransformation(other.m_iCurTransformation),
1356 22 : m_options(other.m_options), m_recordDifferentOperationsUsed(false),
1357 66 : m_lastPjUsedPROJString(std::string()), m_differentOperationsUsed(false)
1358 : {
1359 22 : }
1360 :
1361 : /************************************************************************/
1362 : /* ~OGRProjCT() */
1363 : /************************************************************************/
1364 :
1365 10436 : OGRProjCT::~OGRProjCT()
1366 :
1367 : {
1368 5218 : if (poSRSSource != nullptr)
1369 : {
1370 5207 : poSRSSource->Release();
1371 : }
1372 :
1373 5218 : if (poSRSTarget != nullptr)
1374 : {
1375 5205 : poSRSTarget->Release();
1376 : }
1377 10436 : }
1378 :
1379 : /************************************************************************/
1380 : /* ComputeThreshold() */
1381 : /************************************************************************/
1382 :
1383 5660 : void OGRProjCT::ComputeThreshold()
1384 : {
1385 : // The threshold is experimental. Works well with the cases of ticket #2305.
1386 5660 : if (bSourceLatLong)
1387 : {
1388 : // coverity[tainted_data]
1389 3974 : dfThreshold = CPLAtof(CPLGetConfigOption("THRESHOLD", ".1"));
1390 : }
1391 : else
1392 : {
1393 : // 1 works well for most projections, except for +proj=aeqd that
1394 : // requires a tolerance of 10000.
1395 : // coverity[tainted_data]
1396 1686 : dfThreshold = CPLAtof(CPLGetConfigOption("THRESHOLD", "10000"));
1397 : }
1398 5660 : }
1399 :
1400 : /************************************************************************/
1401 : /* DetectWebMercatorToWGS84() */
1402 : /************************************************************************/
1403 :
1404 5660 : void OGRProjCT::DetectWebMercatorToWGS84()
1405 : {
1406 : // Detect webmercator to WGS84
1407 11284 : if (m_options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1408 18217 : poSRSSource->IsProjected() && poSRSTarget->IsGeographic() &&
1409 1273 : ((m_eTargetFirstAxisOrient == OAO_North &&
1410 1145 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1411 6805 : std::vector<int>{2, 1}) ||
1412 151 : (m_eTargetFirstAxisOrient == OAO_East &&
1413 128 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1414 5788 : std::vector<int>{1, 2})))
1415 : {
1416 : // Examine SRS ID before going to Proj4 string for faster execution
1417 : // This assumes that the SRS definition is "not lying", that is, it
1418 : // is equivalent to the resolution of the official EPSG code.
1419 1242 : const char *pszSourceAuth = poSRSSource->GetAuthorityName();
1420 1242 : const char *pszSourceCode = poSRSSource->GetAuthorityCode();
1421 1242 : const char *pszTargetAuth = poSRSTarget->GetAuthorityName();
1422 1242 : const char *pszTargetCode = poSRSTarget->GetAuthorityCode();
1423 1242 : if (pszSourceAuth && pszSourceCode && pszTargetAuth && pszTargetCode &&
1424 855 : EQUAL(pszSourceAuth, "EPSG") && EQUAL(pszTargetAuth, "EPSG"))
1425 : {
1426 843 : bWebMercatorToWGS84LongLat =
1427 1217 : (EQUAL(pszSourceCode, "3857") ||
1428 374 : EQUAL(pszSourceCode, "3785") || // deprecated
1429 1686 : EQUAL(pszSourceCode, "900913")) && // deprecated
1430 469 : EQUAL(pszTargetCode, "4326");
1431 : }
1432 : else
1433 : {
1434 399 : CPLPushErrorHandler(CPLQuietErrorHandler);
1435 399 : char *pszSrcProj4Defn = nullptr;
1436 399 : poSRSSource->exportToProj4(&pszSrcProj4Defn);
1437 :
1438 399 : char *pszDstProj4Defn = nullptr;
1439 399 : poSRSTarget->exportToProj4(&pszDstProj4Defn);
1440 399 : CPLPopErrorHandler();
1441 :
1442 399 : if (pszSrcProj4Defn && pszDstProj4Defn)
1443 : {
1444 399 : if (pszSrcProj4Defn[0] != '\0' &&
1445 399 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] == ' ')
1446 0 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] = 0;
1447 399 : if (pszDstProj4Defn[0] != '\0' &&
1448 399 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] == ' ')
1449 0 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] = 0;
1450 399 : char *pszNeedle = strstr(pszSrcProj4Defn, " ");
1451 399 : if (pszNeedle)
1452 0 : memmove(pszNeedle, pszNeedle + 1,
1453 0 : strlen(pszNeedle + 1) + 1);
1454 399 : pszNeedle = strstr(pszDstProj4Defn, " ");
1455 399 : if (pszNeedle)
1456 0 : memmove(pszNeedle, pszNeedle + 1,
1457 0 : strlen(pszNeedle + 1) + 1);
1458 :
1459 399 : if ((strstr(pszDstProj4Defn, "+datum=WGS84") != nullptr ||
1460 171 : strstr(pszDstProj4Defn,
1461 : "+ellps=WGS84 +towgs84=0,0,0,0,0,0,0 ") !=
1462 228 : nullptr) &&
1463 228 : strstr(pszSrcProj4Defn, "+nadgrids=@null ") != nullptr &&
1464 76 : strstr(pszSrcProj4Defn, "+towgs84") == nullptr)
1465 : {
1466 : char *pszDst =
1467 76 : strstr(pszDstProj4Defn, "+towgs84=0,0,0,0,0,0,0 ");
1468 76 : if (pszDst != nullptr)
1469 : {
1470 0 : char *pszSrc =
1471 : pszDst + strlen("+towgs84=0,0,0,0,0,0,0 ");
1472 0 : memmove(pszDst, pszSrc, strlen(pszSrc) + 1);
1473 : }
1474 : else
1475 : {
1476 76 : memcpy(strstr(pszDstProj4Defn, "+datum=WGS84"),
1477 : "+ellps", 6);
1478 : }
1479 :
1480 76 : pszDst = strstr(pszSrcProj4Defn, "+nadgrids=@null ");
1481 76 : char *pszSrc = pszDst + strlen("+nadgrids=@null ");
1482 76 : memmove(pszDst, pszSrc, strlen(pszSrc) + 1);
1483 :
1484 76 : pszDst = strstr(pszSrcProj4Defn, "+wktext ");
1485 76 : if (pszDst)
1486 : {
1487 76 : pszSrc = pszDst + strlen("+wktext ");
1488 76 : memmove(pszDst, pszSrc, strlen(pszSrc) + 1);
1489 : }
1490 76 : bWebMercatorToWGS84LongLat =
1491 76 : strcmp(pszDstProj4Defn,
1492 152 : "+proj=longlat +ellps=WGS84 +no_defs") == 0 &&
1493 76 : (strcmp(pszSrcProj4Defn,
1494 : "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 "
1495 : "+lon_0=0.0 "
1496 : "+x_0=0.0 +y_0=0 +k=1.0 +units=m +no_defs") ==
1497 76 : 0 ||
1498 76 : strcmp(pszSrcProj4Defn,
1499 : "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 "
1500 : "+lon_0=0 "
1501 : "+x_0=0 +y_0=0 +k=1 +units=m +no_defs") == 0);
1502 : }
1503 : }
1504 :
1505 399 : CPLFree(pszSrcProj4Defn);
1506 399 : CPLFree(pszDstProj4Defn);
1507 : }
1508 :
1509 1242 : if (bWebMercatorToWGS84LongLat)
1510 : {
1511 545 : CPLDebug("OGRCT", "Using WebMercator to WGS84 optimization");
1512 : }
1513 : }
1514 5660 : }
1515 :
1516 : /************************************************************************/
1517 : /* Initialize() */
1518 : /************************************************************************/
1519 :
1520 3765 : int OGRProjCT::Initialize(const OGRSpatialReference *poSourceIn,
1521 : const char *pszSrcSRS,
1522 : const OGRSpatialReference *poTargetIn,
1523 : const char *pszTargetSRS,
1524 : const OGRCoordinateTransformationOptions &options)
1525 :
1526 : {
1527 3765 : m_options = options;
1528 :
1529 3765 : if (poSourceIn == nullptr || poTargetIn == nullptr)
1530 : {
1531 13 : if (options.d->osCoordOperation.empty())
1532 : {
1533 1 : CPLError(CE_Failure, CPLE_AppDefined,
1534 : "OGRProjCT::Initialize(): if source and/or target CRS "
1535 : "are null, a coordinate operation must be specified");
1536 1 : return FALSE;
1537 : }
1538 : }
1539 :
1540 3764 : if (poSourceIn)
1541 : {
1542 3754 : poSRSSource = poSourceIn->Clone();
1543 3754 : m_osSrcSRS = pszSrcSRS;
1544 : }
1545 3764 : if (poTargetIn)
1546 : {
1547 3753 : poSRSTarget = poTargetIn->Clone();
1548 3753 : m_osTargetSRS = pszTargetSRS;
1549 : }
1550 :
1551 : // To easy quick&dirty compatibility with GDAL < 3.0
1552 3764 : if (CPLTestBool(
1553 : CPLGetConfigOption("OGR_CT_FORCE_TRADITIONAL_GIS_ORDER", "NO")))
1554 : {
1555 0 : if (poSRSSource)
1556 0 : poSRSSource->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1557 0 : if (poSRSTarget)
1558 0 : poSRSTarget->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1559 : }
1560 :
1561 3764 : if (poSRSSource)
1562 : {
1563 3754 : bSourceLatLong = CPL_TO_BOOL(poSRSSource->IsGeographic());
1564 3754 : bSourceIsDynamicCRS = poSRSSource->IsDynamic();
1565 3754 : dfSourceCoordinateEpoch = poSRSSource->GetCoordinateEpoch();
1566 3754 : if (!bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0)
1567 : {
1568 3 : bSourceIsDynamicCRS = poSRSSource->HasPointMotionOperation();
1569 : }
1570 3754 : poSRSSource->GetAxis(nullptr, 0, &m_eSourceFirstAxisOrient);
1571 : }
1572 3764 : if (poSRSTarget)
1573 : {
1574 3753 : bTargetLatLong = CPL_TO_BOOL(poSRSTarget->IsGeographic());
1575 3753 : bTargetIsDynamicCRS = poSRSTarget->IsDynamic();
1576 3753 : dfTargetCoordinateEpoch = poSRSTarget->GetCoordinateEpoch();
1577 3753 : if (!bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
1578 : {
1579 2 : bTargetIsDynamicCRS = poSRSTarget->HasPointMotionOperation();
1580 : }
1581 3753 : poSRSTarget->GetAxis(nullptr, 0, &m_eTargetFirstAxisOrient);
1582 : }
1583 :
1584 : #if PROJ_VERSION_MAJOR < 9 || \
1585 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR < 4)
1586 3764 : if (bSourceIsDynamicCRS && bTargetIsDynamicCRS &&
1587 2363 : dfSourceCoordinateEpoch > 0 && dfTargetCoordinateEpoch > 0 &&
1588 0 : dfSourceCoordinateEpoch != dfTargetCoordinateEpoch)
1589 : {
1590 0 : CPLError(CE_Warning, CPLE_AppDefined,
1591 : "Coordinate transformation between different epochs only"
1592 : "supported since PROJ 9.4");
1593 : }
1594 : #endif
1595 :
1596 : /* -------------------------------------------------------------------- */
1597 : /* Setup source and target translations to radians for lat/long */
1598 : /* systems. */
1599 : /* -------------------------------------------------------------------- */
1600 3764 : bSourceWrap = false;
1601 3764 : dfSourceWrapLong = 0.0;
1602 :
1603 3764 : bTargetWrap = false;
1604 3764 : dfTargetWrapLong = 0.0;
1605 :
1606 : /* -------------------------------------------------------------------- */
1607 : /* Preliminary logic to setup wrapping. */
1608 : /* -------------------------------------------------------------------- */
1609 3764 : if (CPLGetConfigOption("CENTER_LONG", nullptr) != nullptr)
1610 : {
1611 0 : bSourceWrap = true;
1612 0 : bTargetWrap = true;
1613 : // coverity[tainted_data]
1614 0 : dfSourceWrapLong = dfTargetWrapLong =
1615 0 : CPLAtof(CPLGetConfigOption("CENTER_LONG", ""));
1616 0 : CPLDebug("OGRCT", "Wrap at %g.", dfSourceWrapLong);
1617 : }
1618 :
1619 : const char *pszCENTER_LONG;
1620 : {
1621 3764 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1622 3764 : pszCENTER_LONG =
1623 3764 : poSRSSource ? poSRSSource->GetExtension("GEOGCS", "CENTER_LONG")
1624 : : nullptr;
1625 : }
1626 3764 : if (pszCENTER_LONG != nullptr)
1627 : {
1628 0 : dfSourceWrapLong = CPLAtof(pszCENTER_LONG);
1629 0 : bSourceWrap = true;
1630 0 : CPLDebug("OGRCT", "Wrap source at %g.", dfSourceWrapLong);
1631 : }
1632 3764 : else if (bSourceLatLong && options.d->bHasSourceCenterLong)
1633 : {
1634 336 : dfSourceWrapLong = options.d->dfSourceCenterLong;
1635 336 : bSourceWrap = true;
1636 336 : CPLDebug("OGRCT", "Wrap source at %g.", dfSourceWrapLong);
1637 : }
1638 :
1639 : {
1640 3764 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1641 3764 : pszCENTER_LONG =
1642 3764 : poSRSTarget ? poSRSTarget->GetExtension("GEOGCS", "CENTER_LONG")
1643 : : nullptr;
1644 : }
1645 3764 : if (pszCENTER_LONG != nullptr)
1646 : {
1647 0 : dfTargetWrapLong = CPLAtof(pszCENTER_LONG);
1648 0 : bTargetWrap = true;
1649 0 : CPLDebug("OGRCT", "Wrap target at %g.", dfTargetWrapLong);
1650 : }
1651 3764 : else if (bTargetLatLong && options.d->bHasTargetCenterLong)
1652 : {
1653 0 : dfTargetWrapLong = options.d->dfTargetCenterLong;
1654 0 : bTargetWrap = true;
1655 0 : CPLDebug("OGRCT", "Wrap target at %g.", dfTargetWrapLong);
1656 : }
1657 :
1658 3764 : ComputeThreshold();
1659 :
1660 3764 : DetectWebMercatorToWGS84();
1661 :
1662 : const char *pszCTOpSelection =
1663 3764 : CPLGetConfigOption("OGR_CT_OP_SELECTION", nullptr);
1664 3764 : if (pszCTOpSelection)
1665 : {
1666 0 : if (EQUAL(pszCTOpSelection, "PROJ"))
1667 0 : m_eStrategy = Strategy::PROJ;
1668 0 : else if (EQUAL(pszCTOpSelection, "BEST_ACCURACY"))
1669 0 : m_eStrategy = Strategy::BEST_ACCURACY;
1670 0 : else if (EQUAL(pszCTOpSelection, "FIRST_MATCHING"))
1671 0 : m_eStrategy = Strategy::FIRST_MATCHING;
1672 : else
1673 0 : CPLError(CE_Warning, CPLE_NotSupported,
1674 : "OGR_CT_OP_SELECTION=%s not supported", pszCTOpSelection);
1675 : }
1676 : #if PROJ_VERSION_MAJOR < 8
1677 : else
1678 : {
1679 3764 : if (options.d->dfAccuracy >= 0 || !options.d->bAllowBallpark)
1680 : {
1681 6 : m_eStrategy = Strategy::BEST_ACCURACY;
1682 : }
1683 : }
1684 : #endif
1685 3764 : if (m_eStrategy == Strategy::PROJ)
1686 : {
1687 : const char *pszUseApproxTMERC =
1688 3758 : CPLGetConfigOption("OSR_USE_APPROX_TMERC", nullptr);
1689 3758 : if (pszUseApproxTMERC && CPLTestBool(pszUseApproxTMERC))
1690 : {
1691 4 : CPLDebug("OSRCT", "Using OGR_CT_OP_SELECTION=BEST_ACCURACY as "
1692 : "OSR_USE_APPROX_TMERC is set");
1693 4 : m_eStrategy = Strategy::BEST_ACCURACY;
1694 : }
1695 : }
1696 :
1697 3764 : if (!options.d->osCoordOperation.empty())
1698 : {
1699 20 : auto ctx = OSRGetProjTLSContext();
1700 20 : m_pj = proj_create(ctx, options.d->osCoordOperation);
1701 20 : if (!m_pj)
1702 : {
1703 0 : CPLError(CE_Failure, CPLE_NotSupported,
1704 : "Cannot instantiate pipeline %s",
1705 0 : options.d->osCoordOperation.c_str());
1706 0 : return FALSE;
1707 : }
1708 20 : m_bReversePj = options.d->bReverseCO;
1709 : #ifdef DEBUG
1710 20 : auto info = proj_pj_info(m_pj);
1711 20 : CPLDebug("OGRCT", "%s %s(user set)", info.definition,
1712 20 : m_bReversePj ? "(reversed) " : "");
1713 : #endif
1714 : }
1715 3744 : else if (!bWebMercatorToWGS84LongLat && poSRSSource && poSRSTarget)
1716 : {
1717 : #ifdef DEBUG_PERF
1718 : struct CPLTimeVal tvStart;
1719 : CPLGettimeofday(&tvStart, nullptr);
1720 : CPLDebug("OGR_CT", "Before proj_create_crs_to_crs()");
1721 : #endif
1722 : #ifdef DEBUG
1723 3529 : CPLDebug("OGR_CT", "Source CRS: '%s'", pszSrcSRS);
1724 3529 : if (dfSourceCoordinateEpoch > 0)
1725 6 : CPLDebug("OGR_CT", "Source coordinate epoch: %.3f",
1726 : dfSourceCoordinateEpoch);
1727 3529 : CPLDebug("OGR_CT", "Target CRS: '%s'", pszTargetSRS);
1728 3529 : if (dfTargetCoordinateEpoch > 0)
1729 10 : CPLDebug("OGR_CT", "Target coordinate epoch: %.3f",
1730 : dfTargetCoordinateEpoch);
1731 : #endif
1732 :
1733 3529 : if (m_eStrategy == Strategy::PROJ)
1734 : {
1735 3519 : PJ_AREA *area = nullptr;
1736 3519 : if (options.d->bHasAreaOfInterest)
1737 : {
1738 473 : area = proj_area_create();
1739 473 : proj_area_set_bbox(area, options.d->dfWestLongitudeDeg,
1740 473 : options.d->dfSouthLatitudeDeg,
1741 473 : options.d->dfEastLongitudeDeg,
1742 473 : options.d->dfNorthLatitudeDeg);
1743 : }
1744 3519 : auto ctx = OSRGetProjTLSContext();
1745 : #if PROJ_VERSION_MAJOR >= 8
1746 : auto srcCRS = proj_create(ctx, pszSrcSRS);
1747 : auto targetCRS = proj_create(ctx, pszTargetSRS);
1748 : if (srcCRS == nullptr || targetCRS == nullptr)
1749 : {
1750 : proj_destroy(srcCRS);
1751 : proj_destroy(targetCRS);
1752 : if (area)
1753 : proj_area_destroy(area);
1754 : return FALSE;
1755 : }
1756 : CPLStringList aosOptions;
1757 : if (options.d->dfAccuracy >= 0)
1758 : aosOptions.SetNameValue(
1759 : "ACCURACY", CPLSPrintf("%.17g", options.d->dfAccuracy));
1760 : if (!options.d->bAllowBallpark)
1761 : aosOptions.SetNameValue("ALLOW_BALLPARK", "NO");
1762 : #if PROJ_VERSION_MAJOR > 9 || \
1763 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 2)
1764 : if (options.d->bOnlyBestOptionSet)
1765 : {
1766 : aosOptions.SetNameValue("ONLY_BEST",
1767 : options.d->bOnlyBest ? "YES" : "NO");
1768 : }
1769 : #endif
1770 :
1771 : #if PROJ_VERSION_MAJOR > 9 || \
1772 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 4)
1773 : if (bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0 &&
1774 : bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
1775 : {
1776 : auto srcCM = proj_coordinate_metadata_create(
1777 : ctx, srcCRS, dfSourceCoordinateEpoch);
1778 : proj_destroy(srcCRS);
1779 : srcCRS = srcCM;
1780 :
1781 : auto targetCM = proj_coordinate_metadata_create(
1782 : ctx, targetCRS, dfTargetCoordinateEpoch);
1783 : proj_destroy(targetCRS);
1784 : targetCRS = targetCM;
1785 : }
1786 : #endif
1787 :
1788 : m_pj = proj_create_crs_to_crs_from_pj(ctx, srcCRS, targetCRS, area,
1789 : aosOptions.List());
1790 : proj_destroy(srcCRS);
1791 : proj_destroy(targetCRS);
1792 : #else
1793 3519 : m_pj = proj_create_crs_to_crs(ctx, pszSrcSRS, pszTargetSRS, area);
1794 : #endif
1795 3519 : if (area)
1796 473 : proj_area_destroy(area);
1797 3519 : if (m_pj == nullptr)
1798 : {
1799 22 : CPLError(CE_Failure, CPLE_NotSupported,
1800 : "Cannot find coordinate operations from `%s' to `%s'",
1801 : pszSrcSRS, pszTargetSRS);
1802 22 : return FALSE;
1803 : }
1804 : }
1805 10 : else if (!ListCoordinateOperations(pszSrcSRS, pszTargetSRS, options))
1806 : {
1807 6 : CPLError(CE_Failure, CPLE_NotSupported,
1808 : "Cannot find coordinate operations from `%s' to `%s'",
1809 : pszSrcSRS, pszTargetSRS);
1810 6 : return FALSE;
1811 : }
1812 : #ifdef DEBUG_PERF
1813 : struct CPLTimeVal tvEnd;
1814 : CPLGettimeofday(&tvEnd, nullptr);
1815 : const double delay = (tvEnd.tv_sec + tvEnd.tv_usec * 1e-6) -
1816 : (tvStart.tv_sec + tvStart.tv_usec * 1e-6);
1817 : g_dfTotalTimeCRStoCRS += delay;
1818 : CPLDebug("OGR_CT", "After proj_create_crs_to_crs(): %d ms",
1819 : static_cast<int>(delay * 1000));
1820 : #endif
1821 : }
1822 :
1823 7452 : if (options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1824 3716 : (dfSourceCoordinateEpoch == 0 || dfTargetCoordinateEpoch == 0 ||
1825 0 : dfSourceCoordinateEpoch == dfTargetCoordinateEpoch))
1826 : {
1827 : // Determine if we can skip the transformation completely.
1828 3716 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
1829 : nullptr};
1830 3716 : bNoTransform =
1831 7098 : !bSourceWrap && !bTargetWrap &&
1832 3382 : CPL_TO_BOOL(poSRSSource->IsSame(poSRSTarget, apszOptionsIsSame));
1833 : }
1834 :
1835 3736 : return TRUE;
1836 : }
1837 :
1838 : /************************************************************************/
1839 : /* op_to_pj() */
1840 : /************************************************************************/
1841 :
1842 4 : static PJ *op_to_pj(PJ_CONTEXT *ctx, PJ *op,
1843 : CPLString *osOutProjString = nullptr)
1844 : {
1845 : // OSR_USE_ETMERC is here just for legacy
1846 4 : bool bForceApproxTMerc = false;
1847 4 : const char *pszUseETMERC = CPLGetConfigOption("OSR_USE_ETMERC", nullptr);
1848 4 : if (pszUseETMERC && pszUseETMERC[0])
1849 : {
1850 0 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
1851 : "OSR_USE_ETMERC is a legacy configuration option, which "
1852 : "now has only effect when set to NO (YES is the default). "
1853 : "Use OSR_USE_APPROX_TMERC=YES instead");
1854 0 : bForceApproxTMerc = !CPLTestBool(pszUseETMERC);
1855 : }
1856 : else
1857 : {
1858 : const char *pszUseApproxTMERC =
1859 4 : CPLGetConfigOption("OSR_USE_APPROX_TMERC", nullptr);
1860 4 : if (pszUseApproxTMERC && pszUseApproxTMERC[0])
1861 : {
1862 4 : bForceApproxTMerc = CPLTestBool(pszUseApproxTMERC);
1863 : }
1864 : }
1865 4 : const char *options[] = {
1866 4 : bForceApproxTMerc ? "USE_APPROX_TMERC=YES" : nullptr, nullptr};
1867 4 : auto proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, options);
1868 4 : if (!proj_string)
1869 : {
1870 0 : return nullptr;
1871 : }
1872 4 : if (osOutProjString)
1873 0 : *osOutProjString = proj_string;
1874 :
1875 4 : if (proj_string[0] == '\0')
1876 : {
1877 : /* Null transform ? */
1878 0 : return proj_create(ctx, "proj=affine");
1879 : }
1880 : else
1881 : {
1882 4 : return proj_create(ctx, proj_string);
1883 : }
1884 : }
1885 :
1886 : /************************************************************************/
1887 : /* ListCoordinateOperations() */
1888 : /************************************************************************/
1889 :
1890 10 : bool OGRProjCT::ListCoordinateOperations(
1891 : const char *pszSrcSRS, const char *pszTargetSRS,
1892 : const OGRCoordinateTransformationOptions &options)
1893 : {
1894 10 : auto ctx = OSRGetProjTLSContext();
1895 :
1896 10 : auto src = proj_create(ctx, pszSrcSRS);
1897 10 : if (!src)
1898 : {
1899 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot instantiate source_crs");
1900 0 : return false;
1901 : }
1902 :
1903 10 : auto dst = proj_create(ctx, pszTargetSRS);
1904 10 : if (!dst)
1905 : {
1906 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot instantiate target_crs");
1907 0 : proj_destroy(src);
1908 0 : return false;
1909 : }
1910 :
1911 10 : auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
1912 10 : if (!operation_ctx)
1913 : {
1914 0 : proj_destroy(src);
1915 0 : proj_destroy(dst);
1916 0 : return false;
1917 : }
1918 :
1919 10 : proj_operation_factory_context_set_spatial_criterion(
1920 : ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
1921 10 : proj_operation_factory_context_set_grid_availability_use(
1922 : ctx, operation_ctx,
1923 : #if PROJ_VERSION_MAJOR >= 7
1924 : proj_context_is_network_enabled(ctx)
1925 : ? PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE
1926 : :
1927 : #endif
1928 : PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
1929 :
1930 10 : if (options.d->bHasAreaOfInterest)
1931 : {
1932 2 : proj_operation_factory_context_set_area_of_interest(
1933 2 : ctx, operation_ctx, options.d->dfWestLongitudeDeg,
1934 2 : options.d->dfSouthLatitudeDeg, options.d->dfEastLongitudeDeg,
1935 2 : options.d->dfNorthLatitudeDeg);
1936 : }
1937 :
1938 10 : if (options.d->dfAccuracy >= 0)
1939 2 : proj_operation_factory_context_set_desired_accuracy(
1940 2 : ctx, operation_ctx, options.d->dfAccuracy);
1941 10 : if (!options.d->bAllowBallpark)
1942 : {
1943 : #if PROJ_VERSION_MAJOR > 7 || \
1944 : (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 1)
1945 : proj_operation_factory_context_set_allow_ballpark_transformations(
1946 : ctx, operation_ctx, FALSE);
1947 : #else
1948 4 : if (options.d->dfAccuracy < 0)
1949 : {
1950 4 : proj_operation_factory_context_set_desired_accuracy(
1951 : ctx, operation_ctx, HUGE_VAL);
1952 : }
1953 : #endif
1954 : }
1955 :
1956 10 : auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
1957 :
1958 10 : if (!op_list)
1959 : {
1960 0 : proj_operation_factory_context_destroy(operation_ctx);
1961 0 : proj_destroy(src);
1962 0 : proj_destroy(dst);
1963 0 : return false;
1964 : }
1965 :
1966 10 : auto op_count = proj_list_get_count(op_list);
1967 10 : if (op_count == 0)
1968 : {
1969 6 : proj_list_destroy(op_list);
1970 6 : proj_operation_factory_context_destroy(operation_ctx);
1971 6 : proj_destroy(src);
1972 6 : proj_destroy(dst);
1973 6 : CPLDebug("OGRCT", "No operation found matching criteria");
1974 6 : return false;
1975 : }
1976 :
1977 0 : if (op_count == 1 || options.d->bHasAreaOfInterest ||
1978 4 : proj_get_type(src) == PJ_TYPE_GEOCENTRIC_CRS ||
1979 0 : proj_get_type(dst) == PJ_TYPE_GEOCENTRIC_CRS)
1980 : {
1981 4 : auto op = proj_list_get(ctx, op_list, 0);
1982 4 : CPLAssert(op);
1983 4 : m_pj = op_to_pj(ctx, op);
1984 8 : CPLString osName;
1985 4 : auto name = proj_get_name(op);
1986 4 : if (name)
1987 4 : osName = name;
1988 4 : proj_destroy(op);
1989 4 : proj_list_destroy(op_list);
1990 4 : proj_operation_factory_context_destroy(operation_ctx);
1991 4 : proj_destroy(src);
1992 4 : proj_destroy(dst);
1993 4 : if (!m_pj)
1994 0 : return false;
1995 : #ifdef DEBUG
1996 4 : auto info = proj_pj_info(m_pj);
1997 4 : CPLDebug("OGRCT", "%s (%s)", info.definition, osName.c_str());
1998 : #endif
1999 4 : return true;
2000 : }
2001 :
2002 : // Create a geographic 2D long-lat degrees CRS that is related to the
2003 : // source CRS
2004 0 : auto geodetic_crs = proj_crs_get_geodetic_crs(ctx, src);
2005 0 : if (!geodetic_crs)
2006 : {
2007 0 : proj_list_destroy(op_list);
2008 0 : proj_operation_factory_context_destroy(operation_ctx);
2009 0 : proj_destroy(src);
2010 0 : proj_destroy(dst);
2011 0 : CPLDebug("OGRCT", "Cannot find geodetic CRS matching source CRS");
2012 0 : return false;
2013 : }
2014 0 : auto geodetic_crs_type = proj_get_type(geodetic_crs);
2015 0 : if (geodetic_crs_type == PJ_TYPE_GEOCENTRIC_CRS ||
2016 0 : geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
2017 : geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS)
2018 : {
2019 0 : auto datum = proj_crs_get_datum(ctx, geodetic_crs);
2020 : #if PROJ_VERSION_MAJOR > 7 || \
2021 : (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 2)
2022 : if (datum == nullptr)
2023 : {
2024 : datum = proj_crs_get_datum_forced(ctx, geodetic_crs);
2025 : }
2026 : #endif
2027 0 : if (datum)
2028 : {
2029 0 : auto ellps = proj_get_ellipsoid(ctx, datum);
2030 0 : proj_destroy(datum);
2031 0 : double semi_major_metre = 0;
2032 0 : double inv_flattening = 0;
2033 0 : proj_ellipsoid_get_parameters(ctx, ellps, &semi_major_metre,
2034 : nullptr, nullptr, &inv_flattening);
2035 0 : auto cs = proj_create_ellipsoidal_2D_cs(
2036 : ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0);
2037 : // It is critical to set the prime meridian to 0
2038 0 : auto temp = proj_create_geographic_crs(
2039 : ctx, "unnamed crs", "unnamed datum", proj_get_name(ellps),
2040 : semi_major_metre, inv_flattening, "Reference prime meridian", 0,
2041 : nullptr, 0, cs);
2042 0 : proj_destroy(ellps);
2043 0 : proj_destroy(cs);
2044 0 : proj_destroy(geodetic_crs);
2045 0 : geodetic_crs = temp;
2046 0 : geodetic_crs_type = proj_get_type(geodetic_crs);
2047 : }
2048 : }
2049 0 : if (geodetic_crs_type != PJ_TYPE_GEOGRAPHIC_2D_CRS)
2050 : {
2051 : // Shouldn't happen
2052 0 : proj_list_destroy(op_list);
2053 0 : proj_operation_factory_context_destroy(operation_ctx);
2054 0 : proj_destroy(src);
2055 0 : proj_destroy(dst);
2056 0 : proj_destroy(geodetic_crs);
2057 0 : CPLDebug("OGRCT", "Cannot find geographic CRS matching source CRS");
2058 0 : return false;
2059 : }
2060 :
2061 : // Create the transformation from this geographic 2D CRS to the source CRS
2062 : auto op_list_to_geodetic =
2063 0 : proj_create_operations(ctx, geodetic_crs, src, operation_ctx);
2064 0 : proj_destroy(geodetic_crs);
2065 :
2066 0 : if (op_list_to_geodetic == nullptr ||
2067 0 : proj_list_get_count(op_list_to_geodetic) == 0)
2068 : {
2069 0 : CPLDebug(
2070 : "OGRCT",
2071 : "Cannot compute transformation from geographic CRS to source CRS");
2072 0 : proj_list_destroy(op_list);
2073 0 : proj_list_destroy(op_list_to_geodetic);
2074 0 : proj_operation_factory_context_destroy(operation_ctx);
2075 0 : proj_destroy(src);
2076 0 : proj_destroy(dst);
2077 0 : return false;
2078 : }
2079 0 : auto opGeogToSrc = proj_list_get(ctx, op_list_to_geodetic, 0);
2080 0 : CPLAssert(opGeogToSrc);
2081 0 : proj_list_destroy(op_list_to_geodetic);
2082 0 : auto pjGeogToSrc = op_to_pj(ctx, opGeogToSrc);
2083 0 : proj_destroy(opGeogToSrc);
2084 0 : if (!pjGeogToSrc)
2085 : {
2086 0 : proj_list_destroy(op_list);
2087 0 : proj_operation_factory_context_destroy(operation_ctx);
2088 0 : proj_destroy(src);
2089 0 : proj_destroy(dst);
2090 0 : return false;
2091 : }
2092 :
2093 : const auto addTransformation =
2094 0 : [this, &pjGeogToSrc, &ctx](PJ *op, double west_lon, double south_lat,
2095 0 : double east_lon, double north_lat)
2096 : {
2097 0 : double minx = -std::numeric_limits<double>::max();
2098 0 : double miny = -std::numeric_limits<double>::max();
2099 0 : double maxx = std::numeric_limits<double>::max();
2100 0 : double maxy = std::numeric_limits<double>::max();
2101 :
2102 0 : if (!(west_lon == -180.0 && east_lon == 180.0 && south_lat == -90.0 &&
2103 : north_lat == 90.0))
2104 : {
2105 0 : minx = -minx;
2106 0 : miny = -miny;
2107 0 : maxx = -maxx;
2108 0 : maxy = -maxy;
2109 :
2110 : double x[21 * 4], y[21 * 4];
2111 0 : for (int j = 0; j <= 20; j++)
2112 : {
2113 0 : x[j] = west_lon + j * (east_lon - west_lon) / 20;
2114 0 : y[j] = south_lat;
2115 0 : x[21 + j] = west_lon + j * (east_lon - west_lon) / 20;
2116 0 : y[21 + j] = north_lat;
2117 0 : x[21 * 2 + j] = west_lon;
2118 0 : y[21 * 2 + j] = south_lat + j * (north_lat - south_lat) / 20;
2119 0 : x[21 * 3 + j] = east_lon;
2120 0 : y[21 * 3 + j] = south_lat + j * (north_lat - south_lat) / 20;
2121 : }
2122 0 : proj_trans_generic(pjGeogToSrc, PJ_FWD, x, sizeof(double), 21 * 4,
2123 : y, sizeof(double), 21 * 4, nullptr, 0, 0,
2124 : nullptr, 0, 0);
2125 0 : for (int j = 0; j < 21 * 4; j++)
2126 : {
2127 0 : if (x[j] != HUGE_VAL && y[j] != HUGE_VAL)
2128 : {
2129 0 : minx = std::min(minx, x[j]);
2130 0 : miny = std::min(miny, y[j]);
2131 0 : maxx = std::max(maxx, x[j]);
2132 0 : maxy = std::max(maxy, y[j]);
2133 : }
2134 : }
2135 : }
2136 :
2137 0 : if (minx <= maxx)
2138 : {
2139 0 : CPLString osProjString;
2140 0 : const double accuracy = proj_coordoperation_get_accuracy(ctx, op);
2141 0 : auto pj = op_to_pj(ctx, op, &osProjString);
2142 0 : CPLString osName;
2143 0 : auto name = proj_get_name(op);
2144 0 : if (name)
2145 0 : osName = name;
2146 0 : proj_destroy(op);
2147 0 : op = nullptr;
2148 0 : if (pj)
2149 : {
2150 : m_oTransformations.emplace_back(minx, miny, maxx, maxy, pj,
2151 0 : osName, osProjString, accuracy);
2152 : }
2153 : }
2154 0 : return op;
2155 0 : };
2156 :
2157 : // Iterate over source->target candidate transformations and reproject
2158 : // their long-lat bounding box into the source CRS.
2159 0 : bool foundWorldTransformation = false;
2160 0 : for (int i = 0; i < op_count; i++)
2161 : {
2162 0 : auto op = proj_list_get(ctx, op_list, i);
2163 0 : CPLAssert(op);
2164 0 : double west_lon = 0.0;
2165 0 : double south_lat = 0.0;
2166 0 : double east_lon = 0.0;
2167 0 : double north_lat = 0.0;
2168 0 : if (proj_get_area_of_use(ctx, op, &west_lon, &south_lat, &east_lon,
2169 0 : &north_lat, nullptr))
2170 : {
2171 0 : if (west_lon <= east_lon)
2172 : {
2173 0 : if (west_lon == -180 && east_lon == 180 && south_lat == -90 &&
2174 0 : north_lat == 90)
2175 : {
2176 0 : foundWorldTransformation = true;
2177 : }
2178 0 : op = addTransformation(op, west_lon, south_lat, east_lon,
2179 : north_lat);
2180 : }
2181 : else
2182 : {
2183 0 : auto op_clone = proj_clone(ctx, op);
2184 :
2185 0 : op = addTransformation(op, west_lon, south_lat, 180, north_lat);
2186 0 : op_clone = addTransformation(op_clone, -180, south_lat,
2187 : east_lon, north_lat);
2188 0 : proj_destroy(op_clone);
2189 : }
2190 : }
2191 :
2192 0 : proj_destroy(op);
2193 : }
2194 :
2195 0 : proj_list_destroy(op_list);
2196 :
2197 : // Sometimes the user will operate even outside the area of use of the
2198 : // source and target CRS, so if no global transformation has been returned
2199 : // previously, trigger the computation of one.
2200 0 : if (!foundWorldTransformation)
2201 : {
2202 0 : proj_operation_factory_context_set_area_of_interest(ctx, operation_ctx,
2203 : -180, -90, 180, 90);
2204 0 : proj_operation_factory_context_set_spatial_criterion(
2205 : ctx, operation_ctx, PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
2206 0 : op_list = proj_create_operations(ctx, src, dst, operation_ctx);
2207 0 : if (op_list)
2208 : {
2209 0 : op_count = proj_list_get_count(op_list);
2210 0 : for (int i = 0; i < op_count; i++)
2211 : {
2212 0 : auto op = proj_list_get(ctx, op_list, i);
2213 0 : CPLAssert(op);
2214 0 : double west_lon = 0.0;
2215 0 : double south_lat = 0.0;
2216 0 : double east_lon = 0.0;
2217 0 : double north_lat = 0.0;
2218 0 : if (proj_get_area_of_use(ctx, op, &west_lon, &south_lat,
2219 0 : &east_lon, &north_lat, nullptr) &&
2220 0 : west_lon == -180 && east_lon == 180 && south_lat == -90 &&
2221 0 : north_lat == 90)
2222 : {
2223 0 : op = addTransformation(op, west_lon, south_lat, east_lon,
2224 : north_lat);
2225 : }
2226 0 : proj_destroy(op);
2227 : }
2228 : }
2229 0 : proj_list_destroy(op_list);
2230 : }
2231 :
2232 0 : proj_operation_factory_context_destroy(operation_ctx);
2233 0 : proj_destroy(src);
2234 0 : proj_destroy(dst);
2235 0 : proj_destroy(pjGeogToSrc);
2236 0 : return !m_oTransformations.empty();
2237 : }
2238 :
2239 : /************************************************************************/
2240 : /* GetSourceCS() */
2241 : /************************************************************************/
2242 :
2243 5788 : const OGRSpatialReference *OGRProjCT::GetSourceCS() const
2244 :
2245 : {
2246 5788 : return poSRSSource;
2247 : }
2248 :
2249 : /************************************************************************/
2250 : /* GetTargetCS() */
2251 : /************************************************************************/
2252 :
2253 3806 : const OGRSpatialReference *OGRProjCT::GetTargetCS() const
2254 :
2255 : {
2256 3806 : return poSRSTarget;
2257 : }
2258 :
2259 : /************************************************************************/
2260 : /* Transform() */
2261 : /************************************************************************/
2262 :
2263 204231 : int OGRCoordinateTransformation::Transform(size_t nCount, double *x, double *y,
2264 : double *z, int *pabSuccessIn)
2265 :
2266 : {
2267 : int *pabSuccess =
2268 : pabSuccessIn
2269 204231 : ? pabSuccessIn
2270 202440 : : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nCount));
2271 204231 : if (!pabSuccess)
2272 0 : return FALSE;
2273 :
2274 204231 : const int bRet = Transform(nCount, x, y, z, nullptr, pabSuccess);
2275 :
2276 204231 : if (pabSuccess != pabSuccessIn)
2277 202440 : CPLFree(pabSuccess);
2278 :
2279 204231 : return bRet;
2280 : }
2281 :
2282 : /************************************************************************/
2283 : /* TransformWithErrorCodes() */
2284 : /************************************************************************/
2285 :
2286 3976 : int OGRCoordinateTransformation::TransformWithErrorCodes(size_t nCount,
2287 : double *x, double *y,
2288 : double *z, double *t,
2289 : int *panErrorCodes)
2290 :
2291 : {
2292 3976 : if (nCount == 1)
2293 : {
2294 3976 : int nSuccess = 0;
2295 3976 : const int bRet = Transform(nCount, x, y, z, t, &nSuccess);
2296 3976 : if (panErrorCodes)
2297 : {
2298 3976 : panErrorCodes[0] = nSuccess ? 0 : -1;
2299 : }
2300 3976 : return bRet;
2301 : }
2302 :
2303 0 : std::vector<int> abSuccess;
2304 : try
2305 : {
2306 0 : abSuccess.resize(nCount);
2307 : }
2308 0 : catch (const std::bad_alloc &)
2309 : {
2310 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
2311 : "Cannot allocate abSuccess[] temporary array");
2312 0 : return FALSE;
2313 : }
2314 :
2315 0 : const int bRet = Transform(nCount, x, y, z, t, abSuccess.data());
2316 :
2317 0 : if (panErrorCodes)
2318 : {
2319 0 : for (size_t i = 0; i < nCount; i++)
2320 : {
2321 0 : panErrorCodes[i] = abSuccess[i] ? 0 : -1;
2322 : }
2323 : }
2324 :
2325 0 : return bRet;
2326 : }
2327 :
2328 : /************************************************************************/
2329 : /* Transform() */
2330 : /************************************************************************/
2331 :
2332 11695800 : int OGRProjCT::Transform(size_t nCount, double *x, double *y, double *z,
2333 : double *t, int *pabSuccess)
2334 :
2335 : {
2336 11695800 : const int bRet = TransformWithErrorCodes(nCount, x, y, z, t, pabSuccess);
2337 :
2338 11695800 : if (pabSuccess)
2339 : {
2340 36328000 : for (size_t i = 0; i < nCount; i++)
2341 : {
2342 24632200 : pabSuccess[i] = (pabSuccess[i] == 0);
2343 : }
2344 : }
2345 :
2346 11695800 : return bRet;
2347 : }
2348 :
2349 : /************************************************************************/
2350 : /* TransformWithErrorCodes() */
2351 : /************************************************************************/
2352 :
2353 : #ifndef PROJ_ERR_COORD_TRANSFM_INVALID_COORD
2354 : #define PROJ_ERR_COORD_TRANSFM_INVALID_COORD 2049
2355 : #define PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN 2050
2356 : #define PROJ_ERR_COORD_TRANSFM_NO_OPERATION 2051
2357 : #endif
2358 :
2359 11710300 : int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
2360 : double *z, double *t, int *panErrorCodes)
2361 :
2362 : {
2363 11710300 : if (nCount == 0)
2364 0 : return TRUE;
2365 :
2366 : // Prevent any coordinate modification when possible
2367 11710300 : if (bNoTransform)
2368 : {
2369 2396 : if (panErrorCodes)
2370 : {
2371 6548 : for (size_t i = 0; i < nCount; i++)
2372 : {
2373 4922 : panErrorCodes[i] = 0;
2374 : }
2375 : }
2376 2396 : return TRUE;
2377 : }
2378 :
2379 : #ifdef DEBUG_VERBOSE
2380 : bool bDebugCT = CPLTestBool(CPLGetConfigOption("OGR_CT_DEBUG", "NO"));
2381 : if (bDebugCT)
2382 : {
2383 : CPLDebug("OGRCT", "count = %u", static_cast<unsigned>(nCount));
2384 : for (size_t i = 0; i < nCount; ++i)
2385 : {
2386 : CPLDebug("OGRCT", " x[%d] = %.16g y[%d] = %.16g", int(i), x[i],
2387 : int(i), y[i]);
2388 : }
2389 : }
2390 : #endif
2391 : #ifdef DEBUG_PERF
2392 : // CPLDebug("OGR_CT", "Begin TransformWithErrorCodes()");
2393 : struct CPLTimeVal tvStart;
2394 : CPLGettimeofday(&tvStart, nullptr);
2395 : #endif
2396 :
2397 : /* -------------------------------------------------------------------- */
2398 : /* Apply data axis to source CRS mapping. */
2399 : /* -------------------------------------------------------------------- */
2400 :
2401 : // Since we may swap the x and y pointers, but cannot tell the caller about this swap,
2402 : // we save the original pointer. The same axis swap code is executed for poSRSTarget.
2403 : // If this nullifies, we save the swap of both axes
2404 11707900 : const auto xOriginal = x;
2405 :
2406 11707900 : if (poSRSSource)
2407 : {
2408 11707000 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
2409 11707000 : if (mapping.size() >= 2)
2410 : {
2411 11707000 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
2412 : {
2413 4881270 : std::swap(x, y);
2414 : }
2415 11707000 : const bool bNegateX = mapping[0] < 0;
2416 11707000 : if (bNegateX)
2417 : {
2418 2 : for (size_t i = 0; i < nCount; i++)
2419 : {
2420 1 : x[i] = -x[i];
2421 : }
2422 : }
2423 11707000 : const bool bNegateY = mapping[1] < 0;
2424 11707000 : if (bNegateY)
2425 : {
2426 0 : for (size_t i = 0; i < nCount; i++)
2427 : {
2428 0 : y[i] = -y[i];
2429 : }
2430 : }
2431 11707000 : if (z && mapping.size() >= 3 && mapping[2] == -3)
2432 : {
2433 2 : for (size_t i = 0; i < nCount; i++)
2434 : {
2435 1 : z[i] = -z[i];
2436 : }
2437 : }
2438 : }
2439 : }
2440 :
2441 : /* -------------------------------------------------------------------- */
2442 : /* Potentially do longitude wrapping. */
2443 : /* -------------------------------------------------------------------- */
2444 11707900 : if (bSourceLatLong && bSourceWrap)
2445 : {
2446 209970 : if (m_eSourceFirstAxisOrient == OAO_East)
2447 : {
2448 14762 : for (size_t i = 0; i < nCount; i++)
2449 : {
2450 13646 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2451 : {
2452 13646 : if (x[i] < dfSourceWrapLong - 180.0)
2453 0 : x[i] += 360.0;
2454 13646 : else if (x[i] > dfSourceWrapLong + 180)
2455 0 : x[i] -= 360.0;
2456 : }
2457 : }
2458 : }
2459 : else
2460 : {
2461 681573 : for (size_t i = 0; i < nCount; i++)
2462 : {
2463 472719 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2464 : {
2465 471284 : if (y[i] < dfSourceWrapLong - 180.0)
2466 35 : y[i] += 360.0;
2467 471249 : else if (y[i] > dfSourceWrapLong + 180)
2468 34 : y[i] -= 360.0;
2469 : }
2470 : }
2471 : }
2472 : }
2473 :
2474 : /* -------------------------------------------------------------------- */
2475 : /* Optimized transform from WebMercator to WGS84 */
2476 : /* -------------------------------------------------------------------- */
2477 11707900 : bool bTransformDone = false;
2478 11707900 : int bRet = TRUE;
2479 11707900 : if (bWebMercatorToWGS84LongLat)
2480 : {
2481 3207150 : constexpr double REVERSE_SPHERE_RADIUS = 1.0 / 6378137.0;
2482 :
2483 3207150 : if (m_eSourceFirstAxisOrient != OAO_East)
2484 : {
2485 0 : std::swap(x, y);
2486 : }
2487 :
2488 3207150 : double y0 = y[0];
2489 9210620 : for (size_t i = 0; i < nCount; i++)
2490 : {
2491 6003470 : if (!std::isfinite(x[i]))
2492 : {
2493 109 : x[i] = HUGE_VAL;
2494 109 : y[i] = HUGE_VAL;
2495 109 : y0 = HUGE_VAL;
2496 109 : bRet = FALSE;
2497 : }
2498 : else
2499 : {
2500 6003360 : x[i] = x[i] * REVERSE_SPHERE_RADIUS;
2501 6003360 : if (x[i] > M_PI)
2502 : {
2503 732 : if (x[i] < M_PI + 1e-14)
2504 : {
2505 34 : x[i] = M_PI;
2506 : }
2507 698 : else if (m_options.d->bCheckWithInvertProj)
2508 : {
2509 0 : x[i] = HUGE_VAL;
2510 0 : y[i] = HUGE_VAL;
2511 0 : y0 = HUGE_VAL;
2512 0 : continue;
2513 : }
2514 : else
2515 : {
2516 698 : x[i] = std::remainder(x[i], 2 * M_PI);
2517 : }
2518 : }
2519 6002630 : else if (x[i] < -M_PI)
2520 : {
2521 1116 : if (x[i] > -M_PI - 1e-14)
2522 : {
2523 8 : x[i] = -M_PI;
2524 : }
2525 1108 : else if (m_options.d->bCheckWithInvertProj)
2526 : {
2527 2 : x[i] = HUGE_VAL;
2528 2 : y[i] = HUGE_VAL;
2529 2 : y0 = HUGE_VAL;
2530 2 : continue;
2531 : }
2532 : else
2533 : {
2534 1106 : x[i] = std::remainder(x[i], 2 * M_PI);
2535 : }
2536 : }
2537 6003360 : constexpr double RAD_TO_DEG = 180. / M_PI;
2538 6003360 : x[i] *= RAD_TO_DEG;
2539 :
2540 : // Optimization for the case where we are provided a whole line
2541 : // of same northing.
2542 6003360 : if (i > 0 && y[i] == y0)
2543 1261920 : y[i] = y[0];
2544 : else
2545 : {
2546 4741440 : y[i] = M_PI / 2.0 -
2547 4741440 : 2.0 * atan(exp(-y[i] * REVERSE_SPHERE_RADIUS));
2548 4741440 : y[i] *= RAD_TO_DEG;
2549 : }
2550 : }
2551 : }
2552 :
2553 3207150 : if (panErrorCodes)
2554 : {
2555 9210600 : for (size_t i = 0; i < nCount; i++)
2556 : {
2557 6003460 : if (x[i] != HUGE_VAL)
2558 6003350 : panErrorCodes[i] = 0;
2559 : else
2560 108 : panErrorCodes[i] =
2561 : PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2562 : }
2563 : }
2564 :
2565 3207150 : if (m_eTargetFirstAxisOrient != OAO_East)
2566 : {
2567 3203240 : std::swap(x, y);
2568 : }
2569 :
2570 3207150 : bTransformDone = true;
2571 : }
2572 :
2573 : // Determine the default coordinate epoch, if not provided in the point to
2574 : // transform.
2575 : // For time-dependent transformations, PROJ can currently only do
2576 : // staticCRS -> dynamicCRS or dynamicCRS -> staticCRS transformations, and
2577 : // in either case, the coordinate epoch of the dynamicCRS must be provided
2578 : // as the input time.
2579 11707900 : double dfDefaultTime = HUGE_VAL;
2580 8367570 : if (bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0 &&
2581 20075500 : !bTargetIsDynamicCRS &&
2582 4 : CPLTestBool(
2583 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2584 : {
2585 4 : dfDefaultTime = dfSourceCoordinateEpoch;
2586 4 : CPLDebug("OGR_CT", "Using coordinate epoch %f from source CRS",
2587 : dfDefaultTime);
2588 : }
2589 11400500 : else if (bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0 &&
2590 23108400 : !bSourceIsDynamicCRS &&
2591 2 : CPLTestBool(
2592 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2593 : {
2594 2 : dfDefaultTime = dfTargetCoordinateEpoch;
2595 2 : CPLDebug("OGR_CT", "Using coordinate epoch %f from target CRS",
2596 : dfDefaultTime);
2597 : }
2598 :
2599 : /* -------------------------------------------------------------------- */
2600 : /* Select dynamically the best transformation for the data, if */
2601 : /* needed. */
2602 : /* -------------------------------------------------------------------- */
2603 11707900 : auto ctx = OSRGetProjTLSContext();
2604 :
2605 11707900 : PJ *pj = m_pj;
2606 11707900 : if (!bTransformDone && !pj)
2607 : {
2608 0 : double avgX = 0.0;
2609 0 : double avgY = 0.0;
2610 0 : size_t nCountValid = 0;
2611 0 : for (size_t i = 0; i < nCount; i++)
2612 : {
2613 0 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2614 : {
2615 0 : avgX += x[i];
2616 0 : avgY += y[i];
2617 0 : nCountValid++;
2618 : }
2619 : }
2620 0 : if (nCountValid != 0)
2621 : {
2622 0 : avgX /= static_cast<double>(nCountValid);
2623 0 : avgY /= static_cast<double>(nCountValid);
2624 : }
2625 :
2626 0 : constexpr int N_MAX_RETRY = 2;
2627 0 : int iExcluded[N_MAX_RETRY] = {-1, -1};
2628 :
2629 0 : const int nOperations = static_cast<int>(m_oTransformations.size());
2630 : PJ_COORD coord;
2631 0 : coord.xyzt.x = avgX;
2632 0 : coord.xyzt.y = avgY;
2633 0 : coord.xyzt.z = z ? z[0] : 0;
2634 0 : coord.xyzt.t = t ? t[0] : dfDefaultTime;
2635 :
2636 : // We may need several attempts. For example the point at
2637 : // lon=-111.5 lat=45.26 falls into the bounding box of the Canadian
2638 : // ntv2_0.gsb grid, except that it is not in any of the subgrids, being
2639 : // in the US. We thus need another retry that will select the conus
2640 : // grid.
2641 0 : for (int iRetry = 0; iRetry <= N_MAX_RETRY; iRetry++)
2642 : {
2643 0 : int iBestTransf = -1;
2644 : // Select transform whose BBOX match our data and has the best
2645 : // accuracy if m_eStrategy == BEST_ACCURACY. Or just the first BBOX
2646 : // matching one, if
2647 : // m_eStrategy == FIRST_MATCHING
2648 0 : double dfBestAccuracy = std::numeric_limits<double>::infinity();
2649 0 : for (int i = 0; i < nOperations; i++)
2650 : {
2651 0 : if (i == iExcluded[0] || i == iExcluded[1])
2652 : {
2653 0 : continue;
2654 : }
2655 0 : const auto &transf = m_oTransformations[i];
2656 0 : if (avgX >= transf.minx && avgX <= transf.maxx &&
2657 0 : avgY >= transf.miny && avgY <= transf.maxy &&
2658 0 : (iBestTransf < 0 || (transf.accuracy >= 0 &&
2659 0 : transf.accuracy < dfBestAccuracy)))
2660 : {
2661 0 : iBestTransf = i;
2662 0 : dfBestAccuracy = transf.accuracy;
2663 0 : if (m_eStrategy == Strategy::FIRST_MATCHING)
2664 0 : break;
2665 : }
2666 : }
2667 0 : if (iBestTransf < 0)
2668 : {
2669 0 : break;
2670 : }
2671 0 : auto &transf = m_oTransformations[iBestTransf];
2672 0 : pj = transf.pj;
2673 0 : proj_assign_context(pj, ctx);
2674 0 : if (iBestTransf != m_iCurTransformation)
2675 : {
2676 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2677 : transf.osProjString.c_str(), transf.osName.c_str());
2678 0 : m_iCurTransformation = iBestTransf;
2679 : }
2680 :
2681 0 : auto res = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2682 0 : if (res.xyzt.x != HUGE_VAL)
2683 : {
2684 0 : break;
2685 : }
2686 0 : pj = nullptr;
2687 0 : CPLDebug("OGRCT", "Did not result in valid result. "
2688 : "Attempting a retry with another operation.");
2689 0 : if (iRetry == N_MAX_RETRY)
2690 : {
2691 0 : break;
2692 : }
2693 0 : iExcluded[iRetry] = iBestTransf;
2694 : }
2695 :
2696 0 : if (!pj)
2697 : {
2698 : // In case we did not find an operation whose area of use is
2699 : // compatible with the input coordinate, then goes through again the
2700 : // list, and use the first operation that does not require grids.
2701 0 : for (int i = 0; i < nOperations; i++)
2702 : {
2703 0 : auto &transf = m_oTransformations[i];
2704 0 : if (proj_coordoperation_get_grid_used_count(ctx, transf.pj) ==
2705 : 0)
2706 : {
2707 0 : pj = transf.pj;
2708 0 : proj_assign_context(pj, ctx);
2709 0 : if (i != m_iCurTransformation)
2710 : {
2711 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2712 : transf.osProjString.c_str(),
2713 : transf.osName.c_str());
2714 0 : m_iCurTransformation = i;
2715 : }
2716 0 : break;
2717 : }
2718 : }
2719 : }
2720 :
2721 0 : if (!pj)
2722 : {
2723 0 : if (m_bEmitErrors && ++nErrorCount < 20)
2724 : {
2725 0 : CPLError(CE_Failure, CPLE_AppDefined,
2726 : "Cannot find transformation for provided coordinates");
2727 : }
2728 0 : else if (nErrorCount == 20)
2729 : {
2730 0 : CPLError(CE_Failure, CPLE_AppDefined,
2731 : "Reprojection failed, further errors will be "
2732 : "suppressed on the transform object.");
2733 : }
2734 :
2735 0 : for (size_t i = 0; i < nCount; i++)
2736 : {
2737 0 : x[i] = HUGE_VAL;
2738 0 : y[i] = HUGE_VAL;
2739 0 : if (panErrorCodes)
2740 0 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_NO_OPERATION;
2741 : }
2742 0 : return FALSE;
2743 : }
2744 : }
2745 11707900 : if (pj)
2746 : {
2747 11706900 : proj_assign_context(pj, ctx);
2748 : }
2749 :
2750 : /* -------------------------------------------------------------------- */
2751 : /* Do the transformation (or not...) using PROJ */
2752 : /* -------------------------------------------------------------------- */
2753 :
2754 11707900 : if (!bTransformDone)
2755 : {
2756 8500740 : const auto nLastErrorCounter = CPLGetErrorCounter();
2757 :
2758 27178100 : for (size_t i = 0; i < nCount; i++)
2759 : {
2760 : PJ_COORD coord;
2761 18677400 : const double xIn = x[i];
2762 18677400 : const double yIn = y[i];
2763 18677400 : if (!std::isfinite(xIn))
2764 : {
2765 10979 : bRet = FALSE;
2766 10979 : x[i] = HUGE_VAL;
2767 10979 : y[i] = HUGE_VAL;
2768 10979 : if (panErrorCodes)
2769 10979 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_INVALID_COORD;
2770 10979 : continue;
2771 : }
2772 18666400 : coord.xyzt.x = x[i];
2773 18666400 : coord.xyzt.y = y[i];
2774 18666400 : coord.xyzt.z = z ? z[i] : 0;
2775 18666400 : coord.xyzt.t = t ? t[i] : dfDefaultTime;
2776 18666400 : proj_errno_reset(pj);
2777 18666400 : coord = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2778 : #if 0
2779 : CPLDebug("OGRCT",
2780 : "Transforming (x=%f,y=%f,z=%f,time=%f) to "
2781 : "(x=%f,y=%f,z=%f,time=%f)",
2782 : x[i], y[i], z ? z[i] : 0, t ? t[i] : dfDefaultTime,
2783 : coord.xyzt.x, coord.xyzt.y, coord.xyzt.z, coord.xyzt.t);
2784 : #endif
2785 18666400 : x[i] = coord.xyzt.x;
2786 18666400 : y[i] = coord.xyzt.y;
2787 18666400 : if (z)
2788 18613800 : z[i] = coord.xyzt.z;
2789 18666400 : if (t)
2790 25 : t[i] = coord.xyzt.t;
2791 18666400 : int err = 0;
2792 18666400 : if (std::isnan(coord.xyzt.x))
2793 : {
2794 136 : bRet = FALSE;
2795 : // This shouldn't normally happen if PROJ projections behave
2796 : // correctly, but e.g inverse laea before PROJ 8.1.1 could
2797 : // do that for points out of domain.
2798 : // See https://github.com/OSGeo/PROJ/pull/2800
2799 136 : x[i] = HUGE_VAL;
2800 136 : y[i] = HUGE_VAL;
2801 136 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2802 :
2803 : #ifdef DEBUG
2804 136 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
2805 : "PROJ returned a NaN value. It should be fixed");
2806 : #else
2807 : CPLDebugOnce("OGR_CT",
2808 : "PROJ returned a NaN value. It should be fixed");
2809 : #endif
2810 : }
2811 18666200 : else if (coord.xyzt.x == HUGE_VAL)
2812 : {
2813 3389550 : bRet = FALSE;
2814 3389550 : err = proj_errno(pj);
2815 : // PROJ should normally emit an error, but in case it does not
2816 : // (e.g PROJ 6.3 with the +ortho projection), synthesize one
2817 3389550 : if (err == 0)
2818 3363810 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2819 : }
2820 : else
2821 : {
2822 15276700 : if (m_recordDifferentOperationsUsed &&
2823 72239 : !m_differentOperationsUsed)
2824 : {
2825 : #if PROJ_VERSION_MAJOR > 9 || \
2826 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 1)
2827 :
2828 : PJ *lastOp = proj_trans_get_last_used_operation(pj);
2829 : if (lastOp)
2830 : {
2831 : const char *projString = proj_as_proj_string(
2832 : ctx, lastOp, PJ_PROJ_5, nullptr);
2833 : if (projString)
2834 : {
2835 : if (m_lastPjUsedPROJString.empty())
2836 : {
2837 : m_lastPjUsedPROJString = projString;
2838 : }
2839 : else if (m_lastPjUsedPROJString != projString)
2840 : {
2841 : m_differentOperationsUsed = true;
2842 : }
2843 : }
2844 : proj_destroy(lastOp);
2845 : }
2846 : #endif
2847 : }
2848 :
2849 15276700 : if (m_options.d->bCheckWithInvertProj)
2850 : {
2851 : // For some projections, we cannot detect if we are trying to
2852 : // reproject coordinates outside the validity area of the
2853 : // projection. So let's do the reverse reprojection and compare
2854 : // with the source coordinates.
2855 : coord =
2856 1852040 : proj_trans(pj, m_bReversePj ? PJ_FWD : PJ_INV, coord);
2857 1852040 : if (fabs(coord.xyzt.x - xIn) > dfThreshold ||
2858 1699480 : fabs(coord.xyzt.y - yIn) > dfThreshold)
2859 : {
2860 155435 : bRet = FALSE;
2861 155435 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2862 155435 : x[i] = HUGE_VAL;
2863 155435 : y[i] = HUGE_VAL;
2864 : }
2865 : }
2866 : }
2867 :
2868 18666400 : if (panErrorCodes)
2869 18664300 : panErrorCodes[i] = err;
2870 :
2871 : /* --------------------------------------------------------------------
2872 : */
2873 : /* Try to report an error through CPL. Get proj error string
2874 : */
2875 : /* if possible. Try to avoid reporting thousands of errors. */
2876 : /* Suppress further error reporting on this OGRProjCT if we */
2877 : /* have already reported 20 errors. */
2878 : /* --------------------------------------------------------------------
2879 : */
2880 18666400 : if (err != 0)
2881 : {
2882 3545120 : if (++nErrorCount < 20)
2883 : {
2884 : #if PROJ_VERSION_MAJOR >= 8
2885 : const char *pszError = proj_context_errno_string(ctx, err);
2886 : #else
2887 4981 : const char *pszError = proj_errno_string(err);
2888 : #endif
2889 4981 : if (m_bEmitErrors
2890 : #ifdef PROJ_ERR_OTHER_NO_INVERSE_OP
2891 : || (i == 0 && err == PROJ_ERR_OTHER_NO_INVERSE_OP)
2892 : #endif
2893 : )
2894 : {
2895 718 : if (nLastErrorCounter != CPLGetErrorCounter() &&
2896 1002 : CPLGetLastErrorType() == CE_Failure &&
2897 284 : strstr(CPLGetLastErrorMsg(), "PROJ:"))
2898 : {
2899 : // do nothing
2900 : }
2901 718 : else if (pszError == nullptr)
2902 0 : CPLError(CE_Failure, CPLE_AppDefined,
2903 : "Reprojection failed, err = %d", err);
2904 : else
2905 718 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
2906 : pszError);
2907 : }
2908 : else
2909 : {
2910 4263 : if (pszError == nullptr)
2911 0 : CPLDebug("OGRCT", "Reprojection failed, err = %d",
2912 : err);
2913 : else
2914 4263 : CPLDebug("OGRCT", "%s for xIn=%.10g yIn=%.10g",
2915 : pszError, xIn, yIn);
2916 : }
2917 : }
2918 3540140 : else if (nErrorCount == 20)
2919 : {
2920 222 : if (m_bEmitErrors)
2921 : {
2922 15 : CPLError(CE_Failure, CPLE_AppDefined,
2923 : "Reprojection failed, err = %d, further "
2924 : "errors will be "
2925 : "suppressed on the transform object.",
2926 : err);
2927 : }
2928 : else
2929 : {
2930 207 : CPLDebug("OGRCT",
2931 : "Reprojection failed, err = %d, further "
2932 : "errors will be "
2933 : "suppressed on the transform object.",
2934 : err);
2935 : }
2936 : }
2937 : }
2938 : }
2939 : }
2940 :
2941 : /* -------------------------------------------------------------------- */
2942 : /* Potentially do longitude wrapping. */
2943 : /* -------------------------------------------------------------------- */
2944 11707900 : if (bTargetLatLong && bTargetWrap)
2945 : {
2946 4500330 : if (m_eTargetFirstAxisOrient == OAO_East)
2947 : {
2948 34837 : for (size_t i = 0; i < nCount; i++)
2949 : {
2950 28825 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2951 : {
2952 28825 : if (x[i] < dfTargetWrapLong - 180.0)
2953 0 : x[i] += 360.0;
2954 28825 : else if (x[i] > dfTargetWrapLong + 180)
2955 174 : x[i] -= 360.0;
2956 : }
2957 : }
2958 : }
2959 : else
2960 : {
2961 19197600 : for (size_t i = 0; i < nCount; i++)
2962 : {
2963 14703300 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2964 : {
2965 11256000 : if (y[i] < dfTargetWrapLong - 180.0)
2966 6662 : y[i] += 360.0;
2967 11249300 : else if (y[i] > dfTargetWrapLong + 180)
2968 43474 : y[i] -= 360.0;
2969 : }
2970 : }
2971 : }
2972 : }
2973 :
2974 : /* -------------------------------------------------------------------- */
2975 : /* Apply data axis to target CRS mapping. */
2976 : /* -------------------------------------------------------------------- */
2977 11707900 : if (poSRSTarget)
2978 : {
2979 11707000 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
2980 11707000 : if (mapping.size() >= 2)
2981 : {
2982 11707000 : const bool bNegateX = mapping[0] < 0;
2983 11707000 : if (bNegateX)
2984 : {
2985 2 : for (size_t i = 0; i < nCount; i++)
2986 : {
2987 1 : x[i] = -x[i];
2988 : }
2989 : }
2990 11707000 : const bool bNegateY = mapping[1] < 0;
2991 11707000 : if (bNegateY)
2992 : {
2993 0 : for (size_t i = 0; i < nCount; i++)
2994 : {
2995 0 : y[i] = -y[i];
2996 : }
2997 : }
2998 11707000 : if (z && mapping.size() >= 3 && mapping[2] == -3)
2999 : {
3000 2 : for (size_t i = 0; i < nCount; i++)
3001 : {
3002 1 : z[i] = -z[i];
3003 : }
3004 : }
3005 :
3006 11707000 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
3007 : {
3008 10457900 : std::swap(x, y);
3009 : }
3010 : }
3011 : }
3012 :
3013 : // Check whether final "genuine" axis swap is really necessary
3014 11707900 : if (x != xOriginal)
3015 : {
3016 4616930 : std::swap_ranges(x, x + nCount, y);
3017 : }
3018 :
3019 : #ifdef DEBUG_VERBOSE
3020 : if (bDebugCT)
3021 : {
3022 : CPLDebug("OGRCT", "Out:");
3023 : for (size_t i = 0; i < nCount; ++i)
3024 : {
3025 : CPLDebug("OGRCT", " x[%d] = %.16g y[%d] = %.16g", int(i), x[i],
3026 : int(i), y[i]);
3027 : }
3028 : }
3029 : #endif
3030 : #ifdef DEBUG_PERF
3031 : struct CPLTimeVal tvEnd;
3032 : CPLGettimeofday(&tvEnd, nullptr);
3033 : const double delay = (tvEnd.tv_sec + tvEnd.tv_usec * 1e-6) -
3034 : (tvStart.tv_sec + tvStart.tv_usec * 1e-6);
3035 : g_dfTotalTimeReprojection += delay;
3036 : // CPLDebug("OGR_CT", "End TransformWithErrorCodes(): %d ms",
3037 : // static_cast<int>(delay * 1000));
3038 : #endif
3039 :
3040 11707900 : return bRet;
3041 : }
3042 :
3043 : /************************************************************************/
3044 : /* TransformBounds() */
3045 : /************************************************************************/
3046 :
3047 : // ---------------------------------------------------------------------------
3048 1229 : static double simple_min(const double *data, const int *panErrorCodes,
3049 : const int arr_len)
3050 : {
3051 1229 : double min_value = HUGE_VAL;
3052 54101 : for (int iii = 0; iii < arr_len; iii++)
3053 : {
3054 52872 : if ((data[iii] < min_value || min_value == HUGE_VAL) &&
3055 14417 : panErrorCodes[iii] == 0)
3056 14324 : min_value = data[iii];
3057 : }
3058 1229 : return min_value;
3059 : }
3060 :
3061 : // ---------------------------------------------------------------------------
3062 1227 : static double simple_max(const double *data, const int *panErrorCodes,
3063 : const int arr_len)
3064 : {
3065 1227 : double max_value = HUGE_VAL;
3066 53847 : for (int iii = 0; iii < arr_len; iii++)
3067 : {
3068 52620 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
3069 8017 : panErrorCodes[iii] == 0)
3070 5926 : max_value = data[iii];
3071 : }
3072 1227 : return max_value;
3073 : }
3074 :
3075 : // ---------------------------------------------------------------------------
3076 53010 : static int _find_previous_index(const int iii, const int *panErrorCodes,
3077 : const int arr_len)
3078 : {
3079 : // find index of nearest valid previous value if exists
3080 53010 : int prev_iii = iii - 1;
3081 53010 : if (prev_iii == -1) // handle wraparound
3082 862 : prev_iii = arr_len - 1;
3083 53032 : while (panErrorCodes[prev_iii] != 0 && prev_iii != iii)
3084 : {
3085 22 : prev_iii--;
3086 22 : if (prev_iii == -1) // handle wraparound
3087 2 : prev_iii = arr_len - 1;
3088 : }
3089 53010 : return prev_iii;
3090 : }
3091 :
3092 : // ---------------------------------------------------------------------------
3093 : /******************************************************************************
3094 : Handles the case when longitude values cross the antimeridian
3095 : when calculating the minimum.
3096 : Note: The data array must be in a linear ring.
3097 : Note: This requires a densified ring with at least 2 additional
3098 : points per edge to correctly handle global extents.
3099 : If only 1 additional point:
3100 : | |
3101 : |RL--x0--|RL--
3102 : | |
3103 : -180 180|-180
3104 : If they are evenly spaced and it crosses the antimeridian:
3105 : x0 - L = 180
3106 : R - x0 = -180
3107 : For example:
3108 : Let R = -179.9, x0 = 0.1, L = -179.89
3109 : x0 - L = 0.1 - -179.9 = 180
3110 : R - x0 = -179.89 - 0.1 ~= -180
3111 : This is the same in the case when it didn't cross the antimeridian.
3112 : If you have 2 additional points:
3113 : | |
3114 : |RL--x0--x1--|RL--
3115 : | |
3116 : -180 180|-180
3117 : If they are evenly spaced and it crosses the antimeridian:
3118 : x0 - L = 120
3119 : x1 - x0 = 120
3120 : R - x1 = -240
3121 : For example:
3122 : Let R = -179.9, x0 = -59.9, x1 = 60.1 L = -179.89
3123 : x0 - L = 59.9 - -179.9 = 120
3124 : x1 - x0 = 60.1 - 59.9 = 120
3125 : R - x1 = -179.89 - 60.1 ~= -240
3126 : However, if they are evenly spaced and it didn't cross the antimeridian:
3127 : x0 - L = 120
3128 : x1 - x0 = 120
3129 : R - x1 = 120
3130 : From this, we have a delta that is guaranteed to be significantly
3131 : large enough to tell the difference reguarless of the direction
3132 : the antimeridian was crossed.
3133 : However, even though the spacing was even in the source projection, it isn't
3134 : guaranteed in the target geographic projection. So, instead of 240, 200 is used
3135 : as it significantly larger than 120 to be sure that the antimeridian was crossed
3136 : but smalller than 240 to account for possible irregularities in distances
3137 : when re-projecting. Also, 200 ensures latitudes are ignored for axis order
3138 : handling.
3139 : ******************************************************************************/
3140 432 : static double antimeridian_min(const double *data, const int *panErrorCodes,
3141 : const int arr_len)
3142 : {
3143 432 : double positive_min = HUGE_VAL;
3144 432 : double min_value = HUGE_VAL;
3145 432 : int crossed_meridian_count = 0;
3146 432 : bool positive_meridian = false;
3147 :
3148 26948 : for (int iii = 0; iii < arr_len; iii++)
3149 : {
3150 26516 : if (panErrorCodes[iii])
3151 11 : continue;
3152 26505 : int prev_iii = _find_previous_index(iii, panErrorCodes, arr_len);
3153 : // check if crossed meridian
3154 26505 : double delta = data[prev_iii] - data[iii];
3155 : // 180 -> -180
3156 26505 : if (delta >= 200 && delta != HUGE_VAL)
3157 : {
3158 9 : if (crossed_meridian_count == 0)
3159 2 : positive_min = min_value;
3160 9 : crossed_meridian_count++;
3161 9 : positive_meridian = false;
3162 : // -180 -> 180
3163 : }
3164 26496 : else if (delta <= -200 && delta != HUGE_VAL)
3165 : {
3166 9 : if (crossed_meridian_count == 0)
3167 6 : positive_min = data[iii];
3168 9 : crossed_meridian_count++;
3169 9 : positive_meridian = true;
3170 : }
3171 : // positive meridian side min
3172 26505 : if (positive_meridian && data[iii] < positive_min)
3173 38 : positive_min = data[iii];
3174 : // track general min value
3175 26505 : if (data[iii] < min_value)
3176 979 : min_value = data[iii];
3177 : }
3178 :
3179 432 : if (crossed_meridian_count == 2)
3180 7 : return positive_min;
3181 425 : else if (crossed_meridian_count == 4)
3182 : // bounds extends beyond -180/180
3183 1 : return -180;
3184 424 : return min_value;
3185 : }
3186 :
3187 : // ---------------------------------------------------------------------------
3188 : // Handles the case when longitude values cross the antimeridian
3189 : // when calculating the minimum.
3190 : // Note: The data array must be in a linear ring.
3191 : // Note: This requires a densified ring with at least 2 additional
3192 : // points per edge to correctly handle global extents.
3193 : // See antimeridian_min docstring for reasoning.
3194 432 : static double antimeridian_max(const double *data, const int *panErrorCodes,
3195 : const int arr_len)
3196 : {
3197 432 : double negative_max = -HUGE_VAL;
3198 432 : double max_value = -HUGE_VAL;
3199 432 : bool negative_meridian = false;
3200 432 : int crossed_meridian_count = 0;
3201 :
3202 26948 : for (int iii = 0; iii < arr_len; iii++)
3203 : {
3204 26516 : if (panErrorCodes[iii])
3205 11 : continue;
3206 26505 : int prev_iii = _find_previous_index(iii, panErrorCodes, arr_len);
3207 : // check if crossed meridian
3208 26505 : double delta = data[prev_iii] - data[iii];
3209 : // 180 -> -180
3210 26505 : if (delta >= 200 && delta != HUGE_VAL)
3211 : {
3212 9 : if (crossed_meridian_count == 0)
3213 2 : negative_max = data[iii];
3214 9 : crossed_meridian_count++;
3215 9 : negative_meridian = true;
3216 : // -180 -> 180
3217 : }
3218 26496 : else if (delta <= -200 && delta != HUGE_VAL)
3219 : {
3220 9 : if (crossed_meridian_count == 0)
3221 6 : negative_max = max_value;
3222 9 : negative_meridian = false;
3223 9 : crossed_meridian_count++;
3224 : }
3225 : // negative meridian side max
3226 26505 : if (negative_meridian &&
3227 142 : (data[iii] > negative_max || negative_max == HUGE_VAL) &&
3228 45 : panErrorCodes[iii] == 0)
3229 45 : negative_max = data[iii];
3230 : // track general max value
3231 26505 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
3232 9372 : panErrorCodes[iii] == 0)
3233 9372 : max_value = data[iii];
3234 : }
3235 432 : if (crossed_meridian_count == 2)
3236 7 : return negative_max;
3237 425 : else if (crossed_meridian_count == 4)
3238 : // bounds extends beyond -180/180
3239 1 : return 180;
3240 424 : return max_value;
3241 : }
3242 :
3243 : // ---------------------------------------------------------------------------
3244 : // Check if the original projected bounds contains
3245 : // the north pole.
3246 : // This assumes that the destination CRS is geographic.
3247 452 : bool OGRProjCT::ContainsNorthPole(const double xmin, const double ymin,
3248 : const double xmax, const double ymax,
3249 : bool lon_lat_order)
3250 : {
3251 452 : double pole_y = 90;
3252 452 : double pole_x = 0;
3253 452 : if (!lon_lat_order)
3254 : {
3255 4 : pole_y = 0;
3256 4 : pole_x = 90;
3257 : }
3258 904 : auto inverseCT = std::unique_ptr<OGRCoordinateTransformation>(GetInverse());
3259 452 : if (!inverseCT)
3260 0 : return false;
3261 452 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3262 904 : const bool success = CPL_TO_BOOL(inverseCT->TransformWithErrorCodes(
3263 452 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr));
3264 461 : return success && xmin < pole_x && pole_x < xmax && ymax > pole_y &&
3265 461 : pole_y > ymin;
3266 : }
3267 :
3268 : // ---------------------------------------------------------------------------
3269 : // Check if the original projected bounds contains
3270 : // the south pole.
3271 : // This assumes that the destination CRS is geographic.
3272 452 : bool OGRProjCT::ContainsSouthPole(const double xmin, const double ymin,
3273 : const double xmax, const double ymax,
3274 : bool lon_lat_order)
3275 : {
3276 452 : double pole_y = -90;
3277 452 : double pole_x = 0;
3278 452 : if (!lon_lat_order)
3279 : {
3280 4 : pole_y = 0;
3281 4 : pole_x = -90;
3282 : }
3283 904 : auto inverseCT = std::unique_ptr<OGRCoordinateTransformation>(GetInverse());
3284 452 : if (!inverseCT)
3285 0 : return false;
3286 452 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3287 904 : const bool success = CPL_TO_BOOL(inverseCT->TransformWithErrorCodes(
3288 452 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr));
3289 517 : return success && xmin < pole_x && pole_x < xmax && ymax > pole_y &&
3290 517 : pole_y > ymin;
3291 : }
3292 :
3293 1122 : int OGRProjCT::TransformBounds(const double xmin, const double ymin,
3294 : const double xmax, const double ymax,
3295 : double *out_xmin, double *out_ymin,
3296 : double *out_xmax, double *out_ymax,
3297 : const int densify_pts)
3298 : {
3299 1122 : if (bNoTransform)
3300 : {
3301 270 : *out_xmin = xmin;
3302 270 : *out_ymin = ymin;
3303 270 : *out_xmax = xmax;
3304 270 : *out_ymax = ymax;
3305 270 : return true;
3306 : }
3307 :
3308 852 : *out_xmin = HUGE_VAL;
3309 852 : *out_ymin = HUGE_VAL;
3310 852 : *out_xmax = HUGE_VAL;
3311 852 : *out_ymax = HUGE_VAL;
3312 :
3313 852 : if (densify_pts < 0 || densify_pts > 10000)
3314 : {
3315 1 : CPLError(CE_Failure, CPLE_AppDefined,
3316 : "densify_pts must be between 0-10000.");
3317 1 : return false;
3318 : }
3319 851 : if (!poSRSSource)
3320 : {
3321 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing source SRS.");
3322 0 : return false;
3323 : }
3324 851 : if (!poSRSTarget)
3325 : {
3326 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing target SRS.");
3327 0 : return false;
3328 : }
3329 :
3330 851 : bool degree_input = false;
3331 851 : bool degree_output = false;
3332 851 : bool input_lon_lat_order = false;
3333 851 : bool output_lon_lat_order = false;
3334 :
3335 851 : if (bSourceLatLong)
3336 : {
3337 384 : degree_input = fabs(poSRSSource->GetAngularUnits(nullptr) -
3338 384 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3339 384 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
3340 702 : if ((mapping[0] == 1 && m_eSourceFirstAxisOrient == OAO_East) ||
3341 318 : (mapping[0] == 2 && m_eSourceFirstAxisOrient != OAO_East))
3342 : {
3343 375 : input_lon_lat_order = true;
3344 : }
3345 : }
3346 851 : if (bTargetLatLong)
3347 : {
3348 453 : degree_output = fabs(poSRSTarget->GetAngularUnits(nullptr) -
3349 453 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3350 453 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
3351 860 : if ((mapping[0] == 1 && m_eTargetFirstAxisOrient == OAO_East) ||
3352 407 : (mapping[0] == 2 && m_eTargetFirstAxisOrient != OAO_East))
3353 : {
3354 448 : output_lon_lat_order = true;
3355 : }
3356 : }
3357 :
3358 851 : if (degree_output && densify_pts < 2)
3359 : {
3360 1 : CPLError(CE_Failure, CPLE_AppDefined,
3361 : "densify_pts must be at least 2 if the output is geograpic.");
3362 1 : return false;
3363 : }
3364 :
3365 850 : int side_pts = densify_pts + 1; // add one because we are densifying
3366 850 : const int boundary_len = side_pts * 4;
3367 1700 : std::vector<double> x_boundary_array;
3368 1700 : std::vector<double> y_boundary_array;
3369 1700 : std::vector<int> anErrorCodes;
3370 : try
3371 : {
3372 850 : x_boundary_array.resize(boundary_len);
3373 850 : y_boundary_array.resize(boundary_len);
3374 850 : anErrorCodes.resize(boundary_len);
3375 : }
3376 0 : catch (const std::exception &e) // memory allocation failure
3377 : {
3378 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
3379 0 : return false;
3380 : }
3381 850 : double delta_x = 0;
3382 850 : double delta_y = 0;
3383 850 : bool north_pole_in_bounds = false;
3384 850 : bool south_pole_in_bounds = false;
3385 850 : if (degree_output)
3386 : {
3387 : north_pole_in_bounds =
3388 452 : ContainsNorthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3389 : south_pole_in_bounds =
3390 452 : ContainsSouthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3391 : }
3392 :
3393 850 : if (degree_input && xmax < xmin)
3394 : {
3395 1 : if (!input_lon_lat_order)
3396 : {
3397 0 : CPLError(CE_Failure, CPLE_AppDefined,
3398 : "latitude max < latitude min.");
3399 0 : return false;
3400 : }
3401 : // handle antimeridian
3402 1 : delta_x = (xmax - xmin + 360.0) / side_pts;
3403 : }
3404 : else
3405 : {
3406 849 : delta_x = (xmax - xmin) / side_pts;
3407 : }
3408 850 : if (degree_input && ymax < ymin)
3409 : {
3410 1 : if (input_lon_lat_order)
3411 : {
3412 0 : CPLError(CE_Failure, CPLE_AppDefined,
3413 : "latitude max < latitude min.");
3414 0 : return false;
3415 : }
3416 : // handle antimeridian
3417 1 : delta_y = (ymax - ymin + 360.0) / side_pts;
3418 : }
3419 : else
3420 : {
3421 849 : delta_y = (ymax - ymin) / side_pts;
3422 : }
3423 :
3424 : // build densified bounding box
3425 : // Note: must be a linear ring for antimeridian logic
3426 11022 : for (int iii = 0; iii < side_pts; iii++)
3427 : {
3428 : // xmin boundary
3429 10172 : y_boundary_array[iii] = ymax - iii * delta_y;
3430 10172 : x_boundary_array[iii] = xmin;
3431 : // ymin boundary
3432 10172 : y_boundary_array[iii + side_pts] = ymin;
3433 10172 : x_boundary_array[iii + side_pts] = xmin + iii * delta_x;
3434 : // xmax boundary
3435 10172 : y_boundary_array[iii + side_pts * 2] = ymin + iii * delta_y;
3436 10172 : x_boundary_array[iii + side_pts * 2] = xmax;
3437 : // ymax boundary
3438 10172 : y_boundary_array[iii + side_pts * 3] = ymax;
3439 10172 : x_boundary_array[iii + side_pts * 3] = xmax - iii * delta_x;
3440 : }
3441 :
3442 : {
3443 850 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3444 1700 : bool success = CPL_TO_BOOL(TransformWithErrorCodes(
3445 850 : boundary_len, &x_boundary_array[0], &y_boundary_array[0], nullptr,
3446 : nullptr, anErrorCodes.data()));
3447 850 : if (!success)
3448 : {
3449 364 : for (int i = 0; i < boundary_len; ++i)
3450 : {
3451 356 : if (anErrorCodes[i] == 0)
3452 : {
3453 135 : success = true;
3454 135 : break;
3455 : }
3456 : }
3457 143 : if (!success)
3458 8 : return false;
3459 : }
3460 : }
3461 :
3462 842 : if (!degree_output)
3463 : {
3464 394 : *out_xmin =
3465 394 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3466 394 : *out_xmax =
3467 394 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3468 394 : *out_ymin =
3469 394 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3470 394 : *out_ymax =
3471 394 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3472 :
3473 394 : if (poSRSTarget->IsProjected())
3474 : {
3475 788 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3476 :
3477 : auto poBaseTarget = std::unique_ptr<OGRSpatialReference>(
3478 788 : poSRSTarget->CloneGeogCS());
3479 394 : poBaseTarget->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3480 :
3481 : auto poCTBaseTargetToSrc =
3482 : std::unique_ptr<OGRCoordinateTransformation>(
3483 394 : OGRCreateCoordinateTransformation(poBaseTarget.get(),
3484 788 : poSRSSource));
3485 394 : if (poCTBaseTargetToSrc)
3486 : {
3487 : const double dfLon0 =
3488 394 : poSRSTarget->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
3489 :
3490 394 : double dfSignedPoleLat = 0;
3491 394 : double dfAbsPoleLat = 90;
3492 394 : bool bIncludesPole = false;
3493 : const char *pszProjection =
3494 394 : poSRSTarget->GetAttrValue("PROJECTION");
3495 394 : if (pszProjection &&
3496 394 : EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) && dfLon0 == 0)
3497 : {
3498 : // This MAX_LAT_MERCATOR values is equivalent to the
3499 : // semi_major_axis * PI easting/northing value only
3500 : // for EPSG:3857, but it is also quite
3501 : // reasonable for other Mercator projections
3502 241 : constexpr double MAX_LAT_MERCATOR = 85.0511287798066;
3503 241 : dfAbsPoleLat = MAX_LAT_MERCATOR;
3504 : }
3505 :
3506 : // Detect if a point at long = central_meridian and
3507 : // lat = +/- 90deg is included in the extent.
3508 : // Helps for example for EPSG:4326 to ESRI:53037
3509 1182 : for (int iSign = -1; iSign <= 1; iSign += 2)
3510 : {
3511 788 : double dfX = dfLon0;
3512 788 : constexpr double EPS = 1e-8;
3513 788 : double dfY = iSign * (dfAbsPoleLat - EPS);
3514 :
3515 788 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3516 788 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3517 788 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3518 1891 : dfY <= ymax &&
3519 315 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr, nullptr,
3520 : nullptr))
3521 : {
3522 315 : dfSignedPoleLat = iSign * dfAbsPoleLat;
3523 315 : bIncludesPole = true;
3524 315 : *out_xmin = std::min(*out_xmin, dfX);
3525 315 : *out_ymin = std::min(*out_ymin, dfY);
3526 315 : *out_xmax = std::max(*out_xmax, dfX);
3527 315 : *out_ymax = std::max(*out_ymax, dfY);
3528 : }
3529 : }
3530 :
3531 : const auto TryAtPlusMinus180 =
3532 253 : [this, dfLon0, xmin, ymin, xmax, ymax, out_xmin, out_ymin,
3533 4524 : out_xmax, out_ymax, &poCTBaseTargetToSrc](double dfLat)
3534 : {
3535 759 : for (int iSign = -1; iSign <= 1; iSign += 2)
3536 : {
3537 506 : constexpr double EPS = 1e-8;
3538 : double dfX =
3539 506 : fmod(dfLon0 + iSign * (180 - EPS) + 180, 360) - 180;
3540 506 : double dfY = dfLat;
3541 :
3542 506 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3543 506 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3544 506 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3545 1335 : dfY <= ymax &&
3546 323 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr,
3547 : nullptr, nullptr))
3548 : {
3549 312 : *out_xmin = std::min(*out_xmin, dfX);
3550 312 : *out_ymin = std::min(*out_ymin, dfY);
3551 312 : *out_xmax = std::max(*out_xmax, dfX);
3552 312 : *out_ymax = std::max(*out_ymax, dfY);
3553 : }
3554 : }
3555 253 : };
3556 :
3557 : // For a projected CRS with a central meridian != 0, try to
3558 : // reproject the points with long = +/- 180deg of the central
3559 : // meridian and at lat = latitude_of_origin.
3560 394 : const double dfLat0 = poSRSTarget->GetNormProjParm(
3561 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
3562 394 : if (dfLon0 != 0)
3563 : {
3564 98 : TryAtPlusMinus180(dfLat0);
3565 : }
3566 :
3567 394 : if (bIncludesPole && dfLat0 != dfSignedPoleLat)
3568 : {
3569 155 : TryAtPlusMinus180(dfSignedPoleLat);
3570 : }
3571 : }
3572 : }
3573 : }
3574 448 : else if (north_pole_in_bounds && output_lon_lat_order)
3575 : {
3576 8 : *out_xmin = -180;
3577 8 : *out_ymin =
3578 8 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3579 8 : *out_xmax = 180;
3580 8 : *out_ymax = 90;
3581 : }
3582 440 : else if (north_pole_in_bounds)
3583 : {
3584 1 : *out_xmin =
3585 1 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3586 1 : *out_ymin = -180;
3587 1 : *out_xmax = 90;
3588 1 : *out_ymax = 180;
3589 : }
3590 439 : else if (south_pole_in_bounds && output_lon_lat_order)
3591 : {
3592 6 : *out_xmin = -180;
3593 6 : *out_ymin = -90;
3594 6 : *out_xmax = 180;
3595 6 : *out_ymax =
3596 6 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3597 : }
3598 433 : else if (south_pole_in_bounds)
3599 : {
3600 1 : *out_xmin = -90;
3601 1 : *out_ymin = -180;
3602 1 : *out_xmax =
3603 1 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3604 1 : *out_ymax = 180;
3605 : }
3606 432 : else if (output_lon_lat_order)
3607 : {
3608 430 : *out_xmin = antimeridian_min(&x_boundary_array[0], anErrorCodes.data(),
3609 : boundary_len);
3610 430 : *out_xmax = antimeridian_max(&x_boundary_array[0], anErrorCodes.data(),
3611 : boundary_len);
3612 430 : *out_ymin =
3613 430 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3614 430 : *out_ymax =
3615 430 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3616 : }
3617 : else
3618 : {
3619 2 : *out_xmin =
3620 2 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3621 2 : *out_xmax =
3622 2 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3623 2 : *out_ymin = antimeridian_min(&y_boundary_array[0], anErrorCodes.data(),
3624 : boundary_len);
3625 2 : *out_ymax = antimeridian_max(&y_boundary_array[0], anErrorCodes.data(),
3626 : boundary_len);
3627 : }
3628 :
3629 842 : return *out_xmin != HUGE_VAL && *out_ymin != HUGE_VAL &&
3630 1684 : *out_xmax != HUGE_VAL && *out_ymax != HUGE_VAL;
3631 : }
3632 :
3633 : /************************************************************************/
3634 : /* Clone() */
3635 : /************************************************************************/
3636 :
3637 22 : OGRCoordinateTransformation *OGRProjCT::Clone() const
3638 : {
3639 44 : auto poNewCT = std::make_unique<OGRProjCT>(*this);
3640 : #if (PROJ_VERSION_MAJOR * 10000 + PROJ_VERSION_MINOR * 100 + \
3641 : PROJ_VERSION_PATCH) < 80001
3642 : // See https://github.com/OSGeo/PROJ/pull/2582
3643 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3644 : // operation being a set of real operations
3645 22 : bool bCloneDone = ((m_pj == nullptr) == (poNewCT->m_pj == nullptr));
3646 22 : if (!bCloneDone)
3647 : {
3648 1 : poNewCT.reset(new OGRProjCT());
3649 : auto ret =
3650 2 : poNewCT->Initialize(poSRSSource, m_osSrcSRS.c_str(), poSRSTarget,
3651 1 : m_osTargetSRS.c_str(), m_options);
3652 1 : if (!ret)
3653 : {
3654 0 : return nullptr;
3655 : }
3656 : }
3657 : #endif // PROJ_VERSION
3658 22 : return poNewCT.release();
3659 : }
3660 :
3661 : /************************************************************************/
3662 : /* GetInverse() */
3663 : /************************************************************************/
3664 :
3665 2253 : OGRCoordinateTransformation *OGRProjCT::GetInverse() const
3666 : {
3667 2253 : PJ *new_pj = nullptr;
3668 : // m_pj can be nullptr if using m_eStrategy != PROJ
3669 2253 : if (m_pj && !bWebMercatorToWGS84LongLat && !bNoTransform)
3670 : {
3671 : // See https://github.com/OSGeo/PROJ/pull/2582
3672 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3673 : // operation being a set of real operations
3674 1944 : new_pj = proj_clone(OSRGetProjTLSContext(), m_pj);
3675 : }
3676 :
3677 4506 : OGRCoordinateTransformationOptions newOptions(m_options);
3678 2253 : std::swap(newOptions.d->bHasSourceCenterLong,
3679 2253 : newOptions.d->bHasTargetCenterLong);
3680 2253 : std::swap(newOptions.d->dfSourceCenterLong,
3681 2253 : newOptions.d->dfTargetCenterLong);
3682 2253 : newOptions.d->bReverseCO = !newOptions.d->bReverseCO;
3683 2253 : newOptions.d->RefreshCheckWithInvertProj();
3684 :
3685 2253 : if (new_pj == nullptr && !bNoTransform)
3686 : {
3687 357 : return OGRCreateCoordinateTransformation(poSRSTarget, poSRSSource,
3688 357 : newOptions);
3689 : }
3690 :
3691 1896 : auto poNewCT = new OGRProjCT();
3692 :
3693 1896 : if (poSRSTarget)
3694 1890 : poNewCT->poSRSSource = poSRSTarget->Clone();
3695 1896 : poNewCT->m_eSourceFirstAxisOrient = m_eTargetFirstAxisOrient;
3696 1896 : poNewCT->bSourceLatLong = bTargetLatLong;
3697 1896 : poNewCT->bSourceWrap = bTargetWrap;
3698 1896 : poNewCT->dfSourceWrapLong = dfTargetWrapLong;
3699 1896 : poNewCT->bSourceIsDynamicCRS = bTargetIsDynamicCRS;
3700 1896 : poNewCT->dfSourceCoordinateEpoch = dfTargetCoordinateEpoch;
3701 1896 : poNewCT->m_osSrcSRS = m_osTargetSRS;
3702 :
3703 1896 : if (poSRSSource)
3704 1890 : poNewCT->poSRSTarget = poSRSSource->Clone();
3705 1896 : poNewCT->m_eTargetFirstAxisOrient = m_eSourceFirstAxisOrient;
3706 1896 : poNewCT->bTargetLatLong = bSourceLatLong;
3707 1896 : poNewCT->bTargetWrap = bSourceWrap;
3708 1896 : poNewCT->dfTargetWrapLong = dfSourceWrapLong;
3709 1896 : poNewCT->bTargetIsDynamicCRS = bSourceIsDynamicCRS;
3710 1896 : poNewCT->dfTargetCoordinateEpoch = dfSourceCoordinateEpoch;
3711 1896 : poNewCT->m_osTargetSRS = m_osSrcSRS;
3712 :
3713 1896 : poNewCT->ComputeThreshold();
3714 :
3715 1896 : poNewCT->m_pj = new_pj;
3716 1896 : poNewCT->m_bReversePj = !m_bReversePj;
3717 1896 : poNewCT->bNoTransform = bNoTransform;
3718 1896 : poNewCT->m_eStrategy = m_eStrategy;
3719 1896 : poNewCT->m_options = newOptions;
3720 :
3721 1896 : poNewCT->DetectWebMercatorToWGS84();
3722 :
3723 1896 : return poNewCT;
3724 : }
3725 :
3726 : /************************************************************************/
3727 : /* OSRCTCleanCache() */
3728 : /************************************************************************/
3729 :
3730 1232 : void OSRCTCleanCache()
3731 : {
3732 1232 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3733 1232 : delete g_poCTCache;
3734 1232 : g_poCTCache = nullptr;
3735 1232 : }
3736 :
3737 : /************************************************************************/
3738 : /* MakeCacheKey() */
3739 : /************************************************************************/
3740 :
3741 10177 : CTCacheKey OGRProjCT::MakeCacheKey(
3742 : const OGRSpatialReference *poSRS1, const char *pszSrcSRS,
3743 : const OGRSpatialReference *poSRS2, const char *pszTargetSRS,
3744 : const OGRCoordinateTransformationOptions &options)
3745 : {
3746 : const auto GetKeyForSRS =
3747 20354 : [](const OGRSpatialReference *poSRS, const char *pszText)
3748 : {
3749 20354 : if (poSRS)
3750 : {
3751 40606 : std::string ret(pszText);
3752 20303 : const auto &mapping = poSRS->GetDataAxisToSRSAxisMapping();
3753 61204 : for (const auto &axis : mapping)
3754 : {
3755 40901 : ret += std::to_string(axis);
3756 : }
3757 20303 : return ret;
3758 : }
3759 : else
3760 : {
3761 51 : return std::string("null");
3762 : }
3763 : };
3764 :
3765 10177 : std::string ret(GetKeyForSRS(poSRS1, pszSrcSRS));
3766 10177 : ret += GetKeyForSRS(poSRS2, pszTargetSRS);
3767 10177 : ret += options.d->GetKey();
3768 20354 : return ret;
3769 : }
3770 :
3771 : /************************************************************************/
3772 : /* InsertIntoCache() */
3773 : /************************************************************************/
3774 :
3775 4525 : void OGRProjCT::InsertIntoCache(OGRProjCT *poCT)
3776 : {
3777 : {
3778 9050 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3779 4525 : if (g_poCTCache == nullptr)
3780 : {
3781 263 : g_poCTCache = new lru11::Cache<CTCacheKey, CTCacheValue>();
3782 : }
3783 : }
3784 4525 : const auto key = MakeCacheKey(poCT->poSRSSource, poCT->m_osSrcSRS.c_str(),
3785 4525 : poCT->poSRSTarget,
3786 4525 : poCT->m_osTargetSRS.c_str(), poCT->m_options);
3787 :
3788 4525 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3789 4525 : if (g_poCTCache->contains(key))
3790 : {
3791 1048 : delete poCT;
3792 1048 : return;
3793 : }
3794 3477 : g_poCTCache->insert(key, std::unique_ptr<OGRProjCT>(poCT));
3795 : }
3796 :
3797 : /************************************************************************/
3798 : /* FindFromCache() */
3799 : /************************************************************************/
3800 :
3801 6054 : OGRProjCT *OGRProjCT::FindFromCache(
3802 : const OGRSpatialReference *poSource, const char *pszSrcSRS,
3803 : const OGRSpatialReference *poTarget, const char *pszTargetSRS,
3804 : const OGRCoordinateTransformationOptions &options)
3805 : {
3806 : {
3807 6054 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3808 6054 : if (g_poCTCache == nullptr || g_poCTCache->empty())
3809 402 : return nullptr;
3810 : }
3811 :
3812 : const auto key =
3813 11304 : MakeCacheKey(poSource, pszSrcSRS, poTarget, pszTargetSRS, options);
3814 : // Get value from cache and remove it
3815 11304 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3816 5652 : CTCacheValue *cachedValue = g_poCTCache->getPtr(key);
3817 5652 : if (cachedValue)
3818 : {
3819 2290 : auto poCT = cachedValue->release();
3820 2290 : g_poCTCache->remove(key);
3821 2290 : return poCT;
3822 : }
3823 3362 : return nullptr;
3824 : }
3825 :
3826 : //! @endcond
3827 :
3828 : /************************************************************************/
3829 : /* OCTTransform() */
3830 : /************************************************************************/
3831 :
3832 : /** Transform an array of points
3833 : *
3834 : * @param hTransform Transformation object
3835 : * @param nCount Number of points
3836 : * @param x Array of nCount x values.
3837 : * @param y Array of nCount y values.
3838 : * @param z Array of nCount z values.
3839 : * @return TRUE if a transformation could be found (but not all points may
3840 : * have necessarily succeed to transform), otherwise FALSE.
3841 : */
3842 1055 : int CPL_STDCALL OCTTransform(OGRCoordinateTransformationH hTransform,
3843 : int nCount, double *x, double *y, double *z)
3844 :
3845 : {
3846 1055 : VALIDATE_POINTER1(hTransform, "OCTTransform", FALSE);
3847 :
3848 : return OGRCoordinateTransformation::FromHandle(hTransform)
3849 1055 : ->Transform(nCount, x, y, z);
3850 : }
3851 :
3852 : /************************************************************************/
3853 : /* OCTTransformEx() */
3854 : /************************************************************************/
3855 :
3856 : /** Transform an array of points
3857 : *
3858 : * @param hTransform Transformation object
3859 : * @param nCount Number of points
3860 : * @param x Array of nCount x values.
3861 : * @param y Array of nCount y values.
3862 : * @param z Array of nCount z values.
3863 : * @param pabSuccess Output array of nCount value that will be set to TRUE/FALSE
3864 : * @return TRUE if a transformation could be found (but not all points may
3865 : * have necessarily succeed to transform), otherwise FALSE.
3866 : */
3867 0 : int CPL_STDCALL OCTTransformEx(OGRCoordinateTransformationH hTransform,
3868 : int nCount, double *x, double *y, double *z,
3869 : int *pabSuccess)
3870 :
3871 : {
3872 0 : VALIDATE_POINTER1(hTransform, "OCTTransformEx", FALSE);
3873 :
3874 : return OGRCoordinateTransformation::FromHandle(hTransform)
3875 0 : ->Transform(nCount, x, y, z, pabSuccess);
3876 : }
3877 :
3878 : /************************************************************************/
3879 : /* OCTTransform4D() */
3880 : /************************************************************************/
3881 :
3882 : /** Transform an array of points
3883 : *
3884 : * @param hTransform Transformation object
3885 : * @param nCount Number of points
3886 : * @param x Array of nCount x values. Should not be NULL
3887 : * @param y Array of nCount y values. Should not be NULL
3888 : * @param z Array of nCount z values. Might be NULL
3889 : * @param t Array of nCount time values. Might be NULL
3890 : * @param pabSuccess Output array of nCount value that will be set to
3891 : * TRUE/FALSE. Might be NULL.
3892 : * @since GDAL 3.0
3893 : * @return TRUE if a transformation could be found (but not all points may
3894 : * have necessarily succeed to transform), otherwise FALSE.
3895 : */
3896 18 : int OCTTransform4D(OGRCoordinateTransformationH hTransform, int nCount,
3897 : double *x, double *y, double *z, double *t, int *pabSuccess)
3898 :
3899 : {
3900 18 : VALIDATE_POINTER1(hTransform, "OCTTransform4D", FALSE);
3901 :
3902 18 : return OGRCoordinateTransformation::FromHandle(hTransform)
3903 18 : ->Transform(nCount, x, y, z, t, pabSuccess);
3904 : }
3905 :
3906 : /************************************************************************/
3907 : /* OCTTransform4DWithErrorCodes() */
3908 : /************************************************************************/
3909 :
3910 : /** Transform an array of points
3911 : *
3912 : * @param hTransform Transformation object
3913 : * @param nCount Number of points
3914 : * @param x Array of nCount x values. Should not be NULL
3915 : * @param y Array of nCount y values. Should not be NULL
3916 : * @param z Array of nCount z values. Might be NULL
3917 : * @param t Array of nCount time values. Might be NULL
3918 : * @param panErrorCodes Output array of nCount value that will be set to 0 for
3919 : * success, or a non-zero value for failure. Refer to
3920 : * PROJ 8 public error codes. Might be NULL
3921 : * @since GDAL 3.3, and PROJ 8 to be able to use PROJ public error codes
3922 : * @return TRUE if a transformation could be found (but not all points may
3923 : * have necessarily succeed to transform), otherwise FALSE.
3924 : */
3925 0 : int OCTTransform4DWithErrorCodes(OGRCoordinateTransformationH hTransform,
3926 : int nCount, double *x, double *y, double *z,
3927 : double *t, int *panErrorCodes)
3928 :
3929 : {
3930 0 : VALIDATE_POINTER1(hTransform, "OCTTransform4DWithErrorCodes", FALSE);
3931 :
3932 0 : return OGRCoordinateTransformation::FromHandle(hTransform)
3933 0 : ->TransformWithErrorCodes(nCount, x, y, z, t, panErrorCodes);
3934 : }
3935 :
3936 : /************************************************************************/
3937 : /* OCTTransformBounds() */
3938 : /************************************************************************/
3939 : /** \brief Transform boundary.
3940 : *
3941 : * Transform boundary densifying the edges to account for nonlinear
3942 : * transformations along these edges and extracting the outermost bounds.
3943 : *
3944 : * If the destination CRS is geographic, the first axis is longitude,
3945 : * and *out_xmax < *out_xmin then the bounds crossed the antimeridian.
3946 : * In this scenario there are two polygons, one on each side of the
3947 : * antimeridian. The first polygon should be constructed with
3948 : * (*out_xmin, *out_ymin, 180, ymax) and the second with
3949 : * (-180, *out_ymin, *out_xmax, *out_ymax).
3950 : *
3951 : * If the destination CRS is geographic, the first axis is latitude,
3952 : * and *out_ymax < *out_ymin then the bounds crossed the antimeridian.
3953 : * In this scenario there are two polygons, one on each side of the
3954 : * antimeridian. The first polygon should be constructed with
3955 : * (*out_ymin, *out_xmin, *out_ymax, 180) and the second with
3956 : * (*out_ymin, -180, *out_ymax, *out_xmax).
3957 : *
3958 : * @param hTransform Transformation object
3959 : * @param xmin Minimum bounding coordinate of the first axis in source CRS.
3960 : * @param ymin Minimum bounding coordinate of the second axis in source CRS.
3961 : * @param xmax Maximum bounding coordinate of the first axis in source CRS.
3962 : * @param ymax Maximum bounding coordinate of the second axis in source CRS.
3963 : * @param out_xmin Minimum bounding coordinate of the first axis in target CRS
3964 : * @param out_ymin Minimum bounding coordinate of the second axis in target CRS.
3965 : * @param out_xmax Maximum bounding coordinate of the first axis in target CRS.
3966 : * @param out_ymax Maximum bounding coordinate of the second axis in target CRS.
3967 : * @param densify_pts Recommended to use 21. This is the number of points
3968 : * to use to densify the bounding polygon in the transformation.
3969 : * @return TRUE if successful. FALSE if failures encountered.
3970 : * @since 3.4
3971 : */
3972 25 : int CPL_STDCALL OCTTransformBounds(OGRCoordinateTransformationH hTransform,
3973 : const double xmin, const double ymin,
3974 : const double xmax, const double ymax,
3975 : double *out_xmin, double *out_ymin,
3976 : double *out_xmax, double *out_ymax,
3977 : int densify_pts)
3978 :
3979 : {
3980 25 : VALIDATE_POINTER1(hTransform, "TransformBounds", FALSE);
3981 :
3982 25 : return OGRProjCT::FromHandle(hTransform)
3983 25 : ->TransformBounds(xmin, ymin, xmax, ymax, out_xmin, out_ymin, out_xmax,
3984 25 : out_ymax, densify_pts);
3985 : }
3986 :
3987 : /************************************************************************/
3988 : /* OGRCTDumpStatistics() */
3989 : /************************************************************************/
3990 :
3991 1232 : void OGRCTDumpStatistics()
3992 : {
3993 : #ifdef DEBUG_PERF
3994 : CPLDebug("OGR_CT", "Total time in proj_create_crs_to_crs(): %d ms",
3995 : static_cast<int>(g_dfTotalTimeCRStoCRS * 1000));
3996 : CPLDebug("OGR_CT", "Total time in coordinate transformation: %d ms",
3997 : static_cast<int>(g_dfTotalTimeReprojection * 1000));
3998 : #endif
3999 1232 : }
|
