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 4282 : 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 1662 : 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 8961 : OGRCoordinateTransformationOptions::Private::Private()
124 : {
125 8960 : RefreshCheckWithInvertProj();
126 8961 : }
127 :
128 : /************************************************************************/
129 : /* GetKey() */
130 : /************************************************************************/
131 :
132 8454 : std::string OGRCoordinateTransformationOptions::Private::GetKey() const
133 : {
134 8454 : std::string ret;
135 8454 : ret += std::to_string(static_cast<int>(bHasAreaOfInterest));
136 8454 : ret += std::to_string(dfWestLongitudeDeg);
137 8454 : ret += std::to_string(dfSouthLatitudeDeg);
138 8454 : ret += std::to_string(dfEastLongitudeDeg);
139 8454 : ret += std::to_string(dfNorthLatitudeDeg);
140 8454 : ret += osCoordOperation;
141 8454 : ret += std::to_string(static_cast<int>(bReverseCO));
142 8454 : ret += std::to_string(static_cast<int>(bAllowBallpark));
143 8454 : ret += std::to_string(dfAccuracy);
144 8454 : ret += std::to_string(static_cast<int>(bOnlyBestOptionSet));
145 8454 : ret += std::to_string(static_cast<int>(bOnlyBest));
146 8454 : ret += std::to_string(static_cast<int>(bHasSourceCenterLong));
147 8454 : ret += std::to_string(dfSourceCenterLong);
148 8454 : ret += std::to_string(static_cast<int>(bHasTargetCenterLong));
149 8454 : ret += std::to_string(dfTargetCenterLong);
150 8454 : ret += std::to_string(static_cast<int>(bCheckWithInvertProj));
151 8454 : return ret;
152 : }
153 :
154 : /************************************************************************/
155 : /* RefreshCheckWithInvertProj() */
156 : /************************************************************************/
157 :
158 10583 : void OGRCoordinateTransformationOptions::Private::RefreshCheckWithInvertProj()
159 : {
160 10585 : bCheckWithInvertProj =
161 10583 : CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO"));
162 10585 : }
163 :
164 : /************************************************************************/
165 : /* GetAsAProjRecognizableString() */
166 : /************************************************************************/
167 :
168 2357 : static char *GetAsAProjRecognizableString(const OGRSpatialReference *poSRS)
169 : {
170 2357 : 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 2357 : char *pszText = nullptr;
174 2357 : if (poSRS->GetExtension(nullptr, "PROJ4", nullptr))
175 : {
176 15 : poSRS->exportToProj4(&pszText);
177 15 : if (strstr(pszText, " +type=crs") == nullptr)
178 : {
179 30 : auto tmpText = std::string(pszText) + " +type=crs";
180 15 : CPLFree(pszText);
181 15 : pszText = CPLStrdup(tmpText.c_str());
182 : }
183 : }
184 2342 : 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 2339 : poSRS->exportToPROJJSON(&pszText, nullptr);
196 : }
197 :
198 4714 : return pszText;
199 : }
200 :
201 : /************************************************************************/
202 : /* GetTextRepresentation() */
203 : /************************************************************************/
204 :
205 9657 : static char *GetTextRepresentation(const OGRSpatialReference *poSRS)
206 : {
207 7302 : const auto CanUseAuthorityDef = [](const OGRSpatialReference *poSRS1,
208 : OGRSpatialReference *poSRSFromAuth,
209 : const char *pszAuth)
210 : {
211 14517 : if (EQUAL(pszAuth, "EPSG") &&
212 7215 : 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 7300 : return true;
236 : };
237 :
238 9657 : 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 9657 : const char *pszAuth = poSRS->GetAuthorityName(nullptr);
244 9657 : const char *pszCode = poSRS->GetAuthorityCode(nullptr);
245 17010 : if (pszAuth && pszCode &&
246 7353 : CPLTestBool(
247 : CPLGetConfigOption("OGR_CT_PREFER_OFFICIAL_SRS_DEF", "YES")))
248 : {
249 14706 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
250 :
251 14706 : CPLString osAuthCode(pszAuth);
252 7353 : osAuthCode += ':';
253 7353 : osAuthCode += pszCode;
254 14706 : OGRSpatialReference oTmpSRS;
255 7353 : if (oTmpSRS.SetFromUserInput(osAuthCode) == OGRERR_NONE)
256 : {
257 7340 : const char *pszAuthAfter = oTmpSRS.GetAuthorityName(nullptr);
258 7340 : const char *pszCodeAfter = oTmpSRS.GetAuthorityCode(nullptr);
259 7340 : if (pszAuthAfter && pszCodeAfter && EQUAL(pszAuthAfter, pszAuth) &&
260 7339 : EQUAL(pszCodeAfter, pszCode))
261 : {
262 7339 : oTmpSRS.SetDataAxisToSRSAxisMapping(
263 : poSRS->GetDataAxisToSRSAxisMapping());
264 7339 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
265 : nullptr};
266 7339 : if (oTmpSRS.IsSame(poSRS, apszOptionsIsSame))
267 : {
268 7302 : if (CanUseAuthorityDef(poSRS, &oTmpSRS, pszAuth))
269 : {
270 7300 : pszText = CPLStrdup(osAuthCode);
271 : }
272 : }
273 : }
274 : }
275 : }
276 9657 : if (pszText == nullptr)
277 : {
278 2357 : pszText = GetAsAProjRecognizableString(poSRS);
279 : }
280 9657 : return pszText;
281 : }
282 :
283 : /************************************************************************/
284 : /* OGRCoordinateTransformationOptions() */
285 : /************************************************************************/
286 :
287 : /** \brief Constructs a new OGRCoordinateTransformationOptions.
288 : *
289 : * @since GDAL 3.0
290 : */
291 8961 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions()
292 8961 : : d(new Private())
293 : {
294 8961 : }
295 :
296 : /************************************************************************/
297 : /* OGRCoordinateTransformationOptions() */
298 : /************************************************************************/
299 :
300 : /** \brief Copy constructor
301 : *
302 : * @since GDAL 3.1
303 : */
304 1663 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions(
305 1663 : const OGRCoordinateTransformationOptions &other)
306 1663 : : d(new Private(*(other.d)))
307 : {
308 1662 : }
309 :
310 : /************************************************************************/
311 : /* operator =() */
312 : /************************************************************************/
313 :
314 : /** \brief Assignment operator
315 : *
316 : * @since GDAL 3.1
317 : */
318 : OGRCoordinateTransformationOptions &
319 4282 : OGRCoordinateTransformationOptions::operator=(
320 : const OGRCoordinateTransformationOptions &other)
321 : {
322 4282 : if (this != &other)
323 : {
324 4282 : *d = *(other.d);
325 : }
326 4281 : return *this;
327 : }
328 :
329 : /************************************************************************/
330 : /* OGRCoordinateTransformationOptions() */
331 : /************************************************************************/
332 :
333 : /** \brief Destroys a OGRCoordinateTransformationOptions.
334 : *
335 : * @since GDAL 3.0
336 : */
337 10032 : OGRCoordinateTransformationOptions::~OGRCoordinateTransformationOptions()
338 : {
339 10033 : }
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 1111 : bool OGRCoordinateTransformationOptions::SetAreaOfInterest(
388 : double dfWestLongitudeDeg, double dfSouthLatitudeDeg,
389 : double dfEastLongitudeDeg, double dfNorthLatitudeDeg)
390 : {
391 1111 : if (std::fabs(dfWestLongitudeDeg) > 180)
392 : {
393 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfWestLongitudeDeg");
394 1 : return false;
395 : }
396 1110 : if (std::fabs(dfSouthLatitudeDeg) > 90)
397 : {
398 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfSouthLatitudeDeg");
399 1 : return false;
400 : }
401 1109 : if (std::fabs(dfEastLongitudeDeg) > 180)
402 : {
403 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfEastLongitudeDeg");
404 1 : return false;
405 : }
406 1108 : if (std::fabs(dfNorthLatitudeDeg) > 90)
407 : {
408 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfNorthLatitudeDeg");
409 1 : return false;
410 : }
411 1107 : if (dfSouthLatitudeDeg > dfNorthLatitudeDeg)
412 : {
413 0 : CPLError(CE_Failure, CPLE_AppDefined,
414 : "dfSouthLatitudeDeg should be lower than dfNorthLatitudeDeg");
415 0 : return false;
416 : }
417 1107 : d->bHasAreaOfInterest = true;
418 1107 : d->dfWestLongitudeDeg = dfWestLongitudeDeg;
419 1107 : d->dfSouthLatitudeDeg = dfSouthLatitudeDeg;
420 1107 : d->dfEastLongitudeDeg = dfEastLongitudeDeg;
421 1107 : d->dfNorthLatitudeDeg = dfNorthLatitudeDeg;
422 1107 : 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 798 : void OGRCoordinateTransformationOptions::SetSourceCenterLong(
480 : double dfCenterLong)
481 : {
482 798 : d->dfSourceCenterLong = dfCenterLong;
483 798 : d->bHasSourceCenterLong = true;
484 798 : }
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 a 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 1172 : bool OGRCoordinateTransformationOptions::SetBallparkAllowed(bool bAllowBallpark)
584 : {
585 1172 : d->bAllowBallpark = bAllowBallpark;
586 1172 : 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 7821 : void reset()
679 : {
680 7821 : if (m_pj)
681 : {
682 3543 : proj_assign_context(m_pj, OSRGetProjTLSContext());
683 3543 : proj_destroy(m_pj);
684 : }
685 7821 : }
686 :
687 : public:
688 4283 : PjPtr() : m_pj(nullptr)
689 : {
690 4283 : }
691 :
692 0 : explicit PjPtr(PJ *pjIn) : m_pj(pjIn)
693 : {
694 0 : }
695 :
696 3722 : ~PjPtr()
697 3722 : {
698 3722 : reset();
699 3722 : }
700 :
701 27 : PjPtr(const PjPtr &other)
702 27 : : m_pj((other.m_pj != nullptr)
703 27 : ? (proj_clone(OSRGetProjTLSContext(), other.m_pj))
704 27 : : (nullptr))
705 : {
706 27 : }
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 4121 : PjPtr &operator=(PJ *pjIn)
726 : {
727 4121 : if (m_pj != pjIn)
728 : {
729 4100 : reset();
730 4099 : m_pj = pjIn;
731 : }
732 4120 : return *this;
733 : }
734 :
735 2098470 : operator PJ *()
736 : {
737 2098470 : return m_pj;
738 : }
739 :
740 3127 : operator const PJ *() const
741 : {
742 3127 : 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(const OGRProjCT &other);
822 : OGRProjCT &operator=(const OGRProjCT &) = delete;
823 :
824 : static CTCacheKey
825 : MakeCacheKey(const OGRSpatialReference *poSRS1, const char *pszSrcSRS,
826 : const OGRSpatialReference *poSRS2, const char *pszTargetSRS,
827 : const OGRCoordinateTransformationOptions &options);
828 : bool ContainsNorthPole(const double xmin, const double ymin,
829 : const double xmax, const double ymax,
830 : bool lon_lat_order);
831 : bool ContainsSouthPole(const double xmin, const double ymin,
832 : const double xmax, const double ymax,
833 : bool lon_lat_order);
834 :
835 : public:
836 : OGRProjCT();
837 : ~OGRProjCT() override;
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 41 : bool GetEmitErrors() const override
859 : {
860 41 : return m_bEmitErrors;
861 : }
862 :
863 4153 : void SetEmitErrors(bool bEmitErrors) override
864 : {
865 4153 : m_bEmitErrors = bEmitErrors;
866 4153 : }
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 576 : void OGRProjCTDifferentOperationsStart(OGRCoordinateTransformation *poCT)
885 : {
886 576 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
887 576 : if (poOGRCT)
888 : {
889 576 : poOGRCT->m_recordDifferentOperationsUsed = true;
890 576 : poOGRCT->m_differentOperationsUsed = false;
891 576 : 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 576 : }
900 :
901 : /************************************************************************/
902 : /* OGRProjCTDifferentOperationsStop() */
903 : /************************************************************************/
904 :
905 576 : void OGRProjCTDifferentOperationsStop(OGRCoordinateTransformation *poCT)
906 : {
907 576 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
908 576 : if (poOGRCT)
909 : {
910 576 : 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 576 : }
919 :
920 : /************************************************************************/
921 : /* OGRProjCTDifferentOperationsUsed() */
922 : /************************************************************************/
923 :
924 576 : bool OGRProjCTDifferentOperationsUsed(OGRCoordinateTransformation *poCT)
925 : {
926 576 : auto poOGRCT = dynamic_cast<OGRProjCT *>(poCT);
927 576 : if (poOGRCT)
928 : {
929 576 : 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 363 : OCTDestroyCoordinateTransformation(OGRCoordinateTransformationH hCT)
962 :
963 : {
964 363 : OGRCoordinateTransformation::DestroyCT(
965 : OGRCoordinateTransformation::FromHandle(hCT));
966 363 : }
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 4001 : void OGRCoordinateTransformation::DestroyCT(OGRCoordinateTransformation *poCT)
989 : {
990 4001 : auto poProjCT = dynamic_cast<OGRProjCT *>(poCT);
991 4001 : if (poProjCT)
992 : {
993 3995 : OGRProjCT::InsertIntoCache(poProjCT);
994 : }
995 : else
996 : {
997 6 : delete poCT;
998 : }
999 4001 : }
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 3365 : OGRCreateCoordinateTransformation(const OGRSpatialReference *poSource,
1028 : const OGRSpatialReference *poTarget)
1029 :
1030 : {
1031 3365 : return OGRCreateCoordinateTransformation(
1032 6730 : 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 4841 : OGRCoordinateTransformation *OGRCreateCoordinateTransformation(
1101 : const OGRSpatialReference *poSource, const OGRSpatialReference *poTarget,
1102 : const OGRCoordinateTransformationOptions &options)
1103 :
1104 : {
1105 4841 : char *pszSrcSRS = poSource ? GetTextRepresentation(poSource) : nullptr;
1106 4841 : char *pszTargetSRS = poTarget ? GetTextRepresentation(poTarget) : nullptr;
1107 : // Try to find if we have a match in the case
1108 4841 : OGRProjCT *poCT = OGRProjCT::FindFromCache(poSource, pszSrcSRS, poTarget,
1109 : pszTargetSRS, options);
1110 4841 : if (poCT == nullptr)
1111 : {
1112 2822 : poCT = new OGRProjCT();
1113 2822 : if (!poCT->Initialize(poSource, pszSrcSRS, poTarget, pszTargetSRS,
1114 : options))
1115 : {
1116 27 : delete poCT;
1117 27 : poCT = nullptr;
1118 : }
1119 : }
1120 4841 : CPLFree(pszSrcSRS);
1121 4840 : CPLFree(pszTargetSRS);
1122 :
1123 4841 : 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 337 : OGRCoordinateTransformationH CPL_STDCALL OCTNewCoordinateTransformation(
1153 : OGRSpatialReferenceH hSourceSRS, OGRSpatialReferenceH hTargetSRS)
1154 :
1155 : {
1156 : return reinterpret_cast<OGRCoordinateTransformationH>(
1157 337 : OGRCreateCoordinateTransformation(
1158 : reinterpret_cast<OGRSpatialReference *>(hSourceSRS),
1159 337 : 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 1 : OGRSpatialReferenceH OCTGetSourceCS(OGRCoordinateTransformationH hTransform)
1260 :
1261 : {
1262 1 : VALIDATE_POINTER1(hTransform, "OCTGetSourceCS", nullptr);
1263 1 : return OGRSpatialReference::ToHandle(const_cast<OGRSpatialReference *>(
1264 1 : 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 1 : OGRSpatialReferenceH OCTGetTargetCS(OGRCoordinateTransformationH hTransform)
1287 :
1288 : {
1289 1 : VALIDATE_POINTER1(hTransform, "OCTGetTargetCS", nullptr);
1290 1 : return OGRSpatialReference::ToHandle(const_cast<OGRSpatialReference *>(
1291 1 : 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 2 : OCTGetInverse(OGRCoordinateTransformationH hTransform)
1312 :
1313 : {
1314 2 : VALIDATE_POINTER1(hTransform, "OCTGetInverse", nullptr);
1315 2 : return OGRCoordinateTransformation::ToHandle(
1316 4 : OGRCoordinateTransformation::FromHandle(hTransform)->GetInverse());
1317 : }
1318 :
1319 : /************************************************************************/
1320 : /* OGRProjCT() */
1321 : /************************************************************************/
1322 :
1323 : //! @cond Doxygen_Suppress
1324 4283 : OGRProjCT::OGRProjCT()
1325 : {
1326 4283 : }
1327 :
1328 : /************************************************************************/
1329 : /* OGRProjCT(const OGRProjCT& other) */
1330 : /************************************************************************/
1331 :
1332 27 : OGRProjCT::OGRProjCT(const OGRProjCT &other)
1333 27 : : poSRSSource((other.poSRSSource != nullptr) ? (other.poSRSSource->Clone())
1334 : : (nullptr)),
1335 27 : m_eSourceFirstAxisOrient(other.m_eSourceFirstAxisOrient),
1336 27 : bSourceLatLong(other.bSourceLatLong), bSourceWrap(other.bSourceWrap),
1337 27 : dfSourceWrapLong(other.dfSourceWrapLong),
1338 27 : bSourceIsDynamicCRS(other.bSourceIsDynamicCRS),
1339 27 : dfSourceCoordinateEpoch(other.dfSourceCoordinateEpoch),
1340 27 : m_osSrcSRS(other.m_osSrcSRS),
1341 27 : poSRSTarget((other.poSRSTarget != nullptr) ? (other.poSRSTarget->Clone())
1342 : : (nullptr)),
1343 27 : m_eTargetFirstAxisOrient(other.m_eTargetFirstAxisOrient),
1344 27 : bTargetLatLong(other.bTargetLatLong), bTargetWrap(other.bTargetWrap),
1345 27 : dfTargetWrapLong(other.dfTargetWrapLong),
1346 27 : bTargetIsDynamicCRS(other.bTargetIsDynamicCRS),
1347 27 : dfTargetCoordinateEpoch(other.dfTargetCoordinateEpoch),
1348 27 : m_osTargetSRS(other.m_osTargetSRS),
1349 27 : bWebMercatorToWGS84LongLat(other.bWebMercatorToWGS84LongLat),
1350 27 : nErrorCount(other.nErrorCount), dfThreshold(other.dfThreshold),
1351 54 : m_pj(other.m_pj), m_bReversePj(other.m_bReversePj),
1352 27 : m_bEmitErrors(other.m_bEmitErrors), bNoTransform(other.bNoTransform),
1353 27 : m_eStrategy(other.m_eStrategy),
1354 27 : m_oTransformations(other.m_oTransformations),
1355 27 : m_iCurTransformation(other.m_iCurTransformation),
1356 27 : m_options(other.m_options), m_recordDifferentOperationsUsed(false),
1357 81 : m_lastPjUsedPROJString(std::string()), m_differentOperationsUsed(false)
1358 : {
1359 27 : }
1360 :
1361 : /************************************************************************/
1362 : /* ~OGRProjCT() */
1363 : /************************************************************************/
1364 :
1365 7444 : OGRProjCT::~OGRProjCT()
1366 :
1367 : {
1368 3722 : if (poSRSSource != nullptr)
1369 : {
1370 3711 : poSRSSource->Release();
1371 : }
1372 :
1373 3722 : if (poSRSTarget != nullptr)
1374 : {
1375 3709 : poSRSTarget->Release();
1376 : }
1377 7444 : }
1378 :
1379 : /************************************************************************/
1380 : /* ComputeThreshold() */
1381 : /************************************************************************/
1382 :
1383 4282 : void OGRProjCT::ComputeThreshold()
1384 : {
1385 : // The threshold is experimental. Works well with the cases of ticket #2305.
1386 4282 : if (bSourceLatLong)
1387 : {
1388 : // coverity[tainted_data]
1389 2948 : 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 1334 : dfThreshold = CPLAtof(CPLGetConfigOption("THRESHOLD", "10000"));
1397 : }
1398 4282 : }
1399 :
1400 : /************************************************************************/
1401 : /* DetectWebMercatorToWGS84() */
1402 : /************************************************************************/
1403 :
1404 4282 : void OGRProjCT::DetectWebMercatorToWGS84()
1405 : {
1406 : // Detect webmercator to WGS84
1407 8527 : if (m_options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1408 13821 : poSRSSource->IsProjected() && poSRSTarget->IsGeographic() &&
1409 1014 : ((m_eTargetFirstAxisOrient == OAO_North &&
1410 896 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1411 5177 : std::vector<int>{2, 1}) ||
1412 139 : (m_eTargetFirstAxisOrient == OAO_East &&
1413 118 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1414 4399 : 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 984 : const char *pszSourceAuth = poSRSSource->GetAuthorityName(nullptr);
1420 984 : const char *pszSourceCode = poSRSSource->GetAuthorityCode(nullptr);
1421 985 : const char *pszTargetAuth = poSRSTarget->GetAuthorityName(nullptr);
1422 985 : const char *pszTargetCode = poSRSTarget->GetAuthorityCode(nullptr);
1423 985 : if (pszSourceAuth && pszSourceCode && pszTargetAuth && pszTargetCode &&
1424 681 : EQUAL(pszSourceAuth, "EPSG") && EQUAL(pszTargetAuth, "EPSG"))
1425 : {
1426 669 : bWebMercatorToWGS84LongLat =
1427 937 : (EQUAL(pszSourceCode, "3857") ||
1428 268 : EQUAL(pszSourceCode, "3785") || // deprecated
1429 1338 : EQUAL(pszSourceCode, "900913")) && // deprecated
1430 402 : EQUAL(pszTargetCode, "4326");
1431 : }
1432 : else
1433 : {
1434 316 : CPLPushErrorHandler(CPLQuietErrorHandler);
1435 316 : char *pszSrcProj4Defn = nullptr;
1436 316 : poSRSSource->exportToProj4(&pszSrcProj4Defn);
1437 :
1438 316 : char *pszDstProj4Defn = nullptr;
1439 316 : poSRSTarget->exportToProj4(&pszDstProj4Defn);
1440 316 : CPLPopErrorHandler();
1441 :
1442 316 : if (pszSrcProj4Defn && pszDstProj4Defn)
1443 : {
1444 316 : if (pszSrcProj4Defn[0] != '\0' &&
1445 316 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] == ' ')
1446 0 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] = 0;
1447 316 : if (pszDstProj4Defn[0] != '\0' &&
1448 316 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] == ' ')
1449 0 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] = 0;
1450 316 : char *pszNeedle = strstr(pszSrcProj4Defn, " ");
1451 316 : if (pszNeedle)
1452 0 : memmove(pszNeedle, pszNeedle + 1,
1453 0 : strlen(pszNeedle + 1) + 1);
1454 316 : pszNeedle = strstr(pszDstProj4Defn, " ");
1455 316 : if (pszNeedle)
1456 0 : memmove(pszNeedle, pszNeedle + 1,
1457 0 : strlen(pszNeedle + 1) + 1);
1458 :
1459 316 : if ((strstr(pszDstProj4Defn, "+datum=WGS84") != nullptr ||
1460 160 : strstr(pszDstProj4Defn,
1461 : "+ellps=WGS84 +towgs84=0,0,0,0,0,0,0 ") !=
1462 156 : nullptr) &&
1463 156 : strstr(pszSrcProj4Defn, "+nadgrids=@null ") != nullptr &&
1464 56 : strstr(pszSrcProj4Defn, "+towgs84") == nullptr)
1465 : {
1466 : char *pszDst =
1467 56 : strstr(pszDstProj4Defn, "+towgs84=0,0,0,0,0,0,0 ");
1468 56 : 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 56 : memcpy(strstr(pszDstProj4Defn, "+datum=WGS84"),
1477 : "+ellps", 6);
1478 : }
1479 :
1480 56 : pszDst = strstr(pszSrcProj4Defn, "+nadgrids=@null ");
1481 56 : char *pszSrc = pszDst + strlen("+nadgrids=@null ");
1482 56 : memmove(pszDst, pszSrc, strlen(pszSrc) + 1);
1483 :
1484 56 : pszDst = strstr(pszSrcProj4Defn, "+wktext ");
1485 56 : if (pszDst)
1486 : {
1487 56 : pszSrc = pszDst + strlen("+wktext ");
1488 56 : memmove(pszDst, pszSrc, strlen(pszSrc) + 1);
1489 : }
1490 56 : bWebMercatorToWGS84LongLat =
1491 56 : strcmp(pszDstProj4Defn,
1492 112 : "+proj=longlat +ellps=WGS84 +no_defs") == 0 &&
1493 56 : (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 56 : 0 ||
1498 56 : 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 316 : CPLFree(pszSrcProj4Defn);
1506 316 : CPLFree(pszDstProj4Defn);
1507 : }
1508 :
1509 985 : if (bWebMercatorToWGS84LongLat)
1510 : {
1511 457 : CPLDebug("OGRCT", "Using WebMercator to WGS84 optimization");
1512 : }
1513 : }
1514 4282 : }
1515 :
1516 : /************************************************************************/
1517 : /* Initialize() */
1518 : /************************************************************************/
1519 :
1520 2823 : 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 2823 : m_options = options;
1528 :
1529 2822 : 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 2821 : if (poSourceIn)
1541 : {
1542 2811 : poSRSSource = poSourceIn->Clone();
1543 2811 : m_osSrcSRS = pszSrcSRS;
1544 : }
1545 2822 : if (poTargetIn)
1546 : {
1547 2811 : poSRSTarget = poTargetIn->Clone();
1548 2811 : m_osTargetSRS = pszTargetSRS;
1549 : }
1550 :
1551 : // To easy quick&dirty compatibility with GDAL < 3.0
1552 2822 : 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 2822 : if (poSRSSource)
1562 : {
1563 2812 : bSourceLatLong = CPL_TO_BOOL(poSRSSource->IsGeographic());
1564 2812 : bSourceIsDynamicCRS = poSRSSource->IsDynamic();
1565 2812 : dfSourceCoordinateEpoch = poSRSSource->GetCoordinateEpoch();
1566 2812 : if (!bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0)
1567 : {
1568 3 : bSourceIsDynamicCRS = poSRSSource->HasPointMotionOperation();
1569 : }
1570 2812 : poSRSSource->GetAxis(nullptr, 0, &m_eSourceFirstAxisOrient);
1571 : }
1572 2822 : if (poSRSTarget)
1573 : {
1574 2810 : bTargetLatLong = CPL_TO_BOOL(poSRSTarget->IsGeographic());
1575 2811 : bTargetIsDynamicCRS = poSRSTarget->IsDynamic();
1576 2811 : dfTargetCoordinateEpoch = poSRSTarget->GetCoordinateEpoch();
1577 2811 : if (!bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
1578 : {
1579 2 : bTargetIsDynamicCRS = poSRSTarget->HasPointMotionOperation();
1580 : }
1581 2811 : poSRSTarget->GetAxis(nullptr, 0, &m_eTargetFirstAxisOrient);
1582 : }
1583 :
1584 : #if PROJ_VERSION_MAJOR < 9 || \
1585 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR < 4)
1586 2823 : if (bSourceIsDynamicCRS && bTargetIsDynamicCRS &&
1587 1757 : 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 2823 : bSourceWrap = false;
1601 2823 : dfSourceWrapLong = 0.0;
1602 :
1603 2823 : bTargetWrap = false;
1604 2823 : dfTargetWrapLong = 0.0;
1605 :
1606 : /* -------------------------------------------------------------------- */
1607 : /* Preliminary logic to setup wrapping. */
1608 : /* -------------------------------------------------------------------- */
1609 2823 : 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 2822 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1622 2822 : pszCENTER_LONG =
1623 2821 : poSRSSource ? poSRSSource->GetExtension("GEOGCS", "CENTER_LONG")
1624 : : nullptr;
1625 : }
1626 2822 : if (pszCENTER_LONG != nullptr)
1627 : {
1628 1 : dfSourceWrapLong = CPLAtof(pszCENTER_LONG);
1629 0 : bSourceWrap = true;
1630 0 : CPLDebug("OGRCT", "Wrap source at %g.", dfSourceWrapLong);
1631 : }
1632 2821 : else if (bSourceLatLong && options.d->bHasSourceCenterLong)
1633 : {
1634 265 : dfSourceWrapLong = options.d->dfSourceCenterLong;
1635 265 : bSourceWrap = true;
1636 265 : CPLDebug("OGRCT", "Wrap source at %g.", dfSourceWrapLong);
1637 : }
1638 :
1639 : {
1640 2821 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1641 2822 : pszCENTER_LONG =
1642 2822 : poSRSTarget ? poSRSTarget->GetExtension("GEOGCS", "CENTER_LONG")
1643 : : nullptr;
1644 : }
1645 2822 : if (pszCENTER_LONG != nullptr)
1646 : {
1647 1 : dfTargetWrapLong = CPLAtof(pszCENTER_LONG);
1648 0 : bTargetWrap = true;
1649 0 : CPLDebug("OGRCT", "Wrap target at %g.", dfTargetWrapLong);
1650 : }
1651 2821 : 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 2821 : ComputeThreshold();
1659 :
1660 2822 : DetectWebMercatorToWGS84();
1661 :
1662 : const char *pszCTOpSelection =
1663 2821 : CPLGetConfigOption("OGR_CT_OP_SELECTION", nullptr);
1664 2822 : 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 2822 : if (options.d->dfAccuracy >= 0 || !options.d->bAllowBallpark)
1680 : {
1681 6 : m_eStrategy = Strategy::BEST_ACCURACY;
1682 : }
1683 : }
1684 : #endif
1685 2822 : if (m_eStrategy == Strategy::PROJ)
1686 : {
1687 : const char *pszUseApproxTMERC =
1688 2816 : CPLGetConfigOption("OSR_USE_APPROX_TMERC", nullptr);
1689 2816 : 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 2822 : 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 2802 : 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 2647 : CPLDebug("OGR_CT", "Source CRS: '%s'", pszSrcSRS);
1724 2647 : if (dfSourceCoordinateEpoch > 0)
1725 6 : CPLDebug("OGR_CT", "Source coordinate epoch: %.3f",
1726 : dfSourceCoordinateEpoch);
1727 2647 : CPLDebug("OGR_CT", "Target CRS: '%s'", pszTargetSRS);
1728 2647 : if (dfTargetCoordinateEpoch > 0)
1729 10 : CPLDebug("OGR_CT", "Target coordinate epoch: %.3f",
1730 : dfTargetCoordinateEpoch);
1731 : #endif
1732 :
1733 2647 : if (m_eStrategy == Strategy::PROJ)
1734 : {
1735 2637 : PJ_AREA *area = nullptr;
1736 2637 : if (options.d->bHasAreaOfInterest)
1737 : {
1738 387 : area = proj_area_create();
1739 387 : proj_area_set_bbox(area, options.d->dfWestLongitudeDeg,
1740 387 : options.d->dfSouthLatitudeDeg,
1741 387 : options.d->dfEastLongitudeDeg,
1742 386 : options.d->dfNorthLatitudeDeg);
1743 : }
1744 2636 : 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 2637 : m_pj = proj_create_crs_to_crs(ctx, pszSrcSRS, pszTargetSRS, area);
1794 : #endif
1795 2636 : if (area)
1796 386 : proj_area_destroy(area);
1797 2637 : if (m_pj == nullptr)
1798 : {
1799 20 : CPLError(CE_Failure, CPLE_NotSupported,
1800 : "Cannot find coordinate operations from `%s' to `%s'",
1801 : pszSrcSRS, pszTargetSRS);
1802 20 : 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 5570 : if (options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1824 2775 : (dfSourceCoordinateEpoch == 0 || dfTargetCoordinateEpoch == 0 ||
1825 0 : dfSourceCoordinateEpoch == dfTargetCoordinateEpoch))
1826 : {
1827 : // Determine if we can skip the transformation completely.
1828 2775 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
1829 : nullptr};
1830 2775 : bNoTransform =
1831 5287 : !bSourceWrap && !bTargetWrap &&
1832 2512 : CPL_TO_BOOL(poSRSSource->IsSame(poSRSTarget, apszOptionsIsSame));
1833 : }
1834 :
1835 2795 : 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 4883 : const OGRSpatialReference *OGRProjCT::GetSourceCS() const
2244 :
2245 : {
2246 4883 : return poSRSSource;
2247 : }
2248 :
2249 : /************************************************************************/
2250 : /* GetTargetCS() */
2251 : /************************************************************************/
2252 :
2253 2796 : const OGRSpatialReference *OGRProjCT::GetTargetCS() const
2254 :
2255 : {
2256 2796 : return poSRSTarget;
2257 : }
2258 :
2259 : /************************************************************************/
2260 : /* Transform() */
2261 : /************************************************************************/
2262 :
2263 200882 : int OGRCoordinateTransformation::Transform(size_t nCount, double *x, double *y,
2264 : double *z, int *pabSuccessIn)
2265 :
2266 : {
2267 : int *pabSuccess =
2268 : pabSuccessIn
2269 200882 : ? pabSuccessIn
2270 199296 : : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nCount));
2271 200882 : if (!pabSuccess)
2272 0 : return FALSE;
2273 :
2274 200882 : const int bRet = Transform(nCount, x, y, z, nullptr, pabSuccess);
2275 :
2276 200882 : if (pabSuccess != pabSuccessIn)
2277 199296 : CPLFree(pabSuccess);
2278 :
2279 200882 : 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 2085410 : int OGRProjCT::Transform(size_t nCount, double *x, double *y, double *z,
2333 : double *t, int *pabSuccess)
2334 :
2335 : {
2336 2085410 : const int bRet = TransformWithErrorCodes(nCount, x, y, z, t, pabSuccess);
2337 :
2338 2085420 : if (pabSuccess)
2339 : {
2340 14844600 : for (size_t i = 0; i < nCount; i++)
2341 : {
2342 12759200 : pabSuccess[i] = (pabSuccess[i] == 0);
2343 : }
2344 : }
2345 :
2346 2085420 : 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 2098690 : int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
2360 : double *z, double *t, int *panErrorCodes)
2361 :
2362 : {
2363 2098690 : if (nCount == 0)
2364 0 : return TRUE;
2365 :
2366 : // Prevent any coordinate modification when possible
2367 2098690 : if (bNoTransform)
2368 : {
2369 2028 : if (panErrorCodes)
2370 : {
2371 5450 : for (size_t i = 0; i < nCount; i++)
2372 : {
2373 4022 : panErrorCodes[i] = 0;
2374 : }
2375 : }
2376 2028 : 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 = %d", 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 2096660 : const auto xOriginal = x;
2405 :
2406 2096660 : if (poSRSSource)
2407 : {
2408 2095730 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
2409 2095740 : if (mapping.size() >= 2)
2410 : {
2411 2095740 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
2412 : {
2413 1592370 : std::swap(x, y);
2414 : }
2415 2095710 : const bool bNegateX = mapping[0] < 0;
2416 2095710 : if (bNegateX)
2417 : {
2418 2 : for (size_t i = 0; i < nCount; i++)
2419 : {
2420 1 : x[i] = -x[i];
2421 : }
2422 : }
2423 2095710 : const bool bNegateY = mapping[1] < 0;
2424 2095720 : if (bNegateY)
2425 : {
2426 0 : for (size_t i = 0; i < nCount; i++)
2427 : {
2428 0 : y[i] = -y[i];
2429 : }
2430 : }
2431 2095720 : 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 2096650 : if (bSourceLatLong && bSourceWrap)
2445 : {
2446 203772 : if (m_eSourceFirstAxisOrient == OAO_East)
2447 : {
2448 11402 : for (size_t i = 0; i < nCount; i++)
2449 : {
2450 10311 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2451 : {
2452 10311 : if (x[i] < dfSourceWrapLong - 180.0)
2453 0 : x[i] += 360.0;
2454 10311 : else if (x[i] > dfSourceWrapLong + 180)
2455 0 : x[i] -= 360.0;
2456 : }
2457 : }
2458 : }
2459 : else
2460 : {
2461 632657 : for (size_t i = 0; i < nCount; i++)
2462 : {
2463 429976 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2464 : {
2465 428646 : if (y[i] < dfSourceWrapLong - 180.0)
2466 35 : y[i] += 360.0;
2467 428611 : 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 2096650 : bool bTransformDone = false;
2478 2096650 : int bRet = TRUE;
2479 2096650 : if (bWebMercatorToWGS84LongLat)
2480 : {
2481 158488 : constexpr double REVERSE_SPHERE_RADIUS = 1.0 / 6378137.0;
2482 :
2483 158488 : if (m_eSourceFirstAxisOrient != OAO_East)
2484 : {
2485 0 : std::swap(x, y);
2486 : }
2487 :
2488 158486 : double y0 = y[0];
2489 1419820 : for (size_t i = 0; i < nCount; i++)
2490 : {
2491 1261330 : if (x[i] == HUGE_VAL)
2492 : {
2493 0 : bRet = FALSE;
2494 : }
2495 : else
2496 : {
2497 1261330 : x[i] = x[i] * REVERSE_SPHERE_RADIUS;
2498 1261330 : if (x[i] > M_PI)
2499 : {
2500 711 : if (x[i] < M_PI + 1e-14)
2501 : {
2502 138 : x[i] = M_PI;
2503 : }
2504 573 : else if (m_options.d->bCheckWithInvertProj)
2505 : {
2506 0 : x[i] = HUGE_VAL;
2507 0 : y[i] = HUGE_VAL;
2508 0 : y0 = HUGE_VAL;
2509 0 : continue;
2510 : }
2511 : else
2512 : {
2513 2 : do
2514 : {
2515 575 : x[i] -= 2 * M_PI;
2516 575 : } while (x[i] > M_PI);
2517 : }
2518 : }
2519 1260620 : else if (x[i] < -M_PI)
2520 : {
2521 835 : if (x[i] > -M_PI - 1e-14)
2522 : {
2523 8 : x[i] = -M_PI;
2524 : }
2525 827 : 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 0 : do
2535 : {
2536 830 : x[i] += 2 * M_PI;
2537 830 : } while (x[i] < -M_PI);
2538 : }
2539 : }
2540 1261330 : constexpr double RAD_TO_DEG = 180. / M_PI;
2541 1261330 : x[i] *= RAD_TO_DEG;
2542 :
2543 : // Optimization for the case where we are provided a whole line
2544 : // of same northing.
2545 1261330 : if (i > 0 && y[i] == y0)
2546 673122 : y[i] = y[0];
2547 : else
2548 : {
2549 588211 : y[i] = M_PI / 2.0 -
2550 588211 : 2.0 * atan(exp(-y[i] * REVERSE_SPHERE_RADIUS));
2551 588211 : y[i] *= RAD_TO_DEG;
2552 : }
2553 : }
2554 : }
2555 :
2556 158491 : if (panErrorCodes)
2557 : {
2558 1419880 : for (size_t i = 0; i < nCount; i++)
2559 : {
2560 1261400 : if (x[i] != HUGE_VAL)
2561 1261340 : panErrorCodes[i] = 0;
2562 : else
2563 57 : panErrorCodes[i] =
2564 : PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2565 : }
2566 : }
2567 :
2568 158491 : if (m_eTargetFirstAxisOrient != OAO_East)
2569 : {
2570 157490 : std::swap(x, y);
2571 : }
2572 :
2573 158460 : bTransformDone = true;
2574 : }
2575 :
2576 : // Determine the default coordinate epoch, if not provided in the point to
2577 : // transform.
2578 : // For time-dependent transformations, PROJ can currently only do
2579 : // staticCRS -> dynamicCRS or dynamicCRS -> staticCRS transformations, and
2580 : // in either case, the coordinate epoch of the dynamicCRS must be provided
2581 : // as the input time.
2582 2096620 : double dfDefaultTime = HUGE_VAL;
2583 2026070 : if (bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0 &&
2584 4122720 : !bTargetIsDynamicCRS &&
2585 4 : CPLTestBool(
2586 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2587 : {
2588 4 : dfDefaultTime = dfSourceCoordinateEpoch;
2589 4 : CPLDebug("OGR_CT", "Using coordinate epoch %f from source CRS",
2590 : dfDefaultTime);
2591 : }
2592 1837290 : else if (bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0 &&
2593 3933930 : !bSourceIsDynamicCRS &&
2594 2 : CPLTestBool(
2595 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2596 : {
2597 2 : dfDefaultTime = dfTargetCoordinateEpoch;
2598 2 : CPLDebug("OGR_CT", "Using coordinate epoch %f from target CRS",
2599 : dfDefaultTime);
2600 : }
2601 :
2602 : /* -------------------------------------------------------------------- */
2603 : /* Select dynamically the best transformation for the data, if */
2604 : /* needed. */
2605 : /* -------------------------------------------------------------------- */
2606 2096640 : auto ctx = OSRGetProjTLSContext();
2607 :
2608 2096680 : PJ *pj = m_pj;
2609 2096680 : if (!bTransformDone && !pj)
2610 : {
2611 0 : double avgX = 0.0;
2612 0 : double avgY = 0.0;
2613 0 : size_t nCountValid = 0;
2614 0 : for (size_t i = 0; i < nCount; i++)
2615 : {
2616 0 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2617 : {
2618 0 : avgX += x[i];
2619 0 : avgY += y[i];
2620 0 : nCountValid++;
2621 : }
2622 : }
2623 0 : if (nCountValid != 0)
2624 : {
2625 0 : avgX /= static_cast<double>(nCountValid);
2626 0 : avgY /= static_cast<double>(nCountValid);
2627 : }
2628 :
2629 0 : constexpr int N_MAX_RETRY = 2;
2630 0 : int iExcluded[N_MAX_RETRY] = {-1, -1};
2631 :
2632 0 : const int nOperations = static_cast<int>(m_oTransformations.size());
2633 : PJ_COORD coord;
2634 0 : coord.xyzt.x = avgX;
2635 0 : coord.xyzt.y = avgY;
2636 0 : coord.xyzt.z = z ? z[0] : 0;
2637 0 : coord.xyzt.t = t ? t[0] : dfDefaultTime;
2638 :
2639 : // We may need several attempts. For example the point at
2640 : // lon=-111.5 lat=45.26 falls into the bounding box of the Canadian
2641 : // ntv2_0.gsb grid, except that it is not in any of the subgrids, being
2642 : // in the US. We thus need another retry that will select the conus
2643 : // grid.
2644 0 : for (int iRetry = 0; iRetry <= N_MAX_RETRY; iRetry++)
2645 : {
2646 0 : int iBestTransf = -1;
2647 : // Select transform whose BBOX match our data and has the best
2648 : // accuracy if m_eStrategy == BEST_ACCURACY. Or just the first BBOX
2649 : // matching one, if
2650 : // m_eStrategy == FIRST_MATCHING
2651 0 : double dfBestAccuracy = std::numeric_limits<double>::infinity();
2652 0 : for (int i = 0; i < nOperations; i++)
2653 : {
2654 0 : if (i == iExcluded[0] || i == iExcluded[1])
2655 : {
2656 0 : continue;
2657 : }
2658 0 : const auto &transf = m_oTransformations[i];
2659 0 : if (avgX >= transf.minx && avgX <= transf.maxx &&
2660 0 : avgY >= transf.miny && avgY <= transf.maxy &&
2661 0 : (iBestTransf < 0 || (transf.accuracy >= 0 &&
2662 0 : transf.accuracy < dfBestAccuracy)))
2663 : {
2664 0 : iBestTransf = i;
2665 0 : dfBestAccuracy = transf.accuracy;
2666 0 : if (m_eStrategy == Strategy::FIRST_MATCHING)
2667 0 : break;
2668 : }
2669 : }
2670 0 : if (iBestTransf < 0)
2671 : {
2672 0 : break;
2673 : }
2674 0 : auto &transf = m_oTransformations[iBestTransf];
2675 0 : pj = transf.pj;
2676 0 : proj_assign_context(pj, ctx);
2677 0 : if (iBestTransf != m_iCurTransformation)
2678 : {
2679 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2680 : transf.osProjString.c_str(), transf.osName.c_str());
2681 0 : m_iCurTransformation = iBestTransf;
2682 : }
2683 :
2684 0 : auto res = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2685 0 : if (res.xyzt.x != HUGE_VAL)
2686 : {
2687 0 : break;
2688 : }
2689 0 : pj = nullptr;
2690 0 : CPLDebug("OGRCT", "Did not result in valid result. "
2691 : "Attempting a retry with another operation.");
2692 0 : if (iRetry == N_MAX_RETRY)
2693 : {
2694 0 : break;
2695 : }
2696 0 : iExcluded[iRetry] = iBestTransf;
2697 : }
2698 :
2699 0 : if (!pj)
2700 : {
2701 : // In case we did not find an operation whose area of use is
2702 : // compatible with the input coordinate, then goes through again the
2703 : // list, and use the first operation that does not require grids.
2704 0 : for (int i = 0; i < nOperations; i++)
2705 : {
2706 0 : auto &transf = m_oTransformations[i];
2707 0 : if (proj_coordoperation_get_grid_used_count(ctx, transf.pj) ==
2708 : 0)
2709 : {
2710 0 : pj = transf.pj;
2711 0 : proj_assign_context(pj, ctx);
2712 0 : if (i != m_iCurTransformation)
2713 : {
2714 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2715 : transf.osProjString.c_str(),
2716 : transf.osName.c_str());
2717 0 : m_iCurTransformation = i;
2718 : }
2719 0 : break;
2720 : }
2721 : }
2722 : }
2723 :
2724 0 : if (!pj)
2725 : {
2726 0 : if (m_bEmitErrors && ++nErrorCount < 20)
2727 : {
2728 0 : CPLError(CE_Failure, CPLE_AppDefined,
2729 : "Cannot find transformation for provided coordinates");
2730 : }
2731 0 : else if (nErrorCount == 20)
2732 : {
2733 0 : CPLError(CE_Failure, CPLE_AppDefined,
2734 : "Reprojection failed, further errors will be "
2735 : "suppressed on the transform object.");
2736 : }
2737 :
2738 0 : for (size_t i = 0; i < nCount; i++)
2739 : {
2740 0 : x[i] = HUGE_VAL;
2741 0 : y[i] = HUGE_VAL;
2742 0 : if (panErrorCodes)
2743 0 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_NO_OPERATION;
2744 : }
2745 0 : return FALSE;
2746 : }
2747 : }
2748 2096680 : if (pj)
2749 : {
2750 2095800 : proj_assign_context(pj, ctx);
2751 : }
2752 :
2753 : /* -------------------------------------------------------------------- */
2754 : /* Do the transformation (or not...) using PROJ */
2755 : /* -------------------------------------------------------------------- */
2756 :
2757 2096520 : if (!bTransformDone)
2758 : {
2759 1938180 : const auto nLastErrorCounter = CPLGetErrorCounter();
2760 :
2761 13461300 : for (size_t i = 0; i < nCount; i++)
2762 : {
2763 : PJ_COORD coord;
2764 11523200 : const double xIn = x[i];
2765 11523200 : const double yIn = y[i];
2766 11523200 : if (!std::isfinite(xIn))
2767 : {
2768 10681 : bRet = FALSE;
2769 10681 : x[i] = HUGE_VAL;
2770 10681 : y[i] = HUGE_VAL;
2771 10681 : if (panErrorCodes)
2772 10681 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_INVALID_COORD;
2773 10681 : continue;
2774 : }
2775 11512500 : coord.xyzt.x = x[i];
2776 11512500 : coord.xyzt.y = y[i];
2777 11512500 : coord.xyzt.z = z ? z[i] : 0;
2778 11512500 : coord.xyzt.t = t ? t[i] : dfDefaultTime;
2779 11512500 : proj_errno_reset(pj);
2780 11512500 : coord = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2781 : #if 0
2782 : CPLDebug("OGRCT",
2783 : "Transforming (x=%f,y=%f,z=%f,time=%f) to "
2784 : "(x=%f,y=%f,z=%f,time=%f)",
2785 : x[i], y[i], z ? z[i] : 0, t ? t[i] : dfDefaultTime,
2786 : coord.xyzt.x, coord.xyzt.y, coord.xyzt.z, coord.xyzt.t);
2787 : #endif
2788 11512500 : x[i] = coord.xyzt.x;
2789 11512500 : y[i] = coord.xyzt.y;
2790 11512500 : if (z)
2791 11474100 : z[i] = coord.xyzt.z;
2792 11512500 : if (t)
2793 25 : t[i] = coord.xyzt.t;
2794 11512500 : int err = 0;
2795 11512500 : if (std::isnan(coord.xyzt.x))
2796 : {
2797 141 : bRet = FALSE;
2798 : // This shouldn't normally happen if PROJ projections behave
2799 : // correctly, but e.g inverse laea before PROJ 8.1.1 could
2800 : // do that for points out of domain.
2801 : // See https://github.com/OSGeo/PROJ/pull/2800
2802 141 : x[i] = HUGE_VAL;
2803 141 : y[i] = HUGE_VAL;
2804 141 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2805 :
2806 : #ifdef DEBUG
2807 141 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
2808 : "PROJ returned a NaN value. It should be fixed");
2809 : #else
2810 : CPLDebugOnce("OGR_CT",
2811 : "PROJ returned a NaN value. It should be fixed");
2812 : #endif
2813 : }
2814 11512300 : else if (coord.xyzt.x == HUGE_VAL)
2815 : {
2816 3388880 : bRet = FALSE;
2817 3388880 : err = proj_errno(pj);
2818 : // PROJ should normally emit an error, but in case it does not
2819 : // (e.g PROJ 6.3 with the +ortho projection), synthesize one
2820 3388880 : if (err == 0)
2821 3363770 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2822 : }
2823 : else
2824 : {
2825 8123460 : if (m_recordDifferentOperationsUsed &&
2826 60989 : !m_differentOperationsUsed)
2827 : {
2828 : #if PROJ_VERSION_MAJOR > 9 || \
2829 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 1)
2830 :
2831 : PJ *lastOp = proj_trans_get_last_used_operation(pj);
2832 : if (lastOp)
2833 : {
2834 : const char *projString = proj_as_proj_string(
2835 : ctx, lastOp, PJ_PROJ_5, nullptr);
2836 : if (projString)
2837 : {
2838 : if (m_lastPjUsedPROJString.empty())
2839 : {
2840 : m_lastPjUsedPROJString = projString;
2841 : }
2842 : else if (m_lastPjUsedPROJString != projString)
2843 : {
2844 : m_differentOperationsUsed = true;
2845 : }
2846 : }
2847 : proj_destroy(lastOp);
2848 : }
2849 : #endif
2850 : }
2851 :
2852 8123460 : if (m_options.d->bCheckWithInvertProj)
2853 : {
2854 : // For some projections, we cannot detect if we are trying to
2855 : // reproject coordinates outside the validity area of the
2856 : // projection. So let's do the reverse reprojection and compare
2857 : // with the source coordinates.
2858 : coord =
2859 1836040 : proj_trans(pj, m_bReversePj ? PJ_FWD : PJ_INV, coord);
2860 1836040 : if (fabs(coord.xyzt.x - xIn) > dfThreshold ||
2861 1684050 : fabs(coord.xyzt.y - yIn) > dfThreshold)
2862 : {
2863 154814 : bRet = FALSE;
2864 154814 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2865 154814 : x[i] = HUGE_VAL;
2866 154814 : y[i] = HUGE_VAL;
2867 : }
2868 : }
2869 : }
2870 :
2871 11512500 : if (panErrorCodes)
2872 11511000 : panErrorCodes[i] = err;
2873 :
2874 : /* --------------------------------------------------------------------
2875 : */
2876 : /* Try to report an error through CPL. Get proj error string
2877 : */
2878 : /* if possible. Try to avoid reporting thousands of errors. */
2879 : /* Suppress further error reporting on this OGRProjCT if we */
2880 : /* have already reported 20 errors. */
2881 : /* --------------------------------------------------------------------
2882 : */
2883 11512500 : if (err != 0)
2884 : {
2885 3543830 : if (++nErrorCount < 20)
2886 : {
2887 : #if PROJ_VERSION_MAJOR >= 8
2888 : const char *pszError = proj_context_errno_string(ctx, err);
2889 : #else
2890 3990 : const char *pszError = proj_errno_string(err);
2891 : #endif
2892 3990 : if (m_bEmitErrors
2893 : #ifdef PROJ_ERR_OTHER_NO_INVERSE_OP
2894 : || (i == 0 && err == PROJ_ERR_OTHER_NO_INVERSE_OP)
2895 : #endif
2896 : )
2897 : {
2898 493 : if (nLastErrorCounter != CPLGetErrorCounter() &&
2899 755 : CPLGetLastErrorType() == CE_Failure &&
2900 262 : strstr(CPLGetLastErrorMsg(), "PROJ:"))
2901 : {
2902 : // do nothing
2903 : }
2904 493 : else if (pszError == nullptr)
2905 0 : CPLError(CE_Failure, CPLE_AppDefined,
2906 : "Reprojection failed, err = %d", err);
2907 : else
2908 493 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
2909 : pszError);
2910 : }
2911 : else
2912 : {
2913 3497 : if (pszError == nullptr)
2914 0 : CPLDebug("OGRCT", "Reprojection failed, err = %d",
2915 : err);
2916 : else
2917 3497 : CPLDebug("OGRCT", "%s", pszError);
2918 : }
2919 : }
2920 3539840 : else if (nErrorCount == 20)
2921 : {
2922 182 : if (m_bEmitErrors)
2923 : {
2924 14 : CPLError(CE_Failure, CPLE_AppDefined,
2925 : "Reprojection failed, err = %d, further "
2926 : "errors will be "
2927 : "suppressed on the transform object.",
2928 : err);
2929 : }
2930 : else
2931 : {
2932 168 : CPLDebug("OGRCT",
2933 : "Reprojection failed, err = %d, further "
2934 : "errors will be "
2935 : "suppressed on the transform object.",
2936 : err);
2937 : }
2938 : }
2939 : }
2940 : }
2941 : }
2942 :
2943 : /* -------------------------------------------------------------------- */
2944 : /* Potentially do longitude wrapping. */
2945 : /* -------------------------------------------------------------------- */
2946 2096480 : if (bTargetLatLong && bTargetWrap)
2947 : {
2948 673357 : if (m_eTargetFirstAxisOrient == OAO_East)
2949 : {
2950 17007 : for (size_t i = 0; i < nCount; i++)
2951 : {
2952 13898 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2953 : {
2954 13898 : if (x[i] < dfTargetWrapLong - 180.0)
2955 0 : x[i] += 360.0;
2956 13898 : else if (x[i] > dfTargetWrapLong + 180)
2957 154 : x[i] -= 360.0;
2958 : }
2959 : }
2960 : }
2961 : else
2962 : {
2963 9175960 : for (size_t i = 0; i < nCount; i++)
2964 : {
2965 8505710 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2966 : {
2967 5059750 : if (y[i] < dfTargetWrapLong - 180.0)
2968 5849 : y[i] += 360.0;
2969 5053900 : else if (y[i] > dfTargetWrapLong + 180)
2970 27852 : y[i] -= 360.0;
2971 : }
2972 : }
2973 : }
2974 : }
2975 :
2976 : /* -------------------------------------------------------------------- */
2977 : /* Apply data axis to target CRS mapping. */
2978 : /* -------------------------------------------------------------------- */
2979 2096480 : if (poSRSTarget)
2980 : {
2981 2095730 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
2982 2095740 : if (mapping.size() >= 2)
2983 : {
2984 2095740 : const bool bNegateX = mapping[0] < 0;
2985 2095730 : if (bNegateX)
2986 : {
2987 2 : for (size_t i = 0; i < nCount; i++)
2988 : {
2989 1 : x[i] = -x[i];
2990 : }
2991 : }
2992 2095730 : const bool bNegateY = mapping[1] < 0;
2993 2095730 : if (bNegateY)
2994 : {
2995 0 : for (size_t i = 0; i < nCount; i++)
2996 : {
2997 0 : y[i] = -y[i];
2998 : }
2999 : }
3000 2095730 : if (z && mapping.size() >= 3 && mapping[2] == -3)
3001 : {
3002 2 : for (size_t i = 0; i < nCount; i++)
3003 : {
3004 1 : z[i] = -z[i];
3005 : }
3006 : }
3007 :
3008 2095730 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
3009 : {
3010 893535 : std::swap(x, y);
3011 : }
3012 : }
3013 : }
3014 :
3015 : // Check whether final "genuine" axis swap is really necessary
3016 2096490 : if (x != xOriginal)
3017 : {
3018 1317560 : std::swap_ranges(x, x + nCount, y);
3019 : }
3020 :
3021 : #ifdef DEBUG_VERBOSE
3022 : if (bDebugCT)
3023 : {
3024 : CPLDebug("OGRCT", "Out:");
3025 : for (size_t i = 0; i < nCount; ++i)
3026 : {
3027 : CPLDebug("OGRCT", " x[%d] = %.16g y[%d] = %.16g", int(i), x[i],
3028 : int(i), y[i]);
3029 : }
3030 : }
3031 : #endif
3032 : #ifdef DEBUG_PERF
3033 : struct CPLTimeVal tvEnd;
3034 : CPLGettimeofday(&tvEnd, nullptr);
3035 : const double delay = (tvEnd.tv_sec + tvEnd.tv_usec * 1e-6) -
3036 : (tvStart.tv_sec + tvStart.tv_usec * 1e-6);
3037 : g_dfTotalTimeReprojection += delay;
3038 : // CPLDebug("OGR_CT", "End TransformWithErrorCodes(): %d ms",
3039 : // static_cast<int>(delay * 1000));
3040 : #endif
3041 :
3042 2096660 : return bRet;
3043 : }
3044 :
3045 : /************************************************************************/
3046 : /* TransformBounds() */
3047 : /************************************************************************/
3048 :
3049 : // ---------------------------------------------------------------------------
3050 801 : static double simple_min(const double *data, const int *panErrorCodes,
3051 : const int arr_len)
3052 : {
3053 801 : double min_value = HUGE_VAL;
3054 25357 : for (int iii = 0; iii < arr_len; iii++)
3055 : {
3056 24556 : if ((data[iii] < min_value || min_value == HUGE_VAL) &&
3057 5945 : panErrorCodes[iii] == 0)
3058 5822 : min_value = data[iii];
3059 : }
3060 801 : return min_value;
3061 : }
3062 :
3063 : // ---------------------------------------------------------------------------
3064 801 : static double simple_max(const double *data, const int *panErrorCodes,
3065 : const int arr_len)
3066 : {
3067 801 : double max_value = HUGE_VAL;
3068 25281 : for (int iii = 0; iii < arr_len; iii++)
3069 : {
3070 24480 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
3071 5163 : panErrorCodes[iii] == 0)
3072 3104 : max_value = data[iii];
3073 : }
3074 801 : return max_value;
3075 : }
3076 :
3077 : // ---------------------------------------------------------------------------
3078 18234 : static int _find_previous_index(const int iii, const int *panErrorCodes,
3079 : const int arr_len)
3080 : {
3081 : // find index of nearest valid previous value if exists
3082 18234 : int prev_iii = iii - 1;
3083 18234 : if (prev_iii == -1) // handle wraparound
3084 396 : prev_iii = arr_len - 1;
3085 18256 : while (panErrorCodes[prev_iii] != 0 && prev_iii != iii)
3086 : {
3087 22 : prev_iii--;
3088 22 : if (prev_iii == -1) // handle wraparound
3089 2 : prev_iii = arr_len - 1;
3090 : }
3091 18234 : return prev_iii;
3092 : }
3093 :
3094 : // ---------------------------------------------------------------------------
3095 : /******************************************************************************
3096 : Handles the case when longitude values cross the antimeridian
3097 : when calculating the minimum.
3098 : Note: The data array must be in a linear ring.
3099 : Note: This requires a densified ring with at least 2 additional
3100 : points per edge to correctly handle global extents.
3101 : If only 1 additional point:
3102 : | |
3103 : |RL--x0--|RL--
3104 : | |
3105 : -180 180|-180
3106 : If they are evenly spaced and it crosses the antimeridian:
3107 : x0 - L = 180
3108 : R - x0 = -180
3109 : For example:
3110 : Let R = -179.9, x0 = 0.1, L = -179.89
3111 : x0 - L = 0.1 - -179.9 = 180
3112 : R - x0 = -179.89 - 0.1 ~= -180
3113 : This is the same in the case when it didn't cross the antimeridian.
3114 : If you have 2 additional points:
3115 : | |
3116 : |RL--x0--x1--|RL--
3117 : | |
3118 : -180 180|-180
3119 : If they are evenly spaced and it crosses the antimeridian:
3120 : x0 - L = 120
3121 : x1 - x0 = 120
3122 : R - x1 = -240
3123 : For example:
3124 : Let R = -179.9, x0 = -59.9, x1 = 60.1 L = -179.89
3125 : x0 - L = 59.9 - -179.9 = 120
3126 : x1 - x0 = 60.1 - 59.9 = 120
3127 : R - x1 = -179.89 - 60.1 ~= -240
3128 : However, if they are evenly spaced and it didn't cross the antimeridian:
3129 : x0 - L = 120
3130 : x1 - x0 = 120
3131 : R - x1 = 120
3132 : From this, we have a delta that is guaranteed to be significantly
3133 : large enough to tell the difference reguarless of the direction
3134 : the antimeridian was crossed.
3135 : However, even though the spacing was even in the source projection, it isn't
3136 : guaranteed in the target geographic projection. So, instead of 240, 200 is used
3137 : as it significantly larger than 120 to be sure that the antimeridian was crossed
3138 : but smalller than 240 to account for possible irregularities in distances
3139 : when re-projecting. Also, 200 ensures latitudes are ignored for axis order
3140 : handling.
3141 : ******************************************************************************/
3142 199 : static double antimeridian_min(const double *data, const int *panErrorCodes,
3143 : const int arr_len)
3144 : {
3145 199 : double positive_min = HUGE_VAL;
3146 199 : double min_value = HUGE_VAL;
3147 199 : int crossed_meridian_count = 0;
3148 199 : bool positive_meridian = false;
3149 :
3150 9327 : for (int iii = 0; iii < arr_len; iii++)
3151 : {
3152 9128 : if (panErrorCodes[iii])
3153 11 : continue;
3154 9117 : int prev_iii = _find_previous_index(iii, panErrorCodes, arr_len);
3155 : // check if crossed meridian
3156 9117 : double delta = data[prev_iii] - data[iii];
3157 : // 180 -> -180
3158 9117 : if (delta >= 200 && delta != HUGE_VAL)
3159 : {
3160 9 : if (crossed_meridian_count == 0)
3161 2 : positive_min = min_value;
3162 9 : crossed_meridian_count++;
3163 9 : positive_meridian = false;
3164 : // -180 -> 180
3165 : }
3166 9108 : else if (delta <= -200 && delta != HUGE_VAL)
3167 : {
3168 9 : if (crossed_meridian_count == 0)
3169 6 : positive_min = data[iii];
3170 9 : crossed_meridian_count++;
3171 9 : positive_meridian = true;
3172 : }
3173 : // positive meridian side min
3174 9117 : if (positive_meridian && data[iii] < positive_min)
3175 38 : positive_min = data[iii];
3176 : // track general min value
3177 9117 : if (data[iii] < min_value)
3178 317 : min_value = data[iii];
3179 : }
3180 :
3181 199 : if (crossed_meridian_count == 2)
3182 7 : return positive_min;
3183 192 : else if (crossed_meridian_count == 4)
3184 : // bounds extends beyond -180/180
3185 1 : return -180;
3186 191 : return min_value;
3187 : }
3188 :
3189 : // ---------------------------------------------------------------------------
3190 : // Handles the case when longitude values cross the antimeridian
3191 : // when calculating the minimum.
3192 : // Note: The data array must be in a linear ring.
3193 : // Note: This requires a densified ring with at least 2 additional
3194 : // points per edge to correctly handle global extents.
3195 : // See antimeridian_min docstring for reasoning.
3196 199 : static double antimeridian_max(const double *data, const int *panErrorCodes,
3197 : const int arr_len)
3198 : {
3199 199 : double negative_max = -HUGE_VAL;
3200 199 : double max_value = -HUGE_VAL;
3201 199 : bool negative_meridian = false;
3202 199 : int crossed_meridian_count = 0;
3203 :
3204 9327 : for (int iii = 0; iii < arr_len; iii++)
3205 : {
3206 9128 : if (panErrorCodes[iii])
3207 11 : continue;
3208 9117 : int prev_iii = _find_previous_index(iii, panErrorCodes, arr_len);
3209 : // check if crossed meridian
3210 9117 : double delta = data[prev_iii] - data[iii];
3211 : // 180 -> -180
3212 9117 : if (delta >= 200 && delta != HUGE_VAL)
3213 : {
3214 9 : if (crossed_meridian_count == 0)
3215 2 : negative_max = data[iii];
3216 9 : crossed_meridian_count++;
3217 9 : negative_meridian = true;
3218 : // -180 -> 180
3219 : }
3220 9108 : else if (delta <= -200 && delta != HUGE_VAL)
3221 : {
3222 9 : if (crossed_meridian_count == 0)
3223 6 : negative_max = max_value;
3224 9 : negative_meridian = false;
3225 9 : crossed_meridian_count++;
3226 : }
3227 : // negative meridian side max
3228 9117 : if (negative_meridian &&
3229 142 : (data[iii] > negative_max || negative_max == HUGE_VAL) &&
3230 45 : panErrorCodes[iii] == 0)
3231 45 : negative_max = data[iii];
3232 : // track general max value
3233 9117 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
3234 2753 : panErrorCodes[iii] == 0)
3235 2753 : max_value = data[iii];
3236 : }
3237 199 : if (crossed_meridian_count == 2)
3238 7 : return negative_max;
3239 192 : else if (crossed_meridian_count == 4)
3240 : // bounds extends beyond -180/180
3241 1 : return 180;
3242 191 : return max_value;
3243 : }
3244 :
3245 : // ---------------------------------------------------------------------------
3246 : // Check if the original projected bounds contains
3247 : // the north pole.
3248 : // This assumes that the destination CRS is geographic.
3249 211 : bool OGRProjCT::ContainsNorthPole(const double xmin, const double ymin,
3250 : const double xmax, const double ymax,
3251 : bool lon_lat_order)
3252 : {
3253 211 : double pole_y = 90;
3254 211 : double pole_x = 0;
3255 211 : if (!lon_lat_order)
3256 : {
3257 4 : pole_y = 0;
3258 4 : pole_x = 90;
3259 : }
3260 422 : auto inverseCT = std::unique_ptr<OGRCoordinateTransformation>(GetInverse());
3261 211 : if (!inverseCT)
3262 0 : return false;
3263 211 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3264 211 : const bool success = inverseCT->TransformWithErrorCodes(
3265 211 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr);
3266 215 : return success && xmin < pole_x && pole_x < xmax && ymax > pole_y &&
3267 215 : pole_y > ymin;
3268 : }
3269 :
3270 : // ---------------------------------------------------------------------------
3271 : // Check if the original projected bounds contains
3272 : // the south pole.
3273 : // This assumes that the destination CRS is geographic.
3274 211 : bool OGRProjCT::ContainsSouthPole(const double xmin, const double ymin,
3275 : const double xmax, const double ymax,
3276 : bool lon_lat_order)
3277 : {
3278 211 : double pole_y = -90;
3279 211 : double pole_x = 0;
3280 211 : if (!lon_lat_order)
3281 : {
3282 4 : pole_y = 0;
3283 4 : pole_x = -90;
3284 : }
3285 422 : auto inverseCT = std::unique_ptr<OGRCoordinateTransformation>(GetInverse());
3286 211 : if (!inverseCT)
3287 0 : return false;
3288 211 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3289 211 : const bool success = inverseCT->TransformWithErrorCodes(
3290 211 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr);
3291 253 : return success && xmin < pole_x && pole_x < xmax && ymax > pole_y &&
3292 253 : pole_y > ymin;
3293 : }
3294 :
3295 544 : int OGRProjCT::TransformBounds(const double xmin, const double ymin,
3296 : const double xmax, const double ymax,
3297 : double *out_xmin, double *out_ymin,
3298 : double *out_xmax, double *out_ymax,
3299 : const int densify_pts)
3300 : {
3301 544 : if (bNoTransform)
3302 : {
3303 29 : *out_xmin = xmin;
3304 29 : *out_ymin = ymin;
3305 29 : *out_xmax = xmax;
3306 29 : *out_ymax = ymax;
3307 29 : return true;
3308 : }
3309 :
3310 515 : *out_xmin = HUGE_VAL;
3311 515 : *out_ymin = HUGE_VAL;
3312 515 : *out_xmax = HUGE_VAL;
3313 515 : *out_ymax = HUGE_VAL;
3314 :
3315 515 : if (densify_pts < 0 || densify_pts > 10000)
3316 : {
3317 1 : CPLError(CE_Failure, CPLE_AppDefined,
3318 : "densify_pts must be between 0-10000.");
3319 1 : return false;
3320 : }
3321 514 : if (!poSRSSource)
3322 : {
3323 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing source SRS.");
3324 0 : return false;
3325 : }
3326 514 : if (!poSRSTarget)
3327 : {
3328 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing target SRS.");
3329 0 : return false;
3330 : }
3331 :
3332 514 : bool degree_input = false;
3333 514 : bool degree_output = false;
3334 514 : bool input_lon_lat_order = false;
3335 514 : bool output_lon_lat_order = false;
3336 :
3337 514 : if (bSourceLatLong)
3338 : {
3339 260 : degree_input = fabs(poSRSSource->GetAngularUnits(nullptr) -
3340 260 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3341 260 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
3342 476 : if ((mapping[0] == 1 && m_eSourceFirstAxisOrient == OAO_East) ||
3343 216 : (mapping[0] == 2 && m_eSourceFirstAxisOrient != OAO_East))
3344 : {
3345 251 : input_lon_lat_order = true;
3346 : }
3347 : }
3348 514 : if (bTargetLatLong)
3349 : {
3350 212 : degree_output = fabs(poSRSTarget->GetAngularUnits(nullptr) -
3351 212 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3352 212 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
3353 394 : if ((mapping[0] == 1 && m_eTargetFirstAxisOrient == OAO_East) ||
3354 182 : (mapping[0] == 2 && m_eTargetFirstAxisOrient != OAO_East))
3355 : {
3356 207 : output_lon_lat_order = true;
3357 : }
3358 : }
3359 :
3360 514 : if (degree_output && densify_pts < 2)
3361 : {
3362 1 : CPLError(CE_Failure, CPLE_AppDefined,
3363 : "densify_pts must be at least 2 if the output is geograpic.");
3364 1 : return false;
3365 : }
3366 :
3367 513 : int side_pts = densify_pts + 1; // add one because we are densifying
3368 513 : const int boundary_len = side_pts * 4;
3369 1026 : std::vector<double> x_boundary_array;
3370 1026 : std::vector<double> y_boundary_array;
3371 1026 : std::vector<int> anErrorCodes;
3372 : try
3373 : {
3374 513 : x_boundary_array.resize(boundary_len);
3375 513 : y_boundary_array.resize(boundary_len);
3376 513 : anErrorCodes.resize(boundary_len);
3377 : }
3378 0 : catch (const std::exception &e) // memory allocation failure
3379 : {
3380 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
3381 0 : return false;
3382 : }
3383 513 : double delta_x = 0;
3384 513 : double delta_y = 0;
3385 513 : bool north_pole_in_bounds = false;
3386 513 : bool south_pole_in_bounds = false;
3387 513 : if (degree_output)
3388 : {
3389 : north_pole_in_bounds =
3390 211 : ContainsNorthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3391 : south_pole_in_bounds =
3392 211 : ContainsSouthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3393 : }
3394 :
3395 513 : if (degree_input && xmax < xmin)
3396 : {
3397 1 : if (!input_lon_lat_order)
3398 : {
3399 0 : CPLError(CE_Failure, CPLE_AppDefined,
3400 : "latitude max < latitude min.");
3401 0 : return false;
3402 : }
3403 : // handle antimeridian
3404 1 : delta_x = (xmax - xmin + 360.0) / side_pts;
3405 : }
3406 : else
3407 : {
3408 512 : delta_x = (xmax - xmin) / side_pts;
3409 : }
3410 513 : if (degree_input && ymax < ymin)
3411 : {
3412 1 : if (input_lon_lat_order)
3413 : {
3414 0 : CPLError(CE_Failure, CPLE_AppDefined,
3415 : "latitude max < latitude min.");
3416 0 : return false;
3417 : }
3418 : // handle antimeridian
3419 1 : delta_y = (ymax - ymin + 360.0) / side_pts;
3420 : }
3421 : else
3422 : {
3423 512 : delta_y = (ymax - ymin) / side_pts;
3424 : }
3425 :
3426 : // build densified bounding box
3427 : // Note: must be a linear ring for antimeridian logic
3428 4829 : for (int iii = 0; iii < side_pts; iii++)
3429 : {
3430 : // xmin boundary
3431 4316 : y_boundary_array[iii] = ymax - iii * delta_y;
3432 4316 : x_boundary_array[iii] = xmin;
3433 : // ymin boundary
3434 4316 : y_boundary_array[iii + side_pts] = ymin;
3435 4316 : x_boundary_array[iii + side_pts] = xmin + iii * delta_x;
3436 : // xmax boundary
3437 4316 : y_boundary_array[iii + side_pts * 2] = ymin + iii * delta_y;
3438 4316 : x_boundary_array[iii + side_pts * 2] = xmax;
3439 : // ymax boundary
3440 4316 : y_boundary_array[iii + side_pts * 3] = ymax;
3441 4316 : x_boundary_array[iii + side_pts * 3] = xmax - iii * delta_x;
3442 : }
3443 :
3444 : {
3445 513 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3446 1026 : bool success = TransformWithErrorCodes(
3447 513 : boundary_len, &x_boundary_array[0], &y_boundary_array[0], nullptr,
3448 513 : nullptr, anErrorCodes.data());
3449 513 : if (!success)
3450 : {
3451 271 : for (int i = 0; i < boundary_len; ++i)
3452 : {
3453 264 : if (anErrorCodes[i] == 0)
3454 : {
3455 116 : success = true;
3456 116 : break;
3457 : }
3458 : }
3459 123 : if (!success)
3460 7 : return false;
3461 : }
3462 : }
3463 :
3464 506 : if (!degree_output)
3465 : {
3466 299 : *out_xmin =
3467 299 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3468 299 : *out_xmax =
3469 299 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3470 299 : *out_ymin =
3471 299 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3472 299 : *out_ymax =
3473 299 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3474 :
3475 299 : if (poSRSTarget->IsProjected())
3476 : {
3477 598 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3478 :
3479 : auto poBaseTarget = std::unique_ptr<OGRSpatialReference>(
3480 598 : poSRSTarget->CloneGeogCS());
3481 299 : poBaseTarget->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3482 :
3483 : auto poCTBaseTargetToSrc =
3484 : std::unique_ptr<OGRCoordinateTransformation>(
3485 299 : OGRCreateCoordinateTransformation(poBaseTarget.get(),
3486 598 : poSRSSource));
3487 299 : if (poCTBaseTargetToSrc)
3488 : {
3489 : const double dfLon0 =
3490 299 : poSRSTarget->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
3491 :
3492 299 : double dfSignedPoleLat = 0;
3493 299 : double dfAbsPoleLat = 90;
3494 299 : bool bIncludesPole = false;
3495 : const char *pszProjection =
3496 299 : poSRSTarget->GetAttrValue("PROJECTION");
3497 299 : if (pszProjection &&
3498 299 : EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) && dfLon0 == 0)
3499 : {
3500 : // This MAX_LAT_MERCATOR values is equivalent to the
3501 : // semi_major_axis * PI easting/northing value only
3502 : // for EPSG:3857, but it is also quite
3503 : // reasonable for other Mercator projections
3504 204 : constexpr double MAX_LAT_MERCATOR = 85.0511287798066;
3505 204 : dfAbsPoleLat = MAX_LAT_MERCATOR;
3506 : }
3507 :
3508 : // Detect if a point at long = central_meridian and
3509 : // lat = +/- 90deg is included in the extent.
3510 : // Helps for example for EPSG:4326 to ESRI:53037
3511 897 : for (int iSign = -1; iSign <= 1; iSign += 2)
3512 : {
3513 598 : double dfX = dfLon0;
3514 598 : constexpr double EPS = 1e-8;
3515 598 : double dfY = iSign * (dfAbsPoleLat - EPS);
3516 :
3517 598 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3518 598 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3519 598 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3520 1470 : dfY <= ymax &&
3521 274 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr, nullptr,
3522 : nullptr))
3523 : {
3524 274 : dfSignedPoleLat = iSign * dfAbsPoleLat;
3525 274 : bIncludesPole = true;
3526 274 : *out_xmin = std::min(*out_xmin, dfX);
3527 274 : *out_ymin = std::min(*out_ymin, dfY);
3528 274 : *out_xmax = std::max(*out_xmax, dfX);
3529 274 : *out_ymax = std::max(*out_ymax, dfY);
3530 : }
3531 : }
3532 :
3533 : const auto TryAtPlusMinus180 =
3534 217 : [this, dfLon0, xmin, ymin, xmax, ymax, out_xmin, out_ymin,
3535 3941 : out_xmax, out_ymax, &poCTBaseTargetToSrc](double dfLat)
3536 : {
3537 651 : for (int iSign = -1; iSign <= 1; iSign += 2)
3538 : {
3539 434 : constexpr double EPS = 1e-8;
3540 : double dfX =
3541 434 : fmod(dfLon0 + iSign * (180 - EPS) + 180, 360) - 180;
3542 434 : double dfY = dfLat;
3543 :
3544 434 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3545 434 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3546 434 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3547 1153 : dfY <= ymax &&
3548 285 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr,
3549 : nullptr, nullptr))
3550 : {
3551 274 : *out_xmin = std::min(*out_xmin, dfX);
3552 274 : *out_ymin = std::min(*out_ymin, dfY);
3553 274 : *out_xmax = std::max(*out_xmax, dfX);
3554 274 : *out_ymax = std::max(*out_ymax, dfY);
3555 : }
3556 : }
3557 217 : };
3558 :
3559 : // For a projected CRS with a central meridian != 0, try to
3560 : // reproject the points with long = +/- 180deg of the central
3561 : // meridian and at lat = latitude_of_origin.
3562 299 : const double dfLat0 = poSRSTarget->GetNormProjParm(
3563 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
3564 299 : if (dfLon0 != 0)
3565 : {
3566 81 : TryAtPlusMinus180(dfLat0);
3567 : }
3568 :
3569 299 : if (bIncludesPole && dfLat0 != dfSignedPoleLat)
3570 : {
3571 136 : TryAtPlusMinus180(dfSignedPoleLat);
3572 : }
3573 : }
3574 : }
3575 : }
3576 207 : else if (north_pole_in_bounds && output_lon_lat_order)
3577 : {
3578 3 : *out_xmin = -180;
3579 3 : *out_ymin =
3580 3 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3581 3 : *out_xmax = 180;
3582 3 : *out_ymax = 90;
3583 : }
3584 204 : else if (north_pole_in_bounds)
3585 : {
3586 1 : *out_xmin =
3587 1 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3588 1 : *out_ymin = -180;
3589 1 : *out_xmax = 90;
3590 1 : *out_ymax = 180;
3591 : }
3592 203 : else if (south_pole_in_bounds && output_lon_lat_order)
3593 : {
3594 3 : *out_xmin = -180;
3595 3 : *out_ymin = -90;
3596 3 : *out_xmax = 180;
3597 3 : *out_ymax =
3598 3 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3599 : }
3600 200 : else if (south_pole_in_bounds)
3601 : {
3602 1 : *out_xmin = -90;
3603 1 : *out_ymin = -180;
3604 1 : *out_xmax =
3605 1 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3606 1 : *out_ymax = 180;
3607 : }
3608 199 : else if (output_lon_lat_order)
3609 : {
3610 197 : *out_xmin = antimeridian_min(&x_boundary_array[0], anErrorCodes.data(),
3611 : boundary_len);
3612 197 : *out_xmax = antimeridian_max(&x_boundary_array[0], anErrorCodes.data(),
3613 : boundary_len);
3614 197 : *out_ymin =
3615 197 : simple_min(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3616 197 : *out_ymax =
3617 197 : simple_max(&y_boundary_array[0], anErrorCodes.data(), boundary_len);
3618 : }
3619 : else
3620 : {
3621 2 : *out_xmin =
3622 2 : simple_min(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3623 2 : *out_xmax =
3624 2 : simple_max(&x_boundary_array[0], anErrorCodes.data(), boundary_len);
3625 2 : *out_ymin = antimeridian_min(&y_boundary_array[0], anErrorCodes.data(),
3626 : boundary_len);
3627 2 : *out_ymax = antimeridian_max(&y_boundary_array[0], anErrorCodes.data(),
3628 : boundary_len);
3629 : }
3630 :
3631 506 : return *out_xmin != HUGE_VAL && *out_ymin != HUGE_VAL &&
3632 1012 : *out_xmax != HUGE_VAL && *out_ymax != HUGE_VAL;
3633 : }
3634 :
3635 : /************************************************************************/
3636 : /* Clone() */
3637 : /************************************************************************/
3638 :
3639 27 : OGRCoordinateTransformation *OGRProjCT::Clone() const
3640 : {
3641 54 : std::unique_ptr<OGRProjCT> poNewCT(new OGRProjCT(*this));
3642 : #if (PROJ_VERSION_MAJOR * 10000 + PROJ_VERSION_MINOR * 100 + \
3643 : PROJ_VERSION_PATCH) < 80001
3644 : // See https://github.com/OSGeo/PROJ/pull/2582
3645 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3646 : // operation being a set of real operations
3647 27 : bool bCloneDone = ((m_pj == nullptr) == (poNewCT->m_pj == nullptr));
3648 27 : if (!bCloneDone)
3649 : {
3650 1 : poNewCT.reset(new OGRProjCT());
3651 : auto ret =
3652 2 : poNewCT->Initialize(poSRSSource, m_osSrcSRS.c_str(), poSRSTarget,
3653 1 : m_osTargetSRS.c_str(), m_options);
3654 1 : if (!ret)
3655 : {
3656 0 : return nullptr;
3657 : }
3658 : }
3659 : #endif // PROJ_VERSION
3660 27 : return poNewCT.release();
3661 : }
3662 :
3663 : /************************************************************************/
3664 : /* GetInverse() */
3665 : /************************************************************************/
3666 :
3667 1624 : OGRCoordinateTransformation *OGRProjCT::GetInverse() const
3668 : {
3669 1624 : PJ *new_pj = nullptr;
3670 : // m_pj can be nullptr if using m_eStrategy != PROJ
3671 1624 : if (m_pj && !bWebMercatorToWGS84LongLat && !bNoTransform)
3672 : {
3673 : // See https://github.com/OSGeo/PROJ/pull/2582
3674 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3675 : // operation being a set of real operations
3676 1476 : new_pj = proj_clone(OSRGetProjTLSContext(), m_pj);
3677 : }
3678 :
3679 3247 : OGRCoordinateTransformationOptions newOptions(m_options);
3680 1624 : std::swap(newOptions.d->bHasSourceCenterLong,
3681 1624 : newOptions.d->bHasTargetCenterLong);
3682 1624 : std::swap(newOptions.d->dfSourceCenterLong,
3683 1624 : newOptions.d->dfTargetCenterLong);
3684 1623 : newOptions.d->bReverseCO = !newOptions.d->bReverseCO;
3685 1623 : newOptions.d->RefreshCheckWithInvertProj();
3686 :
3687 1624 : if (new_pj == nullptr && !bNoTransform)
3688 : {
3689 164 : return OGRCreateCoordinateTransformation(poSRSTarget, poSRSSource,
3690 164 : newOptions);
3691 : }
3692 :
3693 1460 : auto poNewCT = new OGRProjCT();
3694 :
3695 1460 : if (poSRSTarget)
3696 1454 : poNewCT->poSRSSource = poSRSTarget->Clone();
3697 1460 : poNewCT->m_eSourceFirstAxisOrient = m_eTargetFirstAxisOrient;
3698 1460 : poNewCT->bSourceLatLong = bTargetLatLong;
3699 1460 : poNewCT->bSourceWrap = bTargetWrap;
3700 1460 : poNewCT->dfSourceWrapLong = dfTargetWrapLong;
3701 1460 : poNewCT->bSourceIsDynamicCRS = bTargetIsDynamicCRS;
3702 1460 : poNewCT->dfSourceCoordinateEpoch = dfTargetCoordinateEpoch;
3703 1460 : poNewCT->m_osSrcSRS = m_osTargetSRS;
3704 :
3705 1459 : if (poSRSSource)
3706 1454 : poNewCT->poSRSTarget = poSRSSource->Clone();
3707 1458 : poNewCT->m_eTargetFirstAxisOrient = m_eSourceFirstAxisOrient;
3708 1458 : poNewCT->bTargetLatLong = bSourceLatLong;
3709 1458 : poNewCT->bTargetWrap = bSourceWrap;
3710 1458 : poNewCT->dfTargetWrapLong = dfSourceWrapLong;
3711 1458 : poNewCT->bTargetIsDynamicCRS = bSourceIsDynamicCRS;
3712 1458 : poNewCT->dfTargetCoordinateEpoch = dfSourceCoordinateEpoch;
3713 1458 : poNewCT->m_osTargetSRS = m_osSrcSRS;
3714 :
3715 1460 : poNewCT->ComputeThreshold();
3716 :
3717 1459 : poNewCT->m_pj = new_pj;
3718 1459 : poNewCT->m_bReversePj = !m_bReversePj;
3719 1459 : poNewCT->bNoTransform = bNoTransform;
3720 1459 : poNewCT->m_eStrategy = m_eStrategy;
3721 1459 : poNewCT->m_options = newOptions;
3722 :
3723 1459 : poNewCT->DetectWebMercatorToWGS84();
3724 :
3725 1460 : return poNewCT;
3726 : }
3727 :
3728 : /************************************************************************/
3729 : /* OSRCTCleanCache() */
3730 : /************************************************************************/
3731 :
3732 1128 : void OSRCTCleanCache()
3733 : {
3734 1128 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3735 1128 : delete g_poCTCache;
3736 1128 : g_poCTCache = nullptr;
3737 1128 : }
3738 :
3739 : /************************************************************************/
3740 : /* MakeCacheKey() */
3741 : /************************************************************************/
3742 :
3743 8454 : CTCacheKey OGRProjCT::MakeCacheKey(
3744 : const OGRSpatialReference *poSRS1, const char *pszSrcSRS,
3745 : const OGRSpatialReference *poSRS2, const char *pszTargetSRS,
3746 : const OGRCoordinateTransformationOptions &options)
3747 : {
3748 : const auto GetKeyForSRS =
3749 16908 : [](const OGRSpatialReference *poSRS, const char *pszText)
3750 : {
3751 16908 : if (poSRS)
3752 : {
3753 33714 : std::string ret(pszText);
3754 16857 : const auto &mapping = poSRS->GetDataAxisToSRSAxisMapping();
3755 50847 : for (const auto &axis : mapping)
3756 : {
3757 33990 : ret += std::to_string(axis);
3758 : }
3759 16857 : return ret;
3760 : }
3761 : else
3762 : {
3763 51 : return std::string("null");
3764 : }
3765 : };
3766 :
3767 8454 : std::string ret(GetKeyForSRS(poSRS1, pszSrcSRS));
3768 8454 : ret += GetKeyForSRS(poSRS2, pszTargetSRS);
3769 8454 : ret += options.d->GetKey();
3770 16908 : return ret;
3771 : }
3772 :
3773 : /************************************************************************/
3774 : /* InsertIntoCache() */
3775 : /************************************************************************/
3776 :
3777 3995 : void OGRProjCT::InsertIntoCache(OGRProjCT *poCT)
3778 : {
3779 : {
3780 7990 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3781 3995 : if (g_poCTCache == nullptr)
3782 : {
3783 222 : g_poCTCache = new lru11::Cache<CTCacheKey, CTCacheValue>();
3784 : }
3785 : }
3786 3995 : const auto key = MakeCacheKey(poCT->poSRSSource, poCT->m_osSrcSRS.c_str(),
3787 3995 : poCT->poSRSTarget,
3788 3995 : poCT->m_osTargetSRS.c_str(), poCT->m_options);
3789 :
3790 3995 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3791 3995 : if (g_poCTCache->contains(key))
3792 : {
3793 998 : delete poCT;
3794 998 : return;
3795 : }
3796 2997 : g_poCTCache->insert(key, std::unique_ptr<OGRProjCT>(poCT));
3797 : }
3798 :
3799 : /************************************************************************/
3800 : /* FindFromCache() */
3801 : /************************************************************************/
3802 :
3803 4841 : OGRProjCT *OGRProjCT::FindFromCache(
3804 : const OGRSpatialReference *poSource, const char *pszSrcSRS,
3805 : const OGRSpatialReference *poTarget, const char *pszTargetSRS,
3806 : const OGRCoordinateTransformationOptions &options)
3807 : {
3808 : {
3809 4841 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3810 4841 : if (g_poCTCache == nullptr || g_poCTCache->empty())
3811 382 : return nullptr;
3812 : }
3813 :
3814 : const auto key =
3815 8918 : MakeCacheKey(poSource, pszSrcSRS, poTarget, pszTargetSRS, options);
3816 : // Get value from cache and remove it
3817 8918 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3818 4459 : CTCacheValue *cachedValue = g_poCTCache->getPtr(key);
3819 4459 : if (cachedValue)
3820 : {
3821 2019 : auto poCT = cachedValue->release();
3822 2019 : g_poCTCache->remove(key);
3823 2019 : return poCT;
3824 : }
3825 2440 : return nullptr;
3826 : }
3827 :
3828 : //! @endcond
3829 :
3830 : /************************************************************************/
3831 : /* OCTTransform() */
3832 : /************************************************************************/
3833 :
3834 : /** Transform an array of points
3835 : *
3836 : * @param hTransform Transformation object
3837 : * @param nCount Number of points
3838 : * @param x Array of nCount x values.
3839 : * @param y Array of nCount y values.
3840 : * @param z Array of nCount z values.
3841 : * @return TRUE if a transformation could be found (but not all points may
3842 : * have necessarily succeed to transform), otherwise FALSE.
3843 : */
3844 1024 : int CPL_STDCALL OCTTransform(OGRCoordinateTransformationH hTransform,
3845 : int nCount, double *x, double *y, double *z)
3846 :
3847 : {
3848 1024 : VALIDATE_POINTER1(hTransform, "OCTTransform", FALSE);
3849 :
3850 : return OGRCoordinateTransformation::FromHandle(hTransform)
3851 1024 : ->Transform(nCount, x, y, z);
3852 : }
3853 :
3854 : /************************************************************************/
3855 : /* OCTTransformEx() */
3856 : /************************************************************************/
3857 :
3858 : /** Transform an array of points
3859 : *
3860 : * @param hTransform Transformation object
3861 : * @param nCount Number of points
3862 : * @param x Array of nCount x values.
3863 : * @param y Array of nCount y values.
3864 : * @param z Array of nCount z values.
3865 : * @param pabSuccess Output array of nCount value that will be set to TRUE/FALSE
3866 : * @return TRUE if a transformation could be found (but not all points may
3867 : * have necessarily succeed to transform), otherwise FALSE.
3868 : */
3869 0 : int CPL_STDCALL OCTTransformEx(OGRCoordinateTransformationH hTransform,
3870 : int nCount, double *x, double *y, double *z,
3871 : int *pabSuccess)
3872 :
3873 : {
3874 0 : VALIDATE_POINTER1(hTransform, "OCTTransformEx", FALSE);
3875 :
3876 : return OGRCoordinateTransformation::FromHandle(hTransform)
3877 0 : ->Transform(nCount, x, y, z, pabSuccess);
3878 : }
3879 :
3880 : /************************************************************************/
3881 : /* OCTTransform4D() */
3882 : /************************************************************************/
3883 :
3884 : /** Transform an array of points
3885 : *
3886 : * @param hTransform Transformation object
3887 : * @param nCount Number of points
3888 : * @param x Array of nCount x values. Should not be NULL
3889 : * @param y Array of nCount y values. Should not be NULL
3890 : * @param z Array of nCount z values. Might be NULL
3891 : * @param t Array of nCount time values. Might be NULL
3892 : * @param pabSuccess Output array of nCount value that will be set to
3893 : * TRUE/FALSE. Might be NULL.
3894 : * @since GDAL 3.0
3895 : * @return TRUE if a transformation could be found (but not all points may
3896 : * have necessarily succeed to transform), otherwise FALSE.
3897 : */
3898 15 : int OCTTransform4D(OGRCoordinateTransformationH hTransform, int nCount,
3899 : double *x, double *y, double *z, double *t, int *pabSuccess)
3900 :
3901 : {
3902 15 : VALIDATE_POINTER1(hTransform, "OCTTransform4D", FALSE);
3903 :
3904 15 : return OGRCoordinateTransformation::FromHandle(hTransform)
3905 15 : ->Transform(nCount, x, y, z, t, pabSuccess);
3906 : }
3907 :
3908 : /************************************************************************/
3909 : /* OCTTransform4DWithErrorCodes() */
3910 : /************************************************************************/
3911 :
3912 : /** Transform an array of points
3913 : *
3914 : * @param hTransform Transformation object
3915 : * @param nCount Number of points
3916 : * @param x Array of nCount x values. Should not be NULL
3917 : * @param y Array of nCount y values. Should not be NULL
3918 : * @param z Array of nCount z values. Might be NULL
3919 : * @param t Array of nCount time values. Might be NULL
3920 : * @param panErrorCodes Output array of nCount value that will be set to 0 for
3921 : * success, or a non-zero value for failure. Refer to
3922 : * PROJ 8 public error codes. Might be NULL
3923 : * @since GDAL 3.3, and PROJ 8 to be able to use PROJ public error codes
3924 : * @return TRUE if a transformation could be found (but not all points may
3925 : * have necessarily succeed to transform), otherwise FALSE.
3926 : */
3927 0 : int OCTTransform4DWithErrorCodes(OGRCoordinateTransformationH hTransform,
3928 : int nCount, double *x, double *y, double *z,
3929 : double *t, int *panErrorCodes)
3930 :
3931 : {
3932 0 : VALIDATE_POINTER1(hTransform, "OCTTransform4DWithErrorCodes", FALSE);
3933 :
3934 0 : return OGRCoordinateTransformation::FromHandle(hTransform)
3935 0 : ->TransformWithErrorCodes(nCount, x, y, z, t, panErrorCodes);
3936 : }
3937 :
3938 : /************************************************************************/
3939 : /* OCTTransformBounds() */
3940 : /************************************************************************/
3941 : /** \brief Transform boundary.
3942 : *
3943 : * Transform boundary densifying the edges to account for nonlinear
3944 : * transformations along these edges and extracting the outermost bounds.
3945 : *
3946 : * If the destination CRS is geographic, the first axis is longitude,
3947 : * and *out_xmax < *out_xmin then the bounds crossed the antimeridian.
3948 : * In this scenario there are two polygons, one on each side of the
3949 : * antimeridian. The first polygon should be constructed with
3950 : * (*out_xmin, *out_ymin, 180, ymax) and the second with
3951 : * (-180, *out_ymin, *out_xmax, *out_ymax).
3952 : *
3953 : * If the destination CRS is geographic, the first axis is latitude,
3954 : * and *out_ymax < *out_ymin then the bounds crossed the antimeridian.
3955 : * In this scenario there are two polygons, one on each side of the
3956 : * antimeridian. The first polygon should be constructed with
3957 : * (*out_ymin, *out_xmin, *out_ymax, 180) and the second with
3958 : * (*out_ymin, -180, *out_ymax, *out_xmax).
3959 : *
3960 : * @param hTransform Transformation object
3961 : * @param xmin Minimum bounding coordinate of the first axis in source CRS.
3962 : * @param ymin Minimum bounding coordinate of the second axis in source CRS.
3963 : * @param xmax Maximum bounding coordinate of the first axis in source CRS.
3964 : * @param ymax Maximum bounding coordinate of the second axis in source CRS.
3965 : * @param out_xmin Minimum bounding coordinate of the first axis in target CRS
3966 : * @param out_ymin Minimum bounding coordinate of the second axis in target CRS.
3967 : * @param out_xmax Maximum bounding coordinate of the first axis in target CRS.
3968 : * @param out_ymax Maximum bounding coordinate of the second axis in target CRS.
3969 : * @param densify_pts Recommended to use 21. This is the number of points
3970 : * to use to densify the bounding polygon in the transformation.
3971 : * @return TRUE if successful. FALSE if failures encountered.
3972 : * @since 3.4
3973 : */
3974 25 : int CPL_STDCALL OCTTransformBounds(OGRCoordinateTransformationH hTransform,
3975 : const double xmin, const double ymin,
3976 : const double xmax, const double ymax,
3977 : double *out_xmin, double *out_ymin,
3978 : double *out_xmax, double *out_ymax,
3979 : int densify_pts)
3980 :
3981 : {
3982 25 : VALIDATE_POINTER1(hTransform, "TransformBounds", FALSE);
3983 :
3984 25 : return OGRProjCT::FromHandle(hTransform)
3985 25 : ->TransformBounds(xmin, ymin, xmax, ymax, out_xmin, out_ymin, out_xmax,
3986 25 : out_ymax, densify_pts);
3987 : }
3988 :
3989 : /************************************************************************/
3990 : /* OGRCTDumpStatistics() */
3991 : /************************************************************************/
3992 :
3993 1128 : void OGRCTDumpStatistics()
3994 : {
3995 : #ifdef DEBUG_PERF
3996 : CPLDebug("OGR_CT", "Total time in proj_create_crs_to_crs(): %d ms",
3997 : static_cast<int>(g_dfTotalTimeCRStoCRS * 1000));
3998 : CPLDebug("OGR_CT", "Total time in coordinate transformation: %d ms",
3999 : static_cast<int>(g_dfTotalTimeReprojection * 1000));
4000 : #endif
4001 1128 : }
|
