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