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