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 3235 : 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 1185 : 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 6987 : OGRCoordinateTransformationOptions::Private::Private()
124 : {
125 6987 : RefreshCheckWithInvertProj();
126 6987 : }
127 :
128 : /************************************************************************/
129 : /* GetKey() */
130 : /************************************************************************/
131 :
132 6609 : std::string OGRCoordinateTransformationOptions::Private::GetKey() const
133 : {
134 6609 : std::string ret;
135 6609 : ret += std::to_string(static_cast<int>(bHasAreaOfInterest));
136 6609 : ret += std::to_string(dfWestLongitudeDeg);
137 6609 : ret += std::to_string(dfSouthLatitudeDeg);
138 6609 : ret += std::to_string(dfEastLongitudeDeg);
139 6609 : ret += std::to_string(dfNorthLatitudeDeg);
140 6609 : ret += osCoordOperation;
141 6609 : ret += std::to_string(static_cast<int>(bReverseCO));
142 6609 : ret += std::to_string(static_cast<int>(bAllowBallpark));
143 6609 : ret += std::to_string(dfAccuracy);
144 6609 : ret += std::to_string(static_cast<int>(bOnlyBestOptionSet));
145 6609 : ret += std::to_string(static_cast<int>(bOnlyBest));
146 6609 : ret += std::to_string(static_cast<int>(bHasSourceCenterLong));
147 6609 : ret += std::to_string(dfSourceCenterLong);
148 6609 : ret += std::to_string(static_cast<int>(bHasTargetCenterLong));
149 6609 : ret += std::to_string(dfTargetCenterLong);
150 6609 : ret += std::to_string(static_cast<int>(bCheckWithInvertProj));
151 6609 : return ret;
152 : }
153 :
154 : /************************************************************************/
155 : /* RefreshCheckWithInvertProj() */
156 : /************************************************************************/
157 :
158 8137 : void OGRCoordinateTransformationOptions::Private::RefreshCheckWithInvertProj()
159 : {
160 8137 : bCheckWithInvertProj =
161 8137 : CPLTestBool(CPLGetConfigOption("CHECK_WITH_INVERT_PROJ", "NO"));
162 8137 : }
163 :
164 : /************************************************************************/
165 : /* GetAsAProjRecognizableString() */
166 : /************************************************************************/
167 :
168 2049 : static char *GetAsAProjRecognizableString(const OGRSpatialReference *poSRS)
169 : {
170 2049 : 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 2049 : char *pszText = nullptr;
174 2049 : 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 2034 : 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 2033 : poSRS->exportToPROJJSON(&pszText, nullptr);
196 : }
197 :
198 4098 : return pszText;
199 : }
200 :
201 : /************************************************************************/
202 : /* GetTextRepresentation() */
203 : /************************************************************************/
204 :
205 7495 : static char *GetTextRepresentation(const OGRSpatialReference *poSRS)
206 : {
207 5448 : const auto CanUseAuthorityDef = [](const OGRSpatialReference *poSRS1,
208 : OGRSpatialReference *poSRSFromAuth,
209 : const char *pszAuth)
210 : {
211 10872 : if (EQUAL(pszAuth, "EPSG") &&
212 5424 : 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 5446 : return true;
236 : };
237 :
238 7495 : 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 7495 : const char *pszAuth = poSRS->GetAuthorityName(nullptr);
244 7495 : const char *pszCode = poSRS->GetAuthorityCode(nullptr);
245 12980 : if (pszAuth && pszCode &&
246 5485 : CPLTestBool(
247 : CPLGetConfigOption("OGR_CT_PREFER_OFFICIAL_SRS_DEF", "YES")))
248 : {
249 10970 : CPLString osAuthCode(pszAuth);
250 5485 : osAuthCode += ':';
251 5485 : osAuthCode += pszCode;
252 10970 : OGRSpatialReference oTmpSRS;
253 5485 : oTmpSRS.SetFromUserInput(osAuthCode);
254 5485 : oTmpSRS.SetDataAxisToSRSAxisMapping(
255 : poSRS->GetDataAxisToSRSAxisMapping());
256 5485 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
257 : nullptr};
258 5485 : if (oTmpSRS.IsSame(poSRS, apszOptionsIsSame))
259 : {
260 5448 : if (CanUseAuthorityDef(poSRS, &oTmpSRS, pszAuth))
261 : {
262 5446 : pszText = CPLStrdup(osAuthCode);
263 : }
264 : }
265 : }
266 7495 : if (pszText == nullptr)
267 : {
268 2049 : pszText = GetAsAProjRecognizableString(poSRS);
269 : }
270 7495 : return pszText;
271 : }
272 :
273 : /************************************************************************/
274 : /* OGRCoordinateTransformationOptions() */
275 : /************************************************************************/
276 :
277 : /** \brief Constructs a new OGRCoordinateTransformationOptions.
278 : *
279 : * @since GDAL 3.0
280 : */
281 6987 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions()
282 6987 : : d(new Private())
283 : {
284 6987 : }
285 :
286 : /************************************************************************/
287 : /* OGRCoordinateTransformationOptions() */
288 : /************************************************************************/
289 :
290 : /** \brief Copy constructor
291 : *
292 : * @since GDAL 3.1
293 : */
294 1185 : OGRCoordinateTransformationOptions::OGRCoordinateTransformationOptions(
295 1185 : const OGRCoordinateTransformationOptions &other)
296 1185 : : d(new Private(*(other.d)))
297 : {
298 1185 : }
299 :
300 : /************************************************************************/
301 : /* operator =() */
302 : /************************************************************************/
303 :
304 : /** \brief Assignment operator
305 : *
306 : * @since GDAL 3.1
307 : */
308 : OGRCoordinateTransformationOptions &
309 3235 : OGRCoordinateTransformationOptions::operator=(
310 : const OGRCoordinateTransformationOptions &other)
311 : {
312 3235 : if (this != &other)
313 : {
314 3235 : *d = *(other.d);
315 : }
316 3235 : return *this;
317 : }
318 :
319 : /************************************************************************/
320 : /* OGRCoordinateTransformationOptions() */
321 : /************************************************************************/
322 :
323 : /** \brief Destroys a OGRCoordinateTransformationOptions.
324 : *
325 : * @since GDAL 3.0
326 : */
327 7783 : OGRCoordinateTransformationOptions::~OGRCoordinateTransformationOptions()
328 : {
329 7783 : }
330 :
331 : /************************************************************************/
332 : /* OCTNewCoordinateTransformationOptions() */
333 : /************************************************************************/
334 :
335 : /** \brief Create coordinate transformation options.
336 : *
337 : * To be freed with OCTDestroyCoordinateTransformationOptions()
338 : *
339 : * @since GDAL 3.0
340 : */
341 13 : OGRCoordinateTransformationOptionsH OCTNewCoordinateTransformationOptions(void)
342 : {
343 13 : return new OGRCoordinateTransformationOptions();
344 : }
345 :
346 : /************************************************************************/
347 : /* OCTDestroyCoordinateTransformationOptions() */
348 : /************************************************************************/
349 :
350 : /** \brief Destroy coordinate transformation options.
351 : *
352 : * @since GDAL 3.0
353 : */
354 11 : void OCTDestroyCoordinateTransformationOptions(
355 : OGRCoordinateTransformationOptionsH hOptions)
356 : {
357 11 : delete hOptions;
358 11 : }
359 :
360 : /************************************************************************/
361 : /* SetAreaOfInterest() */
362 : /************************************************************************/
363 :
364 : /** \brief Sets an area of interest.
365 : *
366 : * The west longitude is generally lower than the east longitude, except for
367 : * areas of interest that go across the anti-meridian.
368 : *
369 : * @param dfWestLongitudeDeg West longitude (in degree). Must be in [-180,180]
370 : * @param dfSouthLatitudeDeg South latitude (in degree). Must be in [-90,90]
371 : * @param dfEastLongitudeDeg East longitude (in degree). Must be in [-180,180]
372 : * @param dfNorthLatitudeDeg North latitude (in degree). Must be in [-90,90]
373 : * @return true in case of success.
374 : *
375 : * @since GDAL 3.0
376 : */
377 852 : bool OGRCoordinateTransformationOptions::SetAreaOfInterest(
378 : double dfWestLongitudeDeg, double dfSouthLatitudeDeg,
379 : double dfEastLongitudeDeg, double dfNorthLatitudeDeg)
380 : {
381 852 : if (std::fabs(dfWestLongitudeDeg) > 180)
382 : {
383 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfWestLongitudeDeg");
384 1 : return false;
385 : }
386 851 : if (std::fabs(dfSouthLatitudeDeg) > 90)
387 : {
388 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfSouthLatitudeDeg");
389 1 : return false;
390 : }
391 850 : if (std::fabs(dfEastLongitudeDeg) > 180)
392 : {
393 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfEastLongitudeDeg");
394 1 : return false;
395 : }
396 849 : if (std::fabs(dfNorthLatitudeDeg) > 90)
397 : {
398 1 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dfNorthLatitudeDeg");
399 1 : return false;
400 : }
401 848 : if (dfSouthLatitudeDeg > dfNorthLatitudeDeg)
402 : {
403 0 : CPLError(CE_Failure, CPLE_AppDefined,
404 : "dfSouthLatitudeDeg should be lower than dfNorthLatitudeDeg");
405 0 : return false;
406 : }
407 848 : d->bHasAreaOfInterest = true;
408 848 : d->dfWestLongitudeDeg = dfWestLongitudeDeg;
409 848 : d->dfSouthLatitudeDeg = dfSouthLatitudeDeg;
410 848 : d->dfEastLongitudeDeg = dfEastLongitudeDeg;
411 848 : d->dfNorthLatitudeDeg = dfNorthLatitudeDeg;
412 848 : return true;
413 : }
414 :
415 : /************************************************************************/
416 : /* OCTCoordinateTransformationOptionsSetAreaOfInterest() */
417 : /************************************************************************/
418 :
419 : /** \brief Sets an area of interest.
420 : *
421 : * See OGRCoordinateTransformationOptions::SetAreaOfInterest()
422 : *
423 : * @since GDAL 3.0
424 : */
425 5 : int OCTCoordinateTransformationOptionsSetAreaOfInterest(
426 : OGRCoordinateTransformationOptionsH hOptions, double dfWestLongitudeDeg,
427 : double dfSouthLatitudeDeg, double dfEastLongitudeDeg,
428 : double dfNorthLatitudeDeg)
429 : {
430 5 : return hOptions->SetAreaOfInterest(dfWestLongitudeDeg, dfSouthLatitudeDeg,
431 5 : dfEastLongitudeDeg, dfNorthLatitudeDeg);
432 : }
433 :
434 : /************************************************************************/
435 : /* SetCoordinateOperation() */
436 : /************************************************************************/
437 :
438 : /** \brief Sets a coordinate operation.
439 : *
440 : * This is a user override to be used instead of the normally computed pipeline.
441 : *
442 : * The pipeline must take into account the axis order of the source and target
443 : * SRS.
444 : *
445 : * The pipeline may be provided as a PROJ string (single step operation or
446 : * multiple step string starting with +proj=pipeline), a WKT2 string describing
447 : * a CoordinateOperation, or a "urn:ogc:def:coordinateOperation:EPSG::XXXX" URN
448 : *
449 : * @param pszCO PROJ or WKT string describing a coordinate operation
450 : * @param bReverseCO Whether the PROJ or WKT string should be evaluated in the
451 : * reverse path
452 : * @return true in case of success.
453 : *
454 : * @since GDAL 3.0
455 : */
456 22 : bool OGRCoordinateTransformationOptions::SetCoordinateOperation(
457 : const char *pszCO, bool bReverseCO)
458 : {
459 22 : d->osCoordOperation = pszCO ? pszCO : "";
460 22 : d->bReverseCO = bReverseCO;
461 22 : return true;
462 : }
463 :
464 : /************************************************************************/
465 : /* SetSourceCenterLong() */
466 : /************************************************************************/
467 :
468 : /*! @cond Doxygen_Suppress */
469 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 2 : bool OGRCoordinateTransformationOptions::SetDesiredAccuracy(double dfAccuracy)
533 : {
534 2 : d->dfAccuracy = dfAccuracy;
535 2 : return true;
536 : }
537 :
538 : /************************************************************************/
539 : /* OCTCoordinateTransformationOptionsSetDesiredAccuracy() */
540 : /************************************************************************/
541 :
542 : /** \brief Sets the desired accuracy for coordinate operations.
543 : *
544 : * See OGRCoordinateTransformationOptions::SetDesiredAccuracy()
545 : *
546 : * @since GDAL 3.3
547 : */
548 2 : int OCTCoordinateTransformationOptionsSetDesiredAccuracy(
549 : OGRCoordinateTransformationOptionsH hOptions, double dfAccuracy)
550 : {
551 : // cppcheck-suppress knownConditionTrueFalse
552 2 : return hOptions->SetDesiredAccuracy(dfAccuracy);
553 : }
554 :
555 : /************************************************************************/
556 : /* SetBallparkAllowed() */
557 : /************************************************************************/
558 :
559 : /** \brief Sets whether ballpark transformations are allowed.
560 : *
561 : * By default, PROJ may generate "ballpark transformations" (see
562 : * https://proj.org/glossary.html) when precise datum transformations are
563 : * missing. For high accuracy use cases, such transformations might not be
564 : * allowed.
565 : *
566 : * If this option is specified with PROJ < 8, the OGR_CT_OP_SELECTION
567 : * configuration option will default to BEST_ACCURACY.
568 : *
569 : * @param bAllowBallpark false to disable the user of ballpark transformations
570 : *
571 : * @since GDAL 3.3
572 : */
573 2 : bool OGRCoordinateTransformationOptions::SetBallparkAllowed(bool bAllowBallpark)
574 : {
575 2 : d->bAllowBallpark = bAllowBallpark;
576 2 : return true;
577 : }
578 :
579 : /************************************************************************/
580 : /* OCTCoordinateTransformationOptionsSetBallparkAllowed() */
581 : /************************************************************************/
582 :
583 : /** \brief Sets whether ballpark transformations are allowed.
584 : *
585 : * See OGRCoordinateTransformationOptions::SetDesiredAccuracy()
586 : *
587 : * @since GDAL 3.3 and PROJ 8
588 : */
589 2 : int OCTCoordinateTransformationOptionsSetBallparkAllowed(
590 : OGRCoordinateTransformationOptionsH hOptions, int bAllowBallpark)
591 : {
592 : // cppcheck-suppress knownConditionTrueFalse
593 2 : return hOptions->SetBallparkAllowed(CPL_TO_BOOL(bAllowBallpark));
594 : }
595 :
596 : /************************************************************************/
597 : /* SetOnlyBest() */
598 : /************************************************************************/
599 :
600 : /** \brief Sets whether only the "best" operation should be used.
601 : *
602 : * By default (at least in the PROJ 9.x series), PROJ may use coordinate
603 : * operations that are not the "best" if resources (typically grids) needed
604 : * to use them are missing. It will then fallback to other coordinate operations
605 : * that have a lesser accuracy, for example using Helmert transformations,
606 : * or in the absence of such operations, to ones with potential very 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 6060 : void reset()
662 : {
663 6060 : if (m_pj)
664 : {
665 2828 : proj_assign_context(m_pj, OSRGetProjTLSContext());
666 2828 : proj_destroy(m_pj);
667 : }
668 6060 : }
669 :
670 : public:
671 3235 : PjPtr() : m_pj(nullptr)
672 : {
673 3235 : }
674 :
675 0 : explicit PjPtr(PJ *pjIn) : m_pj(pjIn)
676 : {
677 0 : }
678 :
679 2871 : ~PjPtr()
680 2871 : {
681 2871 : reset();
682 2871 : }
683 :
684 23 : PjPtr(const PjPtr &other)
685 23 : : m_pj((other.m_pj != nullptr)
686 23 : ? (proj_clone(OSRGetProjTLSContext(), other.m_pj))
687 23 : : (nullptr))
688 : {
689 23 : }
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 3192 : PjPtr &operator=(PJ *pjIn)
709 : {
710 3192 : if (m_pj != pjIn)
711 : {
712 3189 : reset();
713 3189 : m_pj = pjIn;
714 : }
715 3192 : return *this;
716 : }
717 :
718 1819400 : operator PJ *()
719 : {
720 1819400 : return m_pj;
721 : }
722 :
723 2320 : operator const PJ *() const
724 : {
725 2320 : 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 3255 : void SetEmitErrors(bool bEmitErrors) override
846 : {
847 3255 : m_bEmitErrors = bEmitErrors;
848 3255 : }
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 373 : OCTDestroyCoordinateTransformation(OGRCoordinateTransformationH hCT)
878 :
879 : {
880 373 : OGRCoordinateTransformation::DestroyCT(
881 : OGRCoordinateTransformation::FromHandle(hCT));
882 373 : }
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 3195 : void OGRCoordinateTransformation::DestroyCT(OGRCoordinateTransformation *poCT)
906 : {
907 3195 : auto poProjCT = dynamic_cast<OGRProjCT *>(poCT);
908 3195 : if (poProjCT)
909 : {
910 3189 : OGRProjCT::InsertIntoCache(poProjCT);
911 : }
912 : else
913 : {
914 6 : delete poCT;
915 : }
916 3195 : }
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 2714 : OGRCreateCoordinateTransformation(const OGRSpatialReference *poSource,
945 : const OGRSpatialReference *poTarget)
946 :
947 : {
948 2714 : return OGRCreateCoordinateTransformation(
949 5428 : 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 3760 : OGRCoordinateTransformation *OGRCreateCoordinateTransformation(
1018 : const OGRSpatialReference *poSource, const OGRSpatialReference *poTarget,
1019 : const OGRCoordinateTransformationOptions &options)
1020 :
1021 : {
1022 3760 : char *pszSrcSRS = poSource ? GetTextRepresentation(poSource) : nullptr;
1023 3760 : char *pszTargetSRS = poTarget ? GetTextRepresentation(poTarget) : nullptr;
1024 : // Try to find if we have a match in the case
1025 3760 : OGRProjCT *poCT = OGRProjCT::FindFromCache(poSource, pszSrcSRS, poTarget,
1026 : pszTargetSRS, options);
1027 3760 : if (poCT == nullptr)
1028 : {
1029 2093 : poCT = new OGRProjCT();
1030 2093 : if (!poCT->Initialize(poSource, pszSrcSRS, poTarget, pszTargetSRS,
1031 : options))
1032 : {
1033 8 : delete poCT;
1034 8 : poCT = nullptr;
1035 : }
1036 : }
1037 3760 : CPLFree(pszSrcSRS);
1038 3760 : CPLFree(pszTargetSRS);
1039 :
1040 3760 : 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 346 : OGRCoordinateTransformationH CPL_STDCALL OCTNewCoordinateTransformation(
1070 : OGRSpatialReferenceH hSourceSRS, OGRSpatialReferenceH hTargetSRS)
1071 :
1072 : {
1073 : return reinterpret_cast<OGRCoordinateTransformationH>(
1074 346 : OGRCreateCoordinateTransformation(
1075 : reinterpret_cast<OGRSpatialReference *>(hSourceSRS),
1076 346 : 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 12 : OCTNewCoordinateTransformationEx(OGRSpatialReferenceH hSourceSRS,
1123 : OGRSpatialReferenceH hTargetSRS,
1124 : OGRCoordinateTransformationOptionsH hOptions)
1125 :
1126 : {
1127 : return reinterpret_cast<OGRCoordinateTransformationH>(
1128 12 : OGRCreateCoordinateTransformation(
1129 : reinterpret_cast<OGRSpatialReference *>(hSourceSRS),
1130 : reinterpret_cast<OGRSpatialReference *>(hTargetSRS),
1131 24 : 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 3235 : OGRProjCT::OGRProjCT()
1242 : {
1243 3235 : }
1244 :
1245 : /************************************************************************/
1246 : /* OGRProjCT(const OGRProjCT& other) */
1247 : /************************************************************************/
1248 :
1249 23 : OGRProjCT::OGRProjCT(const OGRProjCT &other)
1250 23 : : poSRSSource((other.poSRSSource != nullptr) ? (other.poSRSSource->Clone())
1251 : : (nullptr)),
1252 23 : m_eSourceFirstAxisOrient(other.m_eSourceFirstAxisOrient),
1253 23 : bSourceLatLong(other.bSourceLatLong), bSourceWrap(other.bSourceWrap),
1254 23 : dfSourceWrapLong(other.dfSourceWrapLong),
1255 23 : bSourceIsDynamicCRS(other.bSourceIsDynamicCRS),
1256 23 : dfSourceCoordinateEpoch(other.dfSourceCoordinateEpoch),
1257 23 : m_osSrcSRS(other.m_osSrcSRS),
1258 23 : poSRSTarget((other.poSRSTarget != nullptr) ? (other.poSRSTarget->Clone())
1259 : : (nullptr)),
1260 23 : m_eTargetFirstAxisOrient(other.m_eTargetFirstAxisOrient),
1261 23 : bTargetLatLong(other.bTargetLatLong), bTargetWrap(other.bTargetWrap),
1262 23 : dfTargetWrapLong(other.dfTargetWrapLong),
1263 23 : bTargetIsDynamicCRS(other.bTargetIsDynamicCRS),
1264 23 : dfTargetCoordinateEpoch(other.dfTargetCoordinateEpoch),
1265 23 : m_osTargetSRS(other.m_osTargetSRS),
1266 23 : bWebMercatorToWGS84LongLat(other.bWebMercatorToWGS84LongLat),
1267 23 : nErrorCount(other.nErrorCount), dfThreshold(other.dfThreshold),
1268 : m_psLastContext(nullptr),
1269 46 : m_nLastContextThreadId(std::this_thread::get_id()), m_pj(other.m_pj),
1270 23 : m_bReversePj(other.m_bReversePj), m_bEmitErrors(other.m_bEmitErrors),
1271 23 : bNoTransform(other.bNoTransform), m_eStrategy(other.m_eStrategy),
1272 23 : m_oTransformations(other.m_oTransformations),
1273 23 : m_iCurTransformation(other.m_iCurTransformation),
1274 69 : m_options(other.m_options)
1275 : {
1276 23 : }
1277 :
1278 : /************************************************************************/
1279 : /* ~OGRProjCT() */
1280 : /************************************************************************/
1281 :
1282 5742 : OGRProjCT::~OGRProjCT()
1283 :
1284 : {
1285 2871 : if (poSRSSource != nullptr)
1286 : {
1287 2860 : poSRSSource->Release();
1288 : }
1289 :
1290 2871 : if (poSRSTarget != nullptr)
1291 : {
1292 2858 : poSRSTarget->Release();
1293 : }
1294 5742 : }
1295 :
1296 : /************************************************************************/
1297 : /* ComputeThreshold() */
1298 : /************************************************************************/
1299 :
1300 3234 : void OGRProjCT::ComputeThreshold()
1301 : {
1302 : // The threshold is experimental. Works well with the cases of ticket #2305.
1303 3234 : if (bSourceLatLong)
1304 : {
1305 : // coverity[tainted_data]
1306 2347 : 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 887 : dfThreshold = CPLAtof(CPLGetConfigOption("THRESHOLD", "10000"));
1314 : }
1315 3234 : }
1316 :
1317 : /************************************************************************/
1318 : /* DetectWebMercatorToWGS84() */
1319 : /************************************************************************/
1320 :
1321 3234 : void OGRProjCT::DetectWebMercatorToWGS84()
1322 : {
1323 : // Detect webmercator to WGS84
1324 6432 : if (m_options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1325 10347 : poSRSSource->IsProjected() && poSRSTarget->IsGeographic() &&
1326 681 : ((m_eTargetFirstAxisOrient == OAO_North &&
1327 574 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1328 3808 : std::vector<int>{2, 1}) ||
1329 128 : (m_eTargetFirstAxisOrient == OAO_East &&
1330 107 : poSRSTarget->GetDataAxisToSRSAxisMapping() ==
1331 3341 : 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 652 : const char *pszSourceAuth = poSRSSource->GetAuthorityName(nullptr);
1337 652 : const char *pszSourceCode = poSRSSource->GetAuthorityCode(nullptr);
1338 652 : const char *pszTargetAuth = poSRSTarget->GetAuthorityName(nullptr);
1339 652 : const char *pszTargetCode = poSRSTarget->GetAuthorityCode(nullptr);
1340 652 : if (pszSourceAuth && pszSourceCode && pszTargetAuth && pszTargetCode &&
1341 381 : EQUAL(pszSourceAuth, "EPSG") && EQUAL(pszTargetAuth, "EPSG"))
1342 : {
1343 369 : bWebMercatorToWGS84LongLat =
1344 602 : (EQUAL(pszSourceCode, "3857") ||
1345 233 : EQUAL(pszSourceCode, "3785") || // deprecated
1346 738 : EQUAL(pszSourceCode, "900913")) && // deprecated
1347 136 : EQUAL(pszTargetCode, "4326");
1348 : }
1349 : else
1350 : {
1351 283 : CPLPushErrorHandler(CPLQuietErrorHandler);
1352 283 : char *pszSrcProj4Defn = nullptr;
1353 283 : poSRSSource->exportToProj4(&pszSrcProj4Defn);
1354 :
1355 283 : char *pszDstProj4Defn = nullptr;
1356 283 : poSRSTarget->exportToProj4(&pszDstProj4Defn);
1357 283 : CPLPopErrorHandler();
1358 :
1359 283 : if (pszSrcProj4Defn && pszDstProj4Defn)
1360 : {
1361 283 : if (pszSrcProj4Defn[0] != '\0' &&
1362 283 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] == ' ')
1363 0 : pszSrcProj4Defn[strlen(pszSrcProj4Defn) - 1] = 0;
1364 283 : if (pszDstProj4Defn[0] != '\0' &&
1365 283 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] == ' ')
1366 0 : pszDstProj4Defn[strlen(pszDstProj4Defn) - 1] = 0;
1367 283 : char *pszNeedle = strstr(pszSrcProj4Defn, " ");
1368 283 : if (pszNeedle)
1369 0 : memmove(pszNeedle, pszNeedle + 1,
1370 0 : strlen(pszNeedle + 1) + 1);
1371 283 : pszNeedle = strstr(pszDstProj4Defn, " ");
1372 283 : if (pszNeedle)
1373 0 : memmove(pszNeedle, pszNeedle + 1,
1374 0 : strlen(pszNeedle + 1) + 1);
1375 :
1376 283 : if ((strstr(pszDstProj4Defn, "+datum=WGS84") != nullptr ||
1377 152 : strstr(pszDstProj4Defn,
1378 : "+ellps=WGS84 +towgs84=0,0,0,0,0,0,0 ") !=
1379 131 : nullptr) &&
1380 131 : 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 283 : CPLFree(pszSrcProj4Defn);
1423 283 : CPLFree(pszDstProj4Defn);
1424 : }
1425 :
1426 652 : if (bWebMercatorToWGS84LongLat)
1427 : {
1428 184 : CPLDebug("OGRCT", "Using WebMercator to WGS84 optimization");
1429 : }
1430 : }
1431 3234 : }
1432 :
1433 : /************************************************************************/
1434 : /* Initialize() */
1435 : /************************************************************************/
1436 :
1437 2094 : 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 2094 : m_options = options;
1445 :
1446 2094 : if (poSourceIn == nullptr || poTargetIn == nullptr)
1447 : {
1448 13 : 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 2093 : if (poSourceIn)
1458 : {
1459 2083 : poSRSSource = poSourceIn->Clone();
1460 2083 : m_osSrcSRS = pszSrcSRS;
1461 : }
1462 2093 : if (poTargetIn)
1463 : {
1464 2082 : poSRSTarget = poTargetIn->Clone();
1465 2082 : m_osTargetSRS = pszTargetSRS;
1466 : }
1467 :
1468 : // To easy quick&dirty compatibility with GDAL < 3.0
1469 2093 : 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 2093 : if (poSRSSource)
1479 : {
1480 2083 : bSourceLatLong = CPL_TO_BOOL(poSRSSource->IsGeographic());
1481 2083 : bSourceIsDynamicCRS = poSRSSource->IsDynamic();
1482 2083 : dfSourceCoordinateEpoch = poSRSSource->GetCoordinateEpoch();
1483 2083 : if (!bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0)
1484 : {
1485 3 : bSourceIsDynamicCRS = poSRSSource->HasPointMotionOperation();
1486 : }
1487 2083 : poSRSSource->GetAxis(nullptr, 0, &m_eSourceFirstAxisOrient);
1488 : }
1489 2093 : if (poSRSTarget)
1490 : {
1491 2082 : bTargetLatLong = CPL_TO_BOOL(poSRSTarget->IsGeographic());
1492 2082 : bTargetIsDynamicCRS = poSRSTarget->IsDynamic();
1493 2082 : dfTargetCoordinateEpoch = poSRSTarget->GetCoordinateEpoch();
1494 2082 : if (!bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0)
1495 : {
1496 2 : bTargetIsDynamicCRS = poSRSTarget->HasPointMotionOperation();
1497 : }
1498 2082 : poSRSTarget->GetAxis(nullptr, 0, &m_eTargetFirstAxisOrient);
1499 : }
1500 :
1501 : #if PROJ_VERSION_MAJOR < 9 || \
1502 : (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR < 4)
1503 2093 : if (bSourceIsDynamicCRS && bTargetIsDynamicCRS &&
1504 1089 : 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 2093 : bSourceWrap = false;
1518 2093 : dfSourceWrapLong = 0.0;
1519 :
1520 2093 : bTargetWrap = false;
1521 2093 : dfTargetWrapLong = 0.0;
1522 :
1523 : /* -------------------------------------------------------------------- */
1524 : /* Preliminary logic to setup wrapping. */
1525 : /* -------------------------------------------------------------------- */
1526 2093 : 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 2093 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1539 2093 : pszCENTER_LONG =
1540 2093 : poSRSSource ? poSRSSource->GetExtension("GEOGCS", "CENTER_LONG")
1541 : : nullptr;
1542 : }
1543 2093 : 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 2093 : 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 2093 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
1558 2093 : pszCENTER_LONG =
1559 2093 : poSRSTarget ? poSRSTarget->GetExtension("GEOGCS", "CENTER_LONG")
1560 : : nullptr;
1561 : }
1562 2093 : 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 2093 : 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 2093 : ComputeThreshold();
1576 :
1577 2093 : DetectWebMercatorToWGS84();
1578 :
1579 : const char *pszCTOpSelection =
1580 2093 : CPLGetConfigOption("OGR_CT_OP_SELECTION", nullptr);
1581 2093 : 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 2093 : if (options.d->dfAccuracy >= 0 || !options.d->bAllowBallpark)
1597 : {
1598 4 : m_eStrategy = Strategy::BEST_ACCURACY;
1599 : }
1600 : }
1601 : #endif
1602 2093 : if (m_eStrategy == Strategy::PROJ)
1603 : {
1604 : const char *pszUseApproxTMERC =
1605 2089 : CPLGetConfigOption("OSR_USE_APPROX_TMERC", nullptr);
1606 2089 : 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 2093 : if (!options.d->osCoordOperation.empty())
1615 : {
1616 20 : auto ctx = OSRGetProjTLSContext();
1617 20 : m_pj = proj_create(ctx, options.d->osCoordOperation);
1618 20 : 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 20 : m_bReversePj = options.d->bReverseCO;
1626 : #ifdef DEBUG
1627 20 : auto info = proj_pj_info(m_pj);
1628 20 : CPLDebug("OGRCT", "%s %s(user set)", info.definition,
1629 20 : m_bReversePj ? "(reversed) " : "");
1630 : #endif
1631 : }
1632 2073 : 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 2035 : CPLDebug("OGR_CT", "Source CRS: '%s'", pszSrcSRS);
1641 2035 : if (dfSourceCoordinateEpoch > 0)
1642 6 : CPLDebug("OGR_CT", "Source coordinate epoch: %.3f",
1643 : dfSourceCoordinateEpoch);
1644 2035 : CPLDebug("OGR_CT", "Target CRS: '%s'", pszTargetSRS);
1645 2035 : if (dfTargetCoordinateEpoch > 0)
1646 10 : CPLDebug("OGR_CT", "Target coordinate epoch: %.3f",
1647 : dfTargetCoordinateEpoch);
1648 : #endif
1649 :
1650 2035 : if (m_eStrategy == Strategy::PROJ)
1651 : {
1652 2027 : PJ_AREA *area = nullptr;
1653 2027 : if (options.d->bHasAreaOfInterest)
1654 : {
1655 230 : area = proj_area_create();
1656 230 : proj_area_set_bbox(area, options.d->dfWestLongitudeDeg,
1657 230 : options.d->dfSouthLatitudeDeg,
1658 230 : options.d->dfEastLongitudeDeg,
1659 230 : options.d->dfNorthLatitudeDeg);
1660 : }
1661 2027 : 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 2027 : m_pj = proj_create_crs_to_crs(ctx, pszSrcSRS, pszTargetSRS, area);
1711 : #endif
1712 2027 : if (area)
1713 230 : proj_area_destroy(area);
1714 2027 : 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 8 : else if (!ListCoordinateOperations(pszSrcSRS, pszTargetSRS, options))
1723 : {
1724 4 : CPLError(CE_Failure, CPLE_NotSupported,
1725 : "Cannot find coordinate operations from `%s' to `%s'",
1726 : pszSrcSRS, pszTargetSRS);
1727 4 : 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 4152 : if (options.d->osCoordOperation.empty() && poSRSSource && poSRSTarget &&
1741 2066 : (dfSourceCoordinateEpoch == 0 || dfTargetCoordinateEpoch == 0 ||
1742 0 : dfSourceCoordinateEpoch == dfTargetCoordinateEpoch))
1743 : {
1744 : // Determine if we can skip the transformation completely.
1745 2066 : const char *const apszOptionsIsSame[] = {"CRITERION=EQUIVALENT",
1746 : nullptr};
1747 2066 : bNoTransform =
1748 3989 : !bSourceWrap && !bTargetWrap &&
1749 1923 : CPL_TO_BOOL(poSRSSource->IsSame(poSRSTarget, apszOptionsIsSame));
1750 : }
1751 :
1752 2086 : 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 0 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
1768 : "OSR_USE_ETMERC is a legacy configuration option, which "
1769 : "now has only effect when set to NO (YES is the default). "
1770 : "Use OSR_USE_APPROX_TMERC=YES instead");
1771 0 : bForceApproxTMerc = !CPLTestBool(pszUseETMERC);
1772 : }
1773 : else
1774 : {
1775 : const char *pszUseApproxTMERC =
1776 4 : CPLGetConfigOption("OSR_USE_APPROX_TMERC", nullptr);
1777 4 : if (pszUseApproxTMERC && pszUseApproxTMERC[0])
1778 : {
1779 4 : bForceApproxTMerc = CPLTestBool(pszUseApproxTMERC);
1780 : }
1781 : }
1782 4 : const char *options[] = {
1783 4 : bForceApproxTMerc ? "USE_APPROX_TMERC=YES" : nullptr, nullptr};
1784 4 : auto proj_string = proj_as_proj_string(ctx, op, PJ_PROJ_5, options);
1785 4 : if (!proj_string)
1786 : {
1787 0 : return nullptr;
1788 : }
1789 4 : if (osOutProjString)
1790 0 : *osOutProjString = proj_string;
1791 :
1792 4 : if (proj_string[0] == '\0')
1793 : {
1794 : /* Null transform ? */
1795 0 : return proj_create(ctx, "proj=affine");
1796 : }
1797 : else
1798 : {
1799 4 : return proj_create(ctx, proj_string);
1800 : }
1801 : }
1802 :
1803 : /************************************************************************/
1804 : /* ListCoordinateOperations() */
1805 : /************************************************************************/
1806 :
1807 8 : bool OGRProjCT::ListCoordinateOperations(
1808 : const char *pszSrcSRS, const char *pszTargetSRS,
1809 : const OGRCoordinateTransformationOptions &options)
1810 : {
1811 8 : auto ctx = OSRGetProjTLSContext();
1812 :
1813 8 : auto src = proj_create(ctx, pszSrcSRS);
1814 8 : if (!src)
1815 : {
1816 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot instantiate source_crs");
1817 0 : return false;
1818 : }
1819 :
1820 8 : auto dst = proj_create(ctx, pszTargetSRS);
1821 8 : if (!dst)
1822 : {
1823 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot instantiate target_crs");
1824 0 : proj_destroy(src);
1825 0 : return false;
1826 : }
1827 :
1828 8 : auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
1829 8 : if (!operation_ctx)
1830 : {
1831 0 : proj_destroy(src);
1832 0 : proj_destroy(dst);
1833 0 : return false;
1834 : }
1835 :
1836 8 : proj_operation_factory_context_set_spatial_criterion(
1837 : ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
1838 8 : proj_operation_factory_context_set_grid_availability_use(
1839 : ctx, operation_ctx,
1840 : #if PROJ_VERSION_MAJOR >= 7
1841 : proj_context_is_network_enabled(ctx)
1842 : ? PROJ_GRID_AVAILABILITY_KNOWN_AVAILABLE
1843 : :
1844 : #endif
1845 : PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
1846 :
1847 8 : if (options.d->bHasAreaOfInterest)
1848 : {
1849 0 : proj_operation_factory_context_set_area_of_interest(
1850 0 : ctx, operation_ctx, options.d->dfWestLongitudeDeg,
1851 0 : options.d->dfSouthLatitudeDeg, options.d->dfEastLongitudeDeg,
1852 0 : options.d->dfNorthLatitudeDeg);
1853 : }
1854 :
1855 8 : if (options.d->dfAccuracy >= 0)
1856 2 : proj_operation_factory_context_set_desired_accuracy(
1857 2 : ctx, operation_ctx, options.d->dfAccuracy);
1858 8 : if (!options.d->bAllowBallpark)
1859 : {
1860 : #if PROJ_VERSION_MAJOR > 7 || \
1861 : (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 1)
1862 : proj_operation_factory_context_set_allow_ballpark_transformations(
1863 : ctx, operation_ctx, FALSE);
1864 : #else
1865 2 : if (options.d->dfAccuracy < 0)
1866 : {
1867 2 : proj_operation_factory_context_set_desired_accuracy(
1868 : ctx, operation_ctx, HUGE_VAL);
1869 : }
1870 : #endif
1871 : }
1872 :
1873 8 : auto op_list = proj_create_operations(ctx, src, dst, operation_ctx);
1874 :
1875 8 : if (!op_list)
1876 : {
1877 0 : proj_operation_factory_context_destroy(operation_ctx);
1878 0 : proj_destroy(src);
1879 0 : proj_destroy(dst);
1880 0 : return false;
1881 : }
1882 :
1883 8 : auto op_count = proj_list_get_count(op_list);
1884 8 : if (op_count == 0)
1885 : {
1886 4 : proj_list_destroy(op_list);
1887 4 : proj_operation_factory_context_destroy(operation_ctx);
1888 4 : proj_destroy(src);
1889 4 : proj_destroy(dst);
1890 4 : CPLDebug("OGRCT", "No operation found matching criteria");
1891 4 : return false;
1892 : }
1893 :
1894 0 : if (op_count == 1 || options.d->bHasAreaOfInterest ||
1895 4 : proj_get_type(src) == PJ_TYPE_GEOCENTRIC_CRS ||
1896 0 : proj_get_type(dst) == PJ_TYPE_GEOCENTRIC_CRS)
1897 : {
1898 4 : auto op = proj_list_get(ctx, op_list, 0);
1899 4 : CPLAssert(op);
1900 4 : m_pj = op_to_pj(ctx, op);
1901 8 : CPLString osName;
1902 4 : auto name = proj_get_name(op);
1903 4 : if (name)
1904 4 : osName = name;
1905 4 : proj_destroy(op);
1906 4 : proj_list_destroy(op_list);
1907 4 : proj_operation_factory_context_destroy(operation_ctx);
1908 4 : proj_destroy(src);
1909 4 : proj_destroy(dst);
1910 4 : if (!m_pj)
1911 0 : return false;
1912 : #ifdef DEBUG
1913 4 : auto info = proj_pj_info(m_pj);
1914 4 : CPLDebug("OGRCT", "%s (%s)", info.definition, osName.c_str());
1915 : #endif
1916 4 : return true;
1917 : }
1918 :
1919 : // Create a geographic 2D long-lat degrees CRS that is related to the
1920 : // source CRS
1921 0 : auto geodetic_crs = proj_crs_get_geodetic_crs(ctx, src);
1922 0 : if (!geodetic_crs)
1923 : {
1924 0 : proj_list_destroy(op_list);
1925 0 : proj_operation_factory_context_destroy(operation_ctx);
1926 0 : proj_destroy(src);
1927 0 : proj_destroy(dst);
1928 0 : CPLDebug("OGRCT", "Cannot find geodetic CRS matching source CRS");
1929 0 : return false;
1930 : }
1931 0 : auto geodetic_crs_type = proj_get_type(geodetic_crs);
1932 0 : if (geodetic_crs_type == PJ_TYPE_GEOCENTRIC_CRS ||
1933 0 : geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_2D_CRS ||
1934 : geodetic_crs_type == PJ_TYPE_GEOGRAPHIC_3D_CRS)
1935 : {
1936 0 : auto datum = proj_crs_get_datum(ctx, geodetic_crs);
1937 : #if PROJ_VERSION_MAJOR > 7 || \
1938 : (PROJ_VERSION_MAJOR == 7 && PROJ_VERSION_MINOR >= 2)
1939 : if (datum == nullptr)
1940 : {
1941 : datum = proj_crs_get_datum_forced(ctx, geodetic_crs);
1942 : }
1943 : #endif
1944 0 : if (datum)
1945 : {
1946 0 : auto ellps = proj_get_ellipsoid(ctx, datum);
1947 0 : proj_destroy(datum);
1948 0 : double semi_major_metre = 0;
1949 0 : double inv_flattening = 0;
1950 0 : proj_ellipsoid_get_parameters(ctx, ellps, &semi_major_metre,
1951 : nullptr, nullptr, &inv_flattening);
1952 0 : auto cs = proj_create_ellipsoidal_2D_cs(
1953 : ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE, nullptr, 0);
1954 : // It is critical to set the prime meridian to 0
1955 0 : auto temp = proj_create_geographic_crs(
1956 : ctx, "unnamed crs", "unnamed datum", proj_get_name(ellps),
1957 : semi_major_metre, inv_flattening, "Reference prime meridian", 0,
1958 : nullptr, 0, cs);
1959 0 : proj_destroy(ellps);
1960 0 : proj_destroy(cs);
1961 0 : proj_destroy(geodetic_crs);
1962 0 : geodetic_crs = temp;
1963 0 : geodetic_crs_type = proj_get_type(geodetic_crs);
1964 : }
1965 : }
1966 0 : if (geodetic_crs_type != PJ_TYPE_GEOGRAPHIC_2D_CRS)
1967 : {
1968 : // Shouldn't happen
1969 0 : proj_list_destroy(op_list);
1970 0 : proj_operation_factory_context_destroy(operation_ctx);
1971 0 : proj_destroy(src);
1972 0 : proj_destroy(dst);
1973 0 : proj_destroy(geodetic_crs);
1974 0 : CPLDebug("OGRCT", "Cannot find geographic CRS matching source CRS");
1975 0 : return false;
1976 : }
1977 :
1978 : // Create the transformation from this geographic 2D CRS to the source CRS
1979 : auto op_list_to_geodetic =
1980 0 : proj_create_operations(ctx, geodetic_crs, src, operation_ctx);
1981 0 : proj_destroy(geodetic_crs);
1982 :
1983 0 : if (op_list_to_geodetic == nullptr ||
1984 0 : proj_list_get_count(op_list_to_geodetic) == 0)
1985 : {
1986 0 : CPLDebug(
1987 : "OGRCT",
1988 : "Cannot compute transformation from geographic CRS to source CRS");
1989 0 : proj_list_destroy(op_list);
1990 0 : proj_list_destroy(op_list_to_geodetic);
1991 0 : proj_operation_factory_context_destroy(operation_ctx);
1992 0 : proj_destroy(src);
1993 0 : proj_destroy(dst);
1994 0 : return false;
1995 : }
1996 0 : auto opGeogToSrc = proj_list_get(ctx, op_list_to_geodetic, 0);
1997 0 : CPLAssert(opGeogToSrc);
1998 0 : proj_list_destroy(op_list_to_geodetic);
1999 0 : auto pjGeogToSrc = op_to_pj(ctx, opGeogToSrc);
2000 0 : proj_destroy(opGeogToSrc);
2001 0 : if (!pjGeogToSrc)
2002 : {
2003 0 : proj_list_destroy(op_list);
2004 0 : proj_operation_factory_context_destroy(operation_ctx);
2005 0 : proj_destroy(src);
2006 0 : proj_destroy(dst);
2007 0 : return false;
2008 : }
2009 :
2010 : const auto addTransformation =
2011 0 : [this, &pjGeogToSrc, &ctx](PJ *op, double west_lon, double south_lat,
2012 0 : double east_lon, double north_lat)
2013 : {
2014 0 : double minx = -std::numeric_limits<double>::max();
2015 0 : double miny = -std::numeric_limits<double>::max();
2016 0 : double maxx = std::numeric_limits<double>::max();
2017 0 : double maxy = std::numeric_limits<double>::max();
2018 :
2019 0 : if (!(west_lon == -180.0 && east_lon == 180.0 && south_lat == -90.0 &&
2020 : north_lat == 90.0))
2021 : {
2022 0 : minx = -minx;
2023 0 : miny = -miny;
2024 0 : maxx = -maxx;
2025 0 : maxy = -maxy;
2026 :
2027 : double x[21 * 4], y[21 * 4];
2028 0 : for (int j = 0; j <= 20; j++)
2029 : {
2030 0 : x[j] = west_lon + j * (east_lon - west_lon) / 20;
2031 0 : y[j] = south_lat;
2032 0 : x[21 + j] = west_lon + j * (east_lon - west_lon) / 20;
2033 0 : y[21 + j] = north_lat;
2034 0 : x[21 * 2 + j] = west_lon;
2035 0 : y[21 * 2 + j] = south_lat + j * (north_lat - south_lat) / 20;
2036 0 : x[21 * 3 + j] = east_lon;
2037 0 : y[21 * 3 + j] = south_lat + j * (north_lat - south_lat) / 20;
2038 : }
2039 0 : proj_trans_generic(pjGeogToSrc, PJ_FWD, x, sizeof(double), 21 * 4,
2040 : y, sizeof(double), 21 * 4, nullptr, 0, 0,
2041 : nullptr, 0, 0);
2042 0 : for (int j = 0; j < 21 * 4; j++)
2043 : {
2044 0 : if (x[j] != HUGE_VAL && y[j] != HUGE_VAL)
2045 : {
2046 0 : minx = std::min(minx, x[j]);
2047 0 : miny = std::min(miny, y[j]);
2048 0 : maxx = std::max(maxx, x[j]);
2049 0 : maxy = std::max(maxy, y[j]);
2050 : }
2051 : }
2052 : }
2053 :
2054 0 : if (minx <= maxx)
2055 : {
2056 0 : CPLString osProjString;
2057 0 : const double accuracy = proj_coordoperation_get_accuracy(ctx, op);
2058 0 : auto pj = op_to_pj(ctx, op, &osProjString);
2059 0 : CPLString osName;
2060 0 : auto name = proj_get_name(op);
2061 0 : if (name)
2062 0 : osName = name;
2063 0 : proj_destroy(op);
2064 0 : op = nullptr;
2065 0 : if (pj)
2066 : {
2067 : m_oTransformations.emplace_back(minx, miny, maxx, maxy, pj,
2068 0 : osName, osProjString, accuracy);
2069 : }
2070 : }
2071 0 : return op;
2072 0 : };
2073 :
2074 : // Iterate over source->target candidate transformations and reproject
2075 : // their long-lat bounding box into the source CRS.
2076 0 : bool foundWorldTransformation = false;
2077 0 : for (int i = 0; i < op_count; i++)
2078 : {
2079 0 : auto op = proj_list_get(ctx, op_list, i);
2080 0 : CPLAssert(op);
2081 0 : double west_lon = 0.0;
2082 0 : double south_lat = 0.0;
2083 0 : double east_lon = 0.0;
2084 0 : double north_lat = 0.0;
2085 0 : if (proj_get_area_of_use(ctx, op, &west_lon, &south_lat, &east_lon,
2086 0 : &north_lat, nullptr))
2087 : {
2088 0 : if (west_lon <= east_lon)
2089 : {
2090 0 : if (west_lon == -180 && east_lon == 180 && south_lat == -90 &&
2091 0 : north_lat == 90)
2092 : {
2093 0 : foundWorldTransformation = true;
2094 : }
2095 0 : op = addTransformation(op, west_lon, south_lat, east_lon,
2096 : north_lat);
2097 : }
2098 : else
2099 : {
2100 0 : auto op_clone = proj_clone(ctx, op);
2101 :
2102 0 : op = addTransformation(op, west_lon, south_lat, 180, north_lat);
2103 0 : op_clone = addTransformation(op_clone, -180, south_lat,
2104 : east_lon, north_lat);
2105 0 : proj_destroy(op_clone);
2106 : }
2107 : }
2108 :
2109 0 : proj_destroy(op);
2110 : }
2111 :
2112 0 : proj_list_destroy(op_list);
2113 :
2114 : // Sometimes the user will operate even outside the area of use of the
2115 : // source and target CRS, so if no global transformation has been returned
2116 : // previously, trigger the computation of one.
2117 0 : if (!foundWorldTransformation)
2118 : {
2119 0 : proj_operation_factory_context_set_area_of_interest(ctx, operation_ctx,
2120 : -180, -90, 180, 90);
2121 0 : proj_operation_factory_context_set_spatial_criterion(
2122 : ctx, operation_ctx, PROJ_SPATIAL_CRITERION_STRICT_CONTAINMENT);
2123 0 : op_list = proj_create_operations(ctx, src, dst, operation_ctx);
2124 0 : if (op_list)
2125 : {
2126 0 : op_count = proj_list_get_count(op_list);
2127 0 : for (int i = 0; i < op_count; i++)
2128 : {
2129 0 : auto op = proj_list_get(ctx, op_list, i);
2130 0 : CPLAssert(op);
2131 0 : double west_lon = 0.0;
2132 0 : double south_lat = 0.0;
2133 0 : double east_lon = 0.0;
2134 0 : double north_lat = 0.0;
2135 0 : if (proj_get_area_of_use(ctx, op, &west_lon, &south_lat,
2136 0 : &east_lon, &north_lat, nullptr) &&
2137 0 : west_lon == -180 && east_lon == 180 && south_lat == -90 &&
2138 0 : north_lat == 90)
2139 : {
2140 0 : op = addTransformation(op, west_lon, south_lat, east_lon,
2141 : north_lat);
2142 : }
2143 0 : proj_destroy(op);
2144 : }
2145 : }
2146 0 : proj_list_destroy(op_list);
2147 : }
2148 :
2149 0 : proj_operation_factory_context_destroy(operation_ctx);
2150 0 : proj_destroy(src);
2151 0 : proj_destroy(dst);
2152 0 : proj_destroy(pjGeogToSrc);
2153 0 : return !m_oTransformations.empty();
2154 : }
2155 :
2156 : /************************************************************************/
2157 : /* GetSourceCS() */
2158 : /************************************************************************/
2159 :
2160 3989 : const OGRSpatialReference *OGRProjCT::GetSourceCS() const
2161 :
2162 : {
2163 3989 : return poSRSSource;
2164 : }
2165 :
2166 : /************************************************************************/
2167 : /* GetTargetCS() */
2168 : /************************************************************************/
2169 :
2170 2266 : const OGRSpatialReference *OGRProjCT::GetTargetCS() const
2171 :
2172 : {
2173 2266 : return poSRSTarget;
2174 : }
2175 :
2176 : /************************************************************************/
2177 : /* Transform() */
2178 : /************************************************************************/
2179 :
2180 200233 : int OGRCoordinateTransformation::Transform(size_t nCount, double *x, double *y,
2181 : double *z, int *pabSuccessIn)
2182 :
2183 : {
2184 : int *pabSuccess =
2185 : pabSuccessIn
2186 200233 : ? pabSuccessIn
2187 198988 : : static_cast<int *>(VSI_MALLOC2_VERBOSE(sizeof(int), nCount));
2188 200233 : if (!pabSuccess)
2189 0 : return FALSE;
2190 :
2191 : bool bOverallSuccess =
2192 200233 : CPL_TO_BOOL(Transform(nCount, x, y, z, nullptr, pabSuccess));
2193 :
2194 424104 : for (size_t i = 0; i < nCount; i++)
2195 : {
2196 223939 : if (!pabSuccess[i])
2197 : {
2198 68 : bOverallSuccess = false;
2199 68 : break;
2200 : }
2201 : }
2202 :
2203 200233 : if (pabSuccess != pabSuccessIn)
2204 198988 : CPLFree(pabSuccess);
2205 :
2206 200233 : return bOverallSuccess;
2207 : }
2208 :
2209 : /************************************************************************/
2210 : /* TransformWithErrorCodes() */
2211 : /************************************************************************/
2212 :
2213 3976 : int OGRCoordinateTransformation::TransformWithErrorCodes(size_t nCount,
2214 : double *x, double *y,
2215 : double *z, double *t,
2216 : int *panErrorCodes)
2217 :
2218 : {
2219 3976 : if (nCount == 1)
2220 : {
2221 3976 : int nSuccess = 0;
2222 : const bool bOverallSuccess =
2223 3976 : CPL_TO_BOOL(Transform(nCount, x, y, z, t, &nSuccess));
2224 3976 : if (panErrorCodes)
2225 : {
2226 3976 : panErrorCodes[0] = nSuccess ? 0 : -1;
2227 : }
2228 3976 : return bOverallSuccess;
2229 : }
2230 :
2231 0 : std::vector<int> abSuccess;
2232 : try
2233 : {
2234 0 : abSuccess.resize(nCount);
2235 : }
2236 0 : catch (const std::bad_alloc &)
2237 : {
2238 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
2239 : "Cannot allocate abSuccess[] temporary array");
2240 0 : return FALSE;
2241 : }
2242 :
2243 : const bool bOverallSuccess =
2244 0 : CPL_TO_BOOL(Transform(nCount, x, y, z, t, abSuccess.data()));
2245 :
2246 0 : if (panErrorCodes)
2247 : {
2248 0 : for (size_t i = 0; i < nCount; i++)
2249 : {
2250 0 : panErrorCodes[i] = abSuccess[i] ? 0 : -1;
2251 : }
2252 : }
2253 :
2254 0 : return bOverallSuccess;
2255 : }
2256 :
2257 : /************************************************************************/
2258 : /* Transform() */
2259 : /************************************************************************/
2260 :
2261 1806810 : int OGRProjCT::Transform(size_t nCount, double *x, double *y, double *z,
2262 : double *t, int *pabSuccess)
2263 :
2264 : {
2265 : bool bOverallSuccess =
2266 1806810 : CPL_TO_BOOL(TransformWithErrorCodes(nCount, x, y, z, t, pabSuccess));
2267 :
2268 1806800 : if (pabSuccess)
2269 : {
2270 16237400 : for (size_t i = 0; i < nCount; i++)
2271 : {
2272 14430600 : pabSuccess[i] = (pabSuccess[i] == 0);
2273 : }
2274 : }
2275 :
2276 1806800 : return bOverallSuccess;
2277 : }
2278 :
2279 : /************************************************************************/
2280 : /* TransformWithErrorCodes() */
2281 : /************************************************************************/
2282 :
2283 : #ifndef PROJ_ERR_COORD_TRANSFM_INVALID_COORD
2284 : #define PROJ_ERR_COORD_TRANSFM_INVALID_COORD 2049
2285 : #define PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN 2050
2286 : #define PROJ_ERR_COORD_TRANSFM_NO_OPERATION 2051
2287 : #endif
2288 :
2289 1818940 : int OGRProjCT::TransformWithErrorCodes(size_t nCount, double *x, double *y,
2290 : double *z, double *t, int *panErrorCodes)
2291 :
2292 : {
2293 1818940 : if (nCount == 0)
2294 0 : return TRUE;
2295 :
2296 : // Prevent any coordinate modification when possible
2297 1818940 : if (bNoTransform)
2298 : {
2299 1628 : if (panErrorCodes)
2300 : {
2301 4516 : for (size_t i = 0; i < nCount; i++)
2302 : {
2303 3184 : panErrorCodes[i] = 0;
2304 : }
2305 : }
2306 1628 : return TRUE;
2307 : }
2308 :
2309 : #ifdef DEBUG_VERBOSE
2310 : bool bDebugCT = CPLTestBool(CPLGetConfigOption("OGR_CT_DEBUG", "NO"));
2311 : if (bDebugCT)
2312 : {
2313 : CPLDebug("OGRCT", "count = %d", nCount);
2314 : for (size_t i = 0; i < nCount; ++i)
2315 : {
2316 : CPLDebug("OGRCT", " x[%d] = %.16g y[%d] = %.16g", int(i), x[i],
2317 : int(i), y[i]);
2318 : }
2319 : }
2320 : #endif
2321 : #ifdef DEBUG_PERF
2322 : // CPLDebug("OGR_CT", "Begin TransformWithErrorCodes()");
2323 : struct CPLTimeVal tvStart;
2324 : CPLGettimeofday(&tvStart, nullptr);
2325 : #endif
2326 :
2327 : /* -------------------------------------------------------------------- */
2328 : /* Apply data axis to source CRS mapping. */
2329 : /* -------------------------------------------------------------------- */
2330 :
2331 : // Since we may swap the x and y pointers, but cannot tell the caller about this swap,
2332 : // we save the original pointer. The same axis swap code is executed for poSRSTarget.
2333 : // If this nullifies, we save the swap of both axes
2334 1817310 : const auto xOriginal = x;
2335 :
2336 1817310 : if (poSRSSource)
2337 : {
2338 1816400 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
2339 1816400 : if (mapping.size() >= 2)
2340 : {
2341 1816400 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
2342 : {
2343 1549110 : std::swap(x, y);
2344 : }
2345 1816400 : const bool bNegateX = mapping[0] < 0;
2346 1816400 : if (bNegateX)
2347 : {
2348 2 : for (size_t i = 0; i < nCount; i++)
2349 : {
2350 1 : x[i] = -x[i];
2351 : }
2352 : }
2353 1816400 : const bool bNegateY = mapping[1] < 0;
2354 1816400 : if (bNegateY)
2355 : {
2356 0 : for (size_t i = 0; i < nCount; i++)
2357 : {
2358 0 : y[i] = -y[i];
2359 : }
2360 : }
2361 1816400 : if (z && mapping.size() >= 3 && mapping[2] == -3)
2362 : {
2363 2 : for (size_t i = 0; i < nCount; i++)
2364 : {
2365 1 : z[i] = -z[i];
2366 : }
2367 : }
2368 : }
2369 : }
2370 :
2371 : /* -------------------------------------------------------------------- */
2372 : /* Potentially do longitude wrapping. */
2373 : /* -------------------------------------------------------------------- */
2374 1817310 : if (bSourceLatLong && bSourceWrap)
2375 : {
2376 190211 : if (m_eSourceFirstAxisOrient == OAO_East)
2377 : {
2378 7206 : for (size_t i = 0; i < nCount; i++)
2379 : {
2380 6130 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2381 : {
2382 6130 : if (x[i] < dfSourceWrapLong - 180.0)
2383 0 : x[i] += 360.0;
2384 6130 : else if (x[i] > dfSourceWrapLong + 180)
2385 0 : x[i] -= 360.0;
2386 : }
2387 : }
2388 : }
2389 : else
2390 : {
2391 477143 : for (size_t i = 0; i < nCount; i++)
2392 : {
2393 288008 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2394 : {
2395 286928 : if (y[i] < dfSourceWrapLong - 180.0)
2396 225 : y[i] += 360.0;
2397 286703 : else if (y[i] > dfSourceWrapLong + 180)
2398 71 : y[i] -= 360.0;
2399 : }
2400 : }
2401 : }
2402 : }
2403 :
2404 : /* -------------------------------------------------------------------- */
2405 : /* Optimized transform from WebMercator to WGS84 */
2406 : /* -------------------------------------------------------------------- */
2407 1817310 : bool bTransformDone = false;
2408 1817310 : if (bWebMercatorToWGS84LongLat)
2409 : {
2410 27495 : constexpr double REVERSE_SPHERE_RADIUS = 1.0 / 6378137.0;
2411 :
2412 27495 : if (m_eSourceFirstAxisOrient != OAO_East)
2413 : {
2414 0 : std::swap(x, y);
2415 : }
2416 :
2417 27495 : double y0 = y[0];
2418 625541 : for (size_t i = 0; i < nCount; i++)
2419 : {
2420 598046 : if (x[i] != HUGE_VAL)
2421 : {
2422 598046 : x[i] = x[i] * REVERSE_SPHERE_RADIUS;
2423 598046 : if (x[i] > M_PI)
2424 : {
2425 783 : if (x[i] < M_PI + 1e-14)
2426 : {
2427 138 : x[i] = M_PI;
2428 : }
2429 645 : else if (m_options.d->bCheckWithInvertProj)
2430 : {
2431 0 : x[i] = HUGE_VAL;
2432 0 : y[i] = HUGE_VAL;
2433 0 : y0 = HUGE_VAL;
2434 0 : continue;
2435 : }
2436 : else
2437 : {
2438 2 : do
2439 : {
2440 647 : x[i] -= 2 * M_PI;
2441 647 : } while (x[i] > M_PI);
2442 : }
2443 : }
2444 597263 : else if (x[i] < -M_PI)
2445 : {
2446 697 : if (x[i] > -M_PI - 1e-14)
2447 : {
2448 8 : x[i] = -M_PI;
2449 : }
2450 689 : else if (m_options.d->bCheckWithInvertProj)
2451 : {
2452 0 : x[i] = HUGE_VAL;
2453 0 : y[i] = HUGE_VAL;
2454 0 : y0 = HUGE_VAL;
2455 0 : continue;
2456 : }
2457 : else
2458 : {
2459 0 : do
2460 : {
2461 689 : x[i] += 2 * M_PI;
2462 689 : } while (x[i] < -M_PI);
2463 : }
2464 : }
2465 598046 : constexpr double RAD_TO_DEG = 180. / M_PI;
2466 598046 : x[i] *= RAD_TO_DEG;
2467 :
2468 : // Optimization for the case where we are provided a whole line
2469 : // of same northing.
2470 598046 : if (i > 0 && y[i] == y0)
2471 400171 : y[i] = y[0];
2472 : else
2473 : {
2474 197875 : y[i] = M_PI / 2.0 -
2475 197875 : 2.0 * atan(exp(-y[i] * REVERSE_SPHERE_RADIUS));
2476 197875 : y[i] *= RAD_TO_DEG;
2477 : }
2478 : }
2479 : }
2480 :
2481 27495 : if (panErrorCodes)
2482 : {
2483 625535 : for (size_t i = 0; i < nCount; i++)
2484 : {
2485 598042 : if (x[i] != HUGE_VAL)
2486 598042 : panErrorCodes[i] = 0;
2487 : else
2488 0 : panErrorCodes[i] =
2489 : PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2490 : }
2491 : }
2492 :
2493 27495 : if (m_eTargetFirstAxisOrient != OAO_East)
2494 : {
2495 27008 : std::swap(x, y);
2496 : }
2497 :
2498 27495 : bTransformDone = true;
2499 : }
2500 :
2501 : // Determine the default coordinate epoch, if not provided in the point to
2502 : // transform.
2503 : // For time-dependent transformations, PROJ can currently only do
2504 : // staticCRS -> dynamicCRS or dynamicCRS -> staticCRS transformations, and
2505 : // in either case, the coordinate epoch of the dynamicCRS must be provided
2506 : // as the input time.
2507 1817310 : double dfDefaultTime = HUGE_VAL;
2508 1766710 : if (bSourceIsDynamicCRS && dfSourceCoordinateEpoch > 0 &&
2509 3584030 : !bTargetIsDynamicCRS &&
2510 4 : CPLTestBool(
2511 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2512 : {
2513 4 : dfDefaultTime = dfSourceCoordinateEpoch;
2514 4 : CPLDebug("OGR_CT", "Using coordinate epoch %f from source CRS",
2515 : dfDefaultTime);
2516 : }
2517 1558880 : else if (bTargetIsDynamicCRS && dfTargetCoordinateEpoch > 0 &&
2518 3376180 : !bSourceIsDynamicCRS &&
2519 2 : CPLTestBool(
2520 : CPLGetConfigOption("OGR_CT_USE_SRS_COORDINATE_EPOCH", "YES")))
2521 : {
2522 2 : dfDefaultTime = dfTargetCoordinateEpoch;
2523 2 : CPLDebug("OGR_CT", "Using coordinate epoch %f from target CRS",
2524 : dfDefaultTime);
2525 : }
2526 :
2527 : /* -------------------------------------------------------------------- */
2528 : /* Select dynamically the best transformation for the data, if */
2529 : /* needed. */
2530 : /* -------------------------------------------------------------------- */
2531 1817310 : PJ_CONTEXT *ctx = m_psLastContext;
2532 1817310 : const auto nThisThreadId = std::this_thread::get_id();
2533 1817310 : if (!ctx || nThisThreadId != m_nLastContextThreadId)
2534 : {
2535 2680 : m_nLastContextThreadId = nThisThreadId;
2536 2680 : m_psLastContext = OSRGetProjTLSContext();
2537 2680 : ctx = m_psLastContext;
2538 : }
2539 :
2540 1817300 : PJ *pj = m_pj;
2541 1817300 : if (!bTransformDone && !pj)
2542 : {
2543 0 : double avgX = 0.0;
2544 0 : double avgY = 0.0;
2545 0 : size_t nCountValid = 0;
2546 0 : for (size_t i = 0; i < nCount; i++)
2547 : {
2548 0 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2549 : {
2550 0 : avgX += x[i];
2551 0 : avgY += y[i];
2552 0 : nCountValid++;
2553 : }
2554 : }
2555 0 : if (nCountValid != 0)
2556 : {
2557 0 : avgX /= static_cast<double>(nCountValid);
2558 0 : avgY /= static_cast<double>(nCountValid);
2559 : }
2560 :
2561 0 : constexpr int N_MAX_RETRY = 2;
2562 0 : int iExcluded[N_MAX_RETRY] = {-1, -1};
2563 :
2564 0 : const int nOperations = static_cast<int>(m_oTransformations.size());
2565 : PJ_COORD coord;
2566 0 : coord.xyzt.x = avgX;
2567 0 : coord.xyzt.y = avgY;
2568 0 : coord.xyzt.z = z ? z[0] : 0;
2569 0 : coord.xyzt.t = t ? t[0] : dfDefaultTime;
2570 :
2571 : // We may need several attempts. For example the point at
2572 : // lon=-111.5 lat=45.26 falls into the bounding box of the Canadian
2573 : // ntv2_0.gsb grid, except that it is not in any of the subgrids, being
2574 : // in the US. We thus need another retry that will select the conus
2575 : // grid.
2576 0 : for (int iRetry = 0; iRetry <= N_MAX_RETRY; iRetry++)
2577 : {
2578 0 : int iBestTransf = -1;
2579 : // Select transform whose BBOX match our data and has the best
2580 : // accuracy if m_eStrategy == BEST_ACCURACY. Or just the first BBOX
2581 : // matching one, if
2582 : // m_eStrategy == FIRST_MATCHING
2583 0 : double dfBestAccuracy = std::numeric_limits<double>::infinity();
2584 0 : for (int i = 0; i < nOperations; i++)
2585 : {
2586 0 : if (i == iExcluded[0] || i == iExcluded[1])
2587 : {
2588 0 : continue;
2589 : }
2590 0 : const auto &transf = m_oTransformations[i];
2591 0 : if (avgX >= transf.minx && avgX <= transf.maxx &&
2592 0 : avgY >= transf.miny && avgY <= transf.maxy &&
2593 0 : (iBestTransf < 0 || (transf.accuracy >= 0 &&
2594 0 : transf.accuracy < dfBestAccuracy)))
2595 : {
2596 0 : iBestTransf = i;
2597 0 : dfBestAccuracy = transf.accuracy;
2598 0 : if (m_eStrategy == Strategy::FIRST_MATCHING)
2599 0 : break;
2600 : }
2601 : }
2602 0 : if (iBestTransf < 0)
2603 : {
2604 0 : break;
2605 : }
2606 0 : auto &transf = m_oTransformations[iBestTransf];
2607 0 : pj = transf.pj;
2608 0 : proj_assign_context(pj, ctx);
2609 0 : if (iBestTransf != m_iCurTransformation)
2610 : {
2611 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2612 : transf.osProjString.c_str(), transf.osName.c_str());
2613 0 : m_iCurTransformation = iBestTransf;
2614 : }
2615 :
2616 0 : auto res = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2617 0 : if (res.xyzt.x != HUGE_VAL)
2618 : {
2619 0 : break;
2620 : }
2621 0 : pj = nullptr;
2622 0 : CPLDebug("OGRCT", "Did not result in valid result. "
2623 : "Attempting a retry with another operation.");
2624 0 : if (iRetry == N_MAX_RETRY)
2625 : {
2626 0 : break;
2627 : }
2628 0 : iExcluded[iRetry] = iBestTransf;
2629 : }
2630 :
2631 0 : if (!pj)
2632 : {
2633 : // In case we did not find an operation whose area of use is
2634 : // compatible with the input coordinate, then goes through again the
2635 : // list, and use the first operation that does not require grids.
2636 0 : for (int i = 0; i < nOperations; i++)
2637 : {
2638 0 : auto &transf = m_oTransformations[i];
2639 0 : if (proj_coordoperation_get_grid_used_count(ctx, transf.pj) ==
2640 : 0)
2641 : {
2642 0 : pj = transf.pj;
2643 0 : proj_assign_context(pj, ctx);
2644 0 : if (i != m_iCurTransformation)
2645 : {
2646 0 : CPLDebug("OGRCT", "Selecting transformation %s (%s)",
2647 : transf.osProjString.c_str(),
2648 : transf.osName.c_str());
2649 0 : m_iCurTransformation = i;
2650 : }
2651 0 : break;
2652 : }
2653 : }
2654 : }
2655 :
2656 0 : if (!pj)
2657 : {
2658 0 : if (m_bEmitErrors && ++nErrorCount < 20)
2659 : {
2660 0 : CPLError(CE_Failure, CPLE_AppDefined,
2661 : "Cannot find transformation for provided coordinates");
2662 : }
2663 0 : else if (nErrorCount == 20)
2664 : {
2665 0 : CPLError(CE_Failure, CPLE_AppDefined,
2666 : "Reprojection failed, further errors will be "
2667 : "suppressed on the transform object.");
2668 : }
2669 :
2670 0 : for (size_t i = 0; i < nCount; i++)
2671 : {
2672 0 : x[i] = HUGE_VAL;
2673 0 : y[i] = HUGE_VAL;
2674 0 : if (panErrorCodes)
2675 0 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_NO_OPERATION;
2676 : }
2677 0 : return FALSE;
2678 : }
2679 : }
2680 1817300 : if (pj)
2681 : {
2682 1816570 : proj_assign_context(pj, ctx);
2683 : }
2684 :
2685 : /* -------------------------------------------------------------------- */
2686 : /* Do the transformation (or not...) using PROJ */
2687 : /* -------------------------------------------------------------------- */
2688 :
2689 1817290 : if (!bTransformDone)
2690 : {
2691 1789810 : const auto nLastErrorCounter = CPLGetErrorCounter();
2692 :
2693 15638700 : for (size_t i = 0; i < nCount; i++)
2694 : {
2695 : PJ_COORD coord;
2696 13848800 : const double xIn = x[i];
2697 13848800 : const double yIn = y[i];
2698 13848800 : if (!std::isfinite(xIn))
2699 : {
2700 10247 : x[i] = HUGE_VAL;
2701 10247 : y[i] = HUGE_VAL;
2702 10247 : if (panErrorCodes)
2703 10247 : panErrorCodes[i] = PROJ_ERR_COORD_TRANSFM_INVALID_COORD;
2704 10247 : continue;
2705 : }
2706 13838600 : coord.xyzt.x = x[i];
2707 13838600 : coord.xyzt.y = y[i];
2708 13838600 : coord.xyzt.z = z ? z[i] : 0;
2709 13838600 : coord.xyzt.t = t ? t[i] : dfDefaultTime;
2710 13838600 : proj_errno_reset(pj);
2711 13838600 : coord = proj_trans(pj, m_bReversePj ? PJ_INV : PJ_FWD, coord);
2712 : #if 0
2713 : CPLDebug("OGRCT",
2714 : "Transforming (x=%f,y=%f,z=%f,time=%f) to "
2715 : "(x=%f,y=%f,z=%f,time=%f)",
2716 : x[i], y[i], z ? z[i] : 0, t ? t[i] : dfDefaultTime,
2717 : coord.xyzt.x, coord.xyzt.y, coord.xyzt.z, coord.xyzt.t);
2718 : #endif
2719 13838600 : x[i] = coord.xyzt.x;
2720 13838600 : y[i] = coord.xyzt.y;
2721 13838600 : if (z)
2722 13801800 : z[i] = coord.xyzt.z;
2723 13838600 : if (t)
2724 25 : t[i] = coord.xyzt.t;
2725 13838600 : int err = 0;
2726 13838600 : if (std::isnan(coord.xyzt.x))
2727 : {
2728 : // This shouldn't normally happen if PROJ projections behave
2729 : // correctly, but e.g inverse laea before PROJ 8.1.1 could
2730 : // do that for points out of domain.
2731 : // See https://github.com/OSGeo/PROJ/pull/2800
2732 70 : x[i] = HUGE_VAL;
2733 70 : y[i] = HUGE_VAL;
2734 70 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2735 :
2736 : #ifdef DEBUG
2737 70 : CPLErrorOnce(CE_Warning, CPLE_AppDefined,
2738 : "PROJ returned a NaN value. It should be fixed");
2739 : #else
2740 : CPLDebugOnce("OGR_CT",
2741 : "PROJ returned a NaN value. It should be fixed");
2742 : #endif
2743 : }
2744 13838500 : else if (coord.xyzt.x == HUGE_VAL)
2745 : {
2746 3354040 : err = proj_errno(pj);
2747 : // PROJ should normally emit an error, but in case it does not
2748 : // (e.g PROJ 6.3 with the +ortho projection), synthetize one
2749 3354040 : if (err == 0)
2750 3337940 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2751 : }
2752 10484500 : else if (m_options.d->bCheckWithInvertProj)
2753 : {
2754 : // For some projections, we cannot detect if we are trying to
2755 : // reproject coordinates outside the validity area of the
2756 : // projection. So let's do the reverse reprojection and compare
2757 : // with the source coordinates.
2758 1754500 : coord = proj_trans(pj, m_bReversePj ? PJ_FWD : PJ_INV, coord);
2759 1754500 : if (fabs(coord.xyzt.x - xIn) > dfThreshold ||
2760 1602580 : fabs(coord.xyzt.y - yIn) > dfThreshold)
2761 : {
2762 154223 : err = PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN;
2763 154223 : x[i] = HUGE_VAL;
2764 154223 : y[i] = HUGE_VAL;
2765 : }
2766 : }
2767 :
2768 13838600 : if (panErrorCodes)
2769 13829900 : panErrorCodes[i] = err;
2770 :
2771 : /* --------------------------------------------------------------------
2772 : */
2773 : /* Try to report an error through CPL. Get proj error string
2774 : */
2775 : /* if possible. Try to avoid reporting thousands of errors. */
2776 : /* Suppress further error reporting on this OGRProjCT if we */
2777 : /* have already reported 20 errors. */
2778 : /* --------------------------------------------------------------------
2779 : */
2780 13838600 : if (err != 0)
2781 : {
2782 3508340 : if (++nErrorCount < 20)
2783 : {
2784 : #if PROJ_VERSION_MAJOR >= 8
2785 : const char *pszError = proj_context_errno_string(ctx, err);
2786 : #else
2787 1903 : const char *pszError = proj_errno_string(err);
2788 : #endif
2789 1903 : if (m_bEmitErrors
2790 : #ifdef PROJ_ERR_OTHER_NO_INVERSE_OP
2791 : || (i == 0 && err == PROJ_ERR_OTHER_NO_INVERSE_OP)
2792 : #endif
2793 : )
2794 : {
2795 345 : if (nLastErrorCounter != CPLGetErrorCounter() &&
2796 607 : CPLGetLastErrorType() == CE_Failure &&
2797 262 : strstr(CPLGetLastErrorMsg(), "PROJ:"))
2798 : {
2799 : // do nothing
2800 : }
2801 345 : else if (pszError == nullptr)
2802 0 : CPLError(CE_Failure, CPLE_AppDefined,
2803 : "Reprojection failed, err = %d", err);
2804 : else
2805 345 : CPLError(CE_Failure, CPLE_AppDefined, "%s",
2806 : pszError);
2807 : }
2808 : else
2809 : {
2810 1558 : if (pszError == nullptr)
2811 0 : CPLDebug("OGRCT", "Reprojection failed, err = %d",
2812 : err);
2813 : else
2814 1558 : CPLDebug("OGRCT", "%s", pszError);
2815 : }
2816 : }
2817 3506440 : else if (nErrorCount == 20)
2818 : {
2819 86 : if (m_bEmitErrors)
2820 : {
2821 14 : CPLError(CE_Failure, CPLE_AppDefined,
2822 : "Reprojection failed, err = %d, further "
2823 : "errors will be "
2824 : "suppressed on the transform object.",
2825 : err);
2826 : }
2827 : else
2828 : {
2829 72 : CPLDebug("OGRCT",
2830 : "Reprojection failed, err = %d, further "
2831 : "errors will be "
2832 : "suppressed on the transform object.",
2833 : err);
2834 : }
2835 : }
2836 : }
2837 : }
2838 : }
2839 :
2840 : /* -------------------------------------------------------------------- */
2841 : /* Potentially do longitude wrapping. */
2842 : /* -------------------------------------------------------------------- */
2843 1817290 : if (bTargetLatLong && bTargetWrap)
2844 : {
2845 491459 : if (m_eTargetFirstAxisOrient == OAO_East)
2846 : {
2847 1461250 : for (size_t i = 0; i < nCount; i++)
2848 : {
2849 1457210 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2850 : {
2851 1457210 : if (x[i] < dfTargetWrapLong - 180.0)
2852 0 : x[i] += 360.0;
2853 1457210 : else if (x[i] > dfTargetWrapLong + 180)
2854 154 : x[i] -= 360.0;
2855 : }
2856 : }
2857 : }
2858 : else
2859 : {
2860 9536450 : for (size_t i = 0; i < nCount; i++)
2861 : {
2862 9049030 : if (x[i] != HUGE_VAL && y[i] != HUGE_VAL)
2863 : {
2864 5628190 : if (y[i] < dfTargetWrapLong - 180.0)
2865 5823 : y[i] += 360.0;
2866 5622370 : else if (y[i] > dfTargetWrapLong + 180)
2867 23389 : y[i] -= 360.0;
2868 : }
2869 : }
2870 : }
2871 : }
2872 :
2873 : /* -------------------------------------------------------------------- */
2874 : /* Apply data axis to target CRS mapping. */
2875 : /* -------------------------------------------------------------------- */
2876 1817290 : if (poSRSTarget)
2877 : {
2878 1816400 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
2879 1816400 : if (mapping.size() >= 2)
2880 : {
2881 1816390 : const bool bNegateX = mapping[0] < 0;
2882 1816390 : if (bNegateX)
2883 : {
2884 2 : for (size_t i = 0; i < nCount; i++)
2885 : {
2886 1 : x[i] = -x[i];
2887 : }
2888 : }
2889 1816390 : const bool bNegateY = mapping[1] < 0;
2890 1816400 : if (bNegateY)
2891 : {
2892 0 : for (size_t i = 0; i < nCount; i++)
2893 : {
2894 0 : y[i] = -y[i];
2895 : }
2896 : }
2897 1816400 : if (z && mapping.size() >= 3 && mapping[2] == -3)
2898 : {
2899 2 : for (size_t i = 0; i < nCount; i++)
2900 : {
2901 1 : z[i] = -z[i];
2902 : }
2903 : }
2904 :
2905 1816400 : if (std::abs(mapping[0]) == 2 && std::abs(mapping[1]) == 1)
2906 : {
2907 698653 : std::swap(x, y);
2908 : }
2909 : }
2910 : }
2911 :
2912 : // Check whether final "genuine" axis swap is really necessary
2913 1817280 : if (x != xOriginal)
2914 : {
2915 1274260 : std::swap_ranges(x, x + nCount, y);
2916 : }
2917 :
2918 : #ifdef DEBUG_VERBOSE
2919 : if (bDebugCT)
2920 : {
2921 : CPLDebug("OGRCT", "Out:");
2922 : for (size_t i = 0; i < nCount; ++i)
2923 : {
2924 : CPLDebug("OGRCT", " x[%d] = %.16g y[%d] = %.16g", int(i), x[i],
2925 : int(i), y[i]);
2926 : }
2927 : }
2928 : #endif
2929 : #ifdef DEBUG_PERF
2930 : struct CPLTimeVal tvEnd;
2931 : CPLGettimeofday(&tvEnd, nullptr);
2932 : const double delay = (tvEnd.tv_sec + tvEnd.tv_usec * 1e-6) -
2933 : (tvStart.tv_sec + tvStart.tv_usec * 1e-6);
2934 : g_dfTotalTimeReprojection += delay;
2935 : // CPLDebug("OGR_CT", "End TransformWithErrorCodes(): %d ms",
2936 : // static_cast<int>(delay * 1000));
2937 : #endif
2938 :
2939 1817300 : return TRUE;
2940 : }
2941 :
2942 : /************************************************************************/
2943 : /* TransformBounds() */
2944 : /************************************************************************/
2945 :
2946 : // ---------------------------------------------------------------------------
2947 464 : static double simple_min(const double *data, const int arr_len)
2948 : {
2949 464 : double min_value = data[0];
2950 13064 : for (int iii = 1; iii < arr_len; iii++)
2951 : {
2952 12600 : if (data[iii] < min_value)
2953 2455 : min_value = data[iii];
2954 : }
2955 464 : return min_value;
2956 : }
2957 :
2958 : // ---------------------------------------------------------------------------
2959 464 : static double simple_max(const double *data, const int arr_len)
2960 : {
2961 464 : double max_value = data[0];
2962 12988 : for (int iii = 1; iii < arr_len; iii++)
2963 : {
2964 12524 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
2965 3057 : data[iii] != HUGE_VAL)
2966 1599 : max_value = data[iii];
2967 : }
2968 464 : return max_value;
2969 : }
2970 :
2971 : // ---------------------------------------------------------------------------
2972 3858 : static int _find_previous_index(const int iii, const double *data,
2973 : const int arr_len)
2974 : {
2975 : // find index of nearest valid previous value if exists
2976 3858 : int prev_iii = iii - 1;
2977 3858 : if (prev_iii == -1) // handle wraparound
2978 186 : prev_iii = arr_len - 1;
2979 3880 : while (data[prev_iii] == HUGE_VAL && prev_iii != iii)
2980 : {
2981 22 : prev_iii--;
2982 22 : if (prev_iii == -1) // handle wraparound
2983 2 : prev_iii = arr_len - 1;
2984 : }
2985 3858 : return prev_iii;
2986 : }
2987 :
2988 : // ---------------------------------------------------------------------------
2989 : /******************************************************************************
2990 : Handles the case when longitude values cross the antimeridian
2991 : when calculating the minimum.
2992 : Note: The data array must be in a linear ring.
2993 : Note: This requires a densified ring with at least 2 additional
2994 : points per edge to correctly handle global extents.
2995 : If only 1 additional point:
2996 : | |
2997 : |RL--x0--|RL--
2998 : | |
2999 : -180 180|-180
3000 : If they are evenly spaced and it crosses the antimeridian:
3001 : x0 - L = 180
3002 : R - x0 = -180
3003 : For example:
3004 : Let R = -179.9, x0 = 0.1, L = -179.89
3005 : x0 - L = 0.1 - -179.9 = 180
3006 : R - x0 = -179.89 - 0.1 ~= -180
3007 : This is the same in the case when it didn't cross the antimeridian.
3008 : If you have 2 additional points:
3009 : | |
3010 : |RL--x0--x1--|RL--
3011 : | |
3012 : -180 180|-180
3013 : If they are evenly spaced and it crosses the antimeridian:
3014 : x0 - L = 120
3015 : x1 - x0 = 120
3016 : R - x1 = -240
3017 : For example:
3018 : Let R = -179.9, x0 = -59.9, x1 = 60.1 L = -179.89
3019 : x0 - L = 59.9 - -179.9 = 120
3020 : x1 - x0 = 60.1 - 59.9 = 120
3021 : R - x1 = -179.89 - 60.1 ~= -240
3022 : However, if they are evenly spaced and it didn't cross the antimeridian:
3023 : x0 - L = 120
3024 : x1 - x0 = 120
3025 : R - x1 = 120
3026 : From this, we have a delta that is guaranteed to be significantly
3027 : large enough to tell the difference reguarless of the direction
3028 : the antimeridian was crossed.
3029 : However, even though the spacing was even in the source projection, it isn't
3030 : guaranteed in the target geographic projection. So, instead of 240, 200 is used
3031 : as it significantly larger than 120 to be sure that the antimeridian was crossed
3032 : but smalller than 240 to account for possible irregularities in distances
3033 : when re-projecting. Also, 200 ensures latitudes are ignored for axis order
3034 : handling.
3035 : ******************************************************************************/
3036 98 : static double antimeridian_min(const double *data, const int arr_len)
3037 : {
3038 98 : double positive_min = HUGE_VAL;
3039 98 : double min_value = HUGE_VAL;
3040 98 : int crossed_meridian_count = 0;
3041 98 : bool positive_meridian = false;
3042 :
3043 2086 : for (int iii = 0; iii < arr_len; iii++)
3044 : {
3045 1988 : if (data[iii] == HUGE_VAL)
3046 59 : continue;
3047 1929 : int prev_iii = _find_previous_index(iii, data, arr_len);
3048 : // check if crossed meridian
3049 1929 : double delta = data[prev_iii] - data[iii];
3050 : // 180 -> -180
3051 1929 : if (delta >= 200 && delta != HUGE_VAL)
3052 : {
3053 9 : if (crossed_meridian_count == 0)
3054 2 : positive_min = min_value;
3055 9 : crossed_meridian_count++;
3056 9 : positive_meridian = false;
3057 : // -180 -> 180
3058 : }
3059 1920 : else if (delta <= -200 && delta != HUGE_VAL)
3060 : {
3061 9 : if (crossed_meridian_count == 0)
3062 6 : positive_min = data[iii];
3063 9 : crossed_meridian_count++;
3064 9 : positive_meridian = true;
3065 : }
3066 : // positive meridian side min
3067 1929 : if (positive_meridian && data[iii] < positive_min)
3068 38 : positive_min = data[iii];
3069 : // track general min value
3070 1929 : if (data[iii] < min_value)
3071 165 : min_value = data[iii];
3072 : }
3073 :
3074 98 : if (crossed_meridian_count == 2)
3075 7 : return positive_min;
3076 91 : else if (crossed_meridian_count == 4)
3077 : // bounds extends beyond -180/180
3078 1 : return -180;
3079 90 : return min_value;
3080 : }
3081 :
3082 : // ---------------------------------------------------------------------------
3083 : // Handles the case when longitude values cross the antimeridian
3084 : // when calculating the minimum.
3085 : // Note: The data array must be in a linear ring.
3086 : // Note: This requires a densified ring with at least 2 additional
3087 : // points per edge to correctly handle global extents.
3088 : // See antimeridian_min docstring for reasoning.
3089 98 : static double antimeridian_max(const double *data, const int arr_len)
3090 : {
3091 98 : double negative_max = -HUGE_VAL;
3092 98 : double max_value = -HUGE_VAL;
3093 98 : bool negative_meridian = false;
3094 98 : int crossed_meridian_count = 0;
3095 :
3096 2086 : for (int iii = 0; iii < arr_len; iii++)
3097 : {
3098 1988 : if (data[iii] == HUGE_VAL)
3099 59 : continue;
3100 1929 : int prev_iii = _find_previous_index(iii, data, arr_len);
3101 : // check if crossed meridian
3102 1929 : double delta = data[prev_iii] - data[iii];
3103 : // 180 -> -180
3104 1929 : if (delta >= 200 && delta != HUGE_VAL)
3105 : {
3106 9 : if (crossed_meridian_count == 0)
3107 2 : negative_max = data[iii];
3108 9 : crossed_meridian_count++;
3109 9 : negative_meridian = true;
3110 : // -180 -> 180
3111 : }
3112 1920 : else if (delta <= -200 && delta != HUGE_VAL)
3113 : {
3114 9 : if (crossed_meridian_count == 0)
3115 6 : negative_max = max_value;
3116 9 : negative_meridian = false;
3117 9 : crossed_meridian_count++;
3118 : }
3119 : // negative meridian side max
3120 1929 : if (negative_meridian &&
3121 142 : (data[iii] > negative_max || negative_max == HUGE_VAL) &&
3122 45 : data[iii] != HUGE_VAL)
3123 45 : negative_max = data[iii];
3124 : // track general max value
3125 1929 : if ((data[iii] > max_value || max_value == HUGE_VAL) &&
3126 764 : data[iii] != HUGE_VAL)
3127 764 : max_value = data[iii];
3128 : }
3129 98 : if (crossed_meridian_count == 2)
3130 7 : return negative_max;
3131 91 : else if (crossed_meridian_count == 4)
3132 : // bounds extends beyond -180/180
3133 1 : return 180;
3134 90 : return max_value;
3135 : }
3136 :
3137 : // ---------------------------------------------------------------------------
3138 : // Check if the original projected bounds contains
3139 : // the north pole.
3140 : // This assumes that the destination CRS is geographic.
3141 106 : bool OGRProjCT::ContainsNorthPole(const double xmin, const double ymin,
3142 : const double xmax, const double ymax,
3143 : bool lon_lat_order)
3144 : {
3145 106 : double pole_y = 90;
3146 106 : double pole_x = 0;
3147 106 : if (!lon_lat_order)
3148 : {
3149 4 : pole_y = 0;
3150 4 : pole_x = 90;
3151 : }
3152 106 : auto inverseCT = GetInverse();
3153 106 : if (!inverseCT)
3154 0 : return false;
3155 212 : bool success = inverseCT->TransformWithErrorCodes(
3156 106 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr);
3157 106 : if (success && CPLGetLastErrorType() != CE_None)
3158 31 : CPLErrorReset();
3159 106 : delete inverseCT;
3160 106 : if (xmin < pole_x && pole_x < xmax && ymax > pole_y && pole_y > ymin)
3161 4 : return true;
3162 102 : return false;
3163 : }
3164 :
3165 : // ---------------------------------------------------------------------------
3166 : // Check if the original projected bounds contains
3167 : // the south pole.
3168 : // This assumes that the destination CRS is geographic.
3169 106 : bool OGRProjCT::ContainsSouthPole(const double xmin, const double ymin,
3170 : const double xmax, const double ymax,
3171 : bool lon_lat_order)
3172 : {
3173 106 : double pole_y = -90;
3174 106 : double pole_x = 0;
3175 106 : if (!lon_lat_order)
3176 : {
3177 4 : pole_y = 0;
3178 4 : pole_x = -90;
3179 : }
3180 106 : auto inverseCT = GetInverse();
3181 106 : if (!inverseCT)
3182 0 : return false;
3183 212 : bool success = inverseCT->TransformWithErrorCodes(
3184 106 : 1, &pole_x, &pole_y, nullptr, nullptr, nullptr);
3185 106 : if (success && CPLGetLastErrorType() != CE_None)
3186 30 : CPLErrorReset();
3187 106 : delete inverseCT;
3188 106 : if (xmin < pole_x && pole_x < xmax && ymax > pole_y && pole_y > ymin)
3189 5 : return true;
3190 101 : return false;
3191 : }
3192 :
3193 290 : int OGRProjCT::TransformBounds(const double xmin, const double ymin,
3194 : const double xmax, const double ymax,
3195 : double *out_xmin, double *out_ymin,
3196 : double *out_xmax, double *out_ymax,
3197 : const int densify_pts)
3198 : {
3199 290 : CPLErrorReset();
3200 :
3201 290 : if (bNoTransform)
3202 : {
3203 1 : *out_xmin = xmin;
3204 1 : *out_ymin = ymin;
3205 1 : *out_xmax = xmax;
3206 1 : *out_ymax = ymax;
3207 1 : return true;
3208 : }
3209 :
3210 289 : *out_xmin = HUGE_VAL;
3211 289 : *out_ymin = HUGE_VAL;
3212 289 : *out_xmax = HUGE_VAL;
3213 289 : *out_ymax = HUGE_VAL;
3214 :
3215 289 : if (densify_pts < 0 || densify_pts > 10000)
3216 : {
3217 1 : CPLError(CE_Failure, CPLE_AppDefined,
3218 : "densify_pts must be between 0-10000.");
3219 1 : return false;
3220 : }
3221 288 : if (!poSRSSource)
3222 : {
3223 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing source SRS.");
3224 0 : return false;
3225 : }
3226 288 : if (!poSRSTarget)
3227 : {
3228 0 : CPLError(CE_Failure, CPLE_AppDefined, "missing target SRS.");
3229 0 : return false;
3230 : }
3231 :
3232 288 : bool degree_input = false;
3233 288 : bool degree_output = false;
3234 288 : bool input_lon_lat_order = false;
3235 288 : bool output_lon_lat_order = false;
3236 :
3237 288 : if (bSourceLatLong)
3238 : {
3239 149 : degree_input = fabs(poSRSSource->GetAngularUnits(nullptr) -
3240 149 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3241 149 : const auto &mapping = poSRSSource->GetDataAxisToSRSAxisMapping();
3242 261 : if ((mapping[0] == 1 && m_eSourceFirstAxisOrient == OAO_East) ||
3243 112 : (mapping[0] == 2 && m_eSourceFirstAxisOrient != OAO_East))
3244 : {
3245 140 : input_lon_lat_order = true;
3246 : }
3247 : }
3248 288 : if (bTargetLatLong)
3249 : {
3250 107 : degree_output = fabs(poSRSTarget->GetAngularUnits(nullptr) -
3251 107 : CPLAtof(SRS_UA_DEGREE_CONV)) < 1e-8;
3252 107 : const auto &mapping = poSRSTarget->GetDataAxisToSRSAxisMapping();
3253 204 : if ((mapping[0] == 1 && m_eTargetFirstAxisOrient == OAO_East) ||
3254 97 : (mapping[0] == 2 && m_eTargetFirstAxisOrient != OAO_East))
3255 : {
3256 102 : output_lon_lat_order = true;
3257 : }
3258 : }
3259 :
3260 288 : if (degree_output && densify_pts < 2)
3261 : {
3262 1 : CPLError(CE_Failure, CPLE_AppDefined,
3263 : "densify_pts must be at least 2 if the output is geograpic.");
3264 1 : return false;
3265 : }
3266 :
3267 287 : int side_pts = densify_pts + 1; // add one because we are densifying
3268 287 : const int boundary_len = side_pts * 4;
3269 574 : std::vector<double> x_boundary_array;
3270 574 : std::vector<double> y_boundary_array;
3271 : try
3272 : {
3273 287 : x_boundary_array.resize(boundary_len);
3274 287 : y_boundary_array.resize(boundary_len);
3275 : }
3276 0 : catch (const std::exception &e) // memory allocation failure
3277 : {
3278 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s", e.what());
3279 0 : return false;
3280 : }
3281 287 : double delta_x = 0;
3282 287 : double delta_y = 0;
3283 287 : bool north_pole_in_bounds = false;
3284 287 : bool south_pole_in_bounds = false;
3285 287 : if (degree_output)
3286 : {
3287 106 : CPLErrorHandlerPusher oErrorHandlerPusher(CPLQuietErrorHandler);
3288 :
3289 : north_pole_in_bounds =
3290 106 : ContainsNorthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3291 106 : if (CPLGetLastErrorType() != CE_None)
3292 : {
3293 0 : return false;
3294 : }
3295 : south_pole_in_bounds =
3296 106 : ContainsSouthPole(xmin, ymin, xmax, ymax, output_lon_lat_order);
3297 106 : if (CPLGetLastErrorType() != CE_None)
3298 : {
3299 0 : return false;
3300 : }
3301 : }
3302 :
3303 287 : if (degree_input && xmax < xmin)
3304 : {
3305 1 : if (!input_lon_lat_order)
3306 : {
3307 0 : CPLError(CE_Failure, CPLE_AppDefined,
3308 : "latitude max < latitude min.");
3309 0 : return false;
3310 : }
3311 : // handle antimeridian
3312 1 : delta_x = (xmax - xmin + 360.0) / side_pts;
3313 : }
3314 : else
3315 : {
3316 286 : delta_x = (xmax - xmin) / side_pts;
3317 : }
3318 287 : if (degree_input && ymax < ymin)
3319 : {
3320 1 : if (input_lon_lat_order)
3321 : {
3322 0 : CPLError(CE_Failure, CPLE_AppDefined,
3323 : "latitude max < latitude min.");
3324 0 : return false;
3325 : }
3326 : // handle antimeridian
3327 1 : delta_y = (ymax - ymin + 360.0) / side_pts;
3328 : }
3329 : else
3330 : {
3331 286 : delta_y = (ymax - ymin) / side_pts;
3332 : }
3333 :
3334 : // build densified bounding box
3335 : // Note: must be a linear ring for antimeridian logic
3336 2253 : for (int iii = 0; iii < side_pts; iii++)
3337 : {
3338 : // xmin boundary
3339 1966 : y_boundary_array[iii] = ymax - iii * delta_y;
3340 1966 : x_boundary_array[iii] = xmin;
3341 : // ymin boundary
3342 1966 : y_boundary_array[iii + side_pts] = ymin;
3343 1966 : x_boundary_array[iii + side_pts] = xmin + iii * delta_x;
3344 : // xmax boundary
3345 1966 : y_boundary_array[iii + side_pts * 2] = ymin + iii * delta_y;
3346 1966 : x_boundary_array[iii + side_pts * 2] = xmax;
3347 : // ymax boundary
3348 1966 : y_boundary_array[iii + side_pts * 3] = ymax;
3349 1966 : x_boundary_array[iii + side_pts * 3] = xmax - iii * delta_x;
3350 : }
3351 :
3352 : {
3353 287 : CPLErrorHandlerPusher oErrorHandlerPusher(CPLQuietErrorHandler);
3354 287 : bool success = TransformWithErrorCodes(
3355 287 : boundary_len, &x_boundary_array[0], &y_boundary_array[0], nullptr,
3356 287 : nullptr, nullptr);
3357 287 : if (success && CPLGetLastErrorType() != CE_None)
3358 : {
3359 12 : CPLErrorReset();
3360 : }
3361 275 : else if (!success)
3362 : {
3363 0 : return false;
3364 : }
3365 : }
3366 :
3367 287 : if (!degree_output)
3368 : {
3369 181 : *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
3370 181 : *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
3371 181 : *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
3372 181 : *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
3373 :
3374 181 : if (poSRSTarget->IsProjected())
3375 : {
3376 362 : CPLErrorStateBackuper oBackuper(CPLQuietErrorHandler);
3377 :
3378 : auto poBaseTarget = std::unique_ptr<OGRSpatialReference>(
3379 362 : poSRSTarget->CloneGeogCS());
3380 181 : poBaseTarget->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
3381 :
3382 : auto poCTBaseTargetToSrc =
3383 : std::unique_ptr<OGRCoordinateTransformation>(
3384 181 : OGRCreateCoordinateTransformation(poBaseTarget.get(),
3385 362 : poSRSSource));
3386 181 : if (poCTBaseTargetToSrc)
3387 : {
3388 : const double dfLon0 =
3389 181 : poSRSTarget->GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0);
3390 :
3391 181 : double dfSignedPoleLat = 0;
3392 181 : double dfAbsPoleLat = 90;
3393 181 : bool bIncludesPole = false;
3394 : const char *pszProjection =
3395 181 : poSRSTarget->GetAttrValue("PROJECTION");
3396 181 : if (pszProjection &&
3397 181 : EQUAL(pszProjection, SRS_PT_MERCATOR_1SP) && dfLon0 == 0)
3398 : {
3399 : // This MAX_LAT_MERCATOR values is equivalent to the
3400 : // semi_major_axis * PI easting/northing value only
3401 : // for EPSG:3857, but it is also quite
3402 : // reasonable for other Mercator projections
3403 109 : constexpr double MAX_LAT_MERCATOR = 85.0511287798066;
3404 109 : dfAbsPoleLat = MAX_LAT_MERCATOR;
3405 : }
3406 :
3407 : // Detect if a point at long = central_meridian and
3408 : // lat = +/- 90deg is included in the extent.
3409 : // Helps for example for EPSG:4326 to ESRI:53037
3410 543 : for (int iSign = -1; iSign <= 1; iSign += 2)
3411 : {
3412 362 : double dfX = dfLon0;
3413 362 : constexpr double EPS = 1e-8;
3414 362 : double dfY = iSign * (dfAbsPoleLat - EPS);
3415 :
3416 362 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3417 362 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3418 362 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3419 852 : dfY <= ymax &&
3420 128 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr, nullptr,
3421 128 : nullptr))
3422 : {
3423 128 : dfSignedPoleLat = iSign * dfAbsPoleLat;
3424 128 : bIncludesPole = true;
3425 128 : *out_xmin = std::min(*out_xmin, dfX);
3426 128 : *out_ymin = std::min(*out_ymin, dfY);
3427 128 : *out_xmax = std::max(*out_xmax, dfX);
3428 128 : *out_ymax = std::max(*out_ymax, dfY);
3429 : }
3430 : }
3431 :
3432 : const auto TryAtPlusMinus180 =
3433 121 : [this, dfLon0, xmin, ymin, xmax, ymax, out_xmin, out_ymin,
3434 2039 : out_xmax, out_ymax, &poCTBaseTargetToSrc](double dfLat)
3435 : {
3436 363 : for (int iSign = -1; iSign <= 1; iSign += 2)
3437 : {
3438 242 : constexpr double EPS = 1e-8;
3439 : double dfX =
3440 242 : fmod(dfLon0 + iSign * (180 - EPS) + 180, 360) - 180;
3441 242 : double dfY = dfLat;
3442 :
3443 242 : if (poCTBaseTargetToSrc->TransformWithErrorCodes(
3444 242 : 1, &dfX, &dfY, nullptr, nullptr, nullptr) &&
3445 242 : dfX >= xmin && dfY >= ymin && dfX <= xmax &&
3446 618 : dfY <= ymax &&
3447 134 : TransformWithErrorCodes(1, &dfX, &dfY, nullptr,
3448 134 : nullptr, nullptr))
3449 : {
3450 134 : *out_xmin = std::min(*out_xmin, dfX);
3451 134 : *out_ymin = std::min(*out_ymin, dfY);
3452 134 : *out_xmax = std::max(*out_xmax, dfX);
3453 134 : *out_ymax = std::max(*out_ymax, dfY);
3454 : }
3455 : }
3456 121 : };
3457 :
3458 : // For a projected CRS with a central meridian != 0, try to
3459 : // reproject the points with long = +/- 180deg of the central
3460 : // meridian and at lat = latitude_of_origin.
3461 181 : const double dfLat0 = poSRSTarget->GetNormProjParm(
3462 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0);
3463 181 : if (dfLon0 != 0)
3464 : {
3465 57 : TryAtPlusMinus180(dfLat0);
3466 : }
3467 :
3468 181 : if (bIncludesPole && dfLat0 != dfSignedPoleLat)
3469 : {
3470 64 : TryAtPlusMinus180(dfSignedPoleLat);
3471 : }
3472 : }
3473 : }
3474 : }
3475 106 : else if (north_pole_in_bounds && output_lon_lat_order)
3476 : {
3477 3 : *out_xmin = -180;
3478 3 : *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
3479 3 : *out_xmax = 180;
3480 3 : *out_ymax = 90;
3481 : }
3482 103 : else if (north_pole_in_bounds)
3483 : {
3484 1 : *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
3485 1 : *out_ymin = -180;
3486 1 : *out_xmax = 90;
3487 1 : *out_ymax = 180;
3488 : }
3489 102 : else if (south_pole_in_bounds && output_lon_lat_order)
3490 : {
3491 3 : *out_xmin = -180;
3492 3 : *out_ymin = -90;
3493 3 : *out_xmax = 180;
3494 3 : *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
3495 : }
3496 99 : else if (south_pole_in_bounds)
3497 : {
3498 1 : *out_xmin = -90;
3499 1 : *out_ymin = -180;
3500 1 : *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
3501 1 : *out_ymax = 180;
3502 : }
3503 98 : else if (output_lon_lat_order)
3504 : {
3505 96 : *out_xmin = antimeridian_min(&x_boundary_array[0], boundary_len);
3506 96 : *out_xmax = antimeridian_max(&x_boundary_array[0], boundary_len);
3507 96 : *out_ymin = simple_min(&y_boundary_array[0], boundary_len);
3508 96 : *out_ymax = simple_max(&y_boundary_array[0], boundary_len);
3509 : }
3510 : else
3511 : {
3512 2 : *out_xmin = simple_min(&x_boundary_array[0], boundary_len);
3513 2 : *out_xmax = simple_max(&x_boundary_array[0], boundary_len);
3514 2 : *out_ymin = antimeridian_min(&y_boundary_array[0], boundary_len);
3515 2 : *out_ymax = antimeridian_max(&y_boundary_array[0], boundary_len);
3516 : }
3517 :
3518 283 : return *out_xmin != HUGE_VAL && *out_ymin != HUGE_VAL &&
3519 570 : *out_xmax != HUGE_VAL && *out_ymax != HUGE_VAL;
3520 : }
3521 :
3522 : /************************************************************************/
3523 : /* Clone() */
3524 : /************************************************************************/
3525 :
3526 23 : OGRCoordinateTransformation *OGRProjCT::Clone() const
3527 : {
3528 46 : std::unique_ptr<OGRProjCT> poNewCT(new OGRProjCT(*this));
3529 : #if (PROJ_VERSION_MAJOR * 10000 + PROJ_VERSION_MINOR * 100 + \
3530 : PROJ_VERSION_PATCH) < 80001
3531 : // See https://github.com/OSGeo/PROJ/pull/2582
3532 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3533 : // operation being a set of real operations
3534 23 : bool bCloneDone = ((m_pj == nullptr) == (poNewCT->m_pj == nullptr));
3535 23 : if (!bCloneDone)
3536 : {
3537 1 : poNewCT.reset(new OGRProjCT());
3538 : auto ret =
3539 2 : poNewCT->Initialize(poSRSSource, m_osSrcSRS.c_str(), poSRSTarget,
3540 1 : m_osTargetSRS.c_str(), m_options);
3541 1 : if (!ret)
3542 : {
3543 0 : return nullptr;
3544 : }
3545 : }
3546 : #endif // PROJ_VERSION
3547 23 : return poNewCT.release();
3548 : }
3549 :
3550 : /************************************************************************/
3551 : /* GetInverse() */
3552 : /************************************************************************/
3553 :
3554 1150 : OGRCoordinateTransformation *OGRProjCT::GetInverse() const
3555 : {
3556 1150 : PJ *new_pj = nullptr;
3557 : // m_pj can be nullptr if using m_eStrategy != PROJ
3558 1150 : if (m_pj && !bWebMercatorToWGS84LongLat && !bNoTransform)
3559 : {
3560 : // See https://github.com/OSGeo/PROJ/pull/2582
3561 : // This may fail before PROJ 8.0.1 if the m_pj object is a "meta"
3562 : // operation being a set of real operations
3563 1147 : new_pj = proj_clone(OSRGetProjTLSContext(), m_pj);
3564 : }
3565 :
3566 2300 : OGRCoordinateTransformationOptions newOptions(m_options);
3567 1150 : std::swap(newOptions.d->bHasSourceCenterLong,
3568 1150 : newOptions.d->bHasTargetCenterLong);
3569 1150 : std::swap(newOptions.d->dfSourceCenterLong,
3570 1150 : newOptions.d->dfTargetCenterLong);
3571 1150 : newOptions.d->bReverseCO = !newOptions.d->bReverseCO;
3572 1150 : newOptions.d->RefreshCheckWithInvertProj();
3573 :
3574 1150 : if (new_pj == nullptr && !bNoTransform)
3575 : {
3576 9 : return OGRCreateCoordinateTransformation(poSRSTarget, poSRSSource,
3577 9 : newOptions);
3578 : }
3579 :
3580 1141 : auto poNewCT = new OGRProjCT();
3581 :
3582 1141 : if (poSRSTarget)
3583 1135 : poNewCT->poSRSSource = poSRSTarget->Clone();
3584 1141 : poNewCT->m_eSourceFirstAxisOrient = m_eTargetFirstAxisOrient;
3585 1141 : poNewCT->bSourceLatLong = bTargetLatLong;
3586 1141 : poNewCT->bSourceWrap = bTargetWrap;
3587 1141 : poNewCT->dfSourceWrapLong = dfTargetWrapLong;
3588 1141 : poNewCT->bSourceIsDynamicCRS = bTargetIsDynamicCRS;
3589 1141 : poNewCT->dfSourceCoordinateEpoch = dfTargetCoordinateEpoch;
3590 1141 : poNewCT->m_osSrcSRS = m_osTargetSRS;
3591 :
3592 1141 : if (poSRSSource)
3593 1135 : poNewCT->poSRSTarget = poSRSSource->Clone();
3594 1141 : poNewCT->m_eTargetFirstAxisOrient = m_eSourceFirstAxisOrient;
3595 1141 : poNewCT->bTargetLatLong = bSourceLatLong;
3596 1141 : poNewCT->bTargetWrap = bSourceWrap;
3597 1141 : poNewCT->dfTargetWrapLong = dfSourceWrapLong;
3598 1141 : poNewCT->bTargetIsDynamicCRS = bSourceIsDynamicCRS;
3599 1141 : poNewCT->dfTargetCoordinateEpoch = dfSourceCoordinateEpoch;
3600 1141 : poNewCT->m_osTargetSRS = m_osSrcSRS;
3601 :
3602 1141 : poNewCT->ComputeThreshold();
3603 :
3604 1141 : poNewCT->m_pj = new_pj;
3605 1141 : poNewCT->m_bReversePj = !m_bReversePj;
3606 1141 : poNewCT->bNoTransform = bNoTransform;
3607 1141 : poNewCT->m_eStrategy = m_eStrategy;
3608 1141 : poNewCT->m_options = newOptions;
3609 :
3610 1141 : poNewCT->DetectWebMercatorToWGS84();
3611 :
3612 1141 : return poNewCT;
3613 : }
3614 :
3615 : /************************************************************************/
3616 : /* OSRCTCleanCache() */
3617 : /************************************************************************/
3618 :
3619 946 : void OSRCTCleanCache()
3620 : {
3621 946 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3622 946 : delete g_poCTCache;
3623 946 : g_poCTCache = nullptr;
3624 946 : }
3625 :
3626 : /************************************************************************/
3627 : /* MakeCacheKey() */
3628 : /************************************************************************/
3629 :
3630 6609 : CTCacheKey OGRProjCT::MakeCacheKey(
3631 : const OGRSpatialReference *poSRS1, const char *pszSrcSRS,
3632 : const OGRSpatialReference *poSRS2, const char *pszTargetSRS,
3633 : const OGRCoordinateTransformationOptions &options)
3634 : {
3635 : const auto GetKeyForSRS =
3636 13218 : [](const OGRSpatialReference *poSRS, const char *pszText)
3637 : {
3638 13218 : if (poSRS)
3639 : {
3640 26334 : std::string ret(pszText);
3641 13167 : const auto &mapping = poSRS->GetDataAxisToSRSAxisMapping();
3642 39751 : for (const auto &axis : mapping)
3643 : {
3644 26584 : ret += std::to_string(axis);
3645 : }
3646 13167 : return ret;
3647 : }
3648 : else
3649 : {
3650 51 : return std::string("null");
3651 : }
3652 : };
3653 :
3654 6609 : std::string ret(GetKeyForSRS(poSRS1, pszSrcSRS));
3655 6609 : ret += GetKeyForSRS(poSRS2, pszTargetSRS);
3656 6609 : ret += options.d->GetKey();
3657 13218 : return ret;
3658 : }
3659 :
3660 : /************************************************************************/
3661 : /* InsertIntoCache() */
3662 : /************************************************************************/
3663 :
3664 3189 : void OGRProjCT::InsertIntoCache(OGRProjCT *poCT)
3665 : {
3666 : {
3667 6378 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3668 3189 : if (g_poCTCache == nullptr)
3669 : {
3670 191 : g_poCTCache = new lru11::Cache<CTCacheKey, CTCacheValue>();
3671 : }
3672 : }
3673 3189 : const auto key = MakeCacheKey(poCT->poSRSSource, poCT->m_osSrcSRS.c_str(),
3674 3189 : poCT->poSRSTarget,
3675 3189 : poCT->m_osTargetSRS.c_str(), poCT->m_options);
3676 :
3677 3189 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3678 3189 : if (g_poCTCache->contains(key))
3679 : {
3680 761 : delete poCT;
3681 761 : return;
3682 : }
3683 2428 : g_poCTCache->insert(key, std::unique_ptr<OGRProjCT>(poCT));
3684 : }
3685 :
3686 : /************************************************************************/
3687 : /* FindFromCache() */
3688 : /************************************************************************/
3689 :
3690 3760 : OGRProjCT *OGRProjCT::FindFromCache(
3691 : const OGRSpatialReference *poSource, const char *pszSrcSRS,
3692 : const OGRSpatialReference *poTarget, const char *pszTargetSRS,
3693 : const OGRCoordinateTransformationOptions &options)
3694 : {
3695 : {
3696 3760 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3697 3760 : if (g_poCTCache == nullptr || g_poCTCache->empty())
3698 340 : return nullptr;
3699 : }
3700 :
3701 : const auto key =
3702 6840 : MakeCacheKey(poSource, pszSrcSRS, poTarget, pszTargetSRS, options);
3703 : // Get value from cache and remove it
3704 6840 : std::lock_guard<std::mutex> oGuard(g_oCTCacheMutex);
3705 3420 : CTCacheValue *cachedValue = g_poCTCache->getPtr(key);
3706 3420 : if (cachedValue)
3707 : {
3708 1667 : auto poCT = cachedValue->release();
3709 1667 : g_poCTCache->remove(key);
3710 1667 : return poCT;
3711 : }
3712 1753 : return nullptr;
3713 : }
3714 :
3715 : //! @endcond
3716 :
3717 : /************************************************************************/
3718 : /* OCTTransform() */
3719 : /************************************************************************/
3720 :
3721 : /** Transform an array of points
3722 : *
3723 : * @param hTransform Transformation object
3724 : * @param nCount Number of points
3725 : * @param x Array of nCount x values.
3726 : * @param y Array of nCount y values.
3727 : * @param z Array of nCount z values.
3728 : * @return TRUE if a transformation could be found (but not all points may
3729 : * have necessarily succeed to transform), otherwise FALSE.
3730 : */
3731 1053 : int CPL_STDCALL OCTTransform(OGRCoordinateTransformationH hTransform,
3732 : int nCount, double *x, double *y, double *z)
3733 :
3734 : {
3735 1053 : VALIDATE_POINTER1(hTransform, "OCTTransform", FALSE);
3736 :
3737 : return OGRCoordinateTransformation::FromHandle(hTransform)
3738 1053 : ->Transform(nCount, x, y, z);
3739 : }
3740 :
3741 : /************************************************************************/
3742 : /* OCTTransformEx() */
3743 : /************************************************************************/
3744 :
3745 : /** Transform an array of points
3746 : *
3747 : * @param hTransform Transformation object
3748 : * @param nCount Number of points
3749 : * @param x Array of nCount x values.
3750 : * @param y Array of nCount y values.
3751 : * @param z Array of nCount z values.
3752 : * @param pabSuccess Output array of nCount value that will be set to TRUE/FALSE
3753 : * @return TRUE if a transformation could be found (but not all points may
3754 : * have necessarily succeed to transform), otherwise FALSE.
3755 : */
3756 0 : int CPL_STDCALL OCTTransformEx(OGRCoordinateTransformationH hTransform,
3757 : int nCount, double *x, double *y, double *z,
3758 : int *pabSuccess)
3759 :
3760 : {
3761 0 : VALIDATE_POINTER1(hTransform, "OCTTransformEx", FALSE);
3762 :
3763 : return OGRCoordinateTransformation::FromHandle(hTransform)
3764 0 : ->Transform(nCount, x, y, z, pabSuccess);
3765 : }
3766 :
3767 : /************************************************************************/
3768 : /* OCTTransform4D() */
3769 : /************************************************************************/
3770 :
3771 : /** Transform an array of points
3772 : *
3773 : * @param hTransform Transformation object
3774 : * @param nCount Number of points
3775 : * @param x Array of nCount x values. Should not be NULL
3776 : * @param y Array of nCount y values. Should not be NULL
3777 : * @param z Array of nCount z values. Might be NULL
3778 : * @param t Array of nCount time values. Might be NULL
3779 : * @param pabSuccess Output array of nCount value that will be set to
3780 : * TRUE/FALSE. Might be NULL.
3781 : * @since GDAL 3.0
3782 : * @return TRUE if a transformation could be found (but not all points may
3783 : * have necessarily succeed to transform), otherwise FALSE.
3784 : */
3785 15 : int OCTTransform4D(OGRCoordinateTransformationH hTransform, int nCount,
3786 : double *x, double *y, double *z, double *t, int *pabSuccess)
3787 :
3788 : {
3789 15 : VALIDATE_POINTER1(hTransform, "OCTTransform4D", FALSE);
3790 :
3791 15 : return OGRCoordinateTransformation::FromHandle(hTransform)
3792 15 : ->Transform(nCount, x, y, z, t, pabSuccess);
3793 : }
3794 :
3795 : /************************************************************************/
3796 : /* OCTTransform4DWithErrorCodes() */
3797 : /************************************************************************/
3798 :
3799 : /** Transform an array of points
3800 : *
3801 : * @param hTransform Transformation object
3802 : * @param nCount Number of points
3803 : * @param x Array of nCount x values. Should not be NULL
3804 : * @param y Array of nCount y values. Should not be NULL
3805 : * @param z Array of nCount z values. Might be NULL
3806 : * @param t Array of nCount time values. Might be NULL
3807 : * @param panErrorCodes Output array of nCount value that will be set to 0 for
3808 : * success, or a non-zero value for failure. Refer to
3809 : * PROJ 8 public error codes. Might be NULL
3810 : * @since GDAL 3.3, and PROJ 8 to be able to use PROJ public error codes
3811 : * @return TRUE if a transformation could be found (but not all points may
3812 : * have necessarily succeed to transform), otherwise FALSE.
3813 : */
3814 0 : int OCTTransform4DWithErrorCodes(OGRCoordinateTransformationH hTransform,
3815 : int nCount, double *x, double *y, double *z,
3816 : double *t, int *panErrorCodes)
3817 :
3818 : {
3819 0 : VALIDATE_POINTER1(hTransform, "OCTTransform4DWithErrorCodes", FALSE);
3820 :
3821 0 : return OGRCoordinateTransformation::FromHandle(hTransform)
3822 0 : ->TransformWithErrorCodes(nCount, x, y, z, t, panErrorCodes);
3823 : }
3824 :
3825 : /************************************************************************/
3826 : /* OCTTransformBounds() */
3827 : /************************************************************************/
3828 : /** \brief Transform boundary.
3829 : *
3830 : * Transform boundary densifying the edges to account for nonlinear
3831 : * transformations along these edges and extracting the outermost bounds.
3832 : *
3833 : * If the destination CRS is geographic, the first axis is longitude,
3834 : * and xmax < xmin then the bounds crossed the antimeridian.
3835 : * In this scenario there are two polygons, one on each side of the
3836 : * antimeridian. The first polygon should be constructed with (xmin, ymin, 180,
3837 : * ymax) and the second with (-180, ymin, xmax, ymax).
3838 : *
3839 : * If the destination CRS is geographic, the first axis is latitude,
3840 : * and ymax < ymin then the bounds crossed the antimeridian.
3841 : * In this scenario there are two polygons, one on each side of the
3842 : * antimeridian. The first polygon should be constructed with (ymin, xmin, ymax,
3843 : * 180) and the second with (ymin, -180, ymax, xmax).
3844 : *
3845 : * @param hTransform Transformation object
3846 : * @param xmin Minimum bounding coordinate of the first axis in source CRS.
3847 : * @param ymin Minimum bounding coordinate of the second axis in source CRS.
3848 : * @param xmax Maximum bounding coordinate of the first axis in source CRS.
3849 : * @param ymax Maximum bounding coordinate of the second axis in source CRS.
3850 : * @param out_xmin Minimum bounding coordinate of the first axis in target CRS
3851 : * @param out_ymin Minimum bounding coordinate of the second axis in target CRS.
3852 : * @param out_xmax Maximum bounding coordinate of the first axis in target CRS.
3853 : * @param out_ymax Maximum bounding coordinate of the second axis in target CRS.
3854 : * @param densify_pts Recommended to use 21. This is the number of points
3855 : * to use to densify the bounding polygon in the transformation.
3856 : * @return TRUE if successful. FALSE if failures encountered.
3857 : * @since 3.4
3858 : */
3859 25 : int CPL_STDCALL OCTTransformBounds(OGRCoordinateTransformationH hTransform,
3860 : const double xmin, const double ymin,
3861 : const double xmax, const double ymax,
3862 : double *out_xmin, double *out_ymin,
3863 : double *out_xmax, double *out_ymax,
3864 : int densify_pts)
3865 :
3866 : {
3867 25 : VALIDATE_POINTER1(hTransform, "TransformBounds", FALSE);
3868 :
3869 25 : return OGRProjCT::FromHandle(hTransform)
3870 25 : ->TransformBounds(xmin, ymin, xmax, ymax, out_xmin, out_ymin, out_xmax,
3871 25 : out_ymax, densify_pts);
3872 : }
3873 :
3874 : /************************************************************************/
3875 : /* OGRCTDumpStatistics() */
3876 : /************************************************************************/
3877 :
3878 946 : void OGRCTDumpStatistics()
3879 : {
3880 : #ifdef DEBUG_PERF
3881 : CPLDebug("OGR_CT", "Total time in proj_create_crs_to_crs(): %d ms",
3882 : static_cast<int>(g_dfTotalTimeCRStoCRS * 1000));
3883 : CPLDebug("OGR_CT", "Total time in coordinate transformation: %d ms",
3884 : static_cast<int>(g_dfTotalTimeReprojection * 1000));
3885 : #endif
3886 946 : }
|
