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