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