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