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