1 /******************************************************************************
2 *
3 * Project: PROJ
4 * Purpose: ISO19111:2019 implementation
5 * Author: Even Rouault <even dot rouault at spatialys dot com>
6 *
7 ******************************************************************************
8 * Copyright (c) 2018, Even Rouault <even dot rouault at spatialys dot com>
9 *
10 * Permission is hereby granted, free of charge, to any person obtaining a
11 * copy of this software and associated documentation files (the "Software"),
12 * to deal in the Software without restriction, including without limitation
13 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 * and/or sell copies of the Software, and to permit persons to whom the
15 * Software is furnished to do so, subject to the following conditions:
16 *
17 * The above copyright notice and this permission notice shall be included
18 * in all copies or substantial portions of the Software.
19 *
20 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 * DEALINGS IN THE SOFTWARE.
27 ****************************************************************************/
28
29 #ifndef FROM_PROJ_CPP
30 #define FROM_PROJ_CPP
31 #endif
32
33 #include "proj/coordinateoperation.hpp"
34 #include "proj/common.hpp"
35 #include "proj/crs.hpp"
36 #include "proj/io.hpp"
37 #include "proj/metadata.hpp"
38 #include "proj/util.hpp"
39
40 #include "proj/internal/internal.hpp"
41 #include "proj/internal/io_internal.hpp"
42 #include "proj/internal/tracing.hpp"
43
44 #include "coordinateoperation_internal.hpp"
45 #include "coordinateoperation_private.hpp"
46 #include "oputils.hpp"
47
48 // PROJ include order is sensitive
49 // clang-format off
50 #include "proj.h"
51 #include "proj_internal.h" // M_PI
52 // clang-format on
53
54 #include <algorithm>
55 #include <cassert>
56 #include <cmath>
57 #include <cstring>
58 #include <memory>
59 #include <set>
60 #include <string>
61 #include <vector>
62
63 // #define TRACE_CREATE_OPERATIONS
64 // #define DEBUG_SORT
65 // #define DEBUG_CONCATENATED_OPERATION
66 #if defined(DEBUG_SORT) || defined(DEBUG_CONCATENATED_OPERATION)
67 #include <iostream>
68
69 void dumpWKT(const NS_PROJ::crs::CRS *crs);
dumpWKT(const NS_PROJ::crs::CRS * crs)70 void dumpWKT(const NS_PROJ::crs::CRS *crs) {
71 auto f(NS_PROJ::io::WKTFormatter::create(
72 NS_PROJ::io::WKTFormatter::Convention::WKT2_2019));
73 std::cerr << crs->exportToWKT(f.get()) << std::endl;
74 }
75
76 void dumpWKT(const NS_PROJ::crs::CRSPtr &crs);
dumpWKT(const NS_PROJ::crs::CRSPtr & crs)77 void dumpWKT(const NS_PROJ::crs::CRSPtr &crs) { dumpWKT(crs.get()); }
78
79 void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs);
dumpWKT(const NS_PROJ::crs::CRSNNPtr & crs)80 void dumpWKT(const NS_PROJ::crs::CRSNNPtr &crs) {
81 dumpWKT(crs.as_nullable().get());
82 }
83
84 void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs);
dumpWKT(const NS_PROJ::crs::GeographicCRSPtr & crs)85 void dumpWKT(const NS_PROJ::crs::GeographicCRSPtr &crs) { dumpWKT(crs.get()); }
86
87 void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs);
dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr & crs)88 void dumpWKT(const NS_PROJ::crs::GeographicCRSNNPtr &crs) {
89 dumpWKT(crs.as_nullable().get());
90 }
91
92 #endif
93
94 using namespace NS_PROJ::internal;
95
96 // ---------------------------------------------------------------------------
97
98 NS_PROJ_START
99 namespace operation {
100
101 // ---------------------------------------------------------------------------
102
103 #ifdef TRACE_CREATE_OPERATIONS
104
105 //! @cond Doxygen_Suppress
106
objectAsStr(const common::IdentifiedObject * obj)107 static std::string objectAsStr(const common::IdentifiedObject *obj) {
108 std::string ret(obj->nameStr());
109 const auto &ids = obj->identifiers();
110 if (!ids.empty()) {
111 ret += " (";
112 ret += (*ids[0]->codeSpace()) + ":" + ids[0]->code();
113 ret += ")";
114 }
115 return ret;
116 }
117 //! @endcond
118
119 #endif
120
121 // ---------------------------------------------------------------------------
122
123 //! @cond Doxygen_Suppress
124
getPseudoArea(const metadata::ExtentPtr & extent)125 static double getPseudoArea(const metadata::ExtentPtr &extent) {
126 if (!extent)
127 return 0.0;
128 const auto &geographicElements = extent->geographicElements();
129 if (geographicElements.empty())
130 return 0.0;
131 auto bbox = dynamic_cast<const metadata::GeographicBoundingBox *>(
132 geographicElements[0].get());
133 if (!bbox)
134 return 0;
135 double w = bbox->westBoundLongitude();
136 double s = bbox->southBoundLatitude();
137 double e = bbox->eastBoundLongitude();
138 double n = bbox->northBoundLatitude();
139 if (w > e) {
140 e += 360.0;
141 }
142 // Integrate cos(lat) between south_lat and north_lat
143 return (e - w) * (std::sin(common::Angle(n).getSIValue()) -
144 std::sin(common::Angle(s).getSIValue()));
145 }
146
147 //! @endcond
148
149 // ---------------------------------------------------------------------------
150
151 //! @cond Doxygen_Suppress
152 struct CoordinateOperationContext::Private {
153 io::AuthorityFactoryPtr authorityFactory_{};
154 metadata::ExtentPtr extent_{};
155 double accuracy_ = 0.0;
156 SourceTargetCRSExtentUse sourceAndTargetCRSExtentUse_ =
157 CoordinateOperationContext::SourceTargetCRSExtentUse::SMALLEST;
158 SpatialCriterion spatialCriterion_ =
159 CoordinateOperationContext::SpatialCriterion::STRICT_CONTAINMENT;
160 bool usePROJNames_ = true;
161 GridAvailabilityUse gridAvailabilityUse_ =
162 GridAvailabilityUse::USE_FOR_SORTING;
163 IntermediateCRSUse allowUseIntermediateCRS_ = CoordinateOperationContext::
164 IntermediateCRSUse::IF_NO_DIRECT_TRANSFORMATION;
165 std::vector<std::pair<std::string, std::string>>
166 intermediateCRSAuthCodes_{};
167 bool discardSuperseded_ = true;
168 bool allowBallpark_ = true;
169 };
170 //! @endcond
171
172 // ---------------------------------------------------------------------------
173
174 //! @cond Doxygen_Suppress
175 CoordinateOperationContext::~CoordinateOperationContext() = default;
176 //! @endcond
177
178 // ---------------------------------------------------------------------------
179
CoordinateOperationContext()180 CoordinateOperationContext::CoordinateOperationContext()
181 : d(internal::make_unique<Private>()) {}
182
183 // ---------------------------------------------------------------------------
184
185 /** \brief Return the authority factory, or null */
186 const io::AuthorityFactoryPtr &
getAuthorityFactory() const187 CoordinateOperationContext::getAuthorityFactory() const {
188 return d->authorityFactory_;
189 }
190
191 // ---------------------------------------------------------------------------
192
193 /** \brief Return the desired area of interest, or null */
194 const metadata::ExtentPtr &
getAreaOfInterest() const195 CoordinateOperationContext::getAreaOfInterest() const {
196 return d->extent_;
197 }
198
199 // ---------------------------------------------------------------------------
200
201 /** \brief Set the desired area of interest, or null */
setAreaOfInterest(const metadata::ExtentPtr & extent)202 void CoordinateOperationContext::setAreaOfInterest(
203 const metadata::ExtentPtr &extent) {
204 d->extent_ = extent;
205 }
206
207 // ---------------------------------------------------------------------------
208
209 /** \brief Return the desired accuracy (in metre), or 0 */
getDesiredAccuracy() const210 double CoordinateOperationContext::getDesiredAccuracy() const {
211 return d->accuracy_;
212 }
213
214 // ---------------------------------------------------------------------------
215
216 /** \brief Set the desired accuracy (in metre), or 0 */
setDesiredAccuracy(double accuracy)217 void CoordinateOperationContext::setDesiredAccuracy(double accuracy) {
218 d->accuracy_ = accuracy;
219 }
220
221 // ---------------------------------------------------------------------------
222
223 /** \brief Return whether ballpark transformations are allowed */
getAllowBallparkTransformations() const224 bool CoordinateOperationContext::getAllowBallparkTransformations() const {
225 return d->allowBallpark_;
226 }
227
228 // ---------------------------------------------------------------------------
229
230 /** \brief Set whether ballpark transformations are allowed */
setAllowBallparkTransformations(bool allow)231 void CoordinateOperationContext::setAllowBallparkTransformations(bool allow) {
232 d->allowBallpark_ = allow;
233 }
234
235 // ---------------------------------------------------------------------------
236
237 /** \brief Set how source and target CRS extent should be used
238 * when considering if a transformation can be used (only takes effect if
239 * no area of interest is explicitly defined).
240 *
241 * The default is
242 * CoordinateOperationContext::SourceTargetCRSExtentUse::SMALLEST.
243 */
setSourceAndTargetCRSExtentUse(SourceTargetCRSExtentUse use)244 void CoordinateOperationContext::setSourceAndTargetCRSExtentUse(
245 SourceTargetCRSExtentUse use) {
246 d->sourceAndTargetCRSExtentUse_ = use;
247 }
248
249 // ---------------------------------------------------------------------------
250
251 /** \brief Return how source and target CRS extent should be used
252 * when considering if a transformation can be used (only takes effect if
253 * no area of interest is explicitly defined).
254 *
255 * The default is
256 * CoordinateOperationContext::SourceTargetCRSExtentUse::SMALLEST.
257 */
258 CoordinateOperationContext::SourceTargetCRSExtentUse
getSourceAndTargetCRSExtentUse() const259 CoordinateOperationContext::getSourceAndTargetCRSExtentUse() const {
260 return d->sourceAndTargetCRSExtentUse_;
261 }
262
263 // ---------------------------------------------------------------------------
264
265 /** \brief Set the spatial criterion to use when comparing the area of
266 * validity
267 * of coordinate operations with the area of interest / area of validity of
268 * source and target CRS.
269 *
270 * The default is STRICT_CONTAINMENT.
271 */
setSpatialCriterion(SpatialCriterion criterion)272 void CoordinateOperationContext::setSpatialCriterion(
273 SpatialCriterion criterion) {
274 d->spatialCriterion_ = criterion;
275 }
276
277 // ---------------------------------------------------------------------------
278
279 /** \brief Return the spatial criterion to use when comparing the area of
280 * validity
281 * of coordinate operations with the area of interest / area of validity of
282 * source and target CRS.
283 *
284 * The default is STRICT_CONTAINMENT.
285 */
286 CoordinateOperationContext::SpatialCriterion
getSpatialCriterion() const287 CoordinateOperationContext::getSpatialCriterion() const {
288 return d->spatialCriterion_;
289 }
290
291 // ---------------------------------------------------------------------------
292
293 /** \brief Set whether PROJ alternative grid names should be substituted to
294 * the official authority names.
295 *
296 * This only has effect is an authority factory with a non-null database context
297 * has been attached to this context.
298 *
299 * If set to false, it is still possible to
300 * obtain later the substitution by using io::PROJStringFormatter::create()
301 * with a non-null database context.
302 *
303 * The default is true.
304 */
setUsePROJAlternativeGridNames(bool usePROJNames)305 void CoordinateOperationContext::setUsePROJAlternativeGridNames(
306 bool usePROJNames) {
307 d->usePROJNames_ = usePROJNames;
308 }
309
310 // ---------------------------------------------------------------------------
311
312 /** \brief Return whether PROJ alternative grid names should be substituted to
313 * the official authority names.
314 *
315 * The default is true.
316 */
getUsePROJAlternativeGridNames() const317 bool CoordinateOperationContext::getUsePROJAlternativeGridNames() const {
318 return d->usePROJNames_;
319 }
320
321 // ---------------------------------------------------------------------------
322
323 /** \brief Return whether transformations that are superseded (but not
324 * deprecated)
325 * should be discarded.
326 *
327 * The default is true.
328 */
getDiscardSuperseded() const329 bool CoordinateOperationContext::getDiscardSuperseded() const {
330 return d->discardSuperseded_;
331 }
332
333 // ---------------------------------------------------------------------------
334
335 /** \brief Set whether transformations that are superseded (but not deprecated)
336 * should be discarded.
337 *
338 * The default is true.
339 */
setDiscardSuperseded(bool discard)340 void CoordinateOperationContext::setDiscardSuperseded(bool discard) {
341 d->discardSuperseded_ = discard;
342 }
343
344 // ---------------------------------------------------------------------------
345
346 /** \brief Set how grid availability is used.
347 *
348 * The default is USE_FOR_SORTING.
349 */
setGridAvailabilityUse(GridAvailabilityUse use)350 void CoordinateOperationContext::setGridAvailabilityUse(
351 GridAvailabilityUse use) {
352 d->gridAvailabilityUse_ = use;
353 }
354
355 // ---------------------------------------------------------------------------
356
357 /** \brief Return how grid availability is used.
358 *
359 * The default is USE_FOR_SORTING.
360 */
361 CoordinateOperationContext::GridAvailabilityUse
getGridAvailabilityUse() const362 CoordinateOperationContext::getGridAvailabilityUse() const {
363 return d->gridAvailabilityUse_;
364 }
365
366 // ---------------------------------------------------------------------------
367
368 /** \brief Set whether an intermediate pivot CRS can be used for researching
369 * coordinate operations between a source and target CRS.
370 *
371 * Concretely if in the database there is an operation from A to C
372 * (or C to A), and another one from C to B (or B to C), but no direct
373 * operation between A and B, setting this parameter to
374 * ALWAYS/IF_NO_DIRECT_TRANSFORMATION, allow chaining both operations.
375 *
376 * The current implementation is limited to researching one intermediate
377 * step.
378 *
379 * By default, with the IF_NO_DIRECT_TRANSFORMATION stratgey, all potential
380 * C candidates will be used if there is no direct transformation.
381 */
setAllowUseIntermediateCRS(IntermediateCRSUse use)382 void CoordinateOperationContext::setAllowUseIntermediateCRS(
383 IntermediateCRSUse use) {
384 d->allowUseIntermediateCRS_ = use;
385 }
386
387 // ---------------------------------------------------------------------------
388
389 /** \brief Return whether an intermediate pivot CRS can be used for researching
390 * coordinate operations between a source and target CRS.
391 *
392 * Concretely if in the database there is an operation from A to C
393 * (or C to A), and another one from C to B (or B to C), but no direct
394 * operation between A and B, setting this parameter to
395 * ALWAYS/IF_NO_DIRECT_TRANSFORMATION, allow chaining both operations.
396 *
397 * The default is IF_NO_DIRECT_TRANSFORMATION.
398 */
399 CoordinateOperationContext::IntermediateCRSUse
getAllowUseIntermediateCRS() const400 CoordinateOperationContext::getAllowUseIntermediateCRS() const {
401 return d->allowUseIntermediateCRS_;
402 }
403
404 // ---------------------------------------------------------------------------
405
406 /** \brief Restrict the potential pivot CRSs that can be used when trying to
407 * build a coordinate operation between two CRS that have no direct operation.
408 *
409 * @param intermediateCRSAuthCodes a vector of (auth_name, code) that can be
410 * used as potential pivot RS
411 */
setIntermediateCRS(const std::vector<std::pair<std::string,std::string>> & intermediateCRSAuthCodes)412 void CoordinateOperationContext::setIntermediateCRS(
413 const std::vector<std::pair<std::string, std::string>>
414 &intermediateCRSAuthCodes) {
415 d->intermediateCRSAuthCodes_ = intermediateCRSAuthCodes;
416 }
417
418 // ---------------------------------------------------------------------------
419
420 /** \brief Return the potential pivot CRSs that can be used when trying to
421 * build a coordinate operation between two CRS that have no direct operation.
422 *
423 */
424 const std::vector<std::pair<std::string, std::string>> &
getIntermediateCRS() const425 CoordinateOperationContext::getIntermediateCRS() const {
426 return d->intermediateCRSAuthCodes_;
427 }
428
429 // ---------------------------------------------------------------------------
430
431 /** \brief Creates a context for a coordinate operation.
432 *
433 * If a non null authorityFactory is provided, the resulting context should
434 * not be used simultaneously by more than one thread.
435 *
436 * If authorityFactory->getAuthority() is the empty string, then coordinate
437 * operations from any authority will be searched, with the restrictions set
438 * in the authority_to_authority_preference database table.
439 * If authorityFactory->getAuthority() is set to "any", then coordinate
440 * operations from any authority will be searched
441 * If authorityFactory->getAuthority() is a non-empty string different of "any",
442 * then coordinate operatiosn will be searched only in that authority namespace.
443 *
444 * @param authorityFactory Authority factory, or null if no database lookup
445 * is allowed.
446 * Use io::authorityFactory::create(context, std::string()) to allow all
447 * authorities to be used.
448 * @param extent Area of interest, or null if none is known.
449 * @param accuracy Maximum allowed accuracy in metre, as specified in or
450 * 0 to get best accuracy.
451 * @return a new context.
452 */
create(const io::AuthorityFactoryPtr & authorityFactory,const metadata::ExtentPtr & extent,double accuracy)453 CoordinateOperationContextNNPtr CoordinateOperationContext::create(
454 const io::AuthorityFactoryPtr &authorityFactory,
455 const metadata::ExtentPtr &extent, double accuracy) {
456 auto ctxt = NN_NO_CHECK(
457 CoordinateOperationContext::make_unique<CoordinateOperationContext>());
458 ctxt->d->authorityFactory_ = authorityFactory;
459 ctxt->d->extent_ = extent;
460 ctxt->d->accuracy_ = accuracy;
461 return ctxt;
462 }
463
464 // ---------------------------------------------------------------------------
465
466 //! @cond Doxygen_Suppress
467 struct CoordinateOperationFactory::Private {
468
469 struct Context {
470 // This is the extent of the source CRS and target CRS of the initial
471 // CoordinateOperationFactory::createOperations() public call, not
472 // necessarily the ones of intermediate
473 // CoordinateOperationFactory::Private::createOperations() calls.
474 // This is used to compare transformations area of use against the
475 // area of use of the source & target CRS.
476 const metadata::ExtentPtr &extent1;
477 const metadata::ExtentPtr &extent2;
478 const CoordinateOperationContextNNPtr &context;
479 bool inCreateOperationsWithDatumPivotAntiRecursion = false;
480 bool inCreateOperationsGeogToVertWithAlternativeGeog = false;
481 bool inCreateOperationsGeogToVertWithIntermediateVert = false;
482 bool skipHorizontalTransformation = false;
483 std::map<std::pair<io::AuthorityFactory::ObjectType, std::string>,
484 std::list<std::pair<std::string, std::string>>>
485 cacheNameToCRS{};
486
Contextoperation::CoordinateOperationFactory::Private::Context487 Context(const metadata::ExtentPtr &extent1In,
488 const metadata::ExtentPtr &extent2In,
489 const CoordinateOperationContextNNPtr &contextIn)
490 : extent1(extent1In), extent2(extent2In), context(contextIn) {}
491 };
492
493 static std::vector<CoordinateOperationNNPtr>
494 createOperations(const crs::CRSNNPtr &sourceCRS,
495 const crs::CRSNNPtr &targetCRS, Context &context);
496
497 private:
498 static constexpr bool disallowEmptyIntersection = true;
499
500 static void
501 buildCRSIds(const crs::CRSNNPtr &crs, Private::Context &context,
502 std::list<std::pair<std::string, std::string>> &ids);
503
504 static std::vector<CoordinateOperationNNPtr> findOpsInRegistryDirect(
505 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
506 Private::Context &context, bool &resNonEmptyBeforeFiltering);
507
508 static std::vector<CoordinateOperationNNPtr>
509 findOpsInRegistryDirectTo(const crs::CRSNNPtr &targetCRS,
510 Private::Context &context);
511
512 static std::vector<CoordinateOperationNNPtr>
513 findsOpsInRegistryWithIntermediate(
514 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
515 Private::Context &context,
516 bool useCreateBetweenGeodeticCRSWithDatumBasedIntermediates);
517
518 static void createOperationsFromProj4Ext(
519 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
520 const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst,
521 std::vector<CoordinateOperationNNPtr> &res);
522
523 static bool createOperationsFromDatabase(
524 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
525 Private::Context &context, const crs::GeodeticCRS *geodSrc,
526 const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc,
527 const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
528 const crs::VerticalCRS *vertDst,
529 std::vector<CoordinateOperationNNPtr> &res);
530
531 static std::vector<CoordinateOperationNNPtr>
532 createOperationsGeogToVertFromGeoid(const crs::CRSNNPtr &sourceCRS,
533 const crs::CRSNNPtr &targetCRS,
534 const crs::VerticalCRS *vertDst,
535 Context &context);
536
537 static std::vector<CoordinateOperationNNPtr>
538 createOperationsGeogToVertWithIntermediateVert(
539 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
540 const crs::VerticalCRS *vertDst, Context &context);
541
542 static std::vector<CoordinateOperationNNPtr>
543 createOperationsGeogToVertWithAlternativeGeog(
544 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
545 Context &context);
546
547 static void createOperationsFromDatabaseWithVertCRS(
548 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
549 Private::Context &context, const crs::GeographicCRS *geogSrc,
550 const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
551 const crs::VerticalCRS *vertDst,
552 std::vector<CoordinateOperationNNPtr> &res);
553
554 static void createOperationsGeodToGeod(
555 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
556 Private::Context &context, const crs::GeodeticCRS *geodSrc,
557 const crs::GeodeticCRS *geodDst,
558 std::vector<CoordinateOperationNNPtr> &res);
559
560 static void createOperationsDerivedTo(
561 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
562 Private::Context &context, const crs::DerivedCRS *derivedSrc,
563 std::vector<CoordinateOperationNNPtr> &res);
564
565 static void createOperationsBoundToGeog(
566 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
567 Private::Context &context, const crs::BoundCRS *boundSrc,
568 const crs::GeographicCRS *geogDst,
569 std::vector<CoordinateOperationNNPtr> &res);
570
571 static void createOperationsBoundToVert(
572 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
573 Private::Context &context, const crs::BoundCRS *boundSrc,
574 const crs::VerticalCRS *vertDst,
575 std::vector<CoordinateOperationNNPtr> &res);
576
577 static void createOperationsVertToVert(
578 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
579 Private::Context &context, const crs::VerticalCRS *vertSrc,
580 const crs::VerticalCRS *vertDst,
581 std::vector<CoordinateOperationNNPtr> &res);
582
583 static void createOperationsVertToGeog(
584 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
585 Private::Context &context, const crs::VerticalCRS *vertSrc,
586 const crs::GeographicCRS *geogDst,
587 std::vector<CoordinateOperationNNPtr> &res);
588
589 static void createOperationsVertToGeogBallpark(
590 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
591 Private::Context &context, const crs::VerticalCRS *vertSrc,
592 const crs::GeographicCRS *geogDst,
593 std::vector<CoordinateOperationNNPtr> &res);
594
595 static void createOperationsBoundToBound(
596 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
597 Private::Context &context, const crs::BoundCRS *boundSrc,
598 const crs::BoundCRS *boundDst,
599 std::vector<CoordinateOperationNNPtr> &res);
600
601 static void createOperationsCompoundToGeog(
602 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
603 Private::Context &context, const crs::CompoundCRS *compoundSrc,
604 const crs::GeographicCRS *geogDst,
605 std::vector<CoordinateOperationNNPtr> &res);
606
607 static void createOperationsToGeod(
608 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
609 Private::Context &context, const crs::GeodeticCRS *geodDst,
610 std::vector<CoordinateOperationNNPtr> &res);
611
612 static void createOperationsCompoundToCompound(
613 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
614 Private::Context &context, const crs::CompoundCRS *compoundSrc,
615 const crs::CompoundCRS *compoundDst,
616 std::vector<CoordinateOperationNNPtr> &res);
617
618 static void createOperationsBoundToCompound(
619 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
620 Private::Context &context, const crs::BoundCRS *boundSrc,
621 const crs::CompoundCRS *compoundDst,
622 std::vector<CoordinateOperationNNPtr> &res);
623
624 static std::vector<CoordinateOperationNNPtr> createOperationsGeogToGeog(
625 std::vector<CoordinateOperationNNPtr> &res,
626 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
627 Private::Context &context, const crs::GeographicCRS *geogSrc,
628 const crs::GeographicCRS *geogDst);
629
630 static void createOperationsWithDatumPivot(
631 std::vector<CoordinateOperationNNPtr> &res,
632 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
633 const crs::GeodeticCRS *geodSrc, const crs::GeodeticCRS *geodDst,
634 Context &context);
635
636 static bool
637 hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> &res,
638 const Context &context);
639
640 static void setCRSs(CoordinateOperation *co, const crs::CRSNNPtr &sourceCRS,
641 const crs::CRSNNPtr &targetCRS);
642 };
643 //! @endcond
644
645 // ---------------------------------------------------------------------------
646
647 //! @cond Doxygen_Suppress
648 CoordinateOperationFactory::~CoordinateOperationFactory() = default;
649 //! @endcond
650
651 // ---------------------------------------------------------------------------
652
CoordinateOperationFactory()653 CoordinateOperationFactory::CoordinateOperationFactory() : d(nullptr) {}
654
655 // ---------------------------------------------------------------------------
656
657 /** \brief Find a CoordinateOperation from sourceCRS to targetCRS.
658 *
659 * This is a helper of createOperations(), using a coordinate operation
660 * context
661 * with no authority factory (so no catalog searching is done), no desired
662 * accuracy and no area of interest.
663 * This returns the first operation of the result set of createOperations(),
664 * or null if none found.
665 *
666 * @param sourceCRS source CRS.
667 * @param targetCRS source CRS.
668 * @return a CoordinateOperation or nullptr.
669 */
createOperation(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS) const670 CoordinateOperationPtr CoordinateOperationFactory::createOperation(
671 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS) const {
672 auto res = createOperations(
673 sourceCRS, targetCRS,
674 CoordinateOperationContext::create(nullptr, nullptr, 0.0));
675 if (!res.empty()) {
676 return res[0];
677 }
678 return nullptr;
679 }
680
681 // ---------------------------------------------------------------------------
682
683 //! @cond Doxygen_Suppress
684
685 // ---------------------------------------------------------------------------
686
687 struct PrecomputedOpCharacteristics {
688 double area_{};
689 double accuracy_{};
690 bool isPROJExportable_ = false;
691 bool hasGrids_ = false;
692 bool gridsAvailable_ = false;
693 bool gridsKnown_ = false;
694 size_t stepCount_ = 0;
695 bool isApprox_ = false;
696 bool hasBallparkVertical_ = false;
697 bool isNullTransformation_ = false;
698
699 PrecomputedOpCharacteristics() = default;
PrecomputedOpCharacteristicsoperation::PrecomputedOpCharacteristics700 PrecomputedOpCharacteristics(double area, double accuracy,
701 bool isPROJExportable, bool hasGrids,
702 bool gridsAvailable, bool gridsKnown,
703 size_t stepCount, bool isApprox,
704 bool hasBallparkVertical,
705 bool isNullTransformation)
706 : area_(area), accuracy_(accuracy), isPROJExportable_(isPROJExportable),
707 hasGrids_(hasGrids), gridsAvailable_(gridsAvailable),
708 gridsKnown_(gridsKnown), stepCount_(stepCount), isApprox_(isApprox),
709 hasBallparkVertical_(hasBallparkVertical),
710 isNullTransformation_(isNullTransformation) {}
711 };
712
713 // ---------------------------------------------------------------------------
714
715 // We could have used a lambda instead of this old-school way, but
716 // filterAndSort() is already huge.
717 struct SortFunction {
718
719 const std::map<CoordinateOperation *, PrecomputedOpCharacteristics> ↦
720
SortFunctionoperation::SortFunction721 explicit SortFunction(const std::map<CoordinateOperation *,
722 PrecomputedOpCharacteristics> &mapIn)
723 : map(mapIn) {}
724
725 // Sorting function
726 // Return true if a < b
compareoperation::SortFunction727 bool compare(const CoordinateOperationNNPtr &a,
728 const CoordinateOperationNNPtr &b) const {
729 auto iterA = map.find(a.get());
730 assert(iterA != map.end());
731 auto iterB = map.find(b.get());
732 assert(iterB != map.end());
733
734 // CAUTION: the order of the comparisons is extremely important
735 // to get the intended result.
736
737 if (iterA->second.isPROJExportable_ &&
738 !iterB->second.isPROJExportable_) {
739 return true;
740 }
741 if (!iterA->second.isPROJExportable_ &&
742 iterB->second.isPROJExportable_) {
743 return false;
744 }
745
746 if (!iterA->second.isApprox_ && iterB->second.isApprox_) {
747 return true;
748 }
749 if (iterA->second.isApprox_ && !iterB->second.isApprox_) {
750 return false;
751 }
752
753 if (!iterA->second.hasBallparkVertical_ &&
754 iterB->second.hasBallparkVertical_) {
755 return true;
756 }
757 if (iterA->second.hasBallparkVertical_ &&
758 !iterB->second.hasBallparkVertical_) {
759 return false;
760 }
761
762 if (!iterA->second.isNullTransformation_ &&
763 iterB->second.isNullTransformation_) {
764 return true;
765 }
766 if (iterA->second.isNullTransformation_ &&
767 !iterB->second.isNullTransformation_) {
768 return false;
769 }
770
771 // Operations where grids are all available go before other
772 if (iterA->second.gridsAvailable_ && !iterB->second.gridsAvailable_) {
773 return true;
774 }
775 if (iterB->second.gridsAvailable_ && !iterA->second.gridsAvailable_) {
776 return false;
777 }
778
779 // Operations where grids are all known in our DB go before other
780 if (iterA->second.gridsKnown_ && !iterB->second.gridsKnown_) {
781 return true;
782 }
783 if (iterB->second.gridsKnown_ && !iterA->second.gridsKnown_) {
784 return false;
785 }
786
787 // Operations with known accuracy go before those with unknown accuracy
788 const double accuracyA = iterA->second.accuracy_;
789 const double accuracyB = iterB->second.accuracy_;
790 if (accuracyA >= 0 && accuracyB < 0) {
791 return true;
792 }
793 if (accuracyB >= 0 && accuracyA < 0) {
794 return false;
795 }
796
797 if (accuracyA < 0 && accuracyB < 0) {
798 // unknown accuracy ? then prefer operations with grids, which
799 // are likely to have best practical accuracy
800 if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
801 return true;
802 }
803 if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
804 return false;
805 }
806 }
807
808 // Operations with larger non-zero area of use go before those with
809 // lower one
810 const double areaA = iterA->second.area_;
811 const double areaB = iterB->second.area_;
812 if (areaA > 0) {
813 if (areaA > areaB) {
814 return true;
815 }
816 if (areaA < areaB) {
817 return false;
818 }
819 } else if (areaB > 0) {
820 return false;
821 }
822
823 // Operations with better accuracy go before those with worse one
824 if (accuracyA >= 0 && accuracyA < accuracyB) {
825 return true;
826 }
827 if (accuracyB >= 0 && accuracyB < accuracyA) {
828 return false;
829 }
830
831 if (accuracyA >= 0 && accuracyA == accuracyB) {
832 // same accuracy ? then prefer operations without grids
833 if (!iterA->second.hasGrids_ && iterB->second.hasGrids_) {
834 return true;
835 }
836 if (iterA->second.hasGrids_ && !iterB->second.hasGrids_) {
837 return false;
838 }
839 }
840
841 // The less intermediate steps, the better
842 if (iterA->second.stepCount_ < iterB->second.stepCount_) {
843 return true;
844 }
845 if (iterB->second.stepCount_ < iterA->second.stepCount_) {
846 return false;
847 }
848
849 const auto &a_name = a->nameStr();
850 const auto &b_name = b->nameStr();
851 // The shorter name, the better ?
852 if (a_name.size() < b_name.size()) {
853 return true;
854 }
855 if (b_name.size() < a_name.size()) {
856 return false;
857 }
858
859 // Arbitrary final criterion. We actually return the greater element
860 // first, so that "Amersfoort to WGS 84 (4)" is presented before
861 // "Amersfoort to WGS 84 (3)", which is probably a better guess.
862
863 // Except for French NTF (Paris) to NTF, where the (1) conversion
864 // should be preferred because in the remarks of (2), it is mentioned
865 // OGP prefers value from IGN Paris (code 1467)...
866 if (a_name.find("NTF (Paris) to NTF (1)") != std::string::npos &&
867 b_name.find("NTF (Paris) to NTF (2)") != std::string::npos) {
868 return true;
869 }
870 if (a_name.find("NTF (Paris) to NTF (2)") != std::string::npos &&
871 b_name.find("NTF (Paris) to NTF (1)") != std::string::npos) {
872 return false;
873 }
874 if (a_name.find("NTF (Paris) to RGF93 (1)") != std::string::npos &&
875 b_name.find("NTF (Paris) to RGF93 (2)") != std::string::npos) {
876 return true;
877 }
878 if (a_name.find("NTF (Paris) to RGF93 (2)") != std::string::npos &&
879 b_name.find("NTF (Paris) to RGF93 (1)") != std::string::npos) {
880 return false;
881 }
882
883 return a_name > b_name;
884 }
885
operator ()operation::SortFunction886 bool operator()(const CoordinateOperationNNPtr &a,
887 const CoordinateOperationNNPtr &b) const {
888 const bool ret = compare(a, b);
889 #if 0
890 std::cerr << a->nameStr() << " < " << b->nameStr() << " : " << ret << std::endl;
891 #endif
892 return ret;
893 }
894 };
895
896 // ---------------------------------------------------------------------------
897
getStepCount(const CoordinateOperationNNPtr & op)898 static size_t getStepCount(const CoordinateOperationNNPtr &op) {
899 auto concat = dynamic_cast<const ConcatenatedOperation *>(op.get());
900 size_t stepCount = 1;
901 if (concat) {
902 stepCount = concat->operations().size();
903 }
904 return stepCount;
905 }
906
907 // ---------------------------------------------------------------------------
908
909 // Return number of steps that are transformations (and not conversions)
getTransformationStepCount(const CoordinateOperationNNPtr & op)910 static size_t getTransformationStepCount(const CoordinateOperationNNPtr &op) {
911 auto concat = dynamic_cast<const ConcatenatedOperation *>(op.get());
912 size_t stepCount = 1;
913 if (concat) {
914 stepCount = 0;
915 for (const auto &subOp : concat->operations()) {
916 if (dynamic_cast<const Conversion *>(subOp.get()) == nullptr) {
917 stepCount++;
918 }
919 }
920 }
921 return stepCount;
922 }
923
924 // ---------------------------------------------------------------------------
925
isNullTransformation(const std::string & name)926 static bool isNullTransformation(const std::string &name) {
927 if (name.find(" + ") != std::string::npos)
928 return false;
929 return starts_with(name, BALLPARK_GEOCENTRIC_TRANSLATION) ||
930 starts_with(name, BALLPARK_GEOGRAPHIC_OFFSET) ||
931 starts_with(name, NULL_GEOGRAPHIC_OFFSET) ||
932 starts_with(name, NULL_GEOCENTRIC_TRANSLATION);
933 }
934
935 // ---------------------------------------------------------------------------
936
937 struct FilterResults {
938
FilterResultsoperation::FilterResults939 FilterResults(const std::vector<CoordinateOperationNNPtr> &sourceListIn,
940 const CoordinateOperationContextNNPtr &contextIn,
941 const metadata::ExtentPtr &extent1In,
942 const metadata::ExtentPtr &extent2In,
943 bool forceStrictContainmentTest)
944 : sourceList(sourceListIn), context(contextIn), extent1(extent1In),
945 extent2(extent2In), areaOfInterest(context->getAreaOfInterest()),
946 desiredAccuracy(context->getDesiredAccuracy()),
947 sourceAndTargetCRSExtentUse(
948 context->getSourceAndTargetCRSExtentUse()) {
949
950 computeAreaOfInterest();
951 filterOut(forceStrictContainmentTest);
952 }
953
andSortoperation::FilterResults954 FilterResults &andSort() {
955 sort();
956
957 // And now that we have a sorted list, we can remove uninteresting
958 // results
959 // ...
960 removeSyntheticNullTransforms();
961 removeUninterestingOps();
962 removeDuplicateOps();
963 removeSyntheticNullTransforms();
964 return *this;
965 }
966
967 // ----------------------------------------------------------------------
968
969 // cppcheck-suppress functionStatic
getResoperation::FilterResults970 const std::vector<CoordinateOperationNNPtr> &getRes() { return res; }
971
972 // ----------------------------------------------------------------------
973 private:
974 const std::vector<CoordinateOperationNNPtr> &sourceList;
975 const CoordinateOperationContextNNPtr &context;
976 const metadata::ExtentPtr &extent1;
977 const metadata::ExtentPtr &extent2;
978 metadata::ExtentPtr areaOfInterest;
979 const double desiredAccuracy = context->getDesiredAccuracy();
980 const CoordinateOperationContext::SourceTargetCRSExtentUse
981 sourceAndTargetCRSExtentUse;
982
983 bool hasOpThatContainsAreaOfInterestAndNoGrid = false;
984 std::vector<CoordinateOperationNNPtr> res{};
985
986 // ----------------------------------------------------------------------
computeAreaOfInterestoperation::FilterResults987 void computeAreaOfInterest() {
988
989 // Compute an area of interest from the CRS extent if the user did
990 // not specify one
991 if (!areaOfInterest) {
992 if (sourceAndTargetCRSExtentUse ==
993 CoordinateOperationContext::SourceTargetCRSExtentUse::
994 INTERSECTION) {
995 if (extent1 && extent2) {
996 areaOfInterest =
997 extent1->intersection(NN_NO_CHECK(extent2));
998 }
999 } else if (sourceAndTargetCRSExtentUse ==
1000 CoordinateOperationContext::SourceTargetCRSExtentUse::
1001 SMALLEST) {
1002 if (extent1 && extent2) {
1003 if (getPseudoArea(extent1) < getPseudoArea(extent2)) {
1004 areaOfInterest = extent1;
1005 } else {
1006 areaOfInterest = extent2;
1007 }
1008 } else if (extent1) {
1009 areaOfInterest = extent1;
1010 } else {
1011 areaOfInterest = extent2;
1012 }
1013 }
1014 }
1015 }
1016
1017 // ---------------------------------------------------------------------------
1018
filterOutoperation::FilterResults1019 void filterOut(bool forceStrictContainmentTest) {
1020
1021 // Filter out operations that do not match the expected accuracy
1022 // and area of use.
1023 const auto spatialCriterion =
1024 forceStrictContainmentTest
1025 ? CoordinateOperationContext::SpatialCriterion::
1026 STRICT_CONTAINMENT
1027 : context->getSpatialCriterion();
1028 bool hasOnlyBallpark = true;
1029 bool hasNonBallparkWithoutExtent = false;
1030 bool hasNonBallparkOpWithExtent = false;
1031 const bool allowBallpark = context->getAllowBallparkTransformations();
1032 for (const auto &op : sourceList) {
1033 if (desiredAccuracy != 0) {
1034 const double accuracy = getAccuracy(op);
1035 if (accuracy < 0 || accuracy > desiredAccuracy) {
1036 continue;
1037 }
1038 }
1039 if (!allowBallpark && op->hasBallparkTransformation()) {
1040 continue;
1041 }
1042 if (areaOfInterest) {
1043 bool emptyIntersection = false;
1044 auto extent = getExtent(op, true, emptyIntersection);
1045 if (!extent) {
1046 if (!op->hasBallparkTransformation()) {
1047 hasNonBallparkWithoutExtent = true;
1048 }
1049 continue;
1050 }
1051 if (!op->hasBallparkTransformation()) {
1052 hasNonBallparkOpWithExtent = true;
1053 }
1054 bool extentContains =
1055 extent->contains(NN_NO_CHECK(areaOfInterest));
1056 if (!hasOpThatContainsAreaOfInterestAndNoGrid &&
1057 extentContains) {
1058 if (!op->hasBallparkTransformation() &&
1059 op->gridsNeeded(nullptr, true).empty()) {
1060 hasOpThatContainsAreaOfInterestAndNoGrid = true;
1061 }
1062 }
1063 if (spatialCriterion ==
1064 CoordinateOperationContext::SpatialCriterion::
1065 STRICT_CONTAINMENT &&
1066 !extentContains) {
1067 continue;
1068 }
1069 if (spatialCriterion ==
1070 CoordinateOperationContext::SpatialCriterion::
1071 PARTIAL_INTERSECTION &&
1072 !extent->intersects(NN_NO_CHECK(areaOfInterest))) {
1073 continue;
1074 }
1075 } else if (sourceAndTargetCRSExtentUse ==
1076 CoordinateOperationContext::SourceTargetCRSExtentUse::
1077 BOTH) {
1078 bool emptyIntersection = false;
1079 auto extent = getExtent(op, true, emptyIntersection);
1080 if (!extent) {
1081 if (!op->hasBallparkTransformation()) {
1082 hasNonBallparkWithoutExtent = true;
1083 }
1084 continue;
1085 }
1086 if (!op->hasBallparkTransformation()) {
1087 hasNonBallparkOpWithExtent = true;
1088 }
1089 bool extentContainsExtent1 =
1090 !extent1 || extent->contains(NN_NO_CHECK(extent1));
1091 bool extentContainsExtent2 =
1092 !extent2 || extent->contains(NN_NO_CHECK(extent2));
1093 if (!hasOpThatContainsAreaOfInterestAndNoGrid &&
1094 extentContainsExtent1 && extentContainsExtent2) {
1095 if (!op->hasBallparkTransformation() &&
1096 op->gridsNeeded(nullptr, true).empty()) {
1097 hasOpThatContainsAreaOfInterestAndNoGrid = true;
1098 }
1099 }
1100 if (spatialCriterion ==
1101 CoordinateOperationContext::SpatialCriterion::
1102 STRICT_CONTAINMENT) {
1103 if (!extentContainsExtent1 || !extentContainsExtent2) {
1104 continue;
1105 }
1106 } else if (spatialCriterion ==
1107 CoordinateOperationContext::SpatialCriterion::
1108 PARTIAL_INTERSECTION) {
1109 bool extentIntersectsExtent1 =
1110 !extent1 || extent->intersects(NN_NO_CHECK(extent1));
1111 bool extentIntersectsExtent2 =
1112 extent2 && extent->intersects(NN_NO_CHECK(extent2));
1113 if (!extentIntersectsExtent1 || !extentIntersectsExtent2) {
1114 continue;
1115 }
1116 }
1117 }
1118 if (!op->hasBallparkTransformation()) {
1119 hasOnlyBallpark = false;
1120 }
1121 res.emplace_back(op);
1122 }
1123
1124 // In case no operation has an extent and no result is found,
1125 // retain all initial operations that match accuracy criterion.
1126 if ((res.empty() && !hasNonBallparkOpWithExtent) ||
1127 (hasOnlyBallpark && hasNonBallparkWithoutExtent)) {
1128 for (const auto &op : sourceList) {
1129 if (desiredAccuracy != 0) {
1130 const double accuracy = getAccuracy(op);
1131 if (accuracy < 0 || accuracy > desiredAccuracy) {
1132 continue;
1133 }
1134 }
1135 if (!allowBallpark && op->hasBallparkTransformation()) {
1136 continue;
1137 }
1138 res.emplace_back(op);
1139 }
1140 }
1141 }
1142
1143 // ----------------------------------------------------------------------
1144
sortoperation::FilterResults1145 void sort() {
1146
1147 // Precompute a number of parameters for each operation that will be
1148 // useful for the sorting.
1149 std::map<CoordinateOperation *, PrecomputedOpCharacteristics> map;
1150 const auto gridAvailabilityUse = context->getGridAvailabilityUse();
1151 for (const auto &op : res) {
1152 bool dummy = false;
1153 auto extentOp = getExtent(op, true, dummy);
1154 double area = 0.0;
1155 if (extentOp) {
1156 if (areaOfInterest) {
1157 area = getPseudoArea(
1158 extentOp->intersection(NN_NO_CHECK(areaOfInterest)));
1159 } else if (extent1 && extent2) {
1160 auto x = extentOp->intersection(NN_NO_CHECK(extent1));
1161 auto y = extentOp->intersection(NN_NO_CHECK(extent2));
1162 area = getPseudoArea(x) + getPseudoArea(y) -
1163 ((x && y)
1164 ? getPseudoArea(x->intersection(NN_NO_CHECK(y)))
1165 : 0.0);
1166 } else if (extent1) {
1167 area = getPseudoArea(
1168 extentOp->intersection(NN_NO_CHECK(extent1)));
1169 } else if (extent2) {
1170 area = getPseudoArea(
1171 extentOp->intersection(NN_NO_CHECK(extent2)));
1172 } else {
1173 area = getPseudoArea(extentOp);
1174 }
1175 }
1176
1177 bool hasGrids = false;
1178 bool gridsAvailable = true;
1179 bool gridsKnown = true;
1180 if (context->getAuthorityFactory()) {
1181 const auto gridsNeeded = op->gridsNeeded(
1182 context->getAuthorityFactory()->databaseContext(),
1183 gridAvailabilityUse ==
1184 CoordinateOperationContext::GridAvailabilityUse::
1185 KNOWN_AVAILABLE);
1186 for (const auto &gridDesc : gridsNeeded) {
1187 hasGrids = true;
1188 if (gridAvailabilityUse ==
1189 CoordinateOperationContext::GridAvailabilityUse::
1190 USE_FOR_SORTING &&
1191 !gridDesc.available) {
1192 gridsAvailable = false;
1193 }
1194 if (gridDesc.packageName.empty() &&
1195 !(!gridDesc.url.empty() && gridDesc.openLicense) &&
1196 !gridDesc.available) {
1197 gridsKnown = false;
1198 }
1199 }
1200 }
1201
1202 const auto stepCount = getStepCount(op);
1203
1204 bool isPROJExportable = false;
1205 auto formatter = io::PROJStringFormatter::create();
1206 try {
1207 op->exportToPROJString(formatter.get());
1208 // Grids might be missing, but at least this is something
1209 // PROJ could potentially process
1210 isPROJExportable = true;
1211 } catch (const std::exception &) {
1212 }
1213
1214 #if 0
1215 std::cerr << op->nameStr() << " ";
1216 std::cerr << area << " ";
1217 std::cerr << getAccuracy(op) << " ";
1218 std::cerr << isPROJExportable << " ";
1219 std::cerr << hasGrids << " ";
1220 std::cerr << gridsAvailable << " ";
1221 std::cerr << gridsKnown << " ";
1222 std::cerr << stepCount << " ";
1223 std::cerr << op->hasBallparkTransformation() << " ";
1224 std::cerr << isNullTransformation(op->nameStr()) << " ";
1225 std::cerr << std::endl;
1226 #endif
1227 map[op.get()] = PrecomputedOpCharacteristics(
1228 area, getAccuracy(op), isPROJExportable, hasGrids,
1229 gridsAvailable, gridsKnown, stepCount,
1230 op->hasBallparkTransformation(),
1231 op->nameStr().find("ballpark vertical transformation") !=
1232 std::string::npos,
1233 isNullTransformation(op->nameStr()));
1234 }
1235
1236 // Sort !
1237 SortFunction sortFunc(map);
1238 std::sort(res.begin(), res.end(), sortFunc);
1239
1240 // Debug code to check consistency of the sort function
1241 #ifdef DEBUG_SORT
1242 constexpr bool debugSort = true;
1243 #elif !defined(NDEBUG)
1244 const bool debugSort = getenv("PROJ_DEBUG_SORT_FUNCT") != nullptr;
1245 #endif
1246 #if defined(DEBUG_SORT) || !defined(NDEBUG)
1247 if (debugSort) {
1248 const bool assertIfIssue =
1249 !(getenv("PROJ_DEBUG_SORT_FUNCT_ASSERT") != nullptr);
1250 for (size_t i = 0; i < res.size(); ++i) {
1251 for (size_t j = i + 1; j < res.size(); ++j) {
1252 if (sortFunc(res[j], res[i])) {
1253 #ifdef DEBUG_SORT
1254 std::cerr << "Sorting issue with entry " << i << "("
1255 << res[i]->nameStr() << ") and " << j << "("
1256 << res[j]->nameStr() << ")" << std::endl;
1257 #endif
1258 if (assertIfIssue) {
1259 assert(false);
1260 }
1261 }
1262 }
1263 }
1264 }
1265 #endif
1266 }
1267
1268 // ----------------------------------------------------------------------
1269
removeSyntheticNullTransformsoperation::FilterResults1270 void removeSyntheticNullTransforms() {
1271
1272 // If we have more than one result, and than the last result is the
1273 // default "Ballpark geographic offset" or "Ballpark geocentric
1274 // translation" operations we have synthetized, and that at least one
1275 // operation has the desired area of interest and does not require the
1276 // use of grids, remove it as all previous results are necessarily
1277 // better
1278 if (hasOpThatContainsAreaOfInterestAndNoGrid && res.size() > 1) {
1279 const auto &opLast = res.back();
1280 if (opLast->hasBallparkTransformation() ||
1281 isNullTransformation(opLast->nameStr())) {
1282 std::vector<CoordinateOperationNNPtr> resTemp;
1283 for (size_t i = 0; i < res.size() - 1; i++) {
1284 resTemp.emplace_back(res[i]);
1285 }
1286 res = std::move(resTemp);
1287 }
1288 }
1289 }
1290
1291 // ----------------------------------------------------------------------
1292
removeUninterestingOpsoperation::FilterResults1293 void removeUninterestingOps() {
1294
1295 // Eliminate operations that bring nothing, ie for a given area of use,
1296 // do not keep operations that have similar or worse accuracy, but
1297 // involve more (non conversion) steps
1298 std::vector<CoordinateOperationNNPtr> resTemp;
1299 metadata::ExtentPtr lastExtent;
1300 double lastAccuracy = -1;
1301 size_t lastStepCount = 0;
1302 CoordinateOperationPtr lastOp;
1303
1304 bool first = true;
1305 for (const auto &op : res) {
1306 const auto curAccuracy = getAccuracy(op);
1307 bool dummy = false;
1308 const auto curExtent = getExtent(op, true, dummy);
1309 const auto curStepCount = getTransformationStepCount(op);
1310
1311 if (first) {
1312 resTemp.emplace_back(op);
1313 first = false;
1314 } else {
1315 if (lastOp->_isEquivalentTo(op.get())) {
1316 continue;
1317 }
1318 const bool sameExtent =
1319 ((!curExtent && !lastExtent) ||
1320 (curExtent && lastExtent &&
1321 curExtent->contains(NN_NO_CHECK(lastExtent)) &&
1322 lastExtent->contains(NN_NO_CHECK(curExtent))));
1323 if (((curAccuracy >= lastAccuracy && lastAccuracy >= 0) ||
1324 (curAccuracy < 0 && lastAccuracy >= 0)) &&
1325 sameExtent && curStepCount > lastStepCount) {
1326 continue;
1327 }
1328
1329 resTemp.emplace_back(op);
1330 }
1331
1332 lastOp = op.as_nullable();
1333 lastStepCount = curStepCount;
1334 lastExtent = curExtent;
1335 lastAccuracy = curAccuracy;
1336 }
1337 res = std::move(resTemp);
1338 }
1339
1340 // ----------------------------------------------------------------------
1341
1342 // cppcheck-suppress functionStatic
removeDuplicateOpsoperation::FilterResults1343 void removeDuplicateOps() {
1344
1345 if (res.size() <= 1) {
1346 return;
1347 }
1348
1349 // When going from EPSG:4807 (NTF Paris) to EPSG:4171 (RGC93), we get
1350 // EPSG:7811, NTF (Paris) to RGF93 (2), 1 m
1351 // and unknown id, NTF (Paris) to NTF (1) + Inverse of RGF93 to NTF (2),
1352 // 1 m
1353 // both have same PROJ string and extent
1354 // Do not keep the later (that has more steps) as it adds no value.
1355
1356 std::set<std::string> setPROJPlusExtent;
1357 std::vector<CoordinateOperationNNPtr> resTemp;
1358 for (const auto &op : res) {
1359 auto formatter = io::PROJStringFormatter::create();
1360 try {
1361 std::string key(op->exportToPROJString(formatter.get()));
1362 bool dummy = false;
1363 auto extentOp = getExtent(op, true, dummy);
1364 if (extentOp) {
1365 const auto &geogElts = extentOp->geographicElements();
1366 if (geogElts.size() == 1) {
1367 auto bbox = dynamic_cast<
1368 const metadata::GeographicBoundingBox *>(
1369 geogElts[0].get());
1370 if (bbox) {
1371 double w = bbox->westBoundLongitude();
1372 double s = bbox->southBoundLatitude();
1373 double e = bbox->eastBoundLongitude();
1374 double n = bbox->northBoundLatitude();
1375 key += "-";
1376 key += toString(w);
1377 key += "-";
1378 key += toString(s);
1379 key += "-";
1380 key += toString(e);
1381 key += "-";
1382 key += toString(n);
1383 }
1384 }
1385 }
1386
1387 if (setPROJPlusExtent.find(key) == setPROJPlusExtent.end()) {
1388 resTemp.emplace_back(op);
1389 setPROJPlusExtent.insert(key);
1390 }
1391 } catch (const std::exception &) {
1392 resTemp.emplace_back(op);
1393 }
1394 }
1395 res = std::move(resTemp);
1396 }
1397 };
1398
1399 // ---------------------------------------------------------------------------
1400
1401 /** \brief Filter operations and sort them given context.
1402 *
1403 * If a desired accuracy is specified, only keep operations whose accuracy
1404 * is at least the desired one.
1405 * If an area of interest is specified, only keep operations whose area of
1406 * use include the area of interest.
1407 * Then sort remaining operations by descending area of use, and increasing
1408 * accuracy.
1409 */
1410 static std::vector<CoordinateOperationNNPtr>
filterAndSort(const std::vector<CoordinateOperationNNPtr> & sourceList,const CoordinateOperationContextNNPtr & context,const metadata::ExtentPtr & extent1,const metadata::ExtentPtr & extent2)1411 filterAndSort(const std::vector<CoordinateOperationNNPtr> &sourceList,
1412 const CoordinateOperationContextNNPtr &context,
1413 const metadata::ExtentPtr &extent1,
1414 const metadata::ExtentPtr &extent2) {
1415 #ifdef TRACE_CREATE_OPERATIONS
1416 ENTER_FUNCTION();
1417 logTrace("number of results before filter and sort: " +
1418 toString(static_cast<int>(sourceList.size())));
1419 #endif
1420 auto resFiltered =
1421 FilterResults(sourceList, context, extent1, extent2, false)
1422 .andSort()
1423 .getRes();
1424 #ifdef TRACE_CREATE_OPERATIONS
1425 logTrace("number of results after filter and sort: " +
1426 toString(static_cast<int>(resFiltered.size())));
1427 #endif
1428 return resFiltered;
1429 }
1430 //! @endcond
1431
1432 // ---------------------------------------------------------------------------
1433
1434 //! @cond Doxygen_Suppress
1435 // Apply the inverse() method on all elements of the input list
1436 static std::vector<CoordinateOperationNNPtr>
applyInverse(const std::vector<CoordinateOperationNNPtr> & list)1437 applyInverse(const std::vector<CoordinateOperationNNPtr> &list) {
1438 auto res = list;
1439 for (auto &op : res) {
1440 #ifdef DEBUG
1441 auto opNew = op->inverse();
1442 assert(opNew->targetCRS()->isEquivalentTo(op->sourceCRS().get()));
1443 assert(opNew->sourceCRS()->isEquivalentTo(op->targetCRS().get()));
1444 op = opNew;
1445 #else
1446 op = op->inverse();
1447 #endif
1448 }
1449 return res;
1450 }
1451 //! @endcond
1452
1453 // ---------------------------------------------------------------------------
1454
1455 //! @cond Doxygen_Suppress
1456
buildCRSIds(const crs::CRSNNPtr & crs,Private::Context & context,std::list<std::pair<std::string,std::string>> & ids)1457 void CoordinateOperationFactory::Private::buildCRSIds(
1458 const crs::CRSNNPtr &crs, Private::Context &context,
1459 std::list<std::pair<std::string, std::string>> &ids) {
1460 const auto &authFactory = context.context->getAuthorityFactory();
1461 assert(authFactory);
1462 for (const auto &id : crs->identifiers()) {
1463 const auto &authName = *(id->codeSpace());
1464 const auto &code = id->code();
1465 if (!authName.empty()) {
1466 const auto tmpAuthFactory = io::AuthorityFactory::create(
1467 authFactory->databaseContext(), authName);
1468 try {
1469 // Consistency check for the ID attached to the object.
1470 // See https://github.com/OSGeo/PROJ/issues/1982 where EPSG:4656
1471 // is attached to a GeographicCRS whereas it is a ProjectedCRS
1472 if (tmpAuthFactory->createCoordinateReferenceSystem(code)
1473 ->_isEquivalentTo(
1474 crs.get(),
1475 util::IComparable::Criterion::
1476 EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS)) {
1477 ids.emplace_back(authName, code);
1478 } else {
1479 // TODO? log this inconsistency
1480 }
1481 } catch (const std::exception &) {
1482 // TODO? log this inconsistency
1483 }
1484 }
1485 }
1486 if (ids.empty()) {
1487 std::vector<io::AuthorityFactory::ObjectType> allowedObjects;
1488 auto geogCRS = dynamic_cast<const crs::GeographicCRS *>(crs.get());
1489 if (geogCRS) {
1490 allowedObjects.push_back(
1491 geogCRS->coordinateSystem()->axisList().size() == 2
1492 ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
1493 : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
1494 } else if (dynamic_cast<crs::ProjectedCRS *>(crs.get())) {
1495 allowedObjects.push_back(
1496 io::AuthorityFactory::ObjectType::PROJECTED_CRS);
1497 } else if (dynamic_cast<crs::VerticalCRS *>(crs.get())) {
1498 allowedObjects.push_back(
1499 io::AuthorityFactory::ObjectType::VERTICAL_CRS);
1500 }
1501 if (!allowedObjects.empty()) {
1502
1503 const std::pair<io::AuthorityFactory::ObjectType, std::string> key(
1504 allowedObjects[0], crs->nameStr());
1505 auto iter = context.cacheNameToCRS.find(key);
1506 if (iter != context.cacheNameToCRS.end()) {
1507 ids = iter->second;
1508 return;
1509 }
1510
1511 const auto &authFactoryName = authFactory->getAuthority();
1512 try {
1513 const auto tmpAuthFactory = io::AuthorityFactory::create(
1514 authFactory->databaseContext(),
1515 (authFactoryName.empty() || authFactoryName == "any")
1516 ? std::string()
1517 : authFactoryName);
1518
1519 auto matches = tmpAuthFactory->createObjectsFromName(
1520 crs->nameStr(), allowedObjects, false, 2);
1521 if (matches.size() == 1 &&
1522 crs->_isEquivalentTo(
1523 matches.front().get(),
1524 util::IComparable::Criterion::EQUIVALENT) &&
1525 !matches.front()->identifiers().empty()) {
1526 const auto &tmpIds = matches.front()->identifiers();
1527 ids.emplace_back(*(tmpIds[0]->codeSpace()),
1528 tmpIds[0]->code());
1529 }
1530 } catch (const std::exception &) {
1531 }
1532 context.cacheNameToCRS[key] = ids;
1533 }
1534 }
1535 }
1536
1537 // ---------------------------------------------------------------------------
1538
1539 static std::vector<std::string>
getCandidateAuthorities(const io::AuthorityFactoryPtr & authFactory,const std::string & srcAuthName,const std::string & targetAuthName)1540 getCandidateAuthorities(const io::AuthorityFactoryPtr &authFactory,
1541 const std::string &srcAuthName,
1542 const std::string &targetAuthName) {
1543 const auto &authFactoryName = authFactory->getAuthority();
1544 std::vector<std::string> authorities;
1545 if (authFactoryName == "any") {
1546 authorities.emplace_back();
1547 }
1548 if (authFactoryName.empty()) {
1549 authorities = authFactory->databaseContext()->getAllowedAuthorities(
1550 srcAuthName, targetAuthName);
1551 if (authorities.empty()) {
1552 authorities.emplace_back();
1553 }
1554 } else {
1555 authorities.emplace_back(authFactoryName);
1556 }
1557 return authorities;
1558 }
1559
1560 // ---------------------------------------------------------------------------
1561
1562 // Look in the authority registry for operations from sourceCRS to targetCRS
1563 std::vector<CoordinateOperationNNPtr>
findOpsInRegistryDirect(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,bool & resNonEmptyBeforeFiltering)1564 CoordinateOperationFactory::Private::findOpsInRegistryDirect(
1565 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
1566 Private::Context &context, bool &resNonEmptyBeforeFiltering) {
1567 const auto &authFactory = context.context->getAuthorityFactory();
1568 assert(authFactory);
1569
1570 #ifdef TRACE_CREATE_OPERATIONS
1571 ENTER_BLOCK("findOpsInRegistryDirect(" + objectAsStr(sourceCRS.get()) +
1572 " --> " + objectAsStr(targetCRS.get()) + ")");
1573 #endif
1574
1575 resNonEmptyBeforeFiltering = false;
1576 std::list<std::pair<std::string, std::string>> sourceIds;
1577 std::list<std::pair<std::string, std::string>> targetIds;
1578 buildCRSIds(sourceCRS, context, sourceIds);
1579 buildCRSIds(targetCRS, context, targetIds);
1580
1581 const auto gridAvailabilityUse = context.context->getGridAvailabilityUse();
1582 for (const auto &idSrc : sourceIds) {
1583 const auto &srcAuthName = idSrc.first;
1584 const auto &srcCode = idSrc.second;
1585 for (const auto &idTarget : targetIds) {
1586 const auto &targetAuthName = idTarget.first;
1587 const auto &targetCode = idTarget.second;
1588
1589 const auto authorities(getCandidateAuthorities(
1590 authFactory, srcAuthName, targetAuthName));
1591 std::vector<CoordinateOperationNNPtr> res;
1592 for (const auto &authority : authorities) {
1593 const auto authName =
1594 authority == "any" ? std::string() : authority;
1595 const auto tmpAuthFactory = io::AuthorityFactory::create(
1596 authFactory->databaseContext(), authName);
1597 auto resTmp =
1598 tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
1599 srcAuthName, srcCode, targetAuthName, targetCode,
1600 context.context->getUsePROJAlternativeGridNames(),
1601 gridAvailabilityUse ==
1602 CoordinateOperationContext::
1603 GridAvailabilityUse::
1604 DISCARD_OPERATION_IF_MISSING_GRID ||
1605 gridAvailabilityUse ==
1606 CoordinateOperationContext::
1607 GridAvailabilityUse::KNOWN_AVAILABLE,
1608 gridAvailabilityUse ==
1609 CoordinateOperationContext::GridAvailabilityUse::
1610 KNOWN_AVAILABLE,
1611 context.context->getDiscardSuperseded(), true, false,
1612 context.extent1, context.extent2);
1613 res.insert(res.end(), resTmp.begin(), resTmp.end());
1614 if (authName == "PROJ") {
1615 continue;
1616 }
1617 if (!res.empty()) {
1618 resNonEmptyBeforeFiltering = true;
1619 auto resFiltered =
1620 FilterResults(res, context.context, context.extent1,
1621 context.extent2, false)
1622 .getRes();
1623 #ifdef TRACE_CREATE_OPERATIONS
1624 logTrace("filtering reduced from " +
1625 toString(static_cast<int>(res.size())) + " to " +
1626 toString(static_cast<int>(resFiltered.size())));
1627 #endif
1628 return resFiltered;
1629 }
1630 }
1631 }
1632 }
1633 return std::vector<CoordinateOperationNNPtr>();
1634 }
1635
1636 // ---------------------------------------------------------------------------
1637
1638 // Look in the authority registry for operations to targetCRS
1639 std::vector<CoordinateOperationNNPtr>
findOpsInRegistryDirectTo(const crs::CRSNNPtr & targetCRS,Private::Context & context)1640 CoordinateOperationFactory::Private::findOpsInRegistryDirectTo(
1641 const crs::CRSNNPtr &targetCRS, Private::Context &context) {
1642 #ifdef TRACE_CREATE_OPERATIONS
1643 ENTER_BLOCK("findOpsInRegistryDirectTo({any} -->" +
1644 objectAsStr(targetCRS.get()) + ")");
1645 #endif
1646
1647 const auto &authFactory = context.context->getAuthorityFactory();
1648 assert(authFactory);
1649
1650 std::list<std::pair<std::string, std::string>> ids;
1651 buildCRSIds(targetCRS, context, ids);
1652
1653 const auto gridAvailabilityUse = context.context->getGridAvailabilityUse();
1654 for (const auto &id : ids) {
1655 const auto &targetAuthName = id.first;
1656 const auto &targetCode = id.second;
1657
1658 const auto authorities(getCandidateAuthorities(
1659 authFactory, targetAuthName, targetAuthName));
1660 for (const auto &authority : authorities) {
1661 const auto tmpAuthFactory = io::AuthorityFactory::create(
1662 authFactory->databaseContext(),
1663 authority == "any" ? std::string() : authority);
1664 auto res = tmpAuthFactory->createFromCoordinateReferenceSystemCodes(
1665 std::string(), std::string(), targetAuthName, targetCode,
1666 context.context->getUsePROJAlternativeGridNames(),
1667
1668 gridAvailabilityUse ==
1669 CoordinateOperationContext::GridAvailabilityUse::
1670 DISCARD_OPERATION_IF_MISSING_GRID ||
1671 gridAvailabilityUse ==
1672 CoordinateOperationContext::GridAvailabilityUse::
1673 KNOWN_AVAILABLE,
1674 gridAvailabilityUse == CoordinateOperationContext::
1675 GridAvailabilityUse::KNOWN_AVAILABLE,
1676 context.context->getDiscardSuperseded(), true, true,
1677 context.extent1, context.extent2);
1678 if (!res.empty()) {
1679 auto resFiltered =
1680 FilterResults(res, context.context, context.extent1,
1681 context.extent2, false)
1682 .getRes();
1683 #ifdef TRACE_CREATE_OPERATIONS
1684 logTrace("filtering reduced from " +
1685 toString(static_cast<int>(res.size())) + " to " +
1686 toString(static_cast<int>(resFiltered.size())));
1687 #endif
1688 return resFiltered;
1689 }
1690 }
1691 }
1692 return std::vector<CoordinateOperationNNPtr>();
1693 }
1694
1695 //! @endcond
1696
1697 // ---------------------------------------------------------------------------
1698
1699 //! @cond Doxygen_Suppress
1700
1701 // Look in the authority registry for operations from sourceCRS to targetCRS
1702 // using an intermediate pivot
1703 std::vector<CoordinateOperationNNPtr>
findsOpsInRegistryWithIntermediate(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,bool useCreateBetweenGeodeticCRSWithDatumBasedIntermediates)1704 CoordinateOperationFactory::Private::findsOpsInRegistryWithIntermediate(
1705 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
1706 Private::Context &context,
1707 bool useCreateBetweenGeodeticCRSWithDatumBasedIntermediates) {
1708
1709 #ifdef TRACE_CREATE_OPERATIONS
1710 ENTER_BLOCK("findsOpsInRegistryWithIntermediate(" +
1711 objectAsStr(sourceCRS.get()) + " --> " +
1712 objectAsStr(targetCRS.get()) + ")");
1713 #endif
1714
1715 const auto &authFactory = context.context->getAuthorityFactory();
1716 assert(authFactory);
1717
1718 std::list<std::pair<std::string, std::string>> sourceIds;
1719 std::list<std::pair<std::string, std::string>> targetIds;
1720 buildCRSIds(sourceCRS, context, sourceIds);
1721 buildCRSIds(targetCRS, context, targetIds);
1722
1723 const auto gridAvailabilityUse = context.context->getGridAvailabilityUse();
1724 for (const auto &idSrc : sourceIds) {
1725 const auto &srcAuthName = idSrc.first;
1726 const auto &srcCode = idSrc.second;
1727 for (const auto &idTarget : targetIds) {
1728 const auto &targetAuthName = idTarget.first;
1729 const auto &targetCode = idTarget.second;
1730
1731 const auto authorities(getCandidateAuthorities(
1732 authFactory, srcAuthName, targetAuthName));
1733 assert(!authorities.empty());
1734
1735 const auto tmpAuthFactory = io::AuthorityFactory::create(
1736 authFactory->databaseContext(),
1737 (authFactory->getAuthority() == "any" || authorities.size() > 1)
1738 ? std::string()
1739 : authorities.front());
1740
1741 std::vector<CoordinateOperationNNPtr> res;
1742 if (useCreateBetweenGeodeticCRSWithDatumBasedIntermediates) {
1743 res =
1744 tmpAuthFactory
1745 ->createBetweenGeodeticCRSWithDatumBasedIntermediates(
1746 sourceCRS, srcAuthName, srcCode, targetCRS,
1747 targetAuthName, targetCode,
1748 context.context->getUsePROJAlternativeGridNames(),
1749 gridAvailabilityUse ==
1750 CoordinateOperationContext::
1751 GridAvailabilityUse::
1752 DISCARD_OPERATION_IF_MISSING_GRID ||
1753 gridAvailabilityUse ==
1754 CoordinateOperationContext::
1755 GridAvailabilityUse::KNOWN_AVAILABLE,
1756 gridAvailabilityUse ==
1757 CoordinateOperationContext::
1758 GridAvailabilityUse::KNOWN_AVAILABLE,
1759 context.context->getDiscardSuperseded(),
1760 authFactory->getAuthority() != "any" &&
1761 authorities.size() > 1
1762 ? authorities
1763 : std::vector<std::string>(),
1764 context.extent1, context.extent2);
1765 } else {
1766 io::AuthorityFactory::ObjectType intermediateObjectType =
1767 io::AuthorityFactory::ObjectType::CRS;
1768
1769 // If doing GeogCRS --> GeogCRS, only use GeogCRS as
1770 // intermediate CRS
1771 // Avoid weird behavior when doing NAD83 -> NAD83(2011)
1772 // that would go through NAVD88 otherwise.
1773 if (context.context->getIntermediateCRS().empty() &&
1774 dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get()) &&
1775 dynamic_cast<const crs::GeographicCRS *>(targetCRS.get())) {
1776 intermediateObjectType =
1777 io::AuthorityFactory::ObjectType::GEOGRAPHIC_CRS;
1778 }
1779 res = tmpAuthFactory->createFromCRSCodesWithIntermediates(
1780 srcAuthName, srcCode, targetAuthName, targetCode,
1781 context.context->getUsePROJAlternativeGridNames(),
1782 gridAvailabilityUse ==
1783 CoordinateOperationContext::GridAvailabilityUse::
1784 DISCARD_OPERATION_IF_MISSING_GRID ||
1785 gridAvailabilityUse ==
1786 CoordinateOperationContext::GridAvailabilityUse::
1787 KNOWN_AVAILABLE,
1788 gridAvailabilityUse ==
1789 CoordinateOperationContext::GridAvailabilityUse::
1790 KNOWN_AVAILABLE,
1791 context.context->getDiscardSuperseded(),
1792 context.context->getIntermediateCRS(),
1793 intermediateObjectType,
1794 authFactory->getAuthority() != "any" &&
1795 authorities.size() > 1
1796 ? authorities
1797 : std::vector<std::string>(),
1798 context.extent1, context.extent2);
1799 }
1800 if (!res.empty()) {
1801
1802 auto resFiltered =
1803 FilterResults(res, context.context, context.extent1,
1804 context.extent2, false)
1805 .getRes();
1806 #ifdef TRACE_CREATE_OPERATIONS
1807 logTrace("filtering reduced from " +
1808 toString(static_cast<int>(res.size())) + " to " +
1809 toString(static_cast<int>(resFiltered.size())));
1810 #endif
1811 return resFiltered;
1812 }
1813 }
1814 }
1815 return std::vector<CoordinateOperationNNPtr>();
1816 }
1817 //! @endcond
1818
1819 // ---------------------------------------------------------------------------
1820
1821 //! @cond Doxygen_Suppress
1822 static TransformationNNPtr
createBallparkGeographicOffset(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const io::DatabaseContextPtr & dbContext)1823 createBallparkGeographicOffset(const crs::CRSNNPtr &sourceCRS,
1824 const crs::CRSNNPtr &targetCRS,
1825 const io::DatabaseContextPtr &dbContext) {
1826
1827 const crs::GeographicCRS *geogSrc =
1828 dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
1829 const crs::GeographicCRS *geogDst =
1830 dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
1831 const bool isSameDatum = geogSrc && geogDst &&
1832 geogSrc->datumNonNull(dbContext)->_isEquivalentTo(
1833 geogDst->datumNonNull(dbContext).get(),
1834 util::IComparable::Criterion::EQUIVALENT);
1835
1836 auto name = buildOpName(isSameDatum ? NULL_GEOGRAPHIC_OFFSET
1837 : BALLPARK_GEOGRAPHIC_OFFSET,
1838 sourceCRS, targetCRS);
1839
1840 const auto &sourceCRSExtent = getExtent(sourceCRS);
1841 const auto &targetCRSExtent = getExtent(targetCRS);
1842 const bool sameExtent =
1843 sourceCRSExtent && targetCRSExtent &&
1844 sourceCRSExtent->_isEquivalentTo(
1845 targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);
1846
1847 util::PropertyMap map;
1848 map.set(common::IdentifiedObject::NAME_KEY, name)
1849 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
1850 sameExtent ? NN_NO_CHECK(sourceCRSExtent)
1851 : metadata::Extent::WORLD);
1852 const common::Angle angle0(0);
1853
1854 std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
1855 if (isSameDatum) {
1856 accuracies.emplace_back(metadata::PositionalAccuracy::create("0"));
1857 }
1858
1859 if (dynamic_cast<const crs::SingleCRS *>(sourceCRS.get())
1860 ->coordinateSystem()
1861 ->axisList()
1862 .size() == 3 ||
1863 dynamic_cast<const crs::SingleCRS *>(targetCRS.get())
1864 ->coordinateSystem()
1865 ->axisList()
1866 .size() == 3) {
1867 return Transformation::createGeographic3DOffsets(
1868 map, sourceCRS, targetCRS, angle0, angle0, common::Length(0),
1869 accuracies);
1870 } else {
1871 return Transformation::createGeographic2DOffsets(
1872 map, sourceCRS, targetCRS, angle0, angle0, accuracies);
1873 }
1874 }
1875 //! @endcond
1876
1877 // ---------------------------------------------------------------------------
1878
1879 //! @cond Doxygen_Suppress
1880
1881 // ---------------------------------------------------------------------------
1882
1883 struct MyPROJStringExportableGeodToGeod final
1884 : public io::IPROJStringExportable {
1885 crs::GeodeticCRSPtr geodSrc{};
1886 crs::GeodeticCRSPtr geodDst{};
1887
MyPROJStringExportableGeodToGeodoperation::MyPROJStringExportableGeodToGeod1888 MyPROJStringExportableGeodToGeod(const crs::GeodeticCRSPtr &geodSrcIn,
1889 const crs::GeodeticCRSPtr &geodDstIn)
1890 : geodSrc(geodSrcIn), geodDst(geodDstIn) {}
1891
1892 ~MyPROJStringExportableGeodToGeod() override;
1893
1894 void
1895 // cppcheck-suppress functionStatic
_exportToPROJStringoperation::MyPROJStringExportableGeodToGeod1896 _exportToPROJString(io::PROJStringFormatter *formatter) const override {
1897
1898 formatter->startInversion();
1899 geodSrc->_exportToPROJString(formatter);
1900 formatter->stopInversion();
1901 geodDst->_exportToPROJString(formatter);
1902 }
1903 };
1904
1905 MyPROJStringExportableGeodToGeod::~MyPROJStringExportableGeodToGeod() = default;
1906
1907 // ---------------------------------------------------------------------------
1908
1909 struct MyPROJStringExportableHorizVertical final
1910 : public io::IPROJStringExportable {
1911 CoordinateOperationPtr horizTransform{};
1912 CoordinateOperationPtr verticalTransform{};
1913 crs::GeographicCRSPtr geogDst{};
1914
MyPROJStringExportableHorizVerticaloperation::MyPROJStringExportableHorizVertical1915 MyPROJStringExportableHorizVertical(
1916 const CoordinateOperationPtr &horizTransformIn,
1917 const CoordinateOperationPtr &verticalTransformIn,
1918 const crs::GeographicCRSPtr &geogDstIn)
1919 : horizTransform(horizTransformIn),
1920 verticalTransform(verticalTransformIn), geogDst(geogDstIn) {}
1921
1922 ~MyPROJStringExportableHorizVertical() override;
1923
1924 void
1925 // cppcheck-suppress functionStatic
_exportToPROJStringoperation::MyPROJStringExportableHorizVertical1926 _exportToPROJString(io::PROJStringFormatter *formatter) const override {
1927
1928 formatter->pushOmitZUnitConversion();
1929
1930 horizTransform->_exportToPROJString(formatter);
1931
1932 formatter->startInversion();
1933 geogDst->addAngularUnitConvertAndAxisSwap(formatter);
1934 formatter->stopInversion();
1935
1936 formatter->popOmitZUnitConversion();
1937
1938 formatter->pushOmitHorizontalConversionInVertTransformation();
1939 verticalTransform->_exportToPROJString(formatter);
1940 formatter->popOmitHorizontalConversionInVertTransformation();
1941
1942 formatter->pushOmitZUnitConversion();
1943 geogDst->addAngularUnitConvertAndAxisSwap(formatter);
1944 formatter->popOmitZUnitConversion();
1945 }
1946 };
1947
1948 MyPROJStringExportableHorizVertical::~MyPROJStringExportableHorizVertical() =
1949 default;
1950
1951 // ---------------------------------------------------------------------------
1952
1953 struct MyPROJStringExportableHorizVerticalHorizPROJBased final
1954 : public io::IPROJStringExportable {
1955 CoordinateOperationPtr opSrcCRSToGeogCRS{};
1956 CoordinateOperationPtr verticalTransform{};
1957 CoordinateOperationPtr opGeogCRStoDstCRS{};
1958 crs::GeographicCRSPtr interpolationGeogCRS{};
1959
MyPROJStringExportableHorizVerticalHorizPROJBasedoperation::MyPROJStringExportableHorizVerticalHorizPROJBased1960 MyPROJStringExportableHorizVerticalHorizPROJBased(
1961 const CoordinateOperationPtr &opSrcCRSToGeogCRSIn,
1962 const CoordinateOperationPtr &verticalTransformIn,
1963 const CoordinateOperationPtr &opGeogCRStoDstCRSIn,
1964 const crs::GeographicCRSPtr &interpolationGeogCRSIn)
1965 : opSrcCRSToGeogCRS(opSrcCRSToGeogCRSIn),
1966 verticalTransform(verticalTransformIn),
1967 opGeogCRStoDstCRS(opGeogCRStoDstCRSIn),
1968 interpolationGeogCRS(interpolationGeogCRSIn) {}
1969
1970 ~MyPROJStringExportableHorizVerticalHorizPROJBased() override;
1971
1972 void
1973 // cppcheck-suppress functionStatic
_exportToPROJStringoperation::MyPROJStringExportableHorizVerticalHorizPROJBased1974 _exportToPROJString(io::PROJStringFormatter *formatter) const override {
1975
1976 formatter->pushOmitZUnitConversion();
1977
1978 opSrcCRSToGeogCRS->_exportToPROJString(formatter);
1979
1980 formatter->startInversion();
1981 interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
1982 formatter->stopInversion();
1983
1984 formatter->popOmitZUnitConversion();
1985
1986 formatter->pushOmitHorizontalConversionInVertTransformation();
1987 verticalTransform->_exportToPROJString(formatter);
1988 formatter->popOmitHorizontalConversionInVertTransformation();
1989
1990 formatter->pushOmitZUnitConversion();
1991
1992 interpolationGeogCRS->addAngularUnitConvertAndAxisSwap(formatter);
1993
1994 opGeogCRStoDstCRS->_exportToPROJString(formatter);
1995
1996 formatter->popOmitZUnitConversion();
1997 }
1998 };
1999
2000 MyPROJStringExportableHorizVerticalHorizPROJBased::
2001 ~MyPROJStringExportableHorizVerticalHorizPROJBased() = default;
2002
2003 //! @endcond
2004
2005 } // namespace operation
2006 NS_PROJ_END
2007
2008 #if 0
2009 namespace dropbox{ namespace oxygen {
2010 template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableGeodToGeod>>::~nn() = default;
2011 template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizVertical>>::~nn() = default;
2012 template<> nn<std::shared_ptr<NS_PROJ::operation::MyPROJStringExportableHorizVerticalHorizPROJBased>>::~nn() = default;
2013 }}
2014 #endif
2015
2016 NS_PROJ_START
2017 namespace operation {
2018
2019 //! @cond Doxygen_Suppress
2020
2021 // ---------------------------------------------------------------------------
2022
buildTransfName(const std::string & srcName,const std::string & targetName)2023 static std::string buildTransfName(const std::string &srcName,
2024 const std::string &targetName) {
2025 std::string name("Transformation from ");
2026 name += srcName;
2027 name += " to ";
2028 name += targetName;
2029 return name;
2030 }
2031
2032 // ---------------------------------------------------------------------------
2033
createPROJBased(const util::PropertyMap & properties,const io::IPROJStringExportableNNPtr & projExportable,const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const crs::CRSPtr & interpolationCRS,const std::vector<metadata::PositionalAccuracyNNPtr> & accuracies,bool hasBallparkTransformation)2034 static SingleOperationNNPtr createPROJBased(
2035 const util::PropertyMap &properties,
2036 const io::IPROJStringExportableNNPtr &projExportable,
2037 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
2038 const crs::CRSPtr &interpolationCRS,
2039 const std::vector<metadata::PositionalAccuracyNNPtr> &accuracies,
2040 bool hasBallparkTransformation) {
2041 return util::nn_static_pointer_cast<SingleOperation>(
2042 PROJBasedOperation::create(properties, projExportable, false, sourceCRS,
2043 targetCRS, interpolationCRS, accuracies,
2044 hasBallparkTransformation));
2045 }
2046
2047 // ---------------------------------------------------------------------------
2048
2049 static CoordinateOperationNNPtr
createGeodToGeodPROJBased(const crs::CRSNNPtr & geodSrc,const crs::CRSNNPtr & geodDst)2050 createGeodToGeodPROJBased(const crs::CRSNNPtr &geodSrc,
2051 const crs::CRSNNPtr &geodDst) {
2052
2053 auto exportable = util::nn_make_shared<MyPROJStringExportableGeodToGeod>(
2054 util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(geodSrc),
2055 util::nn_dynamic_pointer_cast<crs::GeodeticCRS>(geodDst));
2056
2057 auto properties = util::PropertyMap().set(
2058 common::IdentifiedObject::NAME_KEY,
2059 buildTransfName(geodSrc->nameStr(), geodDst->nameStr()));
2060 return createPROJBased(properties, exportable, geodSrc, geodDst, nullptr,
2061 {}, false);
2062 }
2063
2064 // ---------------------------------------------------------------------------
2065
2066 static std::string
getRemarks(const std::vector<operation::CoordinateOperationNNPtr> & ops)2067 getRemarks(const std::vector<operation::CoordinateOperationNNPtr> &ops) {
2068 std::string remarks;
2069 for (const auto &op : ops) {
2070 const auto &opRemarks = op->remarks();
2071 if (!opRemarks.empty()) {
2072 if (!remarks.empty()) {
2073 remarks += '\n';
2074 }
2075
2076 std::string opName(op->nameStr());
2077 if (starts_with(opName, INVERSE_OF)) {
2078 opName = opName.substr(INVERSE_OF.size());
2079 }
2080
2081 remarks += "For ";
2082 remarks += opName;
2083
2084 const auto &ids = op->identifiers();
2085 if (!ids.empty()) {
2086 std::string authority(*ids.front()->codeSpace());
2087 if (starts_with(authority, "INVERSE(") &&
2088 authority.back() == ')') {
2089 authority = authority.substr(strlen("INVERSE("),
2090 authority.size() - 1 -
2091 strlen("INVERSE("));
2092 }
2093 if (starts_with(authority, "DERIVED_FROM(") &&
2094 authority.back() == ')') {
2095 authority = authority.substr(strlen("DERIVED_FROM("),
2096 authority.size() - 1 -
2097 strlen("DERIVED_FROM("));
2098 }
2099
2100 remarks += " (";
2101 remarks += authority;
2102 remarks += ':';
2103 remarks += ids.front()->code();
2104 remarks += ')';
2105 }
2106 remarks += ": ";
2107 remarks += opRemarks;
2108 }
2109 }
2110 return remarks;
2111 }
2112
2113 // ---------------------------------------------------------------------------
2114
createHorizVerticalPROJBased(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const operation::CoordinateOperationNNPtr & horizTransform,const operation::CoordinateOperationNNPtr & verticalTransform,bool checkExtent)2115 static CoordinateOperationNNPtr createHorizVerticalPROJBased(
2116 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
2117 const operation::CoordinateOperationNNPtr &horizTransform,
2118 const operation::CoordinateOperationNNPtr &verticalTransform,
2119 bool checkExtent) {
2120
2121 auto geogDst = util::nn_dynamic_pointer_cast<crs::GeographicCRS>(targetCRS);
2122 assert(geogDst);
2123
2124 auto exportable = util::nn_make_shared<MyPROJStringExportableHorizVertical>(
2125 horizTransform, verticalTransform, geogDst);
2126
2127 const bool horizTransformIsNoOp =
2128 starts_with(horizTransform->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
2129 horizTransform->nameStr().find(" + ") == std::string::npos;
2130 if (horizTransformIsNoOp) {
2131 auto properties = util::PropertyMap();
2132 properties.set(common::IdentifiedObject::NAME_KEY,
2133 verticalTransform->nameStr());
2134 bool dummy = false;
2135 auto extent = getExtent(verticalTransform, true, dummy);
2136 if (extent) {
2137 properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2138 NN_NO_CHECK(extent));
2139 }
2140 const auto &remarks = verticalTransform->remarks();
2141 if (!remarks.empty()) {
2142 properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
2143 }
2144 return createPROJBased(
2145 properties, exportable, sourceCRS, targetCRS, nullptr,
2146 verticalTransform->coordinateOperationAccuracies(),
2147 verticalTransform->hasBallparkTransformation());
2148 } else {
2149 bool emptyIntersection = false;
2150 auto ops = std::vector<CoordinateOperationNNPtr>{horizTransform,
2151 verticalTransform};
2152 auto extent = getExtent(ops, true, emptyIntersection);
2153 if (checkExtent && emptyIntersection) {
2154 std::string msg(
2155 "empty intersection of area of validity of concatenated "
2156 "operations");
2157 throw InvalidOperationEmptyIntersection(msg);
2158 }
2159 auto properties = util::PropertyMap();
2160 properties.set(common::IdentifiedObject::NAME_KEY,
2161 computeConcatenatedName(ops));
2162
2163 if (extent) {
2164 properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2165 NN_NO_CHECK(extent));
2166 }
2167
2168 const auto remarks = getRemarks(ops);
2169 if (!remarks.empty()) {
2170 properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
2171 }
2172
2173 std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
2174 const double accuracy = getAccuracy(ops);
2175 if (accuracy >= 0.0) {
2176 accuracies.emplace_back(
2177 metadata::PositionalAccuracy::create(toString(accuracy)));
2178 }
2179
2180 return createPROJBased(
2181 properties, exportable, sourceCRS, targetCRS, nullptr, accuracies,
2182 horizTransform->hasBallparkTransformation() ||
2183 verticalTransform->hasBallparkTransformation());
2184 }
2185 }
2186
2187 // ---------------------------------------------------------------------------
2188
createHorizVerticalHorizPROJBased(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const operation::CoordinateOperationNNPtr & opSrcCRSToGeogCRS,const operation::CoordinateOperationNNPtr & verticalTransform,const operation::CoordinateOperationNNPtr & opGeogCRStoDstCRS,const crs::GeographicCRSPtr & interpolationGeogCRS,bool checkExtent)2189 static CoordinateOperationNNPtr createHorizVerticalHorizPROJBased(
2190 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
2191 const operation::CoordinateOperationNNPtr &opSrcCRSToGeogCRS,
2192 const operation::CoordinateOperationNNPtr &verticalTransform,
2193 const operation::CoordinateOperationNNPtr &opGeogCRStoDstCRS,
2194 const crs::GeographicCRSPtr &interpolationGeogCRS, bool checkExtent) {
2195
2196 auto exportable =
2197 util::nn_make_shared<MyPROJStringExportableHorizVerticalHorizPROJBased>(
2198 opSrcCRSToGeogCRS, verticalTransform, opGeogCRStoDstCRS,
2199 interpolationGeogCRS);
2200
2201 std::vector<CoordinateOperationNNPtr> ops;
2202 if (!(starts_with(opSrcCRSToGeogCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
2203 opSrcCRSToGeogCRS->nameStr().find(" + ") == std::string::npos)) {
2204 ops.emplace_back(opSrcCRSToGeogCRS);
2205 }
2206 ops.emplace_back(verticalTransform);
2207 if (!(starts_with(opGeogCRStoDstCRS->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
2208 opGeogCRStoDstCRS->nameStr().find(" + ") == std::string::npos)) {
2209 ops.emplace_back(opGeogCRStoDstCRS);
2210 }
2211
2212 bool hasBallparkTransformation = false;
2213 for (const auto &op : ops) {
2214 hasBallparkTransformation |= op->hasBallparkTransformation();
2215 }
2216 bool emptyIntersection = false;
2217 auto extent = getExtent(ops, false, emptyIntersection);
2218 if (checkExtent && emptyIntersection) {
2219 std::string msg(
2220 "empty intersection of area of validity of concatenated "
2221 "operations");
2222 throw InvalidOperationEmptyIntersection(msg);
2223 }
2224 auto properties = util::PropertyMap();
2225 properties.set(common::IdentifiedObject::NAME_KEY,
2226 computeConcatenatedName(ops));
2227
2228 if (extent) {
2229 properties.set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2230 NN_NO_CHECK(extent));
2231 }
2232
2233 const auto remarks = getRemarks(ops);
2234 if (!remarks.empty()) {
2235 properties.set(common::IdentifiedObject::REMARKS_KEY, remarks);
2236 }
2237
2238 std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
2239 const double accuracy = getAccuracy(ops);
2240 if (accuracy >= 0.0) {
2241 accuracies.emplace_back(
2242 metadata::PositionalAccuracy::create(toString(accuracy)));
2243 }
2244
2245 return createPROJBased(properties, exportable, sourceCRS, targetCRS,
2246 nullptr, accuracies, hasBallparkTransformation);
2247 }
2248
2249 //! @endcond
2250
2251 // ---------------------------------------------------------------------------
2252
2253 //! @cond Doxygen_Suppress
2254
2255 std::vector<CoordinateOperationNNPtr>
createOperationsGeogToGeog(std::vector<CoordinateOperationNNPtr> & res,const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::GeographicCRS * geogSrc,const crs::GeographicCRS * geogDst)2256 CoordinateOperationFactory::Private::createOperationsGeogToGeog(
2257 std::vector<CoordinateOperationNNPtr> &res, const crs::CRSNNPtr &sourceCRS,
2258 const crs::CRSNNPtr &targetCRS, Private::Context &context,
2259 const crs::GeographicCRS *geogSrc, const crs::GeographicCRS *geogDst) {
2260
2261 assert(sourceCRS.get() == geogSrc);
2262 assert(targetCRS.get() == geogDst);
2263
2264 const auto &src_pm = geogSrc->primeMeridian()->longitude();
2265 const auto &dst_pm = geogDst->primeMeridian()->longitude();
2266 auto offset_pm =
2267 (src_pm.unit() == dst_pm.unit())
2268 ? common::Angle(src_pm.value() - dst_pm.value(), src_pm.unit())
2269 : common::Angle(
2270 src_pm.convertToUnit(common::UnitOfMeasure::DEGREE) -
2271 dst_pm.convertToUnit(common::UnitOfMeasure::DEGREE),
2272 common::UnitOfMeasure::DEGREE);
2273
2274 double vconvSrc = 1.0;
2275 const auto &srcCS = geogSrc->coordinateSystem();
2276 const auto &srcAxisList = srcCS->axisList();
2277 if (srcAxisList.size() == 3) {
2278 vconvSrc = srcAxisList[2]->unit().conversionToSI();
2279 }
2280 double vconvDst = 1.0;
2281 const auto &dstCS = geogDst->coordinateSystem();
2282 const auto &dstAxisList = dstCS->axisList();
2283 if (dstAxisList.size() == 3) {
2284 vconvDst = dstAxisList[2]->unit().conversionToSI();
2285 }
2286
2287 std::string name(buildTransfName(geogSrc->nameStr(), geogDst->nameStr()));
2288
2289 const auto &authFactory = context.context->getAuthorityFactory();
2290 const auto dbContext =
2291 authFactory ? authFactory->databaseContext().as_nullable() : nullptr;
2292
2293 const bool sameDatum = geogSrc->datumNonNull(dbContext)->_isEquivalentTo(
2294 geogDst->datumNonNull(dbContext).get(),
2295 util::IComparable::Criterion::EQUIVALENT);
2296
2297 // Do the CRS differ by their axis order ?
2298 bool axisReversal2D = false;
2299 bool axisReversal3D = false;
2300 if (!srcCS->_isEquivalentTo(dstCS.get(),
2301 util::IComparable::Criterion::EQUIVALENT)) {
2302 auto srcOrder = srcCS->axisOrder();
2303 auto dstOrder = dstCS->axisOrder();
2304 if (((srcOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST ||
2305 srcOrder == cs::EllipsoidalCS::AxisOrder::
2306 LAT_NORTH_LONG_EAST_HEIGHT_UP) &&
2307 (dstOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ||
2308 dstOrder == cs::EllipsoidalCS::AxisOrder::
2309 LONG_EAST_LAT_NORTH_HEIGHT_UP)) ||
2310 ((srcOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ||
2311 srcOrder == cs::EllipsoidalCS::AxisOrder::
2312 LONG_EAST_LAT_NORTH_HEIGHT_UP) &&
2313 (dstOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST ||
2314 dstOrder == cs::EllipsoidalCS::AxisOrder::
2315 LAT_NORTH_LONG_EAST_HEIGHT_UP))) {
2316 if (srcAxisList.size() == 3 || dstAxisList.size() == 3)
2317 axisReversal3D = true;
2318 else
2319 axisReversal2D = true;
2320 }
2321 }
2322
2323 // Do they differ by vertical units ?
2324 if (vconvSrc != vconvDst &&
2325 geogSrc->ellipsoid()->_isEquivalentTo(
2326 geogDst->ellipsoid().get(),
2327 util::IComparable::Criterion::EQUIVALENT)) {
2328 if (offset_pm.value() == 0 && !axisReversal2D && !axisReversal3D) {
2329 // If only by vertical units, use a Change of Vertical
2330 // Unit
2331 // transformation
2332 const double factor = vconvSrc / vconvDst;
2333 auto conv = Conversion::createChangeVerticalUnit(
2334 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
2335 name),
2336 common::Scale(factor));
2337 conv->setCRSs(sourceCRS, targetCRS, nullptr);
2338 conv->setHasBallparkTransformation(!sameDatum);
2339 res.push_back(conv);
2340 return res;
2341 } else {
2342 auto op = createGeodToGeodPROJBased(sourceCRS, targetCRS);
2343 op->setHasBallparkTransformation(!sameDatum);
2344 res.emplace_back(op);
2345 return res;
2346 }
2347 }
2348
2349 // Do the CRS differ only by their axis order ?
2350 if (sameDatum && (axisReversal2D || axisReversal3D)) {
2351 auto conv = Conversion::createAxisOrderReversal(axisReversal3D);
2352 conv->setCRSs(sourceCRS, targetCRS, nullptr);
2353 res.emplace_back(conv);
2354 return res;
2355 }
2356
2357 std::vector<CoordinateOperationNNPtr> steps;
2358 // If both are geographic and only differ by their prime
2359 // meridian,
2360 // apply a longitude rotation transformation.
2361 if (geogSrc->ellipsoid()->_isEquivalentTo(
2362 geogDst->ellipsoid().get(),
2363 util::IComparable::Criterion::EQUIVALENT) &&
2364 src_pm.getSIValue() != dst_pm.getSIValue()) {
2365
2366 steps.emplace_back(Transformation::createLongitudeRotation(
2367 util::PropertyMap()
2368 .set(common::IdentifiedObject::NAME_KEY, name)
2369 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2370 metadata::Extent::WORLD),
2371 sourceCRS, targetCRS, offset_pm));
2372 // If only the target has a non-zero prime meridian, chain a
2373 // null geographic offset and then the longitude rotation
2374 } else if (src_pm.getSIValue() == 0 && dst_pm.getSIValue() != 0) {
2375 auto datum = datum::GeodeticReferenceFrame::create(
2376 util::PropertyMap(), geogDst->ellipsoid(),
2377 util::optional<std::string>(), geogSrc->primeMeridian());
2378 std::string interm_crs_name(geogDst->nameStr());
2379 interm_crs_name += " altered to use prime meridian of ";
2380 interm_crs_name += geogSrc->nameStr();
2381 auto interm_crs =
2382 util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
2383 util::PropertyMap()
2384 .set(common::IdentifiedObject::NAME_KEY, interm_crs_name)
2385 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2386 metadata::Extent::WORLD),
2387 datum, dstCS));
2388
2389 steps.emplace_back(
2390 createBallparkGeographicOffset(sourceCRS, interm_crs, dbContext));
2391
2392 steps.emplace_back(Transformation::createLongitudeRotation(
2393 util::PropertyMap()
2394 .set(common::IdentifiedObject::NAME_KEY,
2395 buildTransfName(geogSrc->nameStr(), interm_crs->nameStr()))
2396 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2397 metadata::Extent::WORLD),
2398 interm_crs, targetCRS, offset_pm));
2399
2400 } else {
2401 // If the prime meridians are different, chain a longitude
2402 // rotation and the null geographic offset.
2403 if (src_pm.getSIValue() != dst_pm.getSIValue()) {
2404 auto datum = datum::GeodeticReferenceFrame::create(
2405 util::PropertyMap(), geogSrc->ellipsoid(),
2406 util::optional<std::string>(), geogDst->primeMeridian());
2407 std::string interm_crs_name(geogSrc->nameStr());
2408 interm_crs_name += " altered to use prime meridian of ";
2409 interm_crs_name += geogDst->nameStr();
2410 auto interm_crs = util::nn_static_pointer_cast<crs::CRS>(
2411 crs::GeographicCRS::create(
2412 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
2413 interm_crs_name),
2414 datum, srcCS));
2415
2416 steps.emplace_back(Transformation::createLongitudeRotation(
2417 util::PropertyMap()
2418 .set(common::IdentifiedObject::NAME_KEY,
2419 buildTransfName(geogSrc->nameStr(),
2420 interm_crs->nameStr()))
2421 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2422 metadata::Extent::WORLD),
2423 sourceCRS, interm_crs, offset_pm));
2424 steps.emplace_back(createBallparkGeographicOffset(
2425 interm_crs, targetCRS, dbContext));
2426 } else {
2427 steps.emplace_back(createBallparkGeographicOffset(
2428 sourceCRS, targetCRS, dbContext));
2429 }
2430 }
2431
2432 auto op = ConcatenatedOperation::createComputeMetadata(
2433 steps, disallowEmptyIntersection);
2434 op->setHasBallparkTransformation(!sameDatum);
2435 res.emplace_back(op);
2436 return res;
2437 }
2438
2439 // ---------------------------------------------------------------------------
2440
hasIdentifiers(const CoordinateOperationNNPtr & op)2441 static bool hasIdentifiers(const CoordinateOperationNNPtr &op) {
2442 if (!op->identifiers().empty()) {
2443 return true;
2444 }
2445 auto concatenated = dynamic_cast<const ConcatenatedOperation *>(op.get());
2446 if (concatenated) {
2447 for (const auto &subOp : concatenated->operations()) {
2448 if (hasIdentifiers(subOp)) {
2449 return true;
2450 }
2451 }
2452 }
2453 return false;
2454 }
2455
2456 // ---------------------------------------------------------------------------
2457
2458 static std::vector<crs::CRSNNPtr>
findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr & authFactory,const crs::GeodeticCRS * crs,const datum::GeodeticReferenceFrame * datum)2459 findCandidateGeodCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
2460 const crs::GeodeticCRS *crs,
2461 const datum::GeodeticReferenceFrame *datum) {
2462 std::vector<crs::CRSNNPtr> candidates;
2463 assert(datum);
2464 const auto &ids = datum->identifiers();
2465 const auto &datumName = datum->nameStr();
2466 if (!ids.empty()) {
2467 for (const auto &id : ids) {
2468 const auto &authName = *(id->codeSpace());
2469 const auto &code = id->code();
2470 if (!authName.empty()) {
2471 const auto crsIds = crs->identifiers();
2472 const auto tmpFactory =
2473 (crsIds.size() == 1 &&
2474 *(crsIds.front()->codeSpace()) == authName)
2475 ? io::AuthorityFactory::create(
2476 authFactory->databaseContext(), authName)
2477 .as_nullable()
2478 : authFactory;
2479 auto l_candidates = tmpFactory->createGeodeticCRSFromDatum(
2480 authName, code, std::string());
2481 for (const auto &candidate : l_candidates) {
2482 candidates.emplace_back(candidate);
2483 }
2484 }
2485 }
2486 } else if (datumName != "unknown" && datumName != "unnamed") {
2487 auto matches = authFactory->createObjectsFromName(
2488 datumName,
2489 {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME}, false,
2490 2);
2491 if (matches.size() == 1) {
2492 const auto &match = matches.front();
2493 if (datum->_isEquivalentTo(
2494 match.get(), util::IComparable::Criterion::EQUIVALENT) &&
2495 !match->identifiers().empty()) {
2496 return findCandidateGeodCRSForDatum(
2497 authFactory, crs,
2498 dynamic_cast<const datum::GeodeticReferenceFrame *>(
2499 match.get()));
2500 }
2501 }
2502 }
2503 return candidates;
2504 }
2505
2506 // ---------------------------------------------------------------------------
2507
setCRSs(CoordinateOperation * co,const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS)2508 void CoordinateOperationFactory::Private::setCRSs(
2509 CoordinateOperation *co, const crs::CRSNNPtr &sourceCRS,
2510 const crs::CRSNNPtr &targetCRS) {
2511 co->setCRSs(sourceCRS, targetCRS, nullptr);
2512
2513 auto invCO = dynamic_cast<InverseCoordinateOperation *>(co);
2514 if (invCO) {
2515 invCO->forwardOperation()->setCRSs(targetCRS, sourceCRS, nullptr);
2516 }
2517
2518 auto transf = dynamic_cast<Transformation *>(co);
2519 if (transf) {
2520 transf->inverseAsTransformation()->setCRSs(targetCRS, sourceCRS,
2521 nullptr);
2522 }
2523
2524 auto concat = dynamic_cast<ConcatenatedOperation *>(co);
2525 if (concat) {
2526 auto first = concat->operations().front().get();
2527 auto &firstTarget(first->targetCRS());
2528 if (firstTarget) {
2529 setCRSs(first, sourceCRS, NN_NO_CHECK(firstTarget));
2530 }
2531 auto last = concat->operations().back().get();
2532 auto &lastSource(last->sourceCRS());
2533 if (lastSource) {
2534 setCRSs(last, NN_NO_CHECK(lastSource), targetCRS);
2535 }
2536 }
2537 }
2538
2539 // ---------------------------------------------------------------------------
2540
hasResultSetOnlyResultsWithPROJStep(const std::vector<CoordinateOperationNNPtr> & res)2541 static bool hasResultSetOnlyResultsWithPROJStep(
2542 const std::vector<CoordinateOperationNNPtr> &res) {
2543 for (const auto &op : res) {
2544 auto concat = dynamic_cast<const ConcatenatedOperation *>(op.get());
2545 if (concat) {
2546 bool hasPROJStep = false;
2547 const auto &steps = concat->operations();
2548 for (const auto &step : steps) {
2549 const auto &ids = step->identifiers();
2550 if (!ids.empty()) {
2551 const auto &opAuthority = *(ids.front()->codeSpace());
2552 if (opAuthority == "PROJ" ||
2553 opAuthority == "INVERSE(PROJ)" ||
2554 opAuthority == "DERIVED_FROM(PROJ)") {
2555 hasPROJStep = true;
2556 break;
2557 }
2558 }
2559 }
2560 if (!hasPROJStep) {
2561 return false;
2562 }
2563 } else {
2564 return false;
2565 }
2566 }
2567 return true;
2568 }
2569
2570 // ---------------------------------------------------------------------------
2571
createOperationsWithDatumPivot(std::vector<CoordinateOperationNNPtr> & res,const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const crs::GeodeticCRS * geodSrc,const crs::GeodeticCRS * geodDst,Private::Context & context)2572 void CoordinateOperationFactory::Private::createOperationsWithDatumPivot(
2573 std::vector<CoordinateOperationNNPtr> &res, const crs::CRSNNPtr &sourceCRS,
2574 const crs::CRSNNPtr &targetCRS, const crs::GeodeticCRS *geodSrc,
2575 const crs::GeodeticCRS *geodDst, Private::Context &context) {
2576
2577 #ifdef TRACE_CREATE_OPERATIONS
2578 ENTER_BLOCK("createOperationsWithDatumPivot(" +
2579 objectAsStr(sourceCRS.get()) + "," +
2580 objectAsStr(targetCRS.get()) + ")");
2581 #endif
2582
2583 struct CreateOperationsWithDatumPivotAntiRecursion {
2584 Context &context;
2585
2586 explicit CreateOperationsWithDatumPivotAntiRecursion(Context &contextIn)
2587 : context(contextIn) {
2588 assert(!context.inCreateOperationsWithDatumPivotAntiRecursion);
2589 context.inCreateOperationsWithDatumPivotAntiRecursion = true;
2590 }
2591
2592 ~CreateOperationsWithDatumPivotAntiRecursion() {
2593 context.inCreateOperationsWithDatumPivotAntiRecursion = false;
2594 }
2595 };
2596 CreateOperationsWithDatumPivotAntiRecursion guard(context);
2597
2598 const auto &authFactory = context.context->getAuthorityFactory();
2599 const auto &dbContext = authFactory->databaseContext();
2600
2601 const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
2602 authFactory, geodSrc,
2603 geodSrc->datumNonNull(dbContext.as_nullable()).get()));
2604 const auto candidatesDstGeod(findCandidateGeodCRSForDatum(
2605 authFactory, geodDst,
2606 geodDst->datumNonNull(dbContext.as_nullable()).get()));
2607
2608 const bool sourceAndTargetAre3D =
2609 geodSrc->coordinateSystem()->axisList().size() == 3 &&
2610 geodDst->coordinateSystem()->axisList().size() == 3;
2611
2612 auto createTransformations = [&](const crs::CRSNNPtr &candidateSrcGeod,
2613 const crs::CRSNNPtr &candidateDstGeod,
2614 const CoordinateOperationNNPtr &opFirst,
2615 bool isNullFirst) {
2616 const auto opsSecond =
2617 createOperations(candidateSrcGeod, candidateDstGeod, context);
2618 const auto opsThird =
2619 createOperations(candidateDstGeod, targetCRS, context);
2620 assert(!opsThird.empty());
2621
2622 for (auto &opSecond : opsSecond) {
2623 // Check that it is not a transformation synthetized by
2624 // ourselves
2625 if (!hasIdentifiers(opSecond)) {
2626 continue;
2627 }
2628 // And even if it is a referenced transformation, check that
2629 // it is not a trivial one
2630 auto so = dynamic_cast<const SingleOperation *>(opSecond.get());
2631 if (so && isAxisOrderReversal(so->method()->getEPSGCode())) {
2632 continue;
2633 }
2634
2635 std::vector<CoordinateOperationNNPtr> subOps;
2636 const bool isNullThird =
2637 isNullTransformation(opsThird[0]->nameStr());
2638 CoordinateOperationNNPtr opSecondCloned(
2639 (isNullFirst || isNullThird || sourceAndTargetAre3D)
2640 ? opSecond->shallowClone()
2641 : opSecond);
2642 if (isNullFirst || isNullThird) {
2643 if (opSecondCloned->identifiers().size() == 1 &&
2644 (*opSecondCloned->identifiers()[0]->codeSpace())
2645 .find("DERIVED_FROM") == std::string::npos) {
2646 {
2647 util::PropertyMap map;
2648 addModifiedIdentifier(map, opSecondCloned.get(), false,
2649 true);
2650 opSecondCloned->setProperties(map);
2651 }
2652 auto invCO = dynamic_cast<InverseCoordinateOperation *>(
2653 opSecondCloned.get());
2654 if (invCO) {
2655 auto invCOForward = invCO->forwardOperation().get();
2656 if (invCOForward->identifiers().size() == 1 &&
2657 (*invCOForward->identifiers()[0]->codeSpace())
2658 .find("DERIVED_FROM") ==
2659 std::string::npos) {
2660 util::PropertyMap map;
2661 addModifiedIdentifier(map, invCOForward, false,
2662 true);
2663 invCOForward->setProperties(map);
2664 }
2665 }
2666 }
2667 }
2668 if (sourceAndTargetAre3D) {
2669 opSecondCloned->getPrivate()->use3DHelmert_ = true;
2670 auto invCO = dynamic_cast<InverseCoordinateOperation *>(
2671 opSecondCloned.get());
2672 if (invCO) {
2673 auto invCOForward = invCO->forwardOperation().get();
2674 invCOForward->getPrivate()->use3DHelmert_ = true;
2675 }
2676 }
2677 if (isNullFirst) {
2678 auto oldTarget(NN_CHECK_ASSERT(opSecondCloned->targetCRS()));
2679 setCRSs(opSecondCloned.get(), sourceCRS, oldTarget);
2680 } else {
2681 subOps.emplace_back(opFirst);
2682 }
2683 if (isNullThird) {
2684 auto oldSource(NN_CHECK_ASSERT(opSecondCloned->sourceCRS()));
2685 setCRSs(opSecondCloned.get(), oldSource, targetCRS);
2686 subOps.emplace_back(opSecondCloned);
2687 } else {
2688 subOps.emplace_back(opSecondCloned);
2689 subOps.emplace_back(opsThird[0]);
2690 }
2691 #ifdef TRACE_CREATE_OPERATIONS
2692 std::string debugStr;
2693 for (const auto &op : subOps) {
2694 if (!debugStr.empty()) {
2695 debugStr += " + ";
2696 }
2697 debugStr += objectAsStr(op.get());
2698 debugStr += " (";
2699 debugStr += objectAsStr(op->sourceCRS().get());
2700 debugStr += "->";
2701 debugStr += objectAsStr(op->targetCRS().get());
2702 debugStr += ")";
2703 }
2704 logTrace("transformation " + debugStr);
2705 #endif
2706 try {
2707 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
2708 subOps, disallowEmptyIntersection));
2709 } catch (const InvalidOperationEmptyIntersection &) {
2710 }
2711 }
2712 };
2713
2714 // Start in priority with candidates that have exactly the same name as
2715 // the sourcCRS and targetCRS. Typically for the case of init=IGNF:XXXX
2716
2717 // Transformation from IGNF:NTFP to IGNF:RGF93G,
2718 // using
2719 // NTF geographiques Paris (gr) vers NTF GEOGRAPHIQUES GREENWICH (DMS) +
2720 // NOUVELLE TRIANGULATION DE LA FRANCE (NTF) vers RGF93 (ETRS89)
2721 // that is using ntf_r93.gsb, is horribly dependent
2722 // of IGNF:RGF93G being returned before IGNF:RGF93GEO in candidatesDstGeod.
2723 // If RGF93GEO is returned before then we go through WGS84 and use
2724 // instead a Helmert transformation.
2725 // The below logic is thus quite fragile, and attempts at changing it
2726 // result in degraded results for other use cases...
2727
2728 for (const auto &candidateSrcGeod : candidatesSrcGeod) {
2729 if (candidateSrcGeod->nameStr() == sourceCRS->nameStr()) {
2730 for (const auto &candidateDstGeod : candidatesDstGeod) {
2731 if (candidateDstGeod->nameStr() == targetCRS->nameStr()) {
2732 #ifdef TRACE_CREATE_OPERATIONS
2733 ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
2734 objectAsStr(candidateSrcGeod.get()) + "->" +
2735 objectAsStr(candidateDstGeod.get()) + "->" +
2736 objectAsStr(targetCRS.get()) + ")");
2737 #endif
2738 const auto opsFirst =
2739 createOperations(sourceCRS, candidateSrcGeod, context);
2740 assert(!opsFirst.empty());
2741 const bool isNullFirst =
2742 isNullTransformation(opsFirst[0]->nameStr());
2743 createTransformations(candidateSrcGeod, candidateDstGeod,
2744 opsFirst[0], isNullFirst);
2745 if (!res.empty()) {
2746 if (hasResultSetOnlyResultsWithPROJStep(res)) {
2747 continue;
2748 }
2749 return;
2750 }
2751 }
2752 }
2753 }
2754 }
2755
2756 for (const auto &candidateSrcGeod : candidatesSrcGeod) {
2757 const bool bSameSrcName =
2758 candidateSrcGeod->nameStr() == sourceCRS->nameStr();
2759 #ifdef TRACE_CREATE_OPERATIONS
2760 ENTER_BLOCK("");
2761 #endif
2762 const auto opsFirst =
2763 createOperations(sourceCRS, candidateSrcGeod, context);
2764 assert(!opsFirst.empty());
2765 const bool isNullFirst = isNullTransformation(opsFirst[0]->nameStr());
2766
2767 for (const auto &candidateDstGeod : candidatesDstGeod) {
2768 if (bSameSrcName &&
2769 candidateDstGeod->nameStr() == targetCRS->nameStr()) {
2770 continue;
2771 }
2772
2773 #ifdef TRACE_CREATE_OPERATIONS
2774 ENTER_BLOCK("try " + objectAsStr(sourceCRS.get()) + "->" +
2775 objectAsStr(candidateSrcGeod.get()) + "->" +
2776 objectAsStr(candidateDstGeod.get()) + "->" +
2777 objectAsStr(targetCRS.get()) + ")");
2778 #endif
2779 createTransformations(candidateSrcGeod, candidateDstGeod,
2780 opsFirst[0], isNullFirst);
2781 if (!res.empty() && !hasResultSetOnlyResultsWithPROJStep(res)) {
2782 return;
2783 }
2784 }
2785 }
2786 }
2787
2788 // ---------------------------------------------------------------------------
2789
2790 static CoordinateOperationNNPtr
createBallparkGeocentricTranslation(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS)2791 createBallparkGeocentricTranslation(const crs::CRSNNPtr &sourceCRS,
2792 const crs::CRSNNPtr &targetCRS) {
2793 std::string name(BALLPARK_GEOCENTRIC_TRANSLATION);
2794 name += " from ";
2795 name += sourceCRS->nameStr();
2796 name += " to ";
2797 name += targetCRS->nameStr();
2798
2799 return util::nn_static_pointer_cast<CoordinateOperation>(
2800 Transformation::createGeocentricTranslations(
2801 util::PropertyMap()
2802 .set(common::IdentifiedObject::NAME_KEY, name)
2803 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
2804 metadata::Extent::WORLD),
2805 sourceCRS, targetCRS, 0.0, 0.0, 0.0, {}));
2806 }
2807
2808 // ---------------------------------------------------------------------------
2809
hasPerfectAccuracyResult(const std::vector<CoordinateOperationNNPtr> & res,const Context & context)2810 bool CoordinateOperationFactory::Private::hasPerfectAccuracyResult(
2811 const std::vector<CoordinateOperationNNPtr> &res, const Context &context) {
2812 auto resTmp = FilterResults(res, context.context, context.extent1,
2813 context.extent2, true)
2814 .getRes();
2815 for (const auto &op : resTmp) {
2816 const double acc = getAccuracy(op);
2817 if (acc == 0.0) {
2818 return true;
2819 }
2820 }
2821 return false;
2822 }
2823
2824 // ---------------------------------------------------------------------------
2825
2826 std::vector<CoordinateOperationNNPtr>
createOperations(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context)2827 CoordinateOperationFactory::Private::createOperations(
2828 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
2829 Private::Context &context) {
2830
2831 #ifdef TRACE_CREATE_OPERATIONS
2832 ENTER_BLOCK("createOperations(" + objectAsStr(sourceCRS.get()) + " --> " +
2833 objectAsStr(targetCRS.get()) + ")");
2834 #endif
2835
2836 std::vector<CoordinateOperationNNPtr> res;
2837
2838 auto boundSrc = dynamic_cast<const crs::BoundCRS *>(sourceCRS.get());
2839 auto boundDst = dynamic_cast<const crs::BoundCRS *>(targetCRS.get());
2840
2841 const auto &sourceProj4Ext = boundSrc
2842 ? boundSrc->baseCRS()->getExtensionProj4()
2843 : sourceCRS->getExtensionProj4();
2844 const auto &targetProj4Ext = boundDst
2845 ? boundDst->baseCRS()->getExtensionProj4()
2846 : targetCRS->getExtensionProj4();
2847 if (!sourceProj4Ext.empty() || !targetProj4Ext.empty()) {
2848 createOperationsFromProj4Ext(sourceCRS, targetCRS, boundSrc, boundDst,
2849 res);
2850 return res;
2851 }
2852
2853 auto geodSrc = dynamic_cast<const crs::GeodeticCRS *>(sourceCRS.get());
2854 auto geodDst = dynamic_cast<const crs::GeodeticCRS *>(targetCRS.get());
2855 auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
2856 auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
2857 auto vertSrc = dynamic_cast<const crs::VerticalCRS *>(sourceCRS.get());
2858 auto vertDst = dynamic_cast<const crs::VerticalCRS *>(targetCRS.get());
2859
2860 // First look-up if the registry provide us with operations.
2861 auto derivedSrc = dynamic_cast<const crs::DerivedCRS *>(sourceCRS.get());
2862 auto derivedDst = dynamic_cast<const crs::DerivedCRS *>(targetCRS.get());
2863 const auto &authFactory = context.context->getAuthorityFactory();
2864 if (authFactory &&
2865 (derivedSrc == nullptr ||
2866 !derivedSrc->baseCRS()->_isEquivalentTo(
2867 targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) &&
2868 (derivedDst == nullptr ||
2869 !derivedDst->baseCRS()->_isEquivalentTo(
2870 sourceCRS.get(), util::IComparable::Criterion::EQUIVALENT))) {
2871
2872 if (createOperationsFromDatabase(sourceCRS, targetCRS, context, geodSrc,
2873 geodDst, geogSrc, geogDst, vertSrc,
2874 vertDst, res)) {
2875 return res;
2876 }
2877 }
2878
2879 // Special case if both CRS are geodetic
2880 if (geodSrc && geodDst && !derivedSrc && !derivedDst) {
2881 createOperationsGeodToGeod(sourceCRS, targetCRS, context, geodSrc,
2882 geodDst, res);
2883 return res;
2884 }
2885
2886 // If the source is a derived CRS, then chain the inverse of its
2887 // deriving conversion, with transforms from its baseCRS to the
2888 // targetCRS
2889 if (derivedSrc) {
2890 createOperationsDerivedTo(sourceCRS, targetCRS, context, derivedSrc,
2891 res);
2892 return res;
2893 }
2894
2895 // reverse of previous case
2896 if (derivedDst) {
2897 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2898 }
2899
2900 // Order of comparison between the geogDst vs geodDst is impotant
2901 if (boundSrc && geogDst) {
2902 createOperationsBoundToGeog(sourceCRS, targetCRS, context, boundSrc,
2903 geogDst, res);
2904 return res;
2905 } else if (boundSrc && geodDst) {
2906 createOperationsToGeod(sourceCRS, targetCRS, context, geodDst, res);
2907 return res;
2908 }
2909
2910 // reverse of previous case
2911 if (geodSrc && boundDst) {
2912 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2913 }
2914
2915 // vertCRS (as boundCRS with transformation to target vertCRS) to
2916 // vertCRS
2917 if (boundSrc && vertDst) {
2918 createOperationsBoundToVert(sourceCRS, targetCRS, context, boundSrc,
2919 vertDst, res);
2920 return res;
2921 }
2922
2923 // reverse of previous case
2924 if (boundDst && vertSrc) {
2925 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2926 }
2927
2928 if (vertSrc && vertDst) {
2929 createOperationsVertToVert(sourceCRS, targetCRS, context, vertSrc,
2930 vertDst, res);
2931 return res;
2932 }
2933
2934 // A bit odd case as we are comparing apples to oranges, but in case
2935 // the vertical unit differ, do something useful.
2936 if (vertSrc && geogDst) {
2937 createOperationsVertToGeog(sourceCRS, targetCRS, context, vertSrc,
2938 geogDst, res);
2939 return res;
2940 }
2941
2942 // reverse of previous case
2943 if (vertDst && geogSrc) {
2944 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2945 }
2946
2947 // boundCRS to boundCRS
2948 if (boundSrc && boundDst) {
2949 createOperationsBoundToBound(sourceCRS, targetCRS, context, boundSrc,
2950 boundDst, res);
2951 return res;
2952 }
2953
2954 auto compoundSrc = dynamic_cast<crs::CompoundCRS *>(sourceCRS.get());
2955 // Order of comparison between the geogDst vs geodDst is impotant
2956 if (compoundSrc && geogDst) {
2957 createOperationsCompoundToGeog(sourceCRS, targetCRS, context,
2958 compoundSrc, geogDst, res);
2959 return res;
2960 } else if (compoundSrc && geodDst) {
2961 createOperationsToGeod(sourceCRS, targetCRS, context, geodDst, res);
2962 return res;
2963 }
2964
2965 // reverse of previous cases
2966 auto compoundDst = dynamic_cast<const crs::CompoundCRS *>(targetCRS.get());
2967 if (geodSrc && compoundDst) {
2968 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2969 }
2970
2971 if (compoundSrc && compoundDst) {
2972 createOperationsCompoundToCompound(sourceCRS, targetCRS, context,
2973 compoundSrc, compoundDst, res);
2974 return res;
2975 }
2976
2977 // '+proj=longlat +ellps=GRS67 +nadgrids=@foo.gsb +type=crs' to
2978 // '+proj=longlat +ellps=GRS80 +nadgrids=@bar.gsb +geoidgrids=@bar.gtx
2979 // +type=crs'
2980 if (boundSrc && compoundDst) {
2981 createOperationsBoundToCompound(sourceCRS, targetCRS, context, boundSrc,
2982 compoundDst, res);
2983 return res;
2984 }
2985
2986 // reverse of previous case
2987 if (boundDst && compoundSrc) {
2988 return applyInverse(createOperations(targetCRS, sourceCRS, context));
2989 }
2990
2991 return res;
2992 }
2993
2994 // ---------------------------------------------------------------------------
2995
createOperationsFromProj4Ext(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const crs::BoundCRS * boundSrc,const crs::BoundCRS * boundDst,std::vector<CoordinateOperationNNPtr> & res)2996 void CoordinateOperationFactory::Private::createOperationsFromProj4Ext(
2997 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
2998 const crs::BoundCRS *boundSrc, const crs::BoundCRS *boundDst,
2999 std::vector<CoordinateOperationNNPtr> &res) {
3000
3001 ENTER_FUNCTION();
3002
3003 auto sourceProjExportable = dynamic_cast<const io::IPROJStringExportable *>(
3004 boundSrc ? boundSrc : sourceCRS.get());
3005 auto targetProjExportable = dynamic_cast<const io::IPROJStringExportable *>(
3006 boundDst ? boundDst : targetCRS.get());
3007 if (!sourceProjExportable) {
3008 throw InvalidOperation("Source CRS is not PROJ exportable");
3009 }
3010 if (!targetProjExportable) {
3011 throw InvalidOperation("Target CRS is not PROJ exportable");
3012 }
3013 auto projFormatter = io::PROJStringFormatter::create();
3014 projFormatter->setCRSExport(true);
3015 projFormatter->setLegacyCRSToCRSContext(true);
3016 projFormatter->startInversion();
3017 sourceProjExportable->_exportToPROJString(projFormatter.get());
3018 auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(sourceCRS.get());
3019 if (geogSrc) {
3020 auto tmpFormatter = io::PROJStringFormatter::create();
3021 geogSrc->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
3022 projFormatter->ingestPROJString(tmpFormatter->toString());
3023 }
3024
3025 projFormatter->stopInversion();
3026
3027 targetProjExportable->_exportToPROJString(projFormatter.get());
3028 auto geogDst = dynamic_cast<const crs::GeographicCRS *>(targetCRS.get());
3029 if (geogDst) {
3030 auto tmpFormatter = io::PROJStringFormatter::create();
3031 geogDst->addAngularUnitConvertAndAxisSwap(tmpFormatter.get());
3032 projFormatter->ingestPROJString(tmpFormatter->toString());
3033 }
3034
3035 const auto PROJString = projFormatter->toString();
3036 auto properties = util::PropertyMap().set(
3037 common::IdentifiedObject::NAME_KEY,
3038 buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()));
3039 res.emplace_back(SingleOperation::createPROJBased(
3040 properties, PROJString, sourceCRS, targetCRS, {}));
3041 }
3042
3043 // ---------------------------------------------------------------------------
3044
createOperationsFromDatabase(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::GeodeticCRS * geodSrc,const crs::GeodeticCRS * geodDst,const crs::GeographicCRS * geogSrc,const crs::GeographicCRS * geogDst,const crs::VerticalCRS * vertSrc,const crs::VerticalCRS * vertDst,std::vector<CoordinateOperationNNPtr> & res)3045 bool CoordinateOperationFactory::Private::createOperationsFromDatabase(
3046 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3047 Private::Context &context, const crs::GeodeticCRS *geodSrc,
3048 const crs::GeodeticCRS *geodDst, const crs::GeographicCRS *geogSrc,
3049 const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
3050 const crs::VerticalCRS *vertDst,
3051 std::vector<CoordinateOperationNNPtr> &res) {
3052
3053 ENTER_FUNCTION();
3054
3055 if (geogSrc && vertDst) {
3056 createOperationsFromDatabase(targetCRS, sourceCRS, context, geodDst,
3057 geodSrc, geogDst, geogSrc, vertDst,
3058 vertSrc, res);
3059 res = applyInverse(res);
3060 } else if (geogDst && vertSrc) {
3061 res = applyInverse(createOperationsGeogToVertFromGeoid(
3062 targetCRS, sourceCRS, vertSrc, context));
3063 if (!res.empty()) {
3064 createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context,
3065 vertSrc, geogDst, res);
3066 }
3067 }
3068
3069 if (!res.empty()) {
3070 return true;
3071 }
3072
3073 bool resFindDirectNonEmptyBeforeFiltering = false;
3074 res = findOpsInRegistryDirect(sourceCRS, targetCRS, context,
3075 resFindDirectNonEmptyBeforeFiltering);
3076
3077 // If we get at least a result with perfect accuracy, do not
3078 // bother generating synthetic transforms.
3079 if (hasPerfectAccuracyResult(res, context)) {
3080 return true;
3081 }
3082
3083 bool doFilterAndCheckPerfectOp = false;
3084
3085 bool sameGeodeticDatum = false;
3086
3087 if (vertSrc || vertDst) {
3088 if (res.empty()) {
3089 if (geogSrc &&
3090 geogSrc->coordinateSystem()->axisList().size() == 2 &&
3091 vertDst) {
3092 auto dbContext =
3093 context.context->getAuthorityFactory()->databaseContext();
3094 auto resTmp = findOpsInRegistryDirect(
3095 sourceCRS->promoteTo3D(std::string(), dbContext), targetCRS,
3096 context, resFindDirectNonEmptyBeforeFiltering);
3097 for (auto &op : resTmp) {
3098 auto newOp = op->shallowClone();
3099 setCRSs(newOp.get(), sourceCRS, targetCRS);
3100 res.emplace_back(newOp);
3101 }
3102 } else if (geogDst &&
3103 geogDst->coordinateSystem()->axisList().size() == 2 &&
3104 vertSrc) {
3105 auto dbContext =
3106 context.context->getAuthorityFactory()->databaseContext();
3107 auto resTmp = findOpsInRegistryDirect(
3108 sourceCRS, targetCRS->promoteTo3D(std::string(), dbContext),
3109 context, resFindDirectNonEmptyBeforeFiltering);
3110 for (auto &op : resTmp) {
3111 auto newOp = op->shallowClone();
3112 setCRSs(newOp.get(), sourceCRS, targetCRS);
3113 res.emplace_back(newOp);
3114 }
3115 }
3116 }
3117 if (res.empty()) {
3118 createOperationsFromDatabaseWithVertCRS(sourceCRS, targetCRS,
3119 context, geogSrc, geogDst,
3120 vertSrc, vertDst, res);
3121 }
3122 } else if (geodSrc && geodDst) {
3123
3124 const auto &authFactory = context.context->getAuthorityFactory();
3125 const auto dbContext = authFactory->databaseContext().as_nullable();
3126
3127 const auto srcDatum = geodSrc->datumNonNull(dbContext);
3128 const auto dstDatum = geodDst->datumNonNull(dbContext);
3129 sameGeodeticDatum = srcDatum->_isEquivalentTo(
3130 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT);
3131
3132 if (res.empty() && !sameGeodeticDatum &&
3133 !context.inCreateOperationsWithDatumPivotAntiRecursion) {
3134 // If we still didn't find a transformation, and that the source
3135 // and target are GeodeticCRS, then go through their underlying
3136 // datum to find potential transformations between other
3137 // GeodeticCRSs
3138 // that are made of those datum
3139 // The typical example is if transforming between two
3140 // GeographicCRS,
3141 // but transformations are only available between their
3142 // corresponding geocentric CRS.
3143 createOperationsWithDatumPivot(res, sourceCRS, targetCRS, geodSrc,
3144 geodDst, context);
3145 doFilterAndCheckPerfectOp = !res.empty();
3146 }
3147 }
3148
3149 bool foundInstantiableOp = false;
3150 // FIXME: the limitation to .size() == 1 is just for the
3151 // -s EPSG:4959+5759 -t "EPSG:4959+7839" case
3152 // finding EPSG:7860 'NZVD2016 height to Auckland 1946
3153 // height (1)', which uses the EPSG:1071 'Vertical Offset by Grid
3154 // Interpolation (NZLVD)' method which is not currently implemented by PROJ
3155 // (cannot deal with .csv files)
3156 // Initially the test was written to iterate over for all operations of a
3157 // non-empty res, but this causes failures in the test suite when no grids
3158 // are installed at all. Ideally we should tweak the test suite to be
3159 // robust to that, or skip some tests.
3160 if (res.size() == 1) {
3161 try {
3162 res.front()->exportToPROJString(
3163 io::PROJStringFormatter::create().get());
3164 foundInstantiableOp = true;
3165 } catch (const std::exception &) {
3166 }
3167 if (!foundInstantiableOp) {
3168 resFindDirectNonEmptyBeforeFiltering = false;
3169 }
3170 } else if (res.size() > 1) {
3171 foundInstantiableOp = true;
3172 }
3173
3174 // NAD27 to NAD83 has tens of results already. No need to look
3175 // for a pivot
3176 if (!sameGeodeticDatum &&
3177 (((res.empty() || !foundInstantiableOp) &&
3178 !resFindDirectNonEmptyBeforeFiltering &&
3179 context.context->getAllowUseIntermediateCRS() ==
3180 CoordinateOperationContext::IntermediateCRSUse::
3181 IF_NO_DIRECT_TRANSFORMATION) ||
3182 context.context->getAllowUseIntermediateCRS() ==
3183 CoordinateOperationContext::IntermediateCRSUse::ALWAYS ||
3184 getenv("PROJ_FORCE_SEARCH_PIVOT"))) {
3185 auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
3186 sourceCRS, targetCRS, context, false);
3187 res.insert(res.end(), resWithIntermediate.begin(),
3188 resWithIntermediate.end());
3189 doFilterAndCheckPerfectOp = !res.empty();
3190
3191 } else if (!context.inCreateOperationsWithDatumPivotAntiRecursion &&
3192 !resFindDirectNonEmptyBeforeFiltering && geodSrc && geodDst &&
3193 !sameGeodeticDatum &&
3194 context.context->getIntermediateCRS().empty() &&
3195 context.context->getAllowUseIntermediateCRS() !=
3196 CoordinateOperationContext::IntermediateCRSUse::NEVER) {
3197
3198 bool tryWithGeodeticDatumIntermediate = res.empty();
3199 if (!tryWithGeodeticDatumIntermediate) {
3200 // This is in particular for the GDA94 to WGS 84 (G1762) case
3201 // As we have a WGS 84 -> WGS 84 (G1762) null-transformation in the
3202 // PROJ authority, previous steps might have use that WGS 84
3203 // intermediate directly. They might also have generated a path
3204 // through ITRF2008, as there is a path
3205 // GDA94 (geoc.) -> ITRF2008 (geoc.) -> WGS84 (G1762) (geoc.)
3206 // But there's a better path using
3207 // GDA94 (geog.) --> GDA2020 (geog.) and
3208 // GDA2020 (geoc.) -> WGS84 (G1762) (geoc.) that requires to
3209 // explore intermediates through their datum, and not directly
3210 // trough the CRS code.
3211 // Do that only if the number of results we got through other
3212 // algorithms is small, or if all results we have go through an
3213 // operation in the PROJ authority.
3214 constexpr size_t ARBITRARY_SMALL_NUMBER = 5U;
3215 tryWithGeodeticDatumIntermediate =
3216 res.size() < ARBITRARY_SMALL_NUMBER ||
3217 hasResultSetOnlyResultsWithPROJStep(res);
3218 }
3219 if (tryWithGeodeticDatumIntermediate) {
3220 auto resWithIntermediate = findsOpsInRegistryWithIntermediate(
3221 sourceCRS, targetCRS, context, true);
3222 res.insert(res.end(), resWithIntermediate.begin(),
3223 resWithIntermediate.end());
3224 doFilterAndCheckPerfectOp = !res.empty();
3225 }
3226 }
3227
3228 if (doFilterAndCheckPerfectOp) {
3229 // If we get at least a result with perfect accuracy, do not bother
3230 // generating synthetic transforms.
3231 if (hasPerfectAccuracyResult(res, context)) {
3232 return true;
3233 }
3234 }
3235 return false;
3236 }
3237
3238 // ---------------------------------------------------------------------------
3239
3240 static std::vector<crs::CRSNNPtr>
findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr & authFactory,const datum::VerticalReferenceFrame * datum)3241 findCandidateVertCRSForDatum(const io::AuthorityFactoryPtr &authFactory,
3242 const datum::VerticalReferenceFrame *datum) {
3243 std::vector<crs::CRSNNPtr> candidates;
3244 assert(datum);
3245 const auto &ids = datum->identifiers();
3246 const auto &datumName = datum->nameStr();
3247 if (!ids.empty()) {
3248 for (const auto &id : ids) {
3249 const auto &authName = *(id->codeSpace());
3250 const auto &code = id->code();
3251 if (!authName.empty()) {
3252 auto l_candidates =
3253 authFactory->createVerticalCRSFromDatum(authName, code);
3254 for (const auto &candidate : l_candidates) {
3255 candidates.emplace_back(candidate);
3256 }
3257 }
3258 }
3259 } else if (datumName != "unknown" && datumName != "unnamed") {
3260 auto matches = authFactory->createObjectsFromName(
3261 datumName,
3262 {io::AuthorityFactory::ObjectType::VERTICAL_REFERENCE_FRAME}, false,
3263 2);
3264 if (matches.size() == 1) {
3265 const auto &match = matches.front();
3266 if (datum->_isEquivalentTo(
3267 match.get(), util::IComparable::Criterion::EQUIVALENT) &&
3268 !match->identifiers().empty()) {
3269 return findCandidateVertCRSForDatum(
3270 authFactory,
3271 dynamic_cast<const datum::VerticalReferenceFrame *>(
3272 match.get()));
3273 }
3274 }
3275 }
3276 return candidates;
3277 }
3278
3279 // ---------------------------------------------------------------------------
3280
3281 std::vector<CoordinateOperationNNPtr>
createOperationsGeogToVertFromGeoid(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const crs::VerticalCRS * vertDst,Private::Context & context)3282 CoordinateOperationFactory::Private::createOperationsGeogToVertFromGeoid(
3283 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3284 const crs::VerticalCRS *vertDst, Private::Context &context) {
3285
3286 ENTER_FUNCTION();
3287
3288 const auto useTransf = [&targetCRS, &context,
3289 vertDst](const CoordinateOperationNNPtr &op) {
3290 const auto targetOp =
3291 dynamic_cast<const crs::VerticalCRS *>(op->targetCRS().get());
3292 assert(targetOp);
3293 if (targetOp->_isEquivalentTo(
3294 vertDst, util::IComparable::Criterion::EQUIVALENT)) {
3295 return op;
3296 }
3297 std::vector<CoordinateOperationNNPtr> tmp;
3298 createOperationsVertToVert(NN_NO_CHECK(op->targetCRS()), targetCRS,
3299 context, targetOp, vertDst, tmp);
3300 assert(!tmp.empty());
3301 auto ret = ConcatenatedOperation::createComputeMetadata(
3302 {op, tmp.front()}, disallowEmptyIntersection);
3303 return ret;
3304 };
3305
3306 const auto getProjGeoidTransformation = [&sourceCRS, &targetCRS, &vertDst,
3307 &context](
3308 const CoordinateOperationNNPtr &model,
3309 const std::string &projFilename) {
3310
3311 const auto getNameVertCRSMetre = [](const std::string &name) {
3312 if (name.empty())
3313 return std::string("unnamed");
3314 auto ret(name);
3315 bool haveOriginalUnit = false;
3316 if (name.back() == ')') {
3317 const auto pos = ret.rfind(" (");
3318 if (pos != std::string::npos) {
3319 haveOriginalUnit = true;
3320 ret = ret.substr(0, pos);
3321 }
3322 }
3323 const auto pos = ret.rfind(" depth");
3324 if (pos != std::string::npos) {
3325 ret = ret.substr(0, pos) + " height";
3326 }
3327 if (!haveOriginalUnit) {
3328 ret += " (metre)";
3329 }
3330 return ret;
3331 };
3332
3333 const auto &axis = vertDst->coordinateSystem()->axisList()[0];
3334 const auto geogSrcCRS =
3335 dynamic_cast<crs::GeographicCRS *>(model->interpolationCRS().get())
3336 ? NN_NO_CHECK(model->interpolationCRS())
3337 : sourceCRS;
3338 const auto vertCRSMetre =
3339 axis->unit() == common::UnitOfMeasure::METRE &&
3340 axis->direction() == cs::AxisDirection::UP
3341 ? targetCRS
3342 : util::nn_static_pointer_cast<crs::CRS>(
3343 crs::VerticalCRS::create(
3344 util::PropertyMap().set(
3345 common::IdentifiedObject::NAME_KEY,
3346 getNameVertCRSMetre(targetCRS->nameStr())),
3347 vertDst->datum(), vertDst->datumEnsemble(),
3348 cs::VerticalCS::createGravityRelatedHeight(
3349 common::UnitOfMeasure::METRE)));
3350 const auto properties = util::PropertyMap().set(
3351 common::IdentifiedObject::NAME_KEY,
3352 buildOpName("Transformation", vertCRSMetre, geogSrcCRS));
3353
3354 // Try to find a representative value for the accuracy of this grid
3355 // from the registered transformations.
3356 std::vector<metadata::PositionalAccuracyNNPtr> accuracies;
3357 const auto &modelAccuracies = model->coordinateOperationAccuracies();
3358 if (modelAccuracies.empty()) {
3359 const auto &authFactory = context.context->getAuthorityFactory();
3360 if (authFactory) {
3361 const auto transformationsForGrid =
3362 io::DatabaseContext::getTransformationsForGridName(
3363 authFactory->databaseContext(), projFilename);
3364 double accuracy = -1;
3365 for (const auto &transf : transformationsForGrid) {
3366 accuracy = std::max(accuracy, getAccuracy(transf));
3367 }
3368 if (accuracy >= 0) {
3369 accuracies.emplace_back(
3370 metadata::PositionalAccuracy::create(
3371 toString(accuracy)));
3372 }
3373 }
3374 }
3375
3376 return Transformation::createGravityRelatedHeightToGeographic3D(
3377 properties, vertCRSMetre, geogSrcCRS, nullptr, projFilename,
3378 !modelAccuracies.empty() ? modelAccuracies : accuracies);
3379 };
3380
3381 std::vector<CoordinateOperationNNPtr> res;
3382 const auto &authFactory = context.context->getAuthorityFactory();
3383 if (authFactory) {
3384 const auto &models = vertDst->geoidModel();
3385 for (const auto &model : models) {
3386 const auto &modelName = model->nameStr();
3387 const auto transformations =
3388 starts_with(modelName, "PROJ ")
3389 ? std::vector<
3390 CoordinateOperationNNPtr>{getProjGeoidTransformation(
3391 model, modelName.substr(strlen("PROJ ")))}
3392 : authFactory->getTransformationsForGeoid(
3393 modelName,
3394 context.context->getUsePROJAlternativeGridNames());
3395 for (const auto &transf : transformations) {
3396 if (dynamic_cast<crs::GeographicCRS *>(
3397 transf->sourceCRS().get()) &&
3398 dynamic_cast<crs::VerticalCRS *>(
3399 transf->targetCRS().get())) {
3400 res.push_back(useTransf(transf));
3401 } else if (dynamic_cast<crs::GeographicCRS *>(
3402 transf->targetCRS().get()) &&
3403 dynamic_cast<crs::VerticalCRS *>(
3404 transf->sourceCRS().get())) {
3405 res.push_back(useTransf(transf->inverse()));
3406 }
3407 }
3408 }
3409 }
3410
3411 return res;
3412 }
3413
3414 // ---------------------------------------------------------------------------
3415
3416 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
createOperationsGeogToVertWithIntermediateVert(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const crs::VerticalCRS * vertDst,Private::Context & context)3417 createOperationsGeogToVertWithIntermediateVert(
3418 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3419 const crs::VerticalCRS *vertDst, Private::Context &context) {
3420
3421 ENTER_FUNCTION();
3422
3423 std::vector<CoordinateOperationNNPtr> res;
3424
3425 struct AntiRecursionGuard {
3426 Context &context;
3427
3428 explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
3429 assert(!context.inCreateOperationsGeogToVertWithIntermediateVert);
3430 context.inCreateOperationsGeogToVertWithIntermediateVert = true;
3431 }
3432
3433 ~AntiRecursionGuard() {
3434 context.inCreateOperationsGeogToVertWithIntermediateVert = false;
3435 }
3436 };
3437 AntiRecursionGuard guard(context);
3438 const auto &authFactory = context.context->getAuthorityFactory();
3439 const auto dbContext = authFactory->databaseContext().as_nullable();
3440
3441 auto candidatesVert = findCandidateVertCRSForDatum(
3442 authFactory, vertDst->datumNonNull(dbContext).get());
3443 for (const auto &candidateVert : candidatesVert) {
3444 auto resTmp = createOperations(sourceCRS, candidateVert, context);
3445 if (!resTmp.empty()) {
3446 const auto opsSecond =
3447 createOperations(candidateVert, targetCRS, context);
3448 if (!opsSecond.empty()) {
3449 // The transformation from candidateVert to targetCRS should
3450 // be just a unit change typically, so take only the first one,
3451 // which is likely/hopefully the only one.
3452 for (const auto &opFirst : resTmp) {
3453 if (hasIdentifiers(opFirst)) {
3454 if (candidateVert->_isEquivalentTo(
3455 targetCRS.get(),
3456 util::IComparable::Criterion::EQUIVALENT)) {
3457 res.emplace_back(opFirst);
3458 } else {
3459 res.emplace_back(
3460 ConcatenatedOperation::createComputeMetadata(
3461 {opFirst, opsSecond.front()},
3462 disallowEmptyIntersection));
3463 }
3464 }
3465 }
3466 if (!res.empty())
3467 break;
3468 }
3469 }
3470 }
3471
3472 return res;
3473 }
3474
3475 // ---------------------------------------------------------------------------
3476
3477 std::vector<CoordinateOperationNNPtr> CoordinateOperationFactory::Private::
createOperationsGeogToVertWithAlternativeGeog(const crs::CRSNNPtr &,const crs::CRSNNPtr & targetCRS,Private::Context & context)3478 createOperationsGeogToVertWithAlternativeGeog(
3479 const crs::CRSNNPtr & /*sourceCRS*/, // geographic CRS
3480 const crs::CRSNNPtr &targetCRS, // vertical CRS
3481 Private::Context &context) {
3482
3483 ENTER_FUNCTION();
3484
3485 std::vector<CoordinateOperationNNPtr> res;
3486
3487 struct AntiRecursionGuard {
3488 Context &context;
3489
3490 explicit AntiRecursionGuard(Context &contextIn) : context(contextIn) {
3491 assert(!context.inCreateOperationsGeogToVertWithAlternativeGeog);
3492 context.inCreateOperationsGeogToVertWithAlternativeGeog = true;
3493 }
3494
3495 ~AntiRecursionGuard() {
3496 context.inCreateOperationsGeogToVertWithAlternativeGeog = false;
3497 }
3498 };
3499 AntiRecursionGuard guard(context);
3500
3501 // Generally EPSG has operations from GeogCrs to VertCRS
3502 auto ops = findOpsInRegistryDirectTo(targetCRS, context);
3503
3504 for (const auto &op : ops) {
3505 const auto tmpCRS = op->sourceCRS();
3506 if (tmpCRS && dynamic_cast<const crs::GeographicCRS *>(tmpCRS.get())) {
3507 res.emplace_back(op);
3508 }
3509 }
3510
3511 return res;
3512 }
3513
3514 // ---------------------------------------------------------------------------
3515
3516 void CoordinateOperationFactory::Private::
createOperationsFromDatabaseWithVertCRS(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::GeographicCRS * geogSrc,const crs::GeographicCRS * geogDst,const crs::VerticalCRS * vertSrc,const crs::VerticalCRS * vertDst,std::vector<CoordinateOperationNNPtr> & res)3517 createOperationsFromDatabaseWithVertCRS(
3518 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3519 Private::Context &context, const crs::GeographicCRS *geogSrc,
3520 const crs::GeographicCRS *geogDst, const crs::VerticalCRS *vertSrc,
3521 const crs::VerticalCRS *vertDst,
3522 std::vector<CoordinateOperationNNPtr> &res) {
3523
3524 // Typically to transform from "NAVD88 height (ftUS)" to a geog CRS
3525 // by using transformations of "NAVD88 height" (metre) to that geog CRS
3526 if (res.empty() &&
3527 !context.inCreateOperationsGeogToVertWithIntermediateVert && geogSrc &&
3528 vertDst) {
3529 res = createOperationsGeogToVertWithIntermediateVert(
3530 sourceCRS, targetCRS, vertDst, context);
3531 } else if (res.empty() &&
3532 !context.inCreateOperationsGeogToVertWithIntermediateVert &&
3533 geogDst && vertSrc) {
3534 res = applyInverse(createOperationsGeogToVertWithIntermediateVert(
3535 targetCRS, sourceCRS, vertSrc, context));
3536 }
3537
3538 // NAD83 only exists in 2D version in EPSG, so if it has been
3539 // promoted to 3D, when researching a vertical to geog
3540 // transformation, try to down cast to 2D.
3541 const auto geog3DToVertTryThroughGeog2D = [&res, &context](
3542 const crs::GeographicCRS *geogSrcIn, const crs::VerticalCRS *vertDstIn,
3543 const crs::CRSNNPtr &targetCRSIn) {
3544 if (res.empty() && geogSrcIn && vertDstIn &&
3545 geogSrcIn->coordinateSystem()->axisList().size() == 3) {
3546 const auto &authFactory = context.context->getAuthorityFactory();
3547 const auto dbContext =
3548 authFactory ? authFactory->databaseContext().as_nullable()
3549 : nullptr;
3550 const auto candidatesSrcGeod(findCandidateGeodCRSForDatum(
3551 authFactory, geogSrcIn,
3552 geogSrcIn->datumNonNull(dbContext).get()));
3553 for (const auto &candidate : candidatesSrcGeod) {
3554 auto geogCandidate =
3555 util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
3556 candidate);
3557 if (geogCandidate &&
3558 geogCandidate->coordinateSystem()->axisList().size() == 2) {
3559 bool ignored;
3560 res =
3561 findOpsInRegistryDirect(NN_NO_CHECK(geogCandidate),
3562 targetCRSIn, context, ignored);
3563 break;
3564 }
3565 }
3566 return true;
3567 }
3568 return false;
3569 };
3570
3571 if (geog3DToVertTryThroughGeog2D(geogSrc, vertDst, targetCRS)) {
3572 // do nothing
3573 } else if (geog3DToVertTryThroughGeog2D(geogDst, vertSrc, sourceCRS)) {
3574 res = applyInverse(res);
3575 }
3576
3577 // There's no direct transformation from NAVD88 height to WGS84,
3578 // so try to research all transformations from NAVD88 to another
3579 // intermediate GeographicCRS.
3580 if (res.empty() &&
3581 !context.inCreateOperationsGeogToVertWithAlternativeGeog && geogSrc &&
3582 vertDst) {
3583 res = createOperationsGeogToVertWithAlternativeGeog(sourceCRS,
3584 targetCRS, context);
3585 } else if (res.empty() &&
3586 !context.inCreateOperationsGeogToVertWithAlternativeGeog &&
3587 geogDst && vertSrc) {
3588 res = applyInverse(createOperationsGeogToVertWithAlternativeGeog(
3589 targetCRS, sourceCRS, context));
3590 }
3591 }
3592
3593 // ---------------------------------------------------------------------------
3594
createOperationsGeodToGeod(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::GeodeticCRS * geodSrc,const crs::GeodeticCRS * geodDst,std::vector<CoordinateOperationNNPtr> & res)3595 void CoordinateOperationFactory::Private::createOperationsGeodToGeod(
3596 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3597 Private::Context &context, const crs::GeodeticCRS *geodSrc,
3598 const crs::GeodeticCRS *geodDst,
3599 std::vector<CoordinateOperationNNPtr> &res) {
3600
3601 ENTER_FUNCTION();
3602
3603 if (geodSrc->ellipsoid()->celestialBody() !=
3604 geodDst->ellipsoid()->celestialBody()) {
3605 throw util::UnsupportedOperationException(
3606 "Source and target ellipsoid do not belong to the same "
3607 "celestial body");
3608 }
3609
3610 auto geogSrc = dynamic_cast<const crs::GeographicCRS *>(geodSrc);
3611 auto geogDst = dynamic_cast<const crs::GeographicCRS *>(geodDst);
3612
3613 if (geogSrc && geogDst) {
3614 createOperationsGeogToGeog(res, sourceCRS, targetCRS, context, geogSrc,
3615 geogDst);
3616 return;
3617 }
3618
3619 const bool isSrcGeocentric = geodSrc->isGeocentric();
3620 const bool isSrcGeographic = geogSrc != nullptr;
3621 const bool isTargetGeocentric = geodDst->isGeocentric();
3622 const bool isTargetGeographic = geogDst != nullptr;
3623
3624 const auto IsSameDatum = [&context, &geodSrc, &geodDst]() {
3625 const auto &authFactory = context.context->getAuthorityFactory();
3626 const auto dbContext =
3627 authFactory ? authFactory->databaseContext().as_nullable()
3628 : nullptr;
3629
3630 return geodSrc->datumNonNull(dbContext)->_isEquivalentTo(
3631 geodDst->datumNonNull(dbContext).get(),
3632 util::IComparable::Criterion::EQUIVALENT);
3633 };
3634
3635 if (((isSrcGeocentric && isTargetGeographic) ||
3636 (isSrcGeographic && isTargetGeocentric))) {
3637
3638 // Same datum ?
3639 if (IsSameDatum()) {
3640 res.emplace_back(
3641 Conversion::createGeographicGeocentric(sourceCRS, targetCRS));
3642 } else if (isSrcGeocentric && geogDst) {
3643 std::string interm_crs_name(geogDst->nameStr());
3644 interm_crs_name += " (geocentric)";
3645 auto interm_crs =
3646 util::nn_static_pointer_cast<crs::CRS>(crs::GeodeticCRS::create(
3647 addDomains(util::PropertyMap().set(
3648 common::IdentifiedObject::NAME_KEY,
3649 interm_crs_name),
3650 geogDst),
3651 geogDst->datum(), geogDst->datumEnsemble(),
3652 NN_CHECK_ASSERT(
3653 util::nn_dynamic_pointer_cast<cs::CartesianCS>(
3654 geodSrc->coordinateSystem()))));
3655 auto opFirst =
3656 createBallparkGeocentricTranslation(sourceCRS, interm_crs);
3657 auto opSecond =
3658 Conversion::createGeographicGeocentric(interm_crs, targetCRS);
3659 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
3660 {opFirst, opSecond}, disallowEmptyIntersection));
3661 } else {
3662 // Apply previous case in reverse way
3663 std::vector<CoordinateOperationNNPtr> resTmp;
3664 createOperationsGeodToGeod(targetCRS, sourceCRS, context, geodDst,
3665 geodSrc, resTmp);
3666 assert(resTmp.size() == 1);
3667 res.emplace_back(resTmp.front()->inverse());
3668 }
3669
3670 return;
3671 }
3672
3673 if (isSrcGeocentric && isTargetGeocentric) {
3674 if (sourceCRS->_isEquivalentTo(
3675 targetCRS.get(), util::IComparable::Criterion::EQUIVALENT) ||
3676 IsSameDatum()) {
3677 std::string name(NULL_GEOCENTRIC_TRANSLATION);
3678 name += " from ";
3679 name += sourceCRS->nameStr();
3680 name += " to ";
3681 name += targetCRS->nameStr();
3682 res.emplace_back(Transformation::createGeocentricTranslations(
3683 util::PropertyMap()
3684 .set(common::IdentifiedObject::NAME_KEY, name)
3685 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
3686 metadata::Extent::WORLD),
3687 sourceCRS, targetCRS, 0.0, 0.0, 0.0,
3688 {metadata::PositionalAccuracy::create("0")}));
3689 } else {
3690 res.emplace_back(
3691 createBallparkGeocentricTranslation(sourceCRS, targetCRS));
3692 }
3693 return;
3694 }
3695
3696 // Transformation between two geodetic systems of unknown type
3697 // This should normally not be triggered with "standard" CRS
3698 res.emplace_back(createGeodToGeodPROJBased(sourceCRS, targetCRS));
3699 }
3700
3701 // ---------------------------------------------------------------------------
3702
createOperationsDerivedTo(const crs::CRSNNPtr &,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::DerivedCRS * derivedSrc,std::vector<CoordinateOperationNNPtr> & res)3703 void CoordinateOperationFactory::Private::createOperationsDerivedTo(
3704 const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
3705 Private::Context &context, const crs::DerivedCRS *derivedSrc,
3706 std::vector<CoordinateOperationNNPtr> &res) {
3707
3708 ENTER_FUNCTION();
3709
3710 auto opFirst = derivedSrc->derivingConversion()->inverse();
3711 // Small optimization if the targetCRS is the baseCRS of the source
3712 // derivedCRS.
3713 if (derivedSrc->baseCRS()->_isEquivalentTo(
3714 targetCRS.get(), util::IComparable::Criterion::EQUIVALENT)) {
3715 res.emplace_back(opFirst);
3716 return;
3717 }
3718 auto opsSecond =
3719 createOperations(derivedSrc->baseCRS(), targetCRS, context);
3720 for (const auto &opSecond : opsSecond) {
3721 try {
3722 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
3723 {opFirst, opSecond}, disallowEmptyIntersection));
3724 } catch (const InvalidOperationEmptyIntersection &) {
3725 }
3726 }
3727 }
3728
3729 // ---------------------------------------------------------------------------
3730
createOperationsBoundToGeog(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::BoundCRS * boundSrc,const crs::GeographicCRS * geogDst,std::vector<CoordinateOperationNNPtr> & res)3731 void CoordinateOperationFactory::Private::createOperationsBoundToGeog(
3732 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
3733 Private::Context &context, const crs::BoundCRS *boundSrc,
3734 const crs::GeographicCRS *geogDst,
3735 std::vector<CoordinateOperationNNPtr> &res) {
3736
3737 ENTER_FUNCTION();
3738
3739 const auto &hubSrc = boundSrc->hubCRS();
3740 auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
3741 auto geogCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractGeographicCRS();
3742 {
3743 // If geogCRSOfBaseOfBoundSrc is a DerivedGeographicCRS, use its base
3744 // instead (if it is a GeographicCRS)
3745 auto derivedGeogCRS =
3746 std::dynamic_pointer_cast<crs::DerivedGeographicCRS>(
3747 geogCRSOfBaseOfBoundSrc);
3748 if (derivedGeogCRS) {
3749 auto baseCRS = std::dynamic_pointer_cast<crs::GeographicCRS>(
3750 derivedGeogCRS->baseCRS().as_nullable());
3751 if (baseCRS) {
3752 geogCRSOfBaseOfBoundSrc = baseCRS;
3753 }
3754 }
3755 }
3756
3757 const auto &authFactory = context.context->getAuthorityFactory();
3758 const auto dbContext =
3759 authFactory ? authFactory->databaseContext().as_nullable() : nullptr;
3760
3761 const auto geogDstDatum = geogDst->datumNonNull(dbContext);
3762
3763 // If the underlying datum of the source is the same as the target, do
3764 // not consider the boundCRS at all, but just its base
3765 if (geogCRSOfBaseOfBoundSrc) {
3766 auto geogCRSOfBaseOfBoundSrcDatum =
3767 geogCRSOfBaseOfBoundSrc->datumNonNull(dbContext);
3768 if (geogCRSOfBaseOfBoundSrcDatum->_isEquivalentTo(
3769 geogDstDatum.get(), util::IComparable::Criterion::EQUIVALENT)) {
3770 res = createOperations(boundSrc->baseCRS(), targetCRS, context);
3771 return;
3772 }
3773 }
3774
3775 bool triedBoundCrsToGeogCRSSameAsHubCRS = false;
3776 // Is it: boundCRS to a geogCRS that is the same as the hubCRS ?
3777 if (hubSrcGeog && geogCRSOfBaseOfBoundSrc &&
3778 (hubSrcGeog->_isEquivalentTo(
3779 geogDst, util::IComparable::Criterion::EQUIVALENT) ||
3780 hubSrcGeog->is2DPartOf3D(NN_NO_CHECK(geogDst), dbContext))) {
3781 triedBoundCrsToGeogCRSSameAsHubCRS = true;
3782
3783 CoordinateOperationPtr opIntermediate;
3784 if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
3785 boundSrc->transformation()->sourceCRS().get(),
3786 util::IComparable::Criterion::EQUIVALENT)) {
3787 auto opsIntermediate = createOperations(
3788 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc),
3789 boundSrc->transformation()->sourceCRS(), context);
3790 assert(!opsIntermediate.empty());
3791 opIntermediate = opsIntermediate.front();
3792 }
3793
3794 if (boundSrc->baseCRS() == geogCRSOfBaseOfBoundSrc) {
3795 if (opIntermediate) {
3796 try {
3797 res.emplace_back(
3798 ConcatenatedOperation::createComputeMetadata(
3799 {NN_NO_CHECK(opIntermediate),
3800 boundSrc->transformation()},
3801 disallowEmptyIntersection));
3802 } catch (const InvalidOperationEmptyIntersection &) {
3803 }
3804 } else {
3805 // Optimization to avoid creating a useless concatenated
3806 // operation
3807 res.emplace_back(boundSrc->transformation());
3808 }
3809 return;
3810 }
3811 auto opsFirst = createOperations(
3812 boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
3813 if (!opsFirst.empty()) {
3814 for (const auto &opFirst : opsFirst) {
3815 try {
3816 std::vector<CoordinateOperationNNPtr> subops;
3817 subops.emplace_back(opFirst);
3818 if (opIntermediate) {
3819 subops.emplace_back(NN_NO_CHECK(opIntermediate));
3820 }
3821 subops.emplace_back(boundSrc->transformation());
3822 res.emplace_back(
3823 ConcatenatedOperation::createComputeMetadata(
3824 subops, disallowEmptyIntersection));
3825 } catch (const InvalidOperationEmptyIntersection &) {
3826 }
3827 }
3828 if (!res.empty()) {
3829 return;
3830 }
3831 }
3832 // If the datum are equivalent, this is also fine
3833 } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog &&
3834 hubSrcGeog->datumNonNull(dbContext)->_isEquivalentTo(
3835 geogDstDatum.get(),
3836 util::IComparable::Criterion::EQUIVALENT)) {
3837 auto opsFirst = createOperations(
3838 boundSrc->baseCRS(), NN_NO_CHECK(geogCRSOfBaseOfBoundSrc), context);
3839 auto opsLast = createOperations(hubSrc, targetCRS, context);
3840 if (!opsFirst.empty() && !opsLast.empty()) {
3841 CoordinateOperationPtr opIntermediate;
3842 if (!geogCRSOfBaseOfBoundSrc->_isEquivalentTo(
3843 boundSrc->transformation()->sourceCRS().get(),
3844 util::IComparable::Criterion::EQUIVALENT)) {
3845 auto opsIntermediate = createOperations(
3846 NN_NO_CHECK(geogCRSOfBaseOfBoundSrc),
3847 boundSrc->transformation()->sourceCRS(), context);
3848 assert(!opsIntermediate.empty());
3849 opIntermediate = opsIntermediate.front();
3850 }
3851 for (const auto &opFirst : opsFirst) {
3852 for (const auto &opLast : opsLast) {
3853 try {
3854 std::vector<CoordinateOperationNNPtr> subops;
3855 subops.emplace_back(opFirst);
3856 if (opIntermediate) {
3857 subops.emplace_back(NN_NO_CHECK(opIntermediate));
3858 }
3859 subops.emplace_back(boundSrc->transformation());
3860 subops.emplace_back(opLast);
3861 res.emplace_back(
3862 ConcatenatedOperation::createComputeMetadata(
3863 subops, disallowEmptyIntersection));
3864 } catch (const InvalidOperationEmptyIntersection &) {
3865 }
3866 }
3867 }
3868 if (!res.empty()) {
3869 return;
3870 }
3871 }
3872 // Consider WGS 84 and NAD83 as equivalent in that context if the
3873 // geogCRSOfBaseOfBoundSrc ellipsoid is Clarke66 (for NAD27)
3874 // Case of "+proj=latlong +ellps=clrk66
3875 // +nadgrids=ntv1_can.dat,conus"
3876 // to "+proj=latlong +datum=NAD83"
3877 } else if (geogCRSOfBaseOfBoundSrc && hubSrcGeog &&
3878 geogCRSOfBaseOfBoundSrc->ellipsoid()->_isEquivalentTo(
3879 datum::Ellipsoid::CLARKE_1866.get(),
3880 util::IComparable::Criterion::EQUIVALENT) &&
3881 hubSrcGeog->datumNonNull(dbContext)->_isEquivalentTo(
3882 datum::GeodeticReferenceFrame::EPSG_6326.get(),
3883 util::IComparable::Criterion::EQUIVALENT) &&
3884 geogDstDatum->_isEquivalentTo(
3885 datum::GeodeticReferenceFrame::EPSG_6269.get(),
3886 util::IComparable::Criterion::EQUIVALENT)) {
3887 auto nnGeogCRSOfBaseOfBoundSrc = NN_NO_CHECK(geogCRSOfBaseOfBoundSrc);
3888 if (boundSrc->baseCRS()->_isEquivalentTo(
3889 nnGeogCRSOfBaseOfBoundSrc.get(),
3890 util::IComparable::Criterion::EQUIVALENT)) {
3891 auto transf = boundSrc->transformation()->shallowClone();
3892 transf->setProperties(util::PropertyMap().set(
3893 common::IdentifiedObject::NAME_KEY,
3894 buildTransfName(boundSrc->baseCRS()->nameStr(),
3895 targetCRS->nameStr())));
3896 transf->setCRSs(boundSrc->baseCRS(), targetCRS, nullptr);
3897 res.emplace_back(transf);
3898 return;
3899 } else {
3900 auto opsFirst = createOperations(
3901 boundSrc->baseCRS(), nnGeogCRSOfBaseOfBoundSrc, context);
3902 auto transf = boundSrc->transformation()->shallowClone();
3903 transf->setProperties(util::PropertyMap().set(
3904 common::IdentifiedObject::NAME_KEY,
3905 buildTransfName(nnGeogCRSOfBaseOfBoundSrc->nameStr(),
3906 targetCRS->nameStr())));
3907 transf->setCRSs(nnGeogCRSOfBaseOfBoundSrc, targetCRS, nullptr);
3908 if (!opsFirst.empty()) {
3909 for (const auto &opFirst : opsFirst) {
3910 try {
3911 res.emplace_back(
3912 ConcatenatedOperation::createComputeMetadata(
3913 {opFirst, transf}, disallowEmptyIntersection));
3914 } catch (const InvalidOperationEmptyIntersection &) {
3915 }
3916 }
3917 if (!res.empty()) {
3918 return;
3919 }
3920 }
3921 }
3922 }
3923
3924 if (hubSrcGeog &&
3925 hubSrcGeog->_isEquivalentTo(geogDst,
3926 util::IComparable::Criterion::EQUIVALENT) &&
3927 dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get())) {
3928 auto transfSrc = boundSrc->transformation()->sourceCRS();
3929 if (dynamic_cast<const crs::VerticalCRS *>(transfSrc.get()) &&
3930 !boundSrc->baseCRS()->_isEquivalentTo(
3931 transfSrc.get(), util::IComparable::Criterion::EQUIVALENT)) {
3932 auto opsFirst =
3933 createOperations(boundSrc->baseCRS(), transfSrc, context);
3934 for (const auto &opFirst : opsFirst) {
3935 try {
3936 res.emplace_back(
3937 ConcatenatedOperation::createComputeMetadata(
3938 {opFirst, boundSrc->transformation()},
3939 disallowEmptyIntersection));
3940 } catch (const InvalidOperationEmptyIntersection &) {
3941 }
3942 }
3943 return;
3944 }
3945
3946 res.emplace_back(boundSrc->transformation());
3947 return;
3948 }
3949
3950 if (!triedBoundCrsToGeogCRSSameAsHubCRS && hubSrcGeog &&
3951 geogCRSOfBaseOfBoundSrc) {
3952 // This one should go to the above 'Is it: boundCRS to a geogCRS
3953 // that is the same as the hubCRS ?' case
3954 auto opsFirst = createOperations(sourceCRS, hubSrc, context);
3955 auto opsLast = createOperations(hubSrc, targetCRS, context);
3956 if (!opsFirst.empty() && !opsLast.empty()) {
3957 for (const auto &opFirst : opsFirst) {
3958 for (const auto &opLast : opsLast) {
3959 // Exclude artificial transformations from the hub
3960 // to the target CRS, if it is the only one.
3961 if (opsLast.size() > 1 ||
3962 !opLast->hasBallparkTransformation()) {
3963 try {
3964 res.emplace_back(
3965 ConcatenatedOperation::createComputeMetadata(
3966 {opFirst, opLast},
3967 disallowEmptyIntersection));
3968 } catch (const InvalidOperationEmptyIntersection &) {
3969 }
3970 } else {
3971 // std::cerr << "excluded " << opLast->nameStr() <<
3972 // std::endl;
3973 }
3974 }
3975 }
3976 if (!res.empty()) {
3977 return;
3978 }
3979 }
3980 }
3981
3982 auto vertCRSOfBaseOfBoundSrc =
3983 dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
3984 if (vertCRSOfBaseOfBoundSrc && hubSrcGeog) {
3985 auto opsFirst = createOperations(sourceCRS, hubSrc, context);
3986 if (context.skipHorizontalTransformation) {
3987 if (!opsFirst.empty()) {
3988 const auto &hubAxisList =
3989 hubSrcGeog->coordinateSystem()->axisList();
3990 const auto &targetAxisList =
3991 geogDst->coordinateSystem()->axisList();
3992 if (hubAxisList.size() == 3 && targetAxisList.size() == 3 &&
3993 !hubAxisList[2]->_isEquivalentTo(
3994 targetAxisList[2].get(),
3995 util::IComparable::Criterion::EQUIVALENT)) {
3996
3997 const auto &srcAxis = hubAxisList[2];
3998 const double convSrc = srcAxis->unit().conversionToSI();
3999 const auto &dstAxis = targetAxisList[2];
4000 const double convDst = dstAxis->unit().conversionToSI();
4001 const bool srcIsUp =
4002 srcAxis->direction() == cs::AxisDirection::UP;
4003 const bool srcIsDown =
4004 srcAxis->direction() == cs::AxisDirection::DOWN;
4005 const bool dstIsUp =
4006 dstAxis->direction() == cs::AxisDirection::UP;
4007 const bool dstIsDown =
4008 dstAxis->direction() == cs::AxisDirection::DOWN;
4009 const bool heightDepthReversal =
4010 ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));
4011
4012 const double factor = convSrc / convDst;
4013 auto conv = Conversion::createChangeVerticalUnit(
4014 util::PropertyMap().set(
4015 common::IdentifiedObject::NAME_KEY,
4016 "Change of vertical unit"),
4017 common::Scale(heightDepthReversal ? -factor : factor));
4018
4019 conv->setCRSs(
4020 hubSrc,
4021 hubSrc->demoteTo2D(std::string(), dbContext)
4022 ->promoteTo3D(std::string(), dbContext, dstAxis),
4023 nullptr);
4024
4025 for (const auto &op : opsFirst) {
4026 try {
4027 res.emplace_back(
4028 ConcatenatedOperation::createComputeMetadata(
4029 {op, conv}, disallowEmptyIntersection));
4030 } catch (const InvalidOperationEmptyIntersection &) {
4031 }
4032 }
4033 } else {
4034 res = opsFirst;
4035 }
4036 }
4037 return;
4038 } else {
4039 auto opsSecond = createOperations(hubSrc, targetCRS, context);
4040 if (!opsFirst.empty() && !opsSecond.empty()) {
4041 for (const auto &opFirst : opsFirst) {
4042 for (const auto &opLast : opsSecond) {
4043 // Exclude artificial transformations from the hub
4044 // to the target CRS
4045 if (!opLast->hasBallparkTransformation()) {
4046 try {
4047 res.emplace_back(
4048 ConcatenatedOperation::
4049 createComputeMetadata(
4050 {opFirst, opLast},
4051 disallowEmptyIntersection));
4052 } catch (
4053 const InvalidOperationEmptyIntersection &) {
4054 }
4055 } else {
4056 // std::cerr << "excluded " << opLast->nameStr() <<
4057 // std::endl;
4058 }
4059 }
4060 }
4061 if (!res.empty()) {
4062 return;
4063 }
4064 }
4065 }
4066 }
4067
4068 res = createOperations(boundSrc->baseCRS(), targetCRS, context);
4069 }
4070
4071 // ---------------------------------------------------------------------------
4072
createOperationsBoundToVert(const crs::CRSNNPtr &,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::BoundCRS * boundSrc,const crs::VerticalCRS * vertDst,std::vector<CoordinateOperationNNPtr> & res)4073 void CoordinateOperationFactory::Private::createOperationsBoundToVert(
4074 const crs::CRSNNPtr & /*sourceCRS*/, const crs::CRSNNPtr &targetCRS,
4075 Private::Context &context, const crs::BoundCRS *boundSrc,
4076 const crs::VerticalCRS *vertDst,
4077 std::vector<CoordinateOperationNNPtr> &res) {
4078
4079 ENTER_FUNCTION();
4080
4081 auto baseSrcVert =
4082 dynamic_cast<const crs::VerticalCRS *>(boundSrc->baseCRS().get());
4083 const auto &hubSrc = boundSrc->hubCRS();
4084 auto hubSrcVert = dynamic_cast<const crs::VerticalCRS *>(hubSrc.get());
4085 if (baseSrcVert && hubSrcVert &&
4086 vertDst->_isEquivalentTo(hubSrcVert,
4087 util::IComparable::Criterion::EQUIVALENT)) {
4088 res.emplace_back(boundSrc->transformation());
4089 return;
4090 }
4091
4092 res = createOperations(boundSrc->baseCRS(), targetCRS, context);
4093 }
4094
4095 // ---------------------------------------------------------------------------
4096
createOperationsVertToVert(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::VerticalCRS * vertSrc,const crs::VerticalCRS * vertDst,std::vector<CoordinateOperationNNPtr> & res)4097 void CoordinateOperationFactory::Private::createOperationsVertToVert(
4098 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4099 Private::Context &context, const crs::VerticalCRS *vertSrc,
4100 const crs::VerticalCRS *vertDst,
4101 std::vector<CoordinateOperationNNPtr> &res) {
4102
4103 ENTER_FUNCTION();
4104
4105 const auto &authFactory = context.context->getAuthorityFactory();
4106 const auto dbContext =
4107 authFactory ? authFactory->databaseContext().as_nullable() : nullptr;
4108
4109 const auto srcDatum = vertSrc->datumNonNull(dbContext);
4110 const auto dstDatum = vertDst->datumNonNull(dbContext);
4111 const bool equivalentVDatum = srcDatum->_isEquivalentTo(
4112 dstDatum.get(), util::IComparable::Criterion::EQUIVALENT);
4113
4114 const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0];
4115 const double convSrc = srcAxis->unit().conversionToSI();
4116 const auto &dstAxis = vertDst->coordinateSystem()->axisList()[0];
4117 const double convDst = dstAxis->unit().conversionToSI();
4118 const bool srcIsUp = srcAxis->direction() == cs::AxisDirection::UP;
4119 const bool srcIsDown = srcAxis->direction() == cs::AxisDirection::DOWN;
4120 const bool dstIsUp = dstAxis->direction() == cs::AxisDirection::UP;
4121 const bool dstIsDown = dstAxis->direction() == cs::AxisDirection::DOWN;
4122 const bool heightDepthReversal =
4123 ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));
4124
4125 const double factor = convSrc / convDst;
4126 auto name = buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr());
4127 if (!equivalentVDatum) {
4128 name += BALLPARK_VERTICAL_TRANSFORMATION;
4129 auto conv = Transformation::createChangeVerticalUnit(
4130 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
4131 sourceCRS, targetCRS,
4132 // In case of a height depth reversal, we should probably have
4133 // 2 steps instead of putting a negative factor...
4134 common::Scale(heightDepthReversal ? -factor : factor), {});
4135 conv->setHasBallparkTransformation(true);
4136 res.push_back(conv);
4137 } else if (convSrc != convDst || !heightDepthReversal) {
4138 auto conv = Conversion::createChangeVerticalUnit(
4139 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name),
4140 // In case of a height depth reversal, we should probably have
4141 // 2 steps instead of putting a negative factor...
4142 common::Scale(heightDepthReversal ? -factor : factor));
4143 conv->setCRSs(sourceCRS, targetCRS, nullptr);
4144 res.push_back(conv);
4145 } else {
4146 auto conv = Conversion::createHeightDepthReversal(
4147 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY, name));
4148 conv->setCRSs(sourceCRS, targetCRS, nullptr);
4149 res.push_back(conv);
4150 }
4151 }
4152
4153 // ---------------------------------------------------------------------------
4154
createOperationsVertToGeog(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::VerticalCRS * vertSrc,const crs::GeographicCRS * geogDst,std::vector<CoordinateOperationNNPtr> & res)4155 void CoordinateOperationFactory::Private::createOperationsVertToGeog(
4156 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4157 Private::Context &context, const crs::VerticalCRS *vertSrc,
4158 const crs::GeographicCRS *geogDst,
4159 std::vector<CoordinateOperationNNPtr> &res) {
4160
4161 ENTER_FUNCTION();
4162
4163 if (vertSrc->identifiers().empty()) {
4164 const auto &vertSrcName = vertSrc->nameStr();
4165 const auto &authFactory = context.context->getAuthorityFactory();
4166 if (authFactory != nullptr && vertSrcName != "unnamed" &&
4167 vertSrcName != "unknown") {
4168 auto matches = authFactory->createObjectsFromName(
4169 vertSrcName, {io::AuthorityFactory::ObjectType::VERTICAL_CRS},
4170 false, 2);
4171 if (matches.size() == 1) {
4172 const auto &match = matches.front();
4173 if (vertSrc->_isEquivalentTo(
4174 match.get(),
4175 util::IComparable::Criterion::EQUIVALENT) &&
4176 !match->identifiers().empty()) {
4177 auto resTmp = createOperations(
4178 NN_NO_CHECK(
4179 util::nn_dynamic_pointer_cast<crs::VerticalCRS>(
4180 match)),
4181 targetCRS, context);
4182 res.insert(res.end(), resTmp.begin(), resTmp.end());
4183 return;
4184 }
4185 }
4186 }
4187 }
4188
4189 createOperationsVertToGeogBallpark(sourceCRS, targetCRS, context, vertSrc,
4190 geogDst, res);
4191 }
4192
4193 // ---------------------------------------------------------------------------
4194
createOperationsVertToGeogBallpark(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context &,const crs::VerticalCRS * vertSrc,const crs::GeographicCRS * geogDst,std::vector<CoordinateOperationNNPtr> & res)4195 void CoordinateOperationFactory::Private::createOperationsVertToGeogBallpark(
4196 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4197 Private::Context &, const crs::VerticalCRS *vertSrc,
4198 const crs::GeographicCRS *geogDst,
4199 std::vector<CoordinateOperationNNPtr> &res) {
4200
4201 ENTER_FUNCTION();
4202
4203 const auto &srcAxis = vertSrc->coordinateSystem()->axisList()[0];
4204 const double convSrc = srcAxis->unit().conversionToSI();
4205 double convDst = 1.0;
4206 const auto &geogAxis = geogDst->coordinateSystem()->axisList();
4207 bool dstIsUp = true;
4208 bool dstIsDown = true;
4209 if (geogAxis.size() == 3) {
4210 const auto &dstAxis = geogAxis[2];
4211 convDst = dstAxis->unit().conversionToSI();
4212 dstIsUp = dstAxis->direction() == cs::AxisDirection::UP;
4213 dstIsDown = dstAxis->direction() == cs::AxisDirection::DOWN;
4214 }
4215 const bool srcIsUp = srcAxis->direction() == cs::AxisDirection::UP;
4216 const bool srcIsDown = srcAxis->direction() == cs::AxisDirection::DOWN;
4217 const bool heightDepthReversal =
4218 ((srcIsUp && dstIsDown) || (srcIsDown && dstIsUp));
4219
4220 const double factor = convSrc / convDst;
4221
4222 const auto &sourceCRSExtent = getExtent(sourceCRS);
4223 const auto &targetCRSExtent = getExtent(targetCRS);
4224 const bool sameExtent =
4225 sourceCRSExtent && targetCRSExtent &&
4226 sourceCRSExtent->_isEquivalentTo(
4227 targetCRSExtent.get(), util::IComparable::Criterion::EQUIVALENT);
4228
4229 util::PropertyMap map;
4230 map.set(common::IdentifiedObject::NAME_KEY,
4231 buildTransfName(sourceCRS->nameStr(), targetCRS->nameStr()) +
4232 BALLPARK_VERTICAL_TRANSFORMATION_NO_ELLIPSOID_VERT_HEIGHT)
4233 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
4234 sameExtent ? NN_NO_CHECK(sourceCRSExtent)
4235 : metadata::Extent::WORLD);
4236
4237 auto conv = Transformation::createChangeVerticalUnit(
4238 map, sourceCRS, targetCRS,
4239 common::Scale(heightDepthReversal ? -factor : factor), {});
4240 conv->setHasBallparkTransformation(true);
4241 res.push_back(conv);
4242 }
4243
4244 // ---------------------------------------------------------------------------
4245
createOperationsBoundToBound(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::BoundCRS * boundSrc,const crs::BoundCRS * boundDst,std::vector<CoordinateOperationNNPtr> & res)4246 void CoordinateOperationFactory::Private::createOperationsBoundToBound(
4247 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4248 Private::Context &context, const crs::BoundCRS *boundSrc,
4249 const crs::BoundCRS *boundDst, std::vector<CoordinateOperationNNPtr> &res) {
4250
4251 ENTER_FUNCTION();
4252
4253 // BoundCRS to BoundCRS of horizontal CRS using the same (geographic) hub
4254 const auto &hubSrc = boundSrc->hubCRS();
4255 auto hubSrcGeog = dynamic_cast<const crs::GeographicCRS *>(hubSrc.get());
4256 const auto &hubDst = boundDst->hubCRS();
4257 auto hubDstGeog = dynamic_cast<const crs::GeographicCRS *>(hubDst.get());
4258 if (hubSrcGeog && hubDstGeog &&
4259 hubSrcGeog->_isEquivalentTo(hubDstGeog,
4260 util::IComparable::Criterion::EQUIVALENT)) {
4261 auto opsFirst = createOperations(sourceCRS, hubSrc, context);
4262 auto opsLast = createOperations(hubSrc, targetCRS, context);
4263 for (const auto &opFirst : opsFirst) {
4264 for (const auto &opLast : opsLast) {
4265 try {
4266 std::vector<CoordinateOperationNNPtr> ops;
4267 ops.push_back(opFirst);
4268 ops.push_back(opLast);
4269 res.emplace_back(
4270 ConcatenatedOperation::createComputeMetadata(
4271 ops, disallowEmptyIntersection));
4272 } catch (const InvalidOperationEmptyIntersection &) {
4273 }
4274 }
4275 }
4276 if (!res.empty()) {
4277 return;
4278 }
4279 }
4280
4281 // BoundCRS to BoundCRS of vertical CRS using the same vertical datum
4282 // ==> ignore the bound transformation
4283 auto baseOfBoundSrcAsVertCRS =
4284 dynamic_cast<crs::VerticalCRS *>(boundSrc->baseCRS().get());
4285 auto baseOfBoundDstAsVertCRS =
4286 dynamic_cast<crs::VerticalCRS *>(boundDst->baseCRS().get());
4287 if (baseOfBoundSrcAsVertCRS && baseOfBoundDstAsVertCRS) {
4288
4289 const auto &authFactory = context.context->getAuthorityFactory();
4290 const auto dbContext =
4291 authFactory ? authFactory->databaseContext().as_nullable()
4292 : nullptr;
4293
4294 const auto datumSrc = baseOfBoundSrcAsVertCRS->datumNonNull(dbContext);
4295 const auto datumDst = baseOfBoundDstAsVertCRS->datumNonNull(dbContext);
4296 if (datumSrc->nameStr() == datumDst->nameStr() &&
4297 (datumSrc->nameStr() != "unknown" ||
4298 boundSrc->transformation()->_isEquivalentTo(
4299 boundDst->transformation().get(),
4300 util::IComparable::Criterion::EQUIVALENT))) {
4301 res = createOperations(boundSrc->baseCRS(), boundDst->baseCRS(),
4302 context);
4303 return;
4304 }
4305 }
4306
4307 // BoundCRS to BoundCRS of vertical CRS
4308 auto vertCRSOfBaseOfBoundSrc = boundSrc->baseCRS()->extractVerticalCRS();
4309 auto vertCRSOfBaseOfBoundDst = boundDst->baseCRS()->extractVerticalCRS();
4310 if (hubSrcGeog && hubDstGeog &&
4311 hubSrcGeog->_isEquivalentTo(hubDstGeog,
4312 util::IComparable::Criterion::EQUIVALENT) &&
4313 vertCRSOfBaseOfBoundSrc && vertCRSOfBaseOfBoundDst) {
4314 auto opsFirst = createOperations(sourceCRS, hubSrc, context);
4315 auto opsLast = createOperations(hubSrc, targetCRS, context);
4316 if (!opsFirst.empty() && !opsLast.empty()) {
4317 for (const auto &opFirst : opsFirst) {
4318 for (const auto &opLast : opsLast) {
4319 try {
4320 res.emplace_back(
4321 ConcatenatedOperation::createComputeMetadata(
4322 {opFirst, opLast}, disallowEmptyIntersection));
4323 } catch (const InvalidOperationEmptyIntersection &) {
4324 }
4325 }
4326 }
4327 if (!res.empty()) {
4328 return;
4329 }
4330 }
4331 }
4332
4333 res = createOperations(boundSrc->baseCRS(), boundDst->baseCRS(), context);
4334 }
4335
4336 // ---------------------------------------------------------------------------
4337
4338 static std::vector<CoordinateOperationNNPtr>
getOps(const CoordinateOperationNNPtr & op)4339 getOps(const CoordinateOperationNNPtr &op) {
4340 auto concatenated = dynamic_cast<const ConcatenatedOperation *>(op.get());
4341 if (concatenated)
4342 return concatenated->operations();
4343 return {op};
4344 }
4345
4346 // ---------------------------------------------------------------------------
4347
useDifferentTransformationsForSameSourceTarget(const CoordinateOperationNNPtr & opA,const CoordinateOperationNNPtr & opB)4348 static bool useDifferentTransformationsForSameSourceTarget(
4349 const CoordinateOperationNNPtr &opA, const CoordinateOperationNNPtr &opB) {
4350 auto subOpsA = getOps(opA);
4351 auto subOpsB = getOps(opB);
4352 for (const auto &subOpA : subOpsA) {
4353 if (!dynamic_cast<const Transformation *>(subOpA.get()))
4354 continue;
4355 if (subOpA->sourceCRS()->nameStr() == "unknown" ||
4356 subOpA->targetCRS()->nameStr() == "unknown")
4357 continue;
4358 for (const auto &subOpB : subOpsB) {
4359 if (!dynamic_cast<const Transformation *>(subOpB.get()))
4360 continue;
4361 if (subOpB->sourceCRS()->nameStr() == "unknown" ||
4362 subOpB->targetCRS()->nameStr() == "unknown")
4363 continue;
4364
4365 if (subOpA->sourceCRS()->nameStr() ==
4366 subOpB->sourceCRS()->nameStr() &&
4367 subOpA->targetCRS()->nameStr() ==
4368 subOpB->targetCRS()->nameStr()) {
4369 if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
4370 starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
4371 continue;
4372 }
4373
4374 if (!subOpA->isEquivalentTo(subOpB.get())) {
4375 return true;
4376 }
4377 } else if (subOpA->sourceCRS()->nameStr() ==
4378 subOpB->targetCRS()->nameStr() &&
4379 subOpA->targetCRS()->nameStr() ==
4380 subOpB->sourceCRS()->nameStr()) {
4381 if (starts_with(subOpA->nameStr(), NULL_GEOGRAPHIC_OFFSET) &&
4382 starts_with(subOpB->nameStr(), NULL_GEOGRAPHIC_OFFSET)) {
4383 continue;
4384 }
4385
4386 if (!subOpA->isEquivalentTo(subOpB->inverse().get())) {
4387 return true;
4388 }
4389 }
4390 }
4391 }
4392 return false;
4393 }
4394
4395 // ---------------------------------------------------------------------------
4396
4397 static crs::GeographicCRSPtr
getInterpolationGeogCRS(const CoordinateOperationNNPtr & verticalTransform,const io::DatabaseContextPtr & dbContext)4398 getInterpolationGeogCRS(const CoordinateOperationNNPtr &verticalTransform,
4399 const io::DatabaseContextPtr &dbContext) {
4400 crs::GeographicCRSPtr interpolationGeogCRS;
4401 auto transformationVerticalTransform =
4402 dynamic_cast<const Transformation *>(verticalTransform.get());
4403 if (transformationVerticalTransform == nullptr) {
4404 const auto concat = dynamic_cast<const ConcatenatedOperation *>(
4405 verticalTransform.get());
4406 if (concat) {
4407 const auto &steps = concat->operations();
4408 // Is this change of unit and/or height depth reversal +
4409 // transformation ?
4410 for (const auto &step : steps) {
4411 const auto transf =
4412 dynamic_cast<const Transformation *>(step.get());
4413 if (transf) {
4414 // Only support a single Transformation in the steps
4415 if (transformationVerticalTransform != nullptr) {
4416 transformationVerticalTransform = nullptr;
4417 break;
4418 }
4419 transformationVerticalTransform = transf;
4420 }
4421 }
4422 }
4423 }
4424 if (transformationVerticalTransform &&
4425 !transformationVerticalTransform->hasBallparkTransformation()) {
4426 auto interpTransformCRS =
4427 transformationVerticalTransform->interpolationCRS();
4428 if (interpTransformCRS) {
4429 interpolationGeogCRS =
4430 std::dynamic_pointer_cast<crs::GeographicCRS>(
4431 interpTransformCRS);
4432 } else {
4433 // If no explicit interpolation CRS, then
4434 // this will be the geographic CRS of the
4435 // vertical to geog transformation
4436 interpolationGeogCRS =
4437 std::dynamic_pointer_cast<crs::GeographicCRS>(
4438 transformationVerticalTransform->targetCRS().as_nullable());
4439 }
4440 }
4441
4442 if (interpolationGeogCRS) {
4443 if (interpolationGeogCRS->coordinateSystem()->axisList().size() == 3) {
4444 // We need to force the interpolation CRS, which
4445 // will
4446 // frequently be 3D, to 2D to avoid transformations
4447 // between source CRS and interpolation CRS to have
4448 // 3D terms.
4449 interpolationGeogCRS =
4450 interpolationGeogCRS->demoteTo2D(std::string(), dbContext)
4451 .as_nullable();
4452 }
4453 }
4454
4455 return interpolationGeogCRS;
4456 }
4457
4458 // ---------------------------------------------------------------------------
4459
createOperationsCompoundToGeog(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::CompoundCRS * compoundSrc,const crs::GeographicCRS * geogDst,std::vector<CoordinateOperationNNPtr> & res)4460 void CoordinateOperationFactory::Private::createOperationsCompoundToGeog(
4461 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4462 Private::Context &context, const crs::CompoundCRS *compoundSrc,
4463 const crs::GeographicCRS *geogDst,
4464 std::vector<CoordinateOperationNNPtr> &res) {
4465
4466 ENTER_FUNCTION();
4467
4468 const auto &authFactory = context.context->getAuthorityFactory();
4469 const auto &componentsSrc = compoundSrc->componentReferenceSystems();
4470 if (!componentsSrc.empty()) {
4471
4472 if (componentsSrc.size() == 2) {
4473 auto derivedHSrc =
4474 dynamic_cast<const crs::DerivedCRS *>(componentsSrc[0].get());
4475 if (derivedHSrc) {
4476 std::vector<crs::CRSNNPtr> intermComponents{
4477 derivedHSrc->baseCRS(), componentsSrc[1]};
4478 auto properties = util::PropertyMap().set(
4479 common::IdentifiedObject::NAME_KEY,
4480 intermComponents[0]->nameStr() + " + " +
4481 intermComponents[1]->nameStr());
4482 auto intermCompound =
4483 crs::CompoundCRS::create(properties, intermComponents);
4484 auto opsFirst =
4485 createOperations(sourceCRS, intermCompound, context);
4486 assert(!opsFirst.empty());
4487 auto opsLast =
4488 createOperations(intermCompound, targetCRS, context);
4489 for (const auto &opLast : opsLast) {
4490 try {
4491 res.emplace_back(
4492 ConcatenatedOperation::createComputeMetadata(
4493 {opsFirst.front(), opLast},
4494 disallowEmptyIntersection));
4495 } catch (const std::exception &) {
4496 }
4497 }
4498 return;
4499 }
4500 }
4501
4502 std::vector<CoordinateOperationNNPtr> horizTransforms;
4503 auto srcGeogCRS = componentsSrc[0]->extractGeographicCRS();
4504 if (srcGeogCRS) {
4505 horizTransforms =
4506 createOperations(componentsSrc[0], targetCRS, context);
4507 }
4508 std::vector<CoordinateOperationNNPtr> verticalTransforms;
4509
4510 const auto dbContext =
4511 authFactory ? authFactory->databaseContext().as_nullable()
4512 : nullptr;
4513 if (componentsSrc.size() >= 2 &&
4514 componentsSrc[1]->extractVerticalCRS()) {
4515
4516 struct SetSkipHorizontalTransform {
4517 Context &context;
4518
4519 explicit SetSkipHorizontalTransform(Context &contextIn)
4520 : context(contextIn) {
4521 assert(!context.skipHorizontalTransformation);
4522 context.skipHorizontalTransformation = true;
4523 }
4524
4525 ~SetSkipHorizontalTransform() {
4526 context.skipHorizontalTransformation = false;
4527 }
4528 };
4529 SetSkipHorizontalTransform setSkipHorizontalTransform(context);
4530
4531 verticalTransforms = createOperations(
4532 componentsSrc[1],
4533 targetCRS->promoteTo3D(std::string(), dbContext), context);
4534 bool foundRegisteredTransformWithAllGridsAvailable = false;
4535 const auto gridAvailabilityUse =
4536 context.context->getGridAvailabilityUse();
4537 const bool ignoreMissingGrids =
4538 gridAvailabilityUse ==
4539 CoordinateOperationContext::GridAvailabilityUse::
4540 IGNORE_GRID_AVAILABILITY;
4541 for (const auto &op : verticalTransforms) {
4542 if (hasIdentifiers(op) && dbContext) {
4543 bool missingGrid = false;
4544 if (!ignoreMissingGrids) {
4545 const auto gridsNeeded = op->gridsNeeded(
4546 dbContext,
4547 gridAvailabilityUse ==
4548 CoordinateOperationContext::
4549 GridAvailabilityUse::KNOWN_AVAILABLE);
4550 for (const auto &gridDesc : gridsNeeded) {
4551 if (!gridDesc.available) {
4552 missingGrid = true;
4553 break;
4554 }
4555 }
4556 }
4557 if (!missingGrid) {
4558 foundRegisteredTransformWithAllGridsAvailable = true;
4559 break;
4560 }
4561 }
4562 }
4563 if (!foundRegisteredTransformWithAllGridsAvailable && srcGeogCRS &&
4564 !srcGeogCRS->_isEquivalentTo(
4565 geogDst, util::IComparable::Criterion::EQUIVALENT) &&
4566 !srcGeogCRS->is2DPartOf3D(NN_NO_CHECK(geogDst), dbContext)) {
4567 auto verticalTransformsTmp = createOperations(
4568 componentsSrc[1],
4569 NN_NO_CHECK(srcGeogCRS)
4570 ->promoteTo3D(std::string(), dbContext),
4571 context);
4572 bool foundRegisteredTransform = false;
4573 foundRegisteredTransformWithAllGridsAvailable = false;
4574 for (const auto &op : verticalTransformsTmp) {
4575 if (hasIdentifiers(op) && dbContext) {
4576 bool missingGrid = false;
4577 if (!ignoreMissingGrids) {
4578 const auto gridsNeeded = op->gridsNeeded(
4579 dbContext,
4580 gridAvailabilityUse ==
4581 CoordinateOperationContext::
4582 GridAvailabilityUse::KNOWN_AVAILABLE);
4583 for (const auto &gridDesc : gridsNeeded) {
4584 if (!gridDesc.available) {
4585 missingGrid = true;
4586 break;
4587 }
4588 }
4589 }
4590 foundRegisteredTransform = true;
4591 if (!missingGrid) {
4592 foundRegisteredTransformWithAllGridsAvailable =
4593 true;
4594 break;
4595 }
4596 }
4597 }
4598 if (foundRegisteredTransformWithAllGridsAvailable) {
4599 verticalTransforms = verticalTransformsTmp;
4600 } else if (foundRegisteredTransform) {
4601 verticalTransforms.insert(verticalTransforms.end(),
4602 verticalTransformsTmp.begin(),
4603 verticalTransformsTmp.end());
4604 }
4605 }
4606 }
4607
4608 if (horizTransforms.empty() || verticalTransforms.empty()) {
4609 res = horizTransforms;
4610 return;
4611 }
4612
4613 typedef std::pair<std::vector<CoordinateOperationNNPtr>,
4614 std::vector<CoordinateOperationNNPtr>>
4615 PairOfTransforms;
4616 std::map<std::string, PairOfTransforms>
4617 cacheHorizToInterpAndInterpToTarget;
4618
4619 for (const auto &verticalTransform : verticalTransforms) {
4620 #ifdef TRACE_CREATE_OPERATIONS
4621 ENTER_BLOCK("Considering vertical transform " +
4622 objectAsStr(verticalTransform.get()));
4623 #endif
4624 crs::GeographicCRSPtr interpolationGeogCRS =
4625 getInterpolationGeogCRS(verticalTransform, dbContext);
4626 if (interpolationGeogCRS) {
4627 #ifdef TRACE_CREATE_OPERATIONS
4628 logTrace("Using " + objectAsStr(interpolationGeogCRS.get()) +
4629 " as interpolation CRS");
4630 #endif
4631 std::vector<CoordinateOperationNNPtr> srcToInterpOps;
4632 std::vector<CoordinateOperationNNPtr> interpToTargetOps;
4633
4634 std::string key;
4635 const auto &ids = interpolationGeogCRS->identifiers();
4636 if (!ids.empty()) {
4637 key =
4638 (*ids.front()->codeSpace()) + ':' + ids.front()->code();
4639 }
4640
4641 const auto computeOpsToInterp =
4642 [&srcToInterpOps, &interpToTargetOps, &componentsSrc,
4643 &interpolationGeogCRS, &targetCRS, &dbContext,
4644 &context]() {
4645 srcToInterpOps = createOperations(
4646 componentsSrc[0], NN_NO_CHECK(interpolationGeogCRS),
4647 context);
4648 auto target2D =
4649 targetCRS->demoteTo2D(std::string(), dbContext);
4650 if (!componentsSrc[0]->isEquivalentTo(
4651 target2D.get(),
4652 util::IComparable::Criterion::EQUIVALENT)) {
4653 // We do the transformation from the
4654 // interpolationCRS
4655 // to the target one in 3D (see #2225)
4656 // But we don't do that between sourceCRS and
4657 // interpolationCRS, as this would mess with an
4658 // orthometric elevation.
4659 auto interp3D = interpolationGeogCRS->promoteTo3D(
4660 std::string(), dbContext);
4661 interpToTargetOps =
4662 createOperations(interp3D, targetCRS, context);
4663 }
4664 };
4665
4666 if (!key.empty()) {
4667 auto iter = cacheHorizToInterpAndInterpToTarget.find(key);
4668 if (iter == cacheHorizToInterpAndInterpToTarget.end()) {
4669 #ifdef TRACE_CREATE_OPERATIONS
4670 ENTER_BLOCK("looking for horizontal transformation "
4671 "from source to interpCRS and interpCRS to "
4672 "target");
4673 #endif
4674 computeOpsToInterp();
4675 cacheHorizToInterpAndInterpToTarget[key] =
4676 PairOfTransforms(srcToInterpOps, interpToTargetOps);
4677 } else {
4678 srcToInterpOps = iter->second.first;
4679 interpToTargetOps = iter->second.second;
4680 }
4681 } else {
4682 #ifdef TRACE_CREATE_OPERATIONS
4683 ENTER_BLOCK("looking for horizontal transformation "
4684 "from source to interpCRS and interpCRS to "
4685 "target");
4686 #endif
4687 computeOpsToInterp();
4688 }
4689
4690 #ifdef TRACE_CREATE_OPERATIONS
4691 ENTER_BLOCK("creating HorizVerticalHorizPROJBased operations");
4692 #endif
4693 for (const auto &srcToInterp : srcToInterpOps) {
4694 if (interpToTargetOps.empty()) {
4695 try {
4696 auto op = createHorizVerticalHorizPROJBased(
4697 sourceCRS, targetCRS, srcToInterp,
4698 verticalTransform, srcToInterp->inverse(),
4699 interpolationGeogCRS, true);
4700 res.emplace_back(op);
4701 } catch (const std::exception &) {
4702 }
4703 } else {
4704 for (const auto &interpToTarget : interpToTargetOps) {
4705
4706 if (useDifferentTransformationsForSameSourceTarget(
4707 srcToInterp, interpToTarget)) {
4708 continue;
4709 }
4710
4711 try {
4712 auto op = createHorizVerticalHorizPROJBased(
4713 sourceCRS, targetCRS, srcToInterp,
4714 verticalTransform, interpToTarget,
4715 interpolationGeogCRS, true);
4716 res.emplace_back(op);
4717 } catch (const std::exception &) {
4718 }
4719 }
4720 }
4721 }
4722 } else {
4723 // This case is probably only correct if
4724 // verticalTransform and horizTransform are independent
4725 // and in particular that verticalTransform does not
4726 // involve a grid, because of the rather arbitrary order
4727 // horizontal then vertical applied
4728 for (const auto &horizTransform : horizTransforms) {
4729 try {
4730 auto op = createHorizVerticalPROJBased(
4731 sourceCRS, targetCRS, horizTransform,
4732 verticalTransform, disallowEmptyIntersection);
4733 res.emplace_back(op);
4734 } catch (const std::exception &) {
4735 }
4736 }
4737 }
4738 }
4739 }
4740 }
4741
4742 // ---------------------------------------------------------------------------
4743
createOperationsToGeod(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::GeodeticCRS * geodDst,std::vector<CoordinateOperationNNPtr> & res)4744 void CoordinateOperationFactory::Private::createOperationsToGeod(
4745 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4746 Private::Context &context, const crs::GeodeticCRS *geodDst,
4747 std::vector<CoordinateOperationNNPtr> &res) {
4748
4749 auto cs = cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
4750 common::UnitOfMeasure::DEGREE, common::UnitOfMeasure::METRE);
4751 auto intermGeog3DCRS =
4752 util::nn_static_pointer_cast<crs::CRS>(crs::GeographicCRS::create(
4753 util::PropertyMap()
4754 .set(common::IdentifiedObject::NAME_KEY, geodDst->nameStr())
4755 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
4756 metadata::Extent::WORLD),
4757 geodDst->datum(), geodDst->datumEnsemble(), cs));
4758 auto sourceToGeog3DOps =
4759 createOperations(sourceCRS, intermGeog3DCRS, context);
4760 auto geog3DToTargetOps =
4761 createOperations(intermGeog3DCRS, targetCRS, context);
4762 if (!geog3DToTargetOps.empty()) {
4763 for (const auto &op : sourceToGeog3DOps) {
4764 auto newOp = op->shallowClone();
4765 setCRSs(newOp.get(), sourceCRS, intermGeog3DCRS);
4766 try {
4767 res.emplace_back(ConcatenatedOperation::createComputeMetadata(
4768 {newOp, geog3DToTargetOps.front()},
4769 disallowEmptyIntersection));
4770 } catch (const InvalidOperationEmptyIntersection &) {
4771 }
4772 }
4773 }
4774 }
4775
4776 // ---------------------------------------------------------------------------
4777
createOperationsCompoundToCompound(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::CompoundCRS * compoundSrc,const crs::CompoundCRS * compoundDst,std::vector<CoordinateOperationNNPtr> & res)4778 void CoordinateOperationFactory::Private::createOperationsCompoundToCompound(
4779 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4780 Private::Context &context, const crs::CompoundCRS *compoundSrc,
4781 const crs::CompoundCRS *compoundDst,
4782 std::vector<CoordinateOperationNNPtr> &res) {
4783
4784 const auto &componentsSrc = compoundSrc->componentReferenceSystems();
4785 const auto &componentsDst = compoundDst->componentReferenceSystems();
4786 if (componentsSrc.empty() || componentsSrc.size() != componentsDst.size()) {
4787 return;
4788 }
4789 const auto srcGeog = componentsSrc[0]->extractGeographicCRS();
4790 const auto dstGeog = componentsDst[0]->extractGeographicCRS();
4791 if (srcGeog == nullptr || dstGeog == nullptr) {
4792 return;
4793 }
4794
4795 std::vector<CoordinateOperationNNPtr> verticalTransforms;
4796 if (componentsSrc.size() >= 2 && componentsSrc[1]->extractVerticalCRS() &&
4797 componentsDst[1]->extractVerticalCRS()) {
4798 if (!componentsSrc[1]->_isEquivalentTo(componentsDst[1].get())) {
4799 verticalTransforms =
4800 createOperations(componentsSrc[1], componentsDst[1], context);
4801 }
4802 }
4803
4804 // If we didn't find a non-ballpark transformation between
4805 // the 2 vertical CRS, then try through intermediate geographic CRS
4806 // For example
4807 // WGS 84 + EGM96 --> ETRS89 + Belfast height where
4808 // there is a geoid model for EGM96 referenced to WGS 84
4809 // and a geoid model for Belfast height referenced to ETRS89
4810 if (verticalTransforms.size() == 1 &&
4811 verticalTransforms.front()->hasBallparkTransformation()) {
4812 auto dbContext =
4813 context.context->getAuthorityFactory()->databaseContext();
4814 const auto intermGeogSrc =
4815 srcGeog->promoteTo3D(std::string(), dbContext);
4816 const bool intermGeogSrcIsSameAsIntermGeogDst =
4817 srcGeog->_isEquivalentTo(dstGeog.get());
4818 const auto intermGeogDst =
4819 intermGeogSrcIsSameAsIntermGeogDst
4820 ? intermGeogSrc
4821 : dstGeog->promoteTo3D(std::string(), dbContext);
4822 const auto opsSrcToGeog =
4823 createOperations(sourceCRS, intermGeogSrc, context);
4824 const auto opsGeogToTarget =
4825 createOperations(intermGeogDst, targetCRS, context);
4826 const bool hasNonTrivalSrcTransf =
4827 !opsSrcToGeog.empty() &&
4828 !opsSrcToGeog.front()->hasBallparkTransformation();
4829 const bool hasNonTrivialTargetTransf =
4830 !opsGeogToTarget.empty() &&
4831 !opsGeogToTarget.front()->hasBallparkTransformation();
4832 if (hasNonTrivalSrcTransf && hasNonTrivialTargetTransf) {
4833 const auto opsGeogSrcToGeogDst =
4834 createOperations(intermGeogSrc, intermGeogDst, context);
4835 for (const auto &op1 : opsSrcToGeog) {
4836 if (op1->hasBallparkTransformation()) {
4837 // std::cerr << "excluded " << op1->nameStr() << std::endl;
4838 continue;
4839 }
4840 for (const auto &op2 : opsGeogSrcToGeogDst) {
4841 for (const auto &op3 : opsGeogToTarget) {
4842 if (op3->hasBallparkTransformation()) {
4843 // std::cerr << "excluded " << op3->nameStr() <<
4844 // std::endl;
4845 continue;
4846 }
4847 try {
4848 res.emplace_back(
4849 ConcatenatedOperation::createComputeMetadata(
4850 intermGeogSrcIsSameAsIntermGeogDst
4851 ? std::vector<
4852 CoordinateOperationNNPtr>{op1,
4853 op3}
4854 : std::vector<
4855 CoordinateOperationNNPtr>{op1,
4856 op2,
4857 op3},
4858 disallowEmptyIntersection));
4859 } catch (const std::exception &) {
4860 }
4861 }
4862 }
4863 }
4864 }
4865 if (!res.empty()) {
4866 return;
4867 }
4868 }
4869
4870 for (const auto &verticalTransform : verticalTransforms) {
4871 auto interpolationGeogCRS = NN_NO_CHECK(srcGeog);
4872 auto interpTransformCRS = verticalTransform->interpolationCRS();
4873 if (interpTransformCRS) {
4874 auto nn_interpTransformCRS = NN_NO_CHECK(interpTransformCRS);
4875 if (dynamic_cast<const crs::GeographicCRS *>(
4876 nn_interpTransformCRS.get())) {
4877 interpolationGeogCRS = NN_NO_CHECK(
4878 util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
4879 nn_interpTransformCRS));
4880 }
4881 } else {
4882 auto compSrc0BoundCrs =
4883 dynamic_cast<crs::BoundCRS *>(componentsSrc[0].get());
4884 auto compDst0BoundCrs =
4885 dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
4886 if (compSrc0BoundCrs && compDst0BoundCrs &&
4887 dynamic_cast<crs::GeographicCRS *>(
4888 compSrc0BoundCrs->hubCRS().get()) &&
4889 compSrc0BoundCrs->hubCRS()->_isEquivalentTo(
4890 compDst0BoundCrs->hubCRS().get())) {
4891 interpolationGeogCRS = NN_NO_CHECK(
4892 util::nn_dynamic_pointer_cast<crs::GeographicCRS>(
4893 compSrc0BoundCrs->hubCRS()));
4894 }
4895 }
4896 auto opSrcCRSToGeogCRS =
4897 createOperations(componentsSrc[0], interpolationGeogCRS, context);
4898 auto opGeogCRStoDstCRS =
4899 createOperations(interpolationGeogCRS, componentsDst[0], context);
4900 for (const auto &opSrc : opSrcCRSToGeogCRS) {
4901 for (const auto &opDst : opGeogCRStoDstCRS) {
4902
4903 try {
4904 auto op = createHorizVerticalHorizPROJBased(
4905 sourceCRS, targetCRS, opSrc, verticalTransform, opDst,
4906 interpolationGeogCRS, true);
4907 res.emplace_back(op);
4908 } catch (const InvalidOperationEmptyIntersection &) {
4909 } catch (const io::FormattingException &) {
4910 }
4911 }
4912 }
4913 }
4914
4915 if (verticalTransforms.empty()) {
4916 auto resTmp =
4917 createOperations(componentsSrc[0], componentsDst[0], context);
4918 for (const auto &op : resTmp) {
4919 auto opClone = op->shallowClone();
4920 setCRSs(opClone.get(), sourceCRS, targetCRS);
4921 res.emplace_back(opClone);
4922 }
4923 }
4924 }
4925
4926 // ---------------------------------------------------------------------------
4927
createOperationsBoundToCompound(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,Private::Context & context,const crs::BoundCRS * boundSrc,const crs::CompoundCRS * compoundDst,std::vector<CoordinateOperationNNPtr> & res)4928 void CoordinateOperationFactory::Private::createOperationsBoundToCompound(
4929 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
4930 Private::Context &context, const crs::BoundCRS *boundSrc,
4931 const crs::CompoundCRS *compoundDst,
4932 std::vector<CoordinateOperationNNPtr> &res) {
4933
4934 const auto &authFactory = context.context->getAuthorityFactory();
4935 const auto dbContext =
4936 authFactory ? authFactory->databaseContext().as_nullable() : nullptr;
4937
4938 const auto &componentsDst = compoundDst->componentReferenceSystems();
4939 if (!componentsDst.empty()) {
4940 auto compDst0BoundCrs =
4941 dynamic_cast<crs::BoundCRS *>(componentsDst[0].get());
4942 if (compDst0BoundCrs) {
4943 auto boundSrcHubAsGeogCRS =
4944 dynamic_cast<crs::GeographicCRS *>(boundSrc->hubCRS().get());
4945 auto compDst0BoundCrsHubAsGeogCRS =
4946 dynamic_cast<crs::GeographicCRS *>(
4947 compDst0BoundCrs->hubCRS().get());
4948 if (boundSrcHubAsGeogCRS && compDst0BoundCrsHubAsGeogCRS) {
4949 const auto boundSrcHubAsGeogCRSDatum =
4950 boundSrcHubAsGeogCRS->datumNonNull(dbContext);
4951 const auto compDst0BoundCrsHubAsGeogCRSDatum =
4952 compDst0BoundCrsHubAsGeogCRS->datumNonNull(dbContext);
4953 if (boundSrcHubAsGeogCRSDatum->_isEquivalentTo(
4954 compDst0BoundCrsHubAsGeogCRSDatum.get())) {
4955 auto cs = cs::EllipsoidalCS::
4956 createLatitudeLongitudeEllipsoidalHeight(
4957 common::UnitOfMeasure::DEGREE,
4958 common::UnitOfMeasure::METRE);
4959 auto intermGeog3DCRS = util::nn_static_pointer_cast<
4960 crs::CRS>(crs::GeographicCRS::create(
4961 util::PropertyMap()
4962 .set(common::IdentifiedObject::NAME_KEY,
4963 boundSrcHubAsGeogCRS->nameStr())
4964 .set(common::ObjectUsage::DOMAIN_OF_VALIDITY_KEY,
4965 metadata::Extent::WORLD),
4966 boundSrcHubAsGeogCRS->datum(),
4967 boundSrcHubAsGeogCRS->datumEnsemble(), cs));
4968 auto sourceToGeog3DOps =
4969 createOperations(sourceCRS, intermGeog3DCRS, context);
4970 auto geog3DToTargetOps =
4971 createOperations(intermGeog3DCRS, targetCRS, context);
4972 for (const auto &opSrc : sourceToGeog3DOps) {
4973 for (const auto &opDst : geog3DToTargetOps) {
4974 if (opSrc->targetCRS() && opDst->sourceCRS() &&
4975 !opSrc->targetCRS()->_isEquivalentTo(
4976 opDst->sourceCRS().get())) {
4977 // Shouldn't happen normally, but typically
4978 // one of them can be 2D and the other 3D
4979 // due to above createOperations() not
4980 // exactly setting the expected source and
4981 // target CRS.
4982 // So create an adapter operation...
4983 auto intermOps = createOperations(
4984 NN_NO_CHECK(opSrc->targetCRS()),
4985 NN_NO_CHECK(opDst->sourceCRS()), context);
4986 if (!intermOps.empty()) {
4987 res.emplace_back(
4988 ConcatenatedOperation::
4989 createComputeMetadata(
4990 {opSrc, intermOps.front(),
4991 opDst},
4992 disallowEmptyIntersection));
4993 }
4994 } else {
4995 res.emplace_back(
4996 ConcatenatedOperation::
4997 createComputeMetadata(
4998 {opSrc, opDst},
4999 disallowEmptyIntersection));
5000 }
5001 }
5002 }
5003 return;
5004 }
5005 }
5006 }
5007 }
5008
5009 // There might be better things to do, but for now just ignore the
5010 // transformation of the bound CRS
5011 res = createOperations(boundSrc->baseCRS(), targetCRS, context);
5012 }
5013 //! @endcond
5014
5015 // ---------------------------------------------------------------------------
5016
5017 /** \brief Find a list of CoordinateOperation from sourceCRS to targetCRS.
5018 *
5019 * The operations are sorted with the most relevant ones first: by
5020 * descending
5021 * area (intersection of the transformation area with the area of interest,
5022 * or intersection of the transformation with the area of use of the CRS),
5023 * and
5024 * by increasing accuracy. Operations with unknown accuracy are sorted last,
5025 * whatever their area.
5026 *
5027 * When one of the source or target CRS has a vertical component but not the
5028 * other one, the one that has no vertical component is automatically promoted
5029 * to a 3D version, where its vertical axis is the ellipsoidal height in metres,
5030 * using the ellipsoid of the base geodetic CRS.
5031 *
5032 * @param sourceCRS source CRS.
5033 * @param targetCRS target CRS.
5034 * @param context Search context.
5035 * @return a list
5036 */
5037 std::vector<CoordinateOperationNNPtr>
createOperations(const crs::CRSNNPtr & sourceCRS,const crs::CRSNNPtr & targetCRS,const CoordinateOperationContextNNPtr & context) const5038 CoordinateOperationFactory::createOperations(
5039 const crs::CRSNNPtr &sourceCRS, const crs::CRSNNPtr &targetCRS,
5040 const CoordinateOperationContextNNPtr &context) const {
5041
5042 #ifdef TRACE_CREATE_OPERATIONS
5043 ENTER_FUNCTION();
5044 #endif
5045 // Look if we are called on CRS that have a link to a 'canonical'
5046 // BoundCRS
5047 // If so, use that one as input
5048 const auto &srcBoundCRS = sourceCRS->canonicalBoundCRS();
5049 const auto &targetBoundCRS = targetCRS->canonicalBoundCRS();
5050 auto l_sourceCRS = srcBoundCRS ? NN_NO_CHECK(srcBoundCRS) : sourceCRS;
5051 auto l_targetCRS = targetBoundCRS ? NN_NO_CHECK(targetBoundCRS) : targetCRS;
5052 const auto &authFactory = context->getAuthorityFactory();
5053
5054 metadata::ExtentPtr sourceCRSExtent;
5055 auto l_resolvedSourceCRS =
5056 crs::CRS::getResolvedCRS(l_sourceCRS, authFactory, sourceCRSExtent);
5057 metadata::ExtentPtr targetCRSExtent;
5058 auto l_resolvedTargetCRS =
5059 crs::CRS::getResolvedCRS(l_targetCRS, authFactory, targetCRSExtent);
5060 Private::Context contextPrivate(sourceCRSExtent, targetCRSExtent, context);
5061
5062 if (context->getSourceAndTargetCRSExtentUse() ==
5063 CoordinateOperationContext::SourceTargetCRSExtentUse::INTERSECTION) {
5064 if (sourceCRSExtent && targetCRSExtent &&
5065 !sourceCRSExtent->intersects(NN_NO_CHECK(targetCRSExtent))) {
5066 return std::vector<CoordinateOperationNNPtr>();
5067 }
5068 }
5069
5070 return filterAndSort(Private::createOperations(l_resolvedSourceCRS,
5071 l_resolvedTargetCRS,
5072 contextPrivate),
5073 context, sourceCRSExtent, targetCRSExtent);
5074 }
5075
5076 // ---------------------------------------------------------------------------
5077
5078 /** \brief Instantiate a CoordinateOperationFactory.
5079 */
create()5080 CoordinateOperationFactoryNNPtr CoordinateOperationFactory::create() {
5081 return NN_NO_CHECK(
5082 CoordinateOperationFactory::make_unique<CoordinateOperationFactory>());
5083 }
5084
5085 // ---------------------------------------------------------------------------
5086
5087 } // namespace operation
5088
5089 namespace crs {
5090 // ---------------------------------------------------------------------------
5091
5092 //! @cond Doxygen_Suppress
5093
getResolvedCRS(const crs::CRSNNPtr & crs,const io::AuthorityFactoryPtr & authFactory,metadata::ExtentPtr & extentOut)5094 crs::CRSNNPtr CRS::getResolvedCRS(const crs::CRSNNPtr &crs,
5095 const io::AuthorityFactoryPtr &authFactory,
5096 metadata::ExtentPtr &extentOut) {
5097 const auto &ids = crs->identifiers();
5098 const auto &name = crs->nameStr();
5099
5100 bool approxExtent;
5101 extentOut = operation::getExtentPossiblySynthetized(crs, approxExtent);
5102
5103 // We try to "identify" the provided CRS with the ones of the database,
5104 // but in a more restricted way that what identify() does.
5105 // If we get a match from id in priority, and from name as a fallback, and
5106 // that they are equivalent to the input CRS, then use the identified CRS.
5107 // Even if they aren't equivalent, we update extentOut with the one of the
5108 // identified CRS if our input one is absent/not reliable.
5109
5110 const auto tryToIdentifyByName = [&crs, &name, &authFactory, approxExtent,
5111 &extentOut](
5112 io::AuthorityFactory::ObjectType objectType) {
5113 if (name != "unknown" && name != "unnamed") {
5114 auto matches = authFactory->createObjectsFromName(
5115 name, {objectType}, false, 2);
5116 if (matches.size() == 1) {
5117 const auto match =
5118 util::nn_static_pointer_cast<crs::CRS>(matches.front());
5119 if (approxExtent || !extentOut) {
5120 extentOut = operation::getExtent(match);
5121 }
5122 if (match->isEquivalentTo(
5123 crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
5124 return match;
5125 }
5126 }
5127 }
5128 return crs;
5129 };
5130
5131 auto geogCRS = dynamic_cast<crs::GeographicCRS *>(crs.get());
5132 if (geogCRS && authFactory) {
5133 if (!ids.empty()) {
5134 const auto tmpAuthFactory = io::AuthorityFactory::create(
5135 authFactory->databaseContext(), *ids.front()->codeSpace());
5136 try {
5137 auto resolvedCrs(
5138 tmpAuthFactory->createGeographicCRS(ids.front()->code()));
5139 if (approxExtent || !extentOut) {
5140 extentOut = operation::getExtent(resolvedCrs);
5141 }
5142 if (resolvedCrs->isEquivalentTo(
5143 crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
5144 return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
5145 }
5146 } catch (const std::exception &) {
5147 }
5148 } else {
5149 return tryToIdentifyByName(
5150 geogCRS->coordinateSystem()->axisList().size() == 2
5151 ? io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS
5152 : io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS);
5153 }
5154 }
5155
5156 auto projectedCrs = dynamic_cast<crs::ProjectedCRS *>(crs.get());
5157 if (projectedCrs && authFactory) {
5158 if (!ids.empty()) {
5159 const auto tmpAuthFactory = io::AuthorityFactory::create(
5160 authFactory->databaseContext(), *ids.front()->codeSpace());
5161 try {
5162 auto resolvedCrs(
5163 tmpAuthFactory->createProjectedCRS(ids.front()->code()));
5164 if (approxExtent || !extentOut) {
5165 extentOut = operation::getExtent(resolvedCrs);
5166 }
5167 if (resolvedCrs->isEquivalentTo(
5168 crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
5169 return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
5170 }
5171 } catch (const std::exception &) {
5172 }
5173 } else {
5174 return tryToIdentifyByName(
5175 io::AuthorityFactory::ObjectType::PROJECTED_CRS);
5176 }
5177 }
5178
5179 auto compoundCrs = dynamic_cast<crs::CompoundCRS *>(crs.get());
5180 if (compoundCrs && authFactory) {
5181 if (!ids.empty()) {
5182 const auto tmpAuthFactory = io::AuthorityFactory::create(
5183 authFactory->databaseContext(), *ids.front()->codeSpace());
5184 try {
5185 auto resolvedCrs(
5186 tmpAuthFactory->createCompoundCRS(ids.front()->code()));
5187 if (approxExtent || !extentOut) {
5188 extentOut = operation::getExtent(resolvedCrs);
5189 }
5190 if (resolvedCrs->isEquivalentTo(
5191 crs.get(), util::IComparable::Criterion::EQUIVALENT)) {
5192 return util::nn_static_pointer_cast<crs::CRS>(resolvedCrs);
5193 }
5194 } catch (const std::exception &) {
5195 }
5196 } else {
5197 auto outCrs = tryToIdentifyByName(
5198 io::AuthorityFactory::ObjectType::COMPOUND_CRS);
5199 const auto &components = compoundCrs->componentReferenceSystems();
5200 if (outCrs.get() != crs.get()) {
5201 bool hasGeoid = false;
5202 if (components.size() == 2) {
5203 auto vertCRS =
5204 dynamic_cast<crs::VerticalCRS *>(components[1].get());
5205 if (vertCRS && !vertCRS->geoidModel().empty()) {
5206 hasGeoid = true;
5207 }
5208 }
5209 if (!hasGeoid) {
5210 return outCrs;
5211 }
5212 }
5213 if (approxExtent || !extentOut) {
5214 // If we still did not get a reliable extent, then try to
5215 // resolve the components of the compoundCRS, and take the
5216 // intersection of their extent.
5217 extentOut = metadata::ExtentPtr();
5218 for (const auto &component : components) {
5219 metadata::ExtentPtr componentExtent;
5220 getResolvedCRS(component, authFactory, componentExtent);
5221 if (extentOut && componentExtent)
5222 extentOut = extentOut->intersection(
5223 NN_NO_CHECK(componentExtent));
5224 else if (componentExtent)
5225 extentOut = componentExtent;
5226 }
5227 }
5228 }
5229 }
5230 return crs;
5231 }
5232
5233 //! @endcond
5234
5235 } // namespace crs
5236 NS_PROJ_END
5237