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 //! @cond Doxygen_Suppress
34 #define DO_NOT_DEFINE_EXTERN_DERIVED_CRS_TEMPLATE
35 //! @endcond
36 
37 #include "proj/crs.hpp"
38 #include "proj/common.hpp"
39 #include "proj/coordinateoperation.hpp"
40 #include "proj/coordinatesystem.hpp"
41 #include "proj/io.hpp"
42 #include "proj/util.hpp"
43 
44 #include "proj/internal/coordinatesystem_internal.hpp"
45 #include "proj/internal/internal.hpp"
46 #include "proj/internal/io_internal.hpp"
47 
48 #include "proj_constants.h"
49 #include "proj_json_streaming_writer.hpp"
50 
51 #include <algorithm>
52 #include <cassert>
53 #include <cmath>
54 #include <cstring>
55 #include <iostream>
56 #include <memory>
57 #include <string>
58 #include <vector>
59 
60 using namespace NS_PROJ::internal;
61 
62 #if 0
63 namespace dropbox{ namespace oxygen {
64 template<> nn<NS_PROJ::crs::CRSPtr>::~nn() = default;
65 template<> nn<NS_PROJ::crs::SingleCRSPtr>::~nn() = default;
66 template<> nn<NS_PROJ::crs::GeodeticCRSPtr>::~nn() = default;
67 template<> nn<NS_PROJ::crs::GeographicCRSPtr>::~nn() = default;
68 template<> nn<NS_PROJ::crs::DerivedCRSPtr>::~nn() = default;
69 template<> nn<NS_PROJ::crs::ProjectedCRSPtr>::~nn() = default;
70 template<> nn<NS_PROJ::crs::VerticalCRSPtr>::~nn() = default;
71 template<> nn<NS_PROJ::crs::CompoundCRSPtr>::~nn() = default;
72 template<> nn<NS_PROJ::crs::TemporalCRSPtr>::~nn() = default;
73 template<> nn<NS_PROJ::crs::EngineeringCRSPtr>::~nn() = default;
74 template<> nn<NS_PROJ::crs::ParametricCRSPtr>::~nn() = default;
75 template<> nn<NS_PROJ::crs::BoundCRSPtr>::~nn() = default;
76 template<> nn<NS_PROJ::crs::DerivedGeodeticCRSPtr>::~nn() = default;
77 template<> nn<NS_PROJ::crs::DerivedGeographicCRSPtr>::~nn() = default;
78 template<> nn<NS_PROJ::crs::DerivedProjectedCRSPtr>::~nn() = default;
79 template<> nn<NS_PROJ::crs::DerivedVerticalCRSPtr>::~nn() = default;
80 template<> nn<NS_PROJ::crs::DerivedTemporalCRSPtr>::~nn() = default;
81 template<> nn<NS_PROJ::crs::DerivedEngineeringCRSPtr>::~nn() = default;
82 template<> nn<NS_PROJ::crs::DerivedParametricCRSPtr>::~nn() = default;
83 }}
84 #endif
85 
86 NS_PROJ_START
87 
88 namespace crs {
89 
90 // ---------------------------------------------------------------------------
91 
92 //! @cond Doxygen_Suppress
93 struct CRS::Private {
94     BoundCRSPtr canonicalBoundCRS_{};
95     std::string extensionProj4_{};
96     bool implicitCS_ = false;
97 
98     bool allowNonConformantWKT1Export_ = false;
99     // for what was initially a COMPD_CS with a VERT_CS with a datum type ==
100     // ellipsoidal height / 2002
101     CompoundCRSPtr originalCompoundCRS_{};
102 
setImplicitCScrs::CRS::Private103     void setImplicitCS(const util::PropertyMap &properties) {
104         const auto pVal = properties.get("IMPLICIT_CS");
105         if (pVal) {
106             if (const auto genVal =
107                     dynamic_cast<const util::BoxedValue *>(pVal->get())) {
108                 if (genVal->type() == util::BoxedValue::Type::BOOLEAN &&
109                     genVal->booleanValue()) {
110                     implicitCS_ = true;
111                 }
112             }
113         }
114     }
115 };
116 //! @endcond
117 
118 // ---------------------------------------------------------------------------
119 
CRS()120 CRS::CRS() : d(internal::make_unique<Private>()) {}
121 
122 // ---------------------------------------------------------------------------
123 
CRS(const CRS & other)124 CRS::CRS(const CRS &other)
125     : ObjectUsage(other), d(internal::make_unique<Private>(*(other.d))) {}
126 
127 // ---------------------------------------------------------------------------
128 
129 //! @cond Doxygen_Suppress
130 CRS::~CRS() = default;
131 //! @endcond
132 
133 // ---------------------------------------------------------------------------
134 
135 //! @cond Doxygen_Suppress
136 
137 /** \brief Return whether the CRS has an implicit coordinate system
138  * (e.g from ESRI WKT) */
hasImplicitCS() const139 bool CRS::hasImplicitCS() const { return d->implicitCS_; }
140 
141 //! @endcond
142 
143 // ---------------------------------------------------------------------------
144 
145 /** \brief Return the BoundCRS potentially attached to this CRS.
146  *
147  * In the case this method is called on a object returned by
148  * BoundCRS::baseCRSWithCanonicalBoundCRS(), this method will return this
149  * BoundCRS
150  *
151  * @return a BoundCRSPtr, that might be null.
152  */
canonicalBoundCRS()153 const BoundCRSPtr &CRS::canonicalBoundCRS() PROJ_PURE_DEFN {
154     return d->canonicalBoundCRS_;
155 }
156 
157 // ---------------------------------------------------------------------------
158 
159 /** \brief Return the GeodeticCRS of the CRS.
160  *
161  * Returns the GeodeticCRS contained in a CRS. This works currently with
162  * input parameters of type GeodeticCRS or derived, ProjectedCRS,
163  * CompoundCRS or BoundCRS.
164  *
165  * @return a GeodeticCRSPtr, that might be null.
166  */
extractGeodeticCRS() const167 GeodeticCRSPtr CRS::extractGeodeticCRS() const {
168     auto raw = extractGeodeticCRSRaw();
169     if (raw) {
170         return std::dynamic_pointer_cast<GeodeticCRS>(
171             raw->shared_from_this().as_nullable());
172     }
173     return nullptr;
174 }
175 
176 // ---------------------------------------------------------------------------
177 
178 //! @cond Doxygen_Suppress
extractGeodeticCRSRaw() const179 const GeodeticCRS *CRS::extractGeodeticCRSRaw() const {
180     auto geodCRS = dynamic_cast<const GeodeticCRS *>(this);
181     if (geodCRS) {
182         return geodCRS;
183     }
184     auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
185     if (projCRS) {
186         return projCRS->baseCRS()->extractGeodeticCRSRaw();
187     }
188     auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
189     if (compoundCRS) {
190         for (const auto &subCrs : compoundCRS->componentReferenceSystems()) {
191             auto retGeogCRS = subCrs->extractGeodeticCRSRaw();
192             if (retGeogCRS) {
193                 return retGeogCRS;
194             }
195         }
196     }
197     auto boundCRS = dynamic_cast<const BoundCRS *>(this);
198     if (boundCRS) {
199         return boundCRS->baseCRS()->extractGeodeticCRSRaw();
200     }
201     return nullptr;
202 }
203 //! @endcond
204 
205 // ---------------------------------------------------------------------------
206 
207 //! @cond Doxygen_Suppress
getExtensionProj4() const208 const std::string &CRS::getExtensionProj4() const noexcept {
209     return d->extensionProj4_;
210 }
211 //! @endcond
212 
213 // ---------------------------------------------------------------------------
214 
215 /** \brief Return the GeographicCRS of the CRS.
216  *
217  * Returns the GeographicCRS contained in a CRS. This works currently with
218  * input parameters of type GeographicCRS or derived, ProjectedCRS,
219  * CompoundCRS or BoundCRS.
220  *
221  * @return a GeographicCRSPtr, that might be null.
222  */
extractGeographicCRS() const223 GeographicCRSPtr CRS::extractGeographicCRS() const {
224     auto raw = extractGeodeticCRSRaw();
225     if (raw) {
226         return std::dynamic_pointer_cast<GeographicCRS>(
227             raw->shared_from_this().as_nullable());
228     }
229     return nullptr;
230 }
231 
232 // ---------------------------------------------------------------------------
233 
234 //! @cond Doxygen_Suppress
235 static util::PropertyMap
createPropertyMap(const common::IdentifiedObject * obj)236 createPropertyMap(const common::IdentifiedObject *obj) {
237     auto props = util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
238                                          obj->nameStr());
239     if (obj->isDeprecated()) {
240         props.set(common::IdentifiedObject::DEPRECATED_KEY, true);
241     }
242     return props;
243 }
244 
245 //! @endcond
246 
247 // ---------------------------------------------------------------------------
248 
249 //! @cond Doxygen_Suppress
alterGeodeticCRS(const GeodeticCRSNNPtr & newGeodCRS) const250 CRSNNPtr CRS::alterGeodeticCRS(const GeodeticCRSNNPtr &newGeodCRS) const {
251     auto geodCRS = dynamic_cast<const GeodeticCRS *>(this);
252     if (geodCRS) {
253         return newGeodCRS;
254     }
255 
256     auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
257     if (projCRS) {
258         return ProjectedCRS::create(createPropertyMap(this), newGeodCRS,
259                                     projCRS->derivingConversion(),
260                                     projCRS->coordinateSystem());
261     }
262 
263     auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
264     if (compoundCRS) {
265         std::vector<CRSNNPtr> components;
266         for (const auto &subCrs : compoundCRS->componentReferenceSystems()) {
267             components.emplace_back(subCrs->alterGeodeticCRS(newGeodCRS));
268         }
269         return CompoundCRS::create(createPropertyMap(this), components);
270     }
271 
272     return NN_NO_CHECK(
273         std::dynamic_pointer_cast<CRS>(shared_from_this().as_nullable()));
274 }
275 //! @endcond
276 
277 // ---------------------------------------------------------------------------
278 
279 //! @cond Doxygen_Suppress
alterCSLinearUnit(const common::UnitOfMeasure & unit) const280 CRSNNPtr CRS::alterCSLinearUnit(const common::UnitOfMeasure &unit) const {
281     {
282         auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
283         if (projCRS) {
284             return ProjectedCRS::create(
285                 createPropertyMap(this), projCRS->baseCRS(),
286                 projCRS->derivingConversion(),
287                 projCRS->coordinateSystem()->alterUnit(unit));
288         }
289     }
290 
291     {
292         auto geodCRS = dynamic_cast<const GeodeticCRS *>(this);
293         if (geodCRS && geodCRS->isGeocentric()) {
294             auto cs = dynamic_cast<const cs::CartesianCS *>(
295                 geodCRS->coordinateSystem().get());
296             assert(cs);
297             return GeodeticCRS::create(
298                 createPropertyMap(this), geodCRS->datum(),
299                 geodCRS->datumEnsemble(), cs->alterUnit(unit));
300         }
301     }
302 
303     {
304         auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
305         if (geogCRS && geogCRS->coordinateSystem()->axisList().size() == 3) {
306             return GeographicCRS::create(
307                 createPropertyMap(this), geogCRS->datum(),
308                 geogCRS->datumEnsemble(),
309                 geogCRS->coordinateSystem()->alterLinearUnit(unit));
310         }
311     }
312 
313     {
314         auto vertCRS = dynamic_cast<const VerticalCRS *>(this);
315         if (vertCRS) {
316             return VerticalCRS::create(
317                 createPropertyMap(this), vertCRS->datum(),
318                 vertCRS->datumEnsemble(),
319                 vertCRS->coordinateSystem()->alterUnit(unit));
320         }
321     }
322 
323     {
324         auto engCRS = dynamic_cast<const EngineeringCRS *>(this);
325         if (engCRS) {
326             auto cartCS = util::nn_dynamic_pointer_cast<cs::CartesianCS>(
327                 engCRS->coordinateSystem());
328             if (cartCS) {
329                 return EngineeringCRS::create(createPropertyMap(this),
330                                               engCRS->datum(),
331                                               cartCS->alterUnit(unit));
332             } else {
333                 auto vertCS = util::nn_dynamic_pointer_cast<cs::VerticalCS>(
334                     engCRS->coordinateSystem());
335                 if (vertCS) {
336                     return EngineeringCRS::create(createPropertyMap(this),
337                                                   engCRS->datum(),
338                                                   vertCS->alterUnit(unit));
339                 }
340             }
341         }
342     }
343 
344     return NN_NO_CHECK(
345         std::dynamic_pointer_cast<CRS>(shared_from_this().as_nullable()));
346 }
347 //! @endcond
348 
349 // ---------------------------------------------------------------------------
350 
351 /** \brief Return the VerticalCRS of the CRS.
352  *
353  * Returns the VerticalCRS contained in a CRS. This works currently with
354  * input parameters of type VerticalCRS or derived, CompoundCRS or BoundCRS.
355  *
356  * @return a VerticalCRSPtr, that might be null.
357  */
extractVerticalCRS() const358 VerticalCRSPtr CRS::extractVerticalCRS() const {
359     auto vertCRS = dynamic_cast<const VerticalCRS *>(this);
360     if (vertCRS) {
361         return std::dynamic_pointer_cast<VerticalCRS>(
362             shared_from_this().as_nullable());
363     }
364     auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
365     if (compoundCRS) {
366         for (const auto &subCrs : compoundCRS->componentReferenceSystems()) {
367             auto retVertCRS = subCrs->extractVerticalCRS();
368             if (retVertCRS) {
369                 return retVertCRS;
370             }
371         }
372     }
373     auto boundCRS = dynamic_cast<const BoundCRS *>(this);
374     if (boundCRS) {
375         return boundCRS->baseCRS()->extractVerticalCRS();
376     }
377     return nullptr;
378 }
379 
380 // ---------------------------------------------------------------------------
381 
382 /** \brief Returns potentially
383  * a BoundCRS, with a transformation to EPSG:4326, wrapping this CRS
384  *
385  * If no such BoundCRS is possible, the object will be returned.
386  *
387  * The purpose of this method is to be able to format a PROJ.4 string with
388  * a +towgs84 parameter or a WKT1:GDAL string with a TOWGS node.
389  *
390  * This method will fetch the GeographicCRS of this CRS and find a
391  * transformation to EPSG:4326 using the domain of the validity of the main CRS,
392  * and there's only one Helmert transformation.
393  *
394  * @return a CRS.
395  */
createBoundCRSToWGS84IfPossible(const io::DatabaseContextPtr & dbContext,operation::CoordinateOperationContext::IntermediateCRSUse allowIntermediateCRSUse) const396 CRSNNPtr CRS::createBoundCRSToWGS84IfPossible(
397     const io::DatabaseContextPtr &dbContext,
398     operation::CoordinateOperationContext::IntermediateCRSUse
399         allowIntermediateCRSUse) const {
400     auto thisAsCRS = NN_NO_CHECK(
401         std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
402     auto boundCRS = util::nn_dynamic_pointer_cast<BoundCRS>(thisAsCRS);
403     if (!boundCRS) {
404         boundCRS = canonicalBoundCRS();
405     }
406     if (boundCRS) {
407         if (boundCRS->hubCRS()->_isEquivalentTo(
408                 GeographicCRS::EPSG_4326.get(),
409                 util::IComparable::Criterion::EQUIVALENT, dbContext)) {
410             return NN_NO_CHECK(boundCRS);
411         }
412     }
413 
414     auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
415     if (compoundCRS) {
416         const auto &comps = compoundCRS->componentReferenceSystems();
417         if (comps.size() == 2) {
418             auto horiz = comps[0]->createBoundCRSToWGS84IfPossible(
419                 dbContext, allowIntermediateCRSUse);
420             auto vert = comps[1]->createBoundCRSToWGS84IfPossible(
421                 dbContext, allowIntermediateCRSUse);
422             if (horiz.get() != comps[0].get() || vert.get() != comps[1].get()) {
423                 return CompoundCRS::create(createPropertyMap(this),
424                                            {horiz, vert});
425             }
426         }
427         return thisAsCRS;
428     }
429 
430     if (!dbContext) {
431         return thisAsCRS;
432     }
433 
434     const auto &l_domains = domains();
435     metadata::ExtentPtr extent;
436     if (!l_domains.empty()) {
437         extent = l_domains[0]->domainOfValidity();
438     }
439 
440     std::string crs_authority;
441     const auto &l_identifiers = identifiers();
442     // If the object has an authority, restrict the transformations to
443     // come from that codespace too. This avoids for example EPSG:4269
444     // (NAD83) to use a (dubious) ESRI transformation.
445     if (!l_identifiers.empty()) {
446         crs_authority = *(l_identifiers[0]->codeSpace());
447     }
448 
449     auto authorities = dbContext->getAllowedAuthorities(crs_authority, "EPSG");
450     if (authorities.empty()) {
451         authorities.emplace_back();
452     }
453 
454     // Vertical CRS ?
455     auto vertCRS = dynamic_cast<const VerticalCRS *>(this);
456     if (vertCRS) {
457         auto hubCRS =
458             util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4979);
459         for (const auto &authority : authorities) {
460             try {
461 
462                 auto authFactory = io::AuthorityFactory::create(
463                     NN_NO_CHECK(dbContext),
464                     authority == "any" ? std::string() : authority);
465                 auto ctxt = operation::CoordinateOperationContext::create(
466                     authFactory, extent, 0.0);
467                 ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse);
468                 // ctxt->setSpatialCriterion(
469                 //    operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
470                 auto list = operation::CoordinateOperationFactory::create()
471                                 ->createOperations(hubCRS, thisAsCRS, ctxt);
472                 CRSPtr candidateBoundCRS;
473                 for (const auto &op : list) {
474                     auto transf = util::nn_dynamic_pointer_cast<
475                         operation::Transformation>(op);
476                     // Only keep transformations that use a known grid
477                     if (transf && !transf->hasBallparkTransformation()) {
478                         auto gridsNeeded = transf->gridsNeeded(dbContext, true);
479                         bool gridsKnown = !gridsNeeded.empty();
480                         for (const auto &gridDesc : gridsNeeded) {
481                             if (gridDesc.packageName.empty() &&
482                                 !(!gridDesc.url.empty() &&
483                                   gridDesc.openLicense) &&
484                                 !gridDesc.available) {
485                                 gridsKnown = false;
486                                 break;
487                             }
488                         }
489                         if (gridsKnown) {
490                             if (candidateBoundCRS) {
491                                 candidateBoundCRS = nullptr;
492                                 break;
493                             }
494                             candidateBoundCRS =
495                                 BoundCRS::create(thisAsCRS, hubCRS,
496                                                  NN_NO_CHECK(transf))
497                                     .as_nullable();
498                         }
499                     }
500                 }
501                 if (candidateBoundCRS) {
502                     return NN_NO_CHECK(candidateBoundCRS);
503                 }
504             } catch (const std::exception &) {
505             }
506         }
507         return thisAsCRS;
508     }
509 
510     // Geodetic/geographic CRS ?
511     auto geodCRS = util::nn_dynamic_pointer_cast<GeodeticCRS>(thisAsCRS);
512     auto geogCRS = extractGeographicCRS();
513     auto hubCRS = util::nn_static_pointer_cast<CRS>(GeographicCRS::EPSG_4326);
514     if (geodCRS && !geogCRS) {
515         if (geodCRS->_isEquivalentTo(GeographicCRS::EPSG_4978.get(),
516                                      util::IComparable::Criterion::EQUIVALENT,
517                                      dbContext)) {
518             return thisAsCRS;
519         }
520         hubCRS = util::nn_static_pointer_cast<CRS>(GeodeticCRS::EPSG_4978);
521     } else if (!geogCRS ||
522                geogCRS->_isEquivalentTo(
523                    GeographicCRS::EPSG_4326.get(),
524                    util::IComparable::Criterion::EQUIVALENT, dbContext)) {
525         return thisAsCRS;
526     } else {
527         geodCRS = geogCRS;
528     }
529 
530     for (const auto &authority : authorities) {
531         try {
532 
533             auto authFactory = io::AuthorityFactory::create(
534                 NN_NO_CHECK(dbContext),
535                 authority == "any" ? std::string() : authority);
536             metadata::ExtentPtr extentResolved(extent);
537             if (!extent) {
538                 getResolvedCRS(thisAsCRS, authFactory, extentResolved);
539             }
540             auto ctxt = operation::CoordinateOperationContext::create(
541                 authFactory, extentResolved, 0.0);
542             ctxt->setAllowUseIntermediateCRS(allowIntermediateCRSUse);
543             // ctxt->setSpatialCriterion(
544             //    operation::CoordinateOperationContext::SpatialCriterion::PARTIAL_INTERSECTION);
545             auto list =
546                 operation::CoordinateOperationFactory::create()
547                     ->createOperations(NN_NO_CHECK(geodCRS), hubCRS, ctxt);
548             CRSPtr candidateBoundCRS;
549             for (const auto &op : list) {
550                 auto transf =
551                     util::nn_dynamic_pointer_cast<operation::Transformation>(
552                         op);
553                 if (transf && !starts_with(transf->nameStr(), "Ballpark geo")) {
554                     try {
555                         transf->getTOWGS84Parameters();
556                     } catch (const std::exception &) {
557                         continue;
558                     }
559                     if (candidateBoundCRS) {
560                         candidateBoundCRS = nullptr;
561                         break;
562                     }
563                     candidateBoundCRS =
564                         BoundCRS::create(thisAsCRS, hubCRS, NN_NO_CHECK(transf))
565                             .as_nullable();
566                 } else {
567                     auto concatenated =
568                         dynamic_cast<const operation::ConcatenatedOperation *>(
569                             op.get());
570                     if (concatenated) {
571                         // Case for EPSG:4807 / "NTF (Paris)" that is made of a
572                         // longitude rotation followed by a Helmert
573                         // The prime meridian shift will be accounted elsewhere
574                         const auto &subops = concatenated->operations();
575                         if (subops.size() == 2) {
576                             auto firstOpIsTransformation =
577                                 dynamic_cast<const operation::Transformation *>(
578                                     subops[0].get());
579                             auto firstOpIsConversion =
580                                 dynamic_cast<const operation::Conversion *>(
581                                     subops[0].get());
582                             if ((firstOpIsTransformation &&
583                                  firstOpIsTransformation
584                                      ->isLongitudeRotation()) ||
585                                 (dynamic_cast<DerivedCRS *>(thisAsCRS.get()) &&
586                                  firstOpIsConversion)) {
587                                 transf = util::nn_dynamic_pointer_cast<
588                                     operation::Transformation>(subops[1]);
589                                 if (transf &&
590                                     !starts_with(transf->nameStr(),
591                                                  "Ballpark geo")) {
592                                     try {
593                                         transf->getTOWGS84Parameters();
594                                     } catch (const std::exception &) {
595                                         continue;
596                                     }
597                                     if (candidateBoundCRS) {
598                                         candidateBoundCRS = nullptr;
599                                         break;
600                                     }
601                                     candidateBoundCRS =
602                                         BoundCRS::create(thisAsCRS, hubCRS,
603                                                          NN_NO_CHECK(transf))
604                                             .as_nullable();
605                                 }
606                             }
607                         }
608                     }
609                 }
610             }
611             if (candidateBoundCRS) {
612                 return NN_NO_CHECK(candidateBoundCRS);
613             }
614         } catch (const std::exception &) {
615         }
616     }
617     return thisAsCRS;
618 }
619 
620 // ---------------------------------------------------------------------------
621 
622 /** \brief Returns a CRS whose coordinate system does not contain a vertical
623  * component
624  *
625  * @return a CRS.
626  */
stripVerticalComponent() const627 CRSNNPtr CRS::stripVerticalComponent() const {
628     auto self = NN_NO_CHECK(
629         std::dynamic_pointer_cast<CRS>(shared_from_this().as_nullable()));
630 
631     auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
632     if (geogCRS) {
633         const auto &axisList = geogCRS->coordinateSystem()->axisList();
634         if (axisList.size() == 3) {
635             auto cs = cs::EllipsoidalCS::create(util::PropertyMap(),
636                                                 axisList[0], axisList[1]);
637             return util::nn_static_pointer_cast<CRS>(GeographicCRS::create(
638                 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
639                                         nameStr()),
640                 geogCRS->datum(), geogCRS->datumEnsemble(), cs));
641         }
642     }
643     auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
644     if (projCRS) {
645         const auto &axisList = projCRS->coordinateSystem()->axisList();
646         if (axisList.size() == 3) {
647             auto cs = cs::CartesianCS::create(util::PropertyMap(), axisList[0],
648                                               axisList[1]);
649             return util::nn_static_pointer_cast<CRS>(ProjectedCRS::create(
650                 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
651                                         nameStr()),
652                 projCRS->baseCRS(), projCRS->derivingConversion(), cs));
653         }
654     }
655     return self;
656 }
657 
658 // ---------------------------------------------------------------------------
659 
660 //! @cond Doxygen_Suppress
661 
662 /** \brief Return a shallow clone of this object. */
shallowClone() const663 CRSNNPtr CRS::shallowClone() const { return _shallowClone(); }
664 
665 //! @endcond
666 
667 // ---------------------------------------------------------------------------
668 
669 //! @cond Doxygen_Suppress
670 
allowNonConformantWKT1Export() const671 CRSNNPtr CRS::allowNonConformantWKT1Export() const {
672     const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
673     if (boundCRS) {
674         return BoundCRS::create(
675             boundCRS->baseCRS()->allowNonConformantWKT1Export(),
676             boundCRS->hubCRS(), boundCRS->transformation());
677     }
678     auto crs(shallowClone());
679     crs->d->allowNonConformantWKT1Export_ = true;
680     return crs;
681 }
682 
683 //! @endcond
684 
685 // ---------------------------------------------------------------------------
686 
687 //! @cond Doxygen_Suppress
688 
689 CRSNNPtr
attachOriginalCompoundCRS(const CompoundCRSNNPtr & compoundCRS) const690 CRS::attachOriginalCompoundCRS(const CompoundCRSNNPtr &compoundCRS) const {
691 
692     const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
693     if (boundCRS) {
694         return BoundCRS::create(
695             boundCRS->baseCRS()->attachOriginalCompoundCRS(compoundCRS),
696             boundCRS->hubCRS(), boundCRS->transformation());
697     }
698 
699     auto crs(shallowClone());
700     crs->d->originalCompoundCRS_ = compoundCRS.as_nullable();
701     return crs;
702 }
703 
704 //! @endcond
705 
706 // ---------------------------------------------------------------------------
707 
708 //! @cond Doxygen_Suppress
709 
alterName(const std::string & newName) const710 CRSNNPtr CRS::alterName(const std::string &newName) const {
711     auto crs = shallowClone();
712     auto newNameMod(newName);
713     auto props = util::PropertyMap();
714     if (ends_with(newNameMod, " (deprecated)")) {
715         newNameMod.resize(newNameMod.size() - strlen(" (deprecated)"));
716         props.set(common::IdentifiedObject::DEPRECATED_KEY, true);
717     }
718     props.set(common::IdentifiedObject::NAME_KEY, newNameMod);
719     crs->setProperties(props);
720     return crs;
721 }
722 
723 //! @endcond
724 
725 // ---------------------------------------------------------------------------
726 
727 //! @cond Doxygen_Suppress
728 
alterId(const std::string & authName,const std::string & code) const729 CRSNNPtr CRS::alterId(const std::string &authName,
730                       const std::string &code) const {
731     auto crs = shallowClone();
732     auto props = util::PropertyMap();
733     props.set(metadata::Identifier::CODESPACE_KEY, authName)
734         .set(metadata::Identifier::CODE_KEY, code);
735     crs->setProperties(props);
736     return crs;
737 }
738 
739 //! @endcond
740 
741 // ---------------------------------------------------------------------------
742 
743 //! @cond Doxygen_Suppress
744 
mustAxisOrderBeSwitchedForVisualizationInternal(const std::vector<cs::CoordinateSystemAxisNNPtr> & axisList)745 static bool mustAxisOrderBeSwitchedForVisualizationInternal(
746     const std::vector<cs::CoordinateSystemAxisNNPtr> &axisList) {
747     const auto &dir0 = axisList[0]->direction();
748     const auto &dir1 = axisList[1]->direction();
749     if (&dir0 == &cs::AxisDirection::NORTH &&
750         &dir1 == &cs::AxisDirection::EAST) {
751         return true;
752     }
753 
754     // Address EPSG:32661 "WGS 84 / UPS North (N,E)"
755     if (&dir0 == &cs::AxisDirection::SOUTH &&
756         &dir1 == &cs::AxisDirection::SOUTH) {
757         const auto &meridian0 = axisList[0]->meridian();
758         const auto &meridian1 = axisList[1]->meridian();
759         return meridian0 != nullptr && meridian1 != nullptr &&
760                std::abs(meridian0->longitude().convertToUnit(
761                             common::UnitOfMeasure::DEGREE) -
762                         180.0) < 1e-10 &&
763                std::abs(meridian1->longitude().convertToUnit(
764                             common::UnitOfMeasure::DEGREE) -
765                         90.0) < 1e-10;
766     }
767 
768     if (&dir0 == &cs::AxisDirection::NORTH &&
769         &dir1 == &cs::AxisDirection::NORTH) {
770         const auto &meridian0 = axisList[0]->meridian();
771         const auto &meridian1 = axisList[1]->meridian();
772         return meridian0 != nullptr && meridian1 != nullptr &&
773                ((
774                     // Address EPSG:32761 "WGS 84 / UPS South (N,E)"
775                     std::abs(meridian0->longitude().convertToUnit(
776                                  common::UnitOfMeasure::DEGREE) -
777                              0.0) < 1e-10 &&
778                     std::abs(meridian1->longitude().convertToUnit(
779                                  common::UnitOfMeasure::DEGREE) -
780                              90.0) < 1e-10) ||
781                 // Address EPSG:5482 "RSRGD2000 / RSPS2000"
782                 (std::abs(meridian0->longitude().convertToUnit(
783                               common::UnitOfMeasure::DEGREE) -
784                           180) < 1e-10 &&
785                  std::abs(meridian1->longitude().convertToUnit(
786                               common::UnitOfMeasure::DEGREE) -
787                           -90.0) < 1e-10));
788     }
789 
790     return false;
791 }
792 // ---------------------------------------------------------------------------
793 
mustAxisOrderBeSwitchedForVisualization() const794 bool CRS::mustAxisOrderBeSwitchedForVisualization() const {
795 
796     const CompoundCRS *compoundCRS = dynamic_cast<const CompoundCRS *>(this);
797     if (compoundCRS) {
798         const auto &comps = compoundCRS->componentReferenceSystems();
799         if (!comps.empty()) {
800             return comps[0]->mustAxisOrderBeSwitchedForVisualization();
801         }
802     }
803 
804     const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this);
805     if (geogCRS) {
806         return mustAxisOrderBeSwitchedForVisualizationInternal(
807             geogCRS->coordinateSystem()->axisList());
808     }
809 
810     const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this);
811     if (projCRS) {
812         return mustAxisOrderBeSwitchedForVisualizationInternal(
813             projCRS->coordinateSystem()->axisList());
814     }
815 
816     return false;
817 }
818 
819 //! @endcond
820 
821 // ---------------------------------------------------------------------------
822 
823 //! @cond Doxygen_Suppress
824 
normalizeForVisualization() const825 CRSNNPtr CRS::normalizeForVisualization() const {
826     auto props = util::PropertyMap().set(
827         common::IdentifiedObject::NAME_KEY,
828         nameStr() + " (with axis order normalized for visualization)");
829 
830     const CompoundCRS *compoundCRS = dynamic_cast<const CompoundCRS *>(this);
831     if (compoundCRS) {
832         const auto &comps = compoundCRS->componentReferenceSystems();
833         if (!comps.empty()) {
834             std::vector<CRSNNPtr> newComps;
835             newComps.emplace_back(comps[0]->normalizeForVisualization());
836             for (size_t i = 1; i < comps.size(); i++) {
837                 newComps.emplace_back(comps[i]);
838             }
839             return util::nn_static_pointer_cast<CRS>(
840                 CompoundCRS::create(props, newComps));
841         }
842     }
843 
844     const GeographicCRS *geogCRS = dynamic_cast<const GeographicCRS *>(this);
845     if (geogCRS) {
846         const auto &axisList = geogCRS->coordinateSystem()->axisList();
847         if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) {
848             auto cs = axisList.size() == 2
849                           ? cs::EllipsoidalCS::create(util::PropertyMap(),
850                                                       axisList[1], axisList[0])
851                           : cs::EllipsoidalCS::create(util::PropertyMap(),
852                                                       axisList[1], axisList[0],
853                                                       axisList[2]);
854             return util::nn_static_pointer_cast<CRS>(GeographicCRS::create(
855                 props, geogCRS->datum(), geogCRS->datumEnsemble(), cs));
856         }
857     }
858 
859     const ProjectedCRS *projCRS = dynamic_cast<const ProjectedCRS *>(this);
860     if (projCRS) {
861         const auto &axisList = projCRS->coordinateSystem()->axisList();
862         if (mustAxisOrderBeSwitchedForVisualizationInternal(axisList)) {
863             auto cs =
864                 axisList.size() == 2
865                     ? cs::CartesianCS::create(util::PropertyMap(), axisList[1],
866                                               axisList[0])
867                     : cs::CartesianCS::create(util::PropertyMap(), axisList[1],
868                                               axisList[0], axisList[2]);
869             return util::nn_static_pointer_cast<CRS>(ProjectedCRS::create(
870                 props, projCRS->baseCRS(), projCRS->derivingConversion(), cs));
871         }
872     }
873 
874     return NN_NO_CHECK(
875         std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
876 }
877 
878 //! @endcond
879 
880 // ---------------------------------------------------------------------------
881 
882 /** \brief Identify the CRS with reference CRSs.
883  *
884  * The candidate CRSs are either hard-coded, or looked in the database when
885  * authorityFactory is not null.
886  *
887  * Note that the implementation uses a set of heuristics to have a good
888  * compromise of successful identifications over execution time. It might miss
889  * legitimate matches in some circumstances.
890  *
891  * The method returns a list of matching reference CRS, and the percentage
892  * (0-100) of confidence in the match. The list is sorted by decreasing
893  * confidence.
894  * <ul>
895  * <li>100% means that the name of the reference entry
896  * perfectly matches the CRS name, and both are equivalent. In which case a
897  * single result is returned.
898  * Note: in the case of a GeographicCRS whose axis
899  * order is implicit in the input definition (for example ESRI WKT), then axis
900  * order is ignored for the purpose of identification. That is the CRS built
901  * from
902  * GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
903  * PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
904  * will be identified to EPSG:4326, but will not pass a
905  * isEquivalentTo(EPSG_4326, util::IComparable::Criterion::EQUIVALENT) test,
906  * but rather isEquivalentTo(EPSG_4326,
907  * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS)
908  * </li>
909  * <li>90% means that CRS are equivalent, but the names are not exactly the
910  * same.</li>
911  * <li>70% means that CRS are equivalent), but the names do not match at
912  * all.</li>
913  * <li>25% means that the CRS are not equivalent, but there is some similarity
914  * in
915  * the names.</li>
916  * </ul>
917  * Other confidence values may be returned by some specialized implementations.
918  *
919  * This is implemented for GeodeticCRS, ProjectedCRS, VerticalCRS and
920  * CompoundCRS.
921  *
922  * @param authorityFactory Authority factory (or null, but degraded
923  * functionality)
924  * @return a list of matching reference CRS, and the percentage (0-100) of
925  * confidence in the match.
926  */
927 std::list<std::pair<CRSNNPtr, int>>
identify(const io::AuthorityFactoryPtr & authorityFactory) const928 CRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
929     return _identify(authorityFactory);
930 }
931 
932 // ---------------------------------------------------------------------------
933 
934 /** \brief Return CRSs that are non-deprecated substitutes for the current CRS.
935  */
936 std::list<CRSNNPtr>
getNonDeprecated(const io::DatabaseContextNNPtr & dbContext) const937 CRS::getNonDeprecated(const io::DatabaseContextNNPtr &dbContext) const {
938     std::list<CRSNNPtr> res;
939     const auto &l_identifiers = identifiers();
940     if (l_identifiers.empty()) {
941         return res;
942     }
943     const char *tableName = nullptr;
944     if (dynamic_cast<const GeodeticCRS *>(this)) {
945         tableName = "geodetic_crs";
946     } else if (dynamic_cast<const ProjectedCRS *>(this)) {
947         tableName = "projected_crs";
948     } else if (dynamic_cast<const VerticalCRS *>(this)) {
949         tableName = "vertical_crs";
950     } else if (dynamic_cast<const CompoundCRS *>(this)) {
951         tableName = "compound_crs";
952     }
953     if (!tableName) {
954         return res;
955     }
956     const auto &id = l_identifiers[0];
957     auto tmpRes =
958         dbContext->getNonDeprecated(tableName, *(id->codeSpace()), id->code());
959     for (const auto &pair : tmpRes) {
960         res.emplace_back(io::AuthorityFactory::create(dbContext, pair.first)
961                              ->createCoordinateReferenceSystem(pair.second));
962     }
963     return res;
964 }
965 
966 // ---------------------------------------------------------------------------
967 
968 /** \brief Return a variant of this CRS "promoted" to a 3D one, if not already
969  * the case.
970  *
971  * The new axis will be ellipsoidal height, oriented upwards, and with metre
972  * units.
973  *
974  * @param newName Name of the new CRS. If empty, nameStr() will be used.
975  * @param dbContext Database context to look for potentially already registered
976  *                  3D CRS. May be nullptr.
977  * @return a new CRS promoted to 3D, or the current one if already 3D or not
978  * applicable.
979  * @since 6.3
980  */
promoteTo3D(const std::string & newName,const io::DatabaseContextPtr & dbContext) const981 CRSNNPtr CRS::promoteTo3D(const std::string &newName,
982                           const io::DatabaseContextPtr &dbContext) const {
983     auto upAxis = cs::CoordinateSystemAxis::create(
984         util::PropertyMap().set(IdentifiedObject::NAME_KEY,
985                                 cs::AxisName::Ellipsoidal_height),
986         cs::AxisAbbreviation::h, cs::AxisDirection::UP,
987         common::UnitOfMeasure::METRE);
988     return promoteTo3D(newName, dbContext, upAxis);
989 }
990 
991 // ---------------------------------------------------------------------------
992 
993 //! @cond Doxygen_Suppress
994 
promoteTo3D(const std::string & newName,const io::DatabaseContextPtr & dbContext,const cs::CoordinateSystemAxisNNPtr & verticalAxisIfNotAlreadyPresent) const995 CRSNNPtr CRS::promoteTo3D(const std::string &newName,
996                           const io::DatabaseContextPtr &dbContext,
997                           const cs::CoordinateSystemAxisNNPtr
998                               &verticalAxisIfNotAlreadyPresent) const {
999 
1000     const auto createProperties = [this, &newName]() {
1001         auto props =
1002             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
1003                                     !newName.empty() ? newName : nameStr());
1004         const auto &l_identifiers = identifiers();
1005         if (l_identifiers.size() == 1) {
1006             std::string remarks("Promoted to 3D from ");
1007             remarks += *(l_identifiers[0]->codeSpace());
1008             remarks += ':';
1009             remarks += l_identifiers[0]->code();
1010             props.set(common::IdentifiedObject::REMARKS_KEY, remarks);
1011         }
1012         return props;
1013     };
1014 
1015     const auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
1016     if (geogCRS) {
1017         const auto &axisList = geogCRS->coordinateSystem()->axisList();
1018         if (axisList.size() == 2) {
1019             const auto &l_identifiers = identifiers();
1020             // First check if there is a Geographic 3D CRS in the database
1021             // of the same name.
1022             // This is the common practice in the EPSG dataset.
1023             if (dbContext && l_identifiers.size() == 1) {
1024                 auto authFactory = io::AuthorityFactory::create(
1025                     NN_NO_CHECK(dbContext), *(l_identifiers[0]->codeSpace()));
1026                 auto res = authFactory->createObjectsFromName(
1027                     nameStr(),
1028                     {io::AuthorityFactory::ObjectType::GEOGRAPHIC_3D_CRS},
1029                     false);
1030                 if (!res.empty()) {
1031                     const auto &firstRes = res.front();
1032                     const auto firstResGeog =
1033                         dynamic_cast<GeographicCRS *>(firstRes.get());
1034                     const auto &firstResAxisList =
1035                         firstResGeog->coordinateSystem()->axisList();
1036                     if (firstResAxisList[2]->_isEquivalentTo(
1037                             verticalAxisIfNotAlreadyPresent.get(),
1038                             util::IComparable::Criterion::EQUIVALENT) &&
1039                         geogCRS->is2DPartOf3D(NN_NO_CHECK(firstResGeog),
1040                                               dbContext)) {
1041                         return NN_NO_CHECK(
1042                             util::nn_dynamic_pointer_cast<CRS>(firstRes));
1043                     }
1044                 }
1045             }
1046 
1047             auto cs = cs::EllipsoidalCS::create(
1048                 util::PropertyMap(), axisList[0], axisList[1],
1049                 verticalAxisIfNotAlreadyPresent);
1050             return util::nn_static_pointer_cast<CRS>(
1051                 GeographicCRS::create(createProperties(), geogCRS->datum(),
1052                                       geogCRS->datumEnsemble(), cs));
1053         }
1054     }
1055 
1056     const auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
1057     if (projCRS) {
1058         const auto &axisList = projCRS->coordinateSystem()->axisList();
1059         if (axisList.size() == 2) {
1060             auto base3DCRS =
1061                 projCRS->baseCRS()->promoteTo3D(std::string(), dbContext);
1062             auto cs = cs::CartesianCS::create(util::PropertyMap(), axisList[0],
1063                                               axisList[1],
1064                                               verticalAxisIfNotAlreadyPresent);
1065             return util::nn_static_pointer_cast<CRS>(ProjectedCRS::create(
1066                 createProperties(),
1067                 NN_NO_CHECK(
1068                     util::nn_dynamic_pointer_cast<GeodeticCRS>(base3DCRS)),
1069                 projCRS->derivingConversion(), cs));
1070         }
1071     }
1072 
1073     const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
1074     if (boundCRS) {
1075         auto base3DCRS = boundCRS->baseCRS()->promoteTo3D(
1076             newName, dbContext, verticalAxisIfNotAlreadyPresent);
1077         auto transf = boundCRS->transformation();
1078         try {
1079             transf->getTOWGS84Parameters();
1080             return BoundCRS::create(
1081                 base3DCRS,
1082                 boundCRS->hubCRS()->promoteTo3D(std::string(), dbContext),
1083                 transf->promoteTo3D(std::string(), dbContext));
1084         } catch (const io::FormattingException &) {
1085             return BoundCRS::create(base3DCRS, boundCRS->hubCRS(), transf);
1086         }
1087     }
1088 
1089     return NN_NO_CHECK(
1090         std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
1091 }
1092 
1093 //! @endcond
1094 
1095 // ---------------------------------------------------------------------------
1096 
1097 /** \brief Return a variant of this CRS "demoted" to a 2D one, if not already
1098  * the case.
1099  *
1100  *
1101  * @param newName Name of the new CRS. If empty, nameStr() will be used.
1102  * @param dbContext Database context to look for potentially already registered
1103  *                  2D CRS. May be nullptr.
1104  * @return a new CRS demoted to 2D, or the current one if already 2D or not
1105  * applicable.
1106  * @since 6.3
1107  */
demoteTo2D(const std::string & newName,const io::DatabaseContextPtr & dbContext) const1108 CRSNNPtr CRS::demoteTo2D(const std::string &newName,
1109                          const io::DatabaseContextPtr &dbContext) const {
1110     const auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
1111     if (geogCRS) {
1112         return geogCRS->demoteTo2D(newName, dbContext);
1113     }
1114 
1115     const auto projCRS = dynamic_cast<const ProjectedCRS *>(this);
1116     if (projCRS) {
1117         return projCRS->demoteTo2D(newName, dbContext);
1118     }
1119 
1120     const auto boundCRS = dynamic_cast<const BoundCRS *>(this);
1121     if (boundCRS) {
1122         auto base2DCRS = boundCRS->baseCRS()->demoteTo2D(newName, dbContext);
1123         auto transf = boundCRS->transformation();
1124         try {
1125             transf->getTOWGS84Parameters();
1126             return BoundCRS::create(
1127                 base2DCRS,
1128                 boundCRS->hubCRS()->demoteTo2D(std::string(), dbContext),
1129                 transf->demoteTo2D(std::string(), dbContext));
1130         } catch (const io::FormattingException &) {
1131             return BoundCRS::create(base2DCRS, boundCRS->hubCRS(), transf);
1132         }
1133     }
1134 
1135     const auto compoundCRS = dynamic_cast<const CompoundCRS *>(this);
1136     if (compoundCRS) {
1137         const auto &components = compoundCRS->componentReferenceSystems();
1138         if (components.size() >= 2) {
1139             return components[0];
1140         }
1141     }
1142 
1143     return NN_NO_CHECK(
1144         std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
1145 }
1146 
1147 // ---------------------------------------------------------------------------
1148 
1149 //! @cond Doxygen_Suppress
1150 
1151 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr &) const1152 CRS::_identify(const io::AuthorityFactoryPtr &) const {
1153     return {};
1154 }
1155 
1156 //! @endcond
1157 
1158 // ---------------------------------------------------------------------------
1159 
1160 //! @cond Doxygen_Suppress
1161 struct SingleCRS::Private {
1162     datum::DatumPtr datum{};
1163     datum::DatumEnsemblePtr datumEnsemble{};
1164     cs::CoordinateSystemNNPtr coordinateSystem;
1165 
Privatecrs::SingleCRS::Private1166     Private(const datum::DatumPtr &datumIn,
1167             const datum::DatumEnsemblePtr &datumEnsembleIn,
1168             const cs::CoordinateSystemNNPtr &csIn)
1169         : datum(datumIn), datumEnsemble(datumEnsembleIn),
1170           coordinateSystem(csIn) {
1171         if ((datum ? 1 : 0) + (datumEnsemble ? 1 : 0) != 1) {
1172             throw util::Exception("datum or datumEnsemble should be set");
1173         }
1174     }
1175 };
1176 //! @endcond
1177 
1178 // ---------------------------------------------------------------------------
1179 
SingleCRS(const datum::DatumPtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::CoordinateSystemNNPtr & csIn)1180 SingleCRS::SingleCRS(const datum::DatumPtr &datumIn,
1181                      const datum::DatumEnsemblePtr &datumEnsembleIn,
1182                      const cs::CoordinateSystemNNPtr &csIn)
1183     : d(internal::make_unique<Private>(datumIn, datumEnsembleIn, csIn)) {}
1184 
1185 // ---------------------------------------------------------------------------
1186 
SingleCRS(const SingleCRS & other)1187 SingleCRS::SingleCRS(const SingleCRS &other)
1188     : CRS(other), d(internal::make_unique<Private>(*other.d)) {}
1189 
1190 // ---------------------------------------------------------------------------
1191 
1192 //! @cond Doxygen_Suppress
1193 SingleCRS::~SingleCRS() = default;
1194 //! @endcond
1195 
1196 // ---------------------------------------------------------------------------
1197 
1198 /** \brief Return the datum::Datum associated with the CRS.
1199  *
1200  * This might be null, in which case datumEnsemble() return will not be null.
1201  *
1202  * @return a Datum that might be null.
1203  */
datum()1204 const datum::DatumPtr &SingleCRS::datum() PROJ_PURE_DEFN { return d->datum; }
1205 
1206 // ---------------------------------------------------------------------------
1207 
1208 /** \brief Return the datum::DatumEnsemble associated with the CRS.
1209  *
1210  * This might be null, in which case datum() return will not be null.
1211  *
1212  * @return a DatumEnsemble that might be null.
1213  */
datumEnsemble()1214 const datum::DatumEnsemblePtr &SingleCRS::datumEnsemble() PROJ_PURE_DEFN {
1215     return d->datumEnsemble;
1216 }
1217 
1218 // ---------------------------------------------------------------------------
1219 
1220 //! @cond Doxygen_Suppress
1221 /** \brief Return the real datum or a synthetized one if a datumEnsemble.
1222  */
1223 const datum::DatumNNPtr
datumNonNull(const io::DatabaseContextPtr & dbContext) const1224 SingleCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
1225     return d->datum ? NN_NO_CHECK(d->datum)
1226                     : d->datumEnsemble->asDatum(dbContext);
1227 }
1228 //! @endcond
1229 
1230 // ---------------------------------------------------------------------------
1231 
1232 /** \brief Return the cs::CoordinateSystem associated with the CRS.
1233  *
1234  * @return a CoordinateSystem.
1235  */
coordinateSystem()1236 const cs::CoordinateSystemNNPtr &SingleCRS::coordinateSystem() PROJ_PURE_DEFN {
1237     return d->coordinateSystem;
1238 }
1239 
1240 // ---------------------------------------------------------------------------
1241 
baseIsEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const1242 bool SingleCRS::baseIsEquivalentTo(
1243     const util::IComparable *other, util::IComparable::Criterion criterion,
1244     const io::DatabaseContextPtr &dbContext) const {
1245     auto otherSingleCRS = dynamic_cast<const SingleCRS *>(other);
1246     if (otherSingleCRS == nullptr ||
1247         (criterion == util::IComparable::Criterion::STRICT &&
1248          !ObjectUsage::_isEquivalentTo(other, criterion, dbContext))) {
1249         return false;
1250     }
1251 
1252     if (criterion == util::IComparable::Criterion::STRICT) {
1253         const auto &thisDatum = d->datum;
1254         const auto &otherDatum = otherSingleCRS->d->datum;
1255         if (thisDatum) {
1256             if (!thisDatum->_isEquivalentTo(otherDatum.get(), criterion,
1257                                             dbContext)) {
1258                 return false;
1259             }
1260         } else {
1261             if (otherDatum) {
1262                 return false;
1263             }
1264         }
1265 
1266         const auto &thisDatumEnsemble = d->datumEnsemble;
1267         const auto &otherDatumEnsemble = otherSingleCRS->d->datumEnsemble;
1268         if (thisDatumEnsemble) {
1269             if (!thisDatumEnsemble->_isEquivalentTo(otherDatumEnsemble.get(),
1270                                                     criterion, dbContext)) {
1271                 return false;
1272             }
1273         } else {
1274             if (otherDatumEnsemble) {
1275                 return false;
1276             }
1277         }
1278     } else {
1279         if (!datumNonNull(dbContext)->_isEquivalentTo(
1280                 otherSingleCRS->datumNonNull(dbContext).get(), criterion,
1281                 dbContext)) {
1282             return false;
1283         }
1284     }
1285 
1286     return d->coordinateSystem->_isEquivalentTo(
1287                otherSingleCRS->d->coordinateSystem.get(), criterion,
1288                dbContext) &&
1289            getExtensionProj4() == otherSingleCRS->getExtensionProj4();
1290 }
1291 
1292 // ---------------------------------------------------------------------------
1293 
1294 //! @cond Doxygen_Suppress
exportDatumOrDatumEnsembleToWkt(io::WKTFormatter * formatter) const1295 void SingleCRS::exportDatumOrDatumEnsembleToWkt(
1296     io::WKTFormatter *formatter) const // throw(io::FormattingException)
1297 {
1298     const auto &l_datum = d->datum;
1299     if (l_datum) {
1300         l_datum->_exportToWKT(formatter);
1301     } else {
1302         const auto &l_datumEnsemble = d->datumEnsemble;
1303         assert(l_datumEnsemble);
1304         l_datumEnsemble->_exportToWKT(formatter);
1305     }
1306 }
1307 //! @endcond
1308 
1309 // ---------------------------------------------------------------------------
1310 
1311 //! @cond Doxygen_Suppress
1312 struct GeodeticCRS::Private {
1313     std::vector<operation::PointMotionOperationNNPtr> velocityModel{};
1314     datum::GeodeticReferenceFramePtr datum_;
1315 
Privatecrs::GeodeticCRS::Private1316     explicit Private(const datum::GeodeticReferenceFramePtr &datumIn)
1317         : datum_(datumIn) {}
1318 };
1319 
1320 // ---------------------------------------------------------------------------
1321 
1322 static const datum::DatumEnsemblePtr &
checkEnsembleForGeodeticCRS(const datum::GeodeticReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & ensemble)1323 checkEnsembleForGeodeticCRS(const datum::GeodeticReferenceFramePtr &datumIn,
1324                             const datum::DatumEnsemblePtr &ensemble) {
1325     const char *msg = "One of Datum or DatumEnsemble should be defined";
1326     if (datumIn) {
1327         if (!ensemble) {
1328             return ensemble;
1329         }
1330         msg = "Datum and DatumEnsemble should not be defined";
1331     } else if (ensemble) {
1332         const auto &datums = ensemble->datums();
1333         assert(!datums.empty());
1334         auto grfFirst =
1335             dynamic_cast<datum::GeodeticReferenceFrame *>(datums[0].get());
1336         if (grfFirst) {
1337             return ensemble;
1338         }
1339         msg = "Ensemble should contain GeodeticReferenceFrame";
1340     }
1341     throw util::Exception(msg);
1342 }
1343 
1344 //! @endcond
1345 
1346 // ---------------------------------------------------------------------------
1347 
GeodeticCRS(const datum::GeodeticReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::EllipsoidalCSNNPtr & csIn)1348 GeodeticCRS::GeodeticCRS(const datum::GeodeticReferenceFramePtr &datumIn,
1349                          const datum::DatumEnsemblePtr &datumEnsembleIn,
1350                          const cs::EllipsoidalCSNNPtr &csIn)
1351     : SingleCRS(datumIn, checkEnsembleForGeodeticCRS(datumIn, datumEnsembleIn),
1352                 csIn),
1353       d(internal::make_unique<Private>(datumIn)) {}
1354 
1355 // ---------------------------------------------------------------------------
1356 
GeodeticCRS(const datum::GeodeticReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::SphericalCSNNPtr & csIn)1357 GeodeticCRS::GeodeticCRS(const datum::GeodeticReferenceFramePtr &datumIn,
1358                          const datum::DatumEnsemblePtr &datumEnsembleIn,
1359                          const cs::SphericalCSNNPtr &csIn)
1360     : SingleCRS(datumIn, checkEnsembleForGeodeticCRS(datumIn, datumEnsembleIn),
1361                 csIn),
1362       d(internal::make_unique<Private>(datumIn)) {}
1363 
1364 // ---------------------------------------------------------------------------
1365 
GeodeticCRS(const datum::GeodeticReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::CartesianCSNNPtr & csIn)1366 GeodeticCRS::GeodeticCRS(const datum::GeodeticReferenceFramePtr &datumIn,
1367                          const datum::DatumEnsemblePtr &datumEnsembleIn,
1368                          const cs::CartesianCSNNPtr &csIn)
1369     : SingleCRS(datumIn, checkEnsembleForGeodeticCRS(datumIn, datumEnsembleIn),
1370                 csIn),
1371       d(internal::make_unique<Private>(datumIn)) {}
1372 
1373 // ---------------------------------------------------------------------------
1374 
GeodeticCRS(const GeodeticCRS & other)1375 GeodeticCRS::GeodeticCRS(const GeodeticCRS &other)
1376     : SingleCRS(other), d(internal::make_unique<Private>(*other.d)) {}
1377 
1378 // ---------------------------------------------------------------------------
1379 
1380 //! @cond Doxygen_Suppress
1381 GeodeticCRS::~GeodeticCRS() = default;
1382 //! @endcond
1383 
1384 // ---------------------------------------------------------------------------
1385 
_shallowClone() const1386 CRSNNPtr GeodeticCRS::_shallowClone() const {
1387     auto crs(GeodeticCRS::nn_make_shared<GeodeticCRS>(*this));
1388     crs->assignSelf(crs);
1389     return crs;
1390 }
1391 
1392 // ---------------------------------------------------------------------------
1393 
1394 /** \brief Return the datum::GeodeticReferenceFrame associated with the CRS.
1395  *
1396  * @return a GeodeticReferenceFrame or null (in which case datumEnsemble()
1397  * should return a non-null pointer.)
1398  */
datum()1399 const datum::GeodeticReferenceFramePtr &GeodeticCRS::datum() PROJ_PURE_DEFN {
1400     return d->datum_;
1401 }
1402 
1403 // ---------------------------------------------------------------------------
1404 
1405 //! @cond Doxygen_Suppress
1406 /** \brief Return the real datum or a synthetized one if a datumEnsemble.
1407  */
1408 const datum::GeodeticReferenceFrameNNPtr
datumNonNull(const io::DatabaseContextPtr & dbContext) const1409 GeodeticCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
1410     return NN_NO_CHECK(
1411         d->datum_
1412             ? d->datum_
1413             : util::nn_dynamic_pointer_cast<datum::GeodeticReferenceFrame>(
1414                   SingleCRS::getPrivate()->datumEnsemble->asDatum(dbContext)));
1415 }
1416 //! @endcond
1417 
1418 // ---------------------------------------------------------------------------
1419 
1420 //! @cond Doxygen_Suppress
oneDatum(const GeodeticCRS * crs)1421 static datum::GeodeticReferenceFrame *oneDatum(const GeodeticCRS *crs) {
1422     const auto &l_datumEnsemble = crs->datumEnsemble();
1423     assert(l_datumEnsemble);
1424     const auto &l_datums = l_datumEnsemble->datums();
1425     return static_cast<datum::GeodeticReferenceFrame *>(l_datums[0].get());
1426 }
1427 //! @endcond
1428 
1429 // ---------------------------------------------------------------------------
1430 
1431 /** \brief Return the PrimeMeridian associated with the GeodeticReferenceFrame
1432  * or with one of the GeodeticReferenceFrame of the datumEnsemble().
1433  *
1434  * @return the PrimeMeridian.
1435  */
primeMeridian()1436 const datum::PrimeMeridianNNPtr &GeodeticCRS::primeMeridian() PROJ_PURE_DEFN {
1437     if (d->datum_) {
1438         return d->datum_->primeMeridian();
1439     }
1440     return oneDatum(this)->primeMeridian();
1441 }
1442 
1443 // ---------------------------------------------------------------------------
1444 
1445 /** \brief Return the ellipsoid associated with the GeodeticReferenceFrame
1446  * or with one of the GeodeticReferenceFrame of the datumEnsemble().
1447  *
1448  * @return the PrimeMeridian.
1449  */
ellipsoid()1450 const datum::EllipsoidNNPtr &GeodeticCRS::ellipsoid() PROJ_PURE_DEFN {
1451     if (d->datum_) {
1452         return d->datum_->ellipsoid();
1453     }
1454     return oneDatum(this)->ellipsoid();
1455 }
1456 
1457 // ---------------------------------------------------------------------------
1458 
1459 /** \brief Return the velocity model associated with the CRS.
1460  *
1461  * @return a velocity model. might be null.
1462  */
1463 const std::vector<operation::PointMotionOperationNNPtr> &
velocityModel()1464 GeodeticCRS::velocityModel() PROJ_PURE_DEFN {
1465     return d->velocityModel;
1466 }
1467 
1468 // ---------------------------------------------------------------------------
1469 
1470 /** \brief Return whether the CRS is a geocentric one.
1471  *
1472  * A geocentric CRS is a geodetic CRS that has a Cartesian coordinate system
1473  * with three axis, whose direction is respectively
1474  * cs::AxisDirection::GEOCENTRIC_X,
1475  * cs::AxisDirection::GEOCENTRIC_Y and cs::AxisDirection::GEOCENTRIC_Z.
1476  *
1477  * @return true if the CRS is a geocentric CRS.
1478  */
isGeocentric()1479 bool GeodeticCRS::isGeocentric() PROJ_PURE_DEFN {
1480     const auto &cs = coordinateSystem();
1481     const auto &axisList = cs->axisList();
1482     return axisList.size() == 3 &&
1483            dynamic_cast<cs::CartesianCS *>(cs.get()) != nullptr &&
1484            &axisList[0]->direction() == &cs::AxisDirection::GEOCENTRIC_X &&
1485            &axisList[1]->direction() == &cs::AxisDirection::GEOCENTRIC_Y &&
1486            &axisList[2]->direction() == &cs::AxisDirection::GEOCENTRIC_Z;
1487 }
1488 
1489 // ---------------------------------------------------------------------------
1490 
1491 /** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame and a
1492  * cs::SphericalCS.
1493  *
1494  * @param properties See \ref general_properties.
1495  * At minimum the name should be defined.
1496  * @param datum The datum of the CRS.
1497  * @param cs a SphericalCS.
1498  * @return new GeodeticCRS.
1499  */
1500 GeodeticCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFrameNNPtr & datum,const cs::SphericalCSNNPtr & cs)1501 GeodeticCRS::create(const util::PropertyMap &properties,
1502                     const datum::GeodeticReferenceFrameNNPtr &datum,
1503                     const cs::SphericalCSNNPtr &cs) {
1504     return create(properties, datum.as_nullable(), nullptr, cs);
1505 }
1506 
1507 // ---------------------------------------------------------------------------
1508 
1509 /** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame or
1510  * datum::DatumEnsemble and a cs::SphericalCS.
1511  *
1512  * One and only one of datum or datumEnsemble should be set to a non-null value.
1513  *
1514  * @param properties See \ref general_properties.
1515  * At minimum the name should be defined.
1516  * @param datum The datum of the CRS, or nullptr
1517  * @param datumEnsemble The datum ensemble of the CRS, or nullptr.
1518  * @param cs a SphericalCS.
1519  * @return new GeodeticCRS.
1520  */
1521 GeodeticCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFramePtr & datum,const datum::DatumEnsemblePtr & datumEnsemble,const cs::SphericalCSNNPtr & cs)1522 GeodeticCRS::create(const util::PropertyMap &properties,
1523                     const datum::GeodeticReferenceFramePtr &datum,
1524                     const datum::DatumEnsemblePtr &datumEnsemble,
1525                     const cs::SphericalCSNNPtr &cs) {
1526     auto crs(
1527         GeodeticCRS::nn_make_shared<GeodeticCRS>(datum, datumEnsemble, cs));
1528     crs->assignSelf(crs);
1529     crs->setProperties(properties);
1530     properties.getStringValue("EXTENSION_PROJ4",
1531                               crs->CRS::getPrivate()->extensionProj4_);
1532 
1533     return crs;
1534 }
1535 
1536 // ---------------------------------------------------------------------------
1537 
1538 /** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame and a
1539  * cs::CartesianCS.
1540  *
1541  * @param properties See \ref general_properties.
1542  * At minimum the name should be defined.
1543  * @param datum The datum of the CRS.
1544  * @param cs a CartesianCS.
1545  * @return new GeodeticCRS.
1546  */
1547 GeodeticCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFrameNNPtr & datum,const cs::CartesianCSNNPtr & cs)1548 GeodeticCRS::create(const util::PropertyMap &properties,
1549                     const datum::GeodeticReferenceFrameNNPtr &datum,
1550                     const cs::CartesianCSNNPtr &cs) {
1551     return create(properties, datum.as_nullable(), nullptr, cs);
1552 }
1553 
1554 // ---------------------------------------------------------------------------
1555 
1556 /** \brief Instantiate a GeodeticCRS from a datum::GeodeticReferenceFrame or
1557  * datum::DatumEnsemble and a cs::CartesianCS.
1558  *
1559  * One and only one of datum or datumEnsemble should be set to a non-null value.
1560  *
1561  * @param properties See \ref general_properties.
1562  * At minimum the name should be defined.
1563  * @param datum The datum of the CRS, or nullptr
1564  * @param datumEnsemble The datum ensemble of the CRS, or nullptr.
1565  * @param cs a CartesianCS
1566  * @return new GeodeticCRS.
1567  */
1568 GeodeticCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFramePtr & datum,const datum::DatumEnsemblePtr & datumEnsemble,const cs::CartesianCSNNPtr & cs)1569 GeodeticCRS::create(const util::PropertyMap &properties,
1570                     const datum::GeodeticReferenceFramePtr &datum,
1571                     const datum::DatumEnsemblePtr &datumEnsemble,
1572                     const cs::CartesianCSNNPtr &cs) {
1573     auto crs(
1574         GeodeticCRS::nn_make_shared<GeodeticCRS>(datum, datumEnsemble, cs));
1575     crs->assignSelf(crs);
1576     crs->setProperties(properties);
1577     properties.getStringValue("EXTENSION_PROJ4",
1578                               crs->CRS::getPrivate()->extensionProj4_);
1579     return crs;
1580 }
1581 
1582 // ---------------------------------------------------------------------------
1583 
1584 //! @cond Doxygen_Suppress
1585 
1586 // Try to format a Geographic/ProjectedCRS 3D CRS as a
1587 // GEOGCS[]/PROJCS[],VERTCS[...,DATUM[],...] if we find corresponding objects
exportAsESRIWktCompoundCRSWithEllipsoidalHeight(const CRS * self,const GeodeticCRS * geodCRS,io::WKTFormatter * formatter)1588 static bool exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
1589     const CRS *self, const GeodeticCRS *geodCRS, io::WKTFormatter *formatter) {
1590     const auto &dbContext = formatter->databaseContext();
1591     if (!dbContext) {
1592         return false;
1593     }
1594     const auto l_datum = geodCRS->datumNonNull(formatter->databaseContext());
1595     auto l_alias = dbContext->getAliasFromOfficialName(
1596         l_datum->nameStr(), "geodetic_datum", "ESRI");
1597     if (l_alias.empty()) {
1598         return false;
1599     }
1600     auto authFactory =
1601         io::AuthorityFactory::create(NN_NO_CHECK(dbContext), std::string());
1602     auto list = authFactory->createObjectsFromName(
1603         l_alias, {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
1604         false /* approximate=false*/);
1605     if (list.empty()) {
1606         return false;
1607     }
1608     auto gdatum = util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
1609     if (gdatum == nullptr || gdatum->identifiers().empty()) {
1610         return false;
1611     }
1612     const auto &gdatum_ids = gdatum->identifiers();
1613     auto vertCRSList = authFactory->createVerticalCRSFromDatum(
1614         "ESRI", "from_geogdatum_" + *gdatum_ids[0]->codeSpace() + '_' +
1615                     gdatum_ids[0]->code());
1616     if (vertCRSList.size() != 1) {
1617         return false;
1618     }
1619     self->demoteTo2D(std::string(), dbContext)->_exportToWKT(formatter);
1620     vertCRSList.front()->_exportToWKT(formatter);
1621     return true;
1622 }
1623 
1624 // ---------------------------------------------------------------------------
1625 
1626 // Try to format a Geographic/ProjectedCRS 3D CRS as a
1627 // GEOGCS[]/PROJCS[],VERTCS["Ellipsoid (metre)",DATUM["Ellipsoid",2002],...]
exportAsWKT1CompoundCRSWithEllipsoidalHeight(const CRSNNPtr & base2DCRS,const cs::CoordinateSystemAxisNNPtr & verticalAxis,io::WKTFormatter * formatter)1628 static bool exportAsWKT1CompoundCRSWithEllipsoidalHeight(
1629     const CRSNNPtr &base2DCRS,
1630     const cs::CoordinateSystemAxisNNPtr &verticalAxis,
1631     io::WKTFormatter *formatter) {
1632     std::string verticalCRSName = "Ellipsoid (";
1633     verticalCRSName += verticalAxis->unit().name();
1634     verticalCRSName += ')';
1635     auto vertDatum = datum::VerticalReferenceFrame::create(
1636         util::PropertyMap()
1637             .set(common::IdentifiedObject::NAME_KEY, "Ellipsoid")
1638             .set("VERT_DATUM_TYPE", "2002"));
1639     auto vertCRS = VerticalCRS::create(
1640         util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
1641                                 verticalCRSName),
1642         vertDatum.as_nullable(), nullptr,
1643         cs::VerticalCS::create(util::PropertyMap(), verticalAxis));
1644     formatter->startNode(io::WKTConstants::COMPD_CS, false);
1645     formatter->addQuotedString(base2DCRS->nameStr() + " + " + verticalCRSName);
1646     base2DCRS->_exportToWKT(formatter);
1647     vertCRS->_exportToWKT(formatter);
1648     formatter->endNode();
1649     return true;
1650 }
1651 //! @endcond
1652 
1653 // ---------------------------------------------------------------------------
1654 
1655 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const1656 void GeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
1657     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
1658     const bool isGeographic =
1659         dynamic_cast<const GeographicCRS *>(this) != nullptr;
1660 
1661     const auto &cs = coordinateSystem();
1662     const auto &axisList = cs->axisList();
1663     const auto oldAxisOutputRule = formatter->outputAxis();
1664     auto l_name = nameStr();
1665     const auto &dbContext = formatter->databaseContext();
1666 
1667     if (!isWKT2 && formatter->useESRIDialect() && axisList.size() == 3) {
1668         if (!isGeographic) {
1669             io::FormattingException::Throw(
1670                 "Geocentric CRS not supported in WKT1_ESRI");
1671         }
1672         // Try to format the Geographic 3D CRS as a GEOGCS[],VERTCS[...,DATUM[]]
1673         // if we find corresponding objects
1674         if (dbContext) {
1675             if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(this, this,
1676                                                                 formatter)) {
1677                 return;
1678             }
1679         }
1680         io::FormattingException::Throw(
1681             "Cannot export this Geographic 3D CRS in WKT1_ESRI");
1682     }
1683 
1684     if (!isWKT2 && formatter->isStrict() && isGeographic &&
1685         axisList.size() != 2 &&
1686         oldAxisOutputRule != io::WKTFormatter::OutputAxisRule::NO) {
1687 
1688         auto geogCRS2D = demoteTo2D(std::string(), dbContext);
1689         if (dbContext) {
1690             const auto res = geogCRS2D->identify(
1691                 io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "EPSG"));
1692             if (res.size() == 1) {
1693                 const auto &front = res.front();
1694                 if (front.second == 100) {
1695                     geogCRS2D = front.first;
1696                 }
1697             }
1698         }
1699 
1700         if (CRS::getPrivate()->allowNonConformantWKT1Export_) {
1701             formatter->startNode(io::WKTConstants::COMPD_CS, false);
1702             formatter->addQuotedString(l_name + " + " + l_name);
1703             geogCRS2D->_exportToWKT(formatter);
1704             const auto oldTOWGSParameters = formatter->getTOWGS84Parameters();
1705             formatter->setTOWGS84Parameters({});
1706             geogCRS2D->_exportToWKT(formatter);
1707             formatter->setTOWGS84Parameters(oldTOWGSParameters);
1708             formatter->endNode();
1709             return;
1710         }
1711 
1712         auto &originalCompoundCRS = CRS::getPrivate()->originalCompoundCRS_;
1713         if (originalCompoundCRS) {
1714             originalCompoundCRS->_exportToWKT(formatter);
1715             return;
1716         }
1717 
1718         if (formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
1719             if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
1720                     geogCRS2D, axisList[2], formatter)) {
1721                 return;
1722             }
1723         }
1724 
1725         io::FormattingException::Throw(
1726             "WKT1 does not support Geographic 3D CRS.");
1727     }
1728 
1729     formatter->startNode(isWKT2
1730                              ? ((formatter->use2019Keywords() && isGeographic)
1731                                     ? io::WKTConstants::GEOGCRS
1732                                     : io::WKTConstants::GEODCRS)
1733                              : isGeocentric() ? io::WKTConstants::GEOCCS
1734                                               : io::WKTConstants::GEOGCS,
1735                          !identifiers().empty());
1736 
1737     if (formatter->useESRIDialect()) {
1738         if (l_name == "WGS 84") {
1739             l_name = "GCS_WGS_1984";
1740         } else {
1741             bool aliasFound = false;
1742             if (dbContext) {
1743                 auto l_alias = dbContext->getAliasFromOfficialName(
1744                     l_name, "geodetic_crs", "ESRI");
1745                 if (!l_alias.empty()) {
1746                     l_name = l_alias;
1747                     aliasFound = true;
1748                 }
1749             }
1750             if (!aliasFound) {
1751                 l_name = io::WKTFormatter::morphNameToESRI(l_name);
1752                 if (!starts_with(l_name, "GCS_")) {
1753                     l_name = "GCS_" + l_name;
1754                 }
1755             }
1756         }
1757     }
1758 
1759     if (!isWKT2 && !formatter->useESRIDialect() && isDeprecated()) {
1760         l_name += " (deprecated)";
1761     }
1762     formatter->addQuotedString(l_name);
1763 
1764     const auto &unit = axisList[0]->unit();
1765     formatter->pushAxisAngularUnit(common::UnitOfMeasure::create(unit));
1766     exportDatumOrDatumEnsembleToWkt(formatter);
1767     primeMeridian()->_exportToWKT(formatter);
1768     formatter->popAxisAngularUnit();
1769     if (!isWKT2) {
1770         unit._exportToWKT(formatter);
1771     }
1772 
1773     if (oldAxisOutputRule ==
1774             io::WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE &&
1775         isGeocentric()) {
1776         formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES);
1777     }
1778     cs->_exportToWKT(formatter);
1779     formatter->setOutputAxis(oldAxisOutputRule);
1780 
1781     ObjectUsage::baseExportToWKT(formatter);
1782 
1783     if (!isWKT2 && !formatter->useESRIDialect()) {
1784         const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
1785         if (!extensionProj4.empty()) {
1786             formatter->startNode(io::WKTConstants::EXTENSION, false);
1787             formatter->addQuotedString("PROJ4");
1788             formatter->addQuotedString(extensionProj4);
1789             formatter->endNode();
1790         }
1791     }
1792 
1793     formatter->endNode();
1794 }
1795 //! @endcond
1796 
1797 // ---------------------------------------------------------------------------
1798 
1799 //! @cond Doxygen_Suppress
addGeocentricUnitConversionIntoPROJString(io::PROJStringFormatter * formatter) const1800 void GeodeticCRS::addGeocentricUnitConversionIntoPROJString(
1801     io::PROJStringFormatter *formatter) const {
1802 
1803     const auto &axisList = coordinateSystem()->axisList();
1804     const auto &unit = axisList[0]->unit();
1805     if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE,
1806                               util::IComparable::Criterion::EQUIVALENT)) {
1807         if (formatter->getCRSExport()) {
1808             io::FormattingException::Throw(
1809                 "GeodeticCRS::exportToPROJString() only "
1810                 "supports metre unit");
1811         }
1812         formatter->addStep("unitconvert");
1813         formatter->addParam("xy_in", "m");
1814         formatter->addParam("z_in", "m");
1815         {
1816             auto projUnit = unit.exportToPROJString();
1817             if (!projUnit.empty()) {
1818                 formatter->addParam("xy_out", projUnit);
1819                 formatter->addParam("z_out", projUnit);
1820                 return;
1821             }
1822         }
1823 
1824         const auto &toSI = unit.conversionToSI();
1825         formatter->addParam("xy_out", toSI);
1826         formatter->addParam("z_out", toSI);
1827     } else if (formatter->getCRSExport()) {
1828         formatter->addParam("units", "m");
1829     }
1830 }
1831 //! @endcond
1832 
1833 // ---------------------------------------------------------------------------
1834 
1835 //! @cond Doxygen_Suppress
_exportToPROJString(io::PROJStringFormatter * formatter) const1836 void GeodeticCRS::_exportToPROJString(
1837     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
1838 {
1839     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
1840     if (!extensionProj4.empty()) {
1841         formatter->ingestPROJString(
1842             replaceAll(extensionProj4, " +type=crs", ""));
1843         formatter->addNoDefs(false);
1844         return;
1845     }
1846 
1847     if (!isGeocentric()) {
1848         io::FormattingException::Throw(
1849             "GeodeticCRS::exportToPROJString() only "
1850             "supports geocentric coordinate systems");
1851     }
1852 
1853     if (!formatter->getCRSExport()) {
1854         formatter->addStep("cart");
1855     } else {
1856         formatter->addStep("geocent");
1857     }
1858     addDatumInfoToPROJString(formatter);
1859     addGeocentricUnitConversionIntoPROJString(formatter);
1860 }
1861 //! @endcond
1862 
1863 // ---------------------------------------------------------------------------
1864 
1865 //! @cond Doxygen_Suppress
addDatumInfoToPROJString(io::PROJStringFormatter * formatter) const1866 void GeodeticCRS::addDatumInfoToPROJString(
1867     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
1868 {
1869     const auto &TOWGS84Params = formatter->getTOWGS84Parameters();
1870     bool datumWritten = false;
1871     const auto &nadgrids = formatter->getHDatumExtension();
1872     const auto l_datum = datumNonNull(formatter->databaseContext());
1873     if (formatter->getCRSExport() && TOWGS84Params.empty() &&
1874         nadgrids.empty()) {
1875         if (l_datum->_isEquivalentTo(
1876                 datum::GeodeticReferenceFrame::EPSG_6326.get(),
1877                 util::IComparable::Criterion::EQUIVALENT)) {
1878             datumWritten = true;
1879             formatter->addParam("datum", "WGS84");
1880         } else if (l_datum->_isEquivalentTo(
1881                        datum::GeodeticReferenceFrame::EPSG_6267.get(),
1882                        util::IComparable::Criterion::EQUIVALENT)) {
1883             datumWritten = true;
1884             formatter->addParam("datum", "NAD27");
1885         } else if (l_datum->_isEquivalentTo(
1886                        datum::GeodeticReferenceFrame::EPSG_6269.get(),
1887                        util::IComparable::Criterion::EQUIVALENT)) {
1888             datumWritten = true;
1889             if (formatter->getLegacyCRSToCRSContext()) {
1890                 // We do not want datum=NAD83 to cause a useless towgs84=0,0,0
1891                 formatter->addParam("ellps", "GRS80");
1892             } else {
1893                 formatter->addParam("datum", "NAD83");
1894             }
1895         }
1896     }
1897     if (!datumWritten) {
1898         ellipsoid()->_exportToPROJString(formatter);
1899         primeMeridian()->_exportToPROJString(formatter);
1900     }
1901     if (TOWGS84Params.size() == 7) {
1902         formatter->addParam("towgs84", TOWGS84Params);
1903     }
1904     if (!nadgrids.empty()) {
1905         formatter->addParam("nadgrids", nadgrids);
1906     }
1907 }
1908 //! @endcond
1909 
1910 // ---------------------------------------------------------------------------
1911 
1912 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const1913 void GeodeticCRS::_exportToJSON(
1914     io::JSONFormatter *formatter) const // throw(io::FormattingException)
1915 {
1916     auto writer = formatter->writer();
1917     auto objectContext(
1918         formatter->MakeObjectContext("GeodeticCRS", !identifiers().empty()));
1919 
1920     writer->AddObjKey("name");
1921     auto l_name = nameStr();
1922     if (l_name.empty()) {
1923         writer->Add("unnamed");
1924     } else {
1925         writer->Add(l_name);
1926     }
1927 
1928     const auto &l_datum(datum());
1929     if (l_datum) {
1930         writer->AddObjKey("datum");
1931         l_datum->_exportToJSON(formatter);
1932     } else {
1933         writer->AddObjKey("datum_ensemble");
1934         formatter->setOmitTypeInImmediateChild();
1935         datumEnsemble()->_exportToJSON(formatter);
1936     }
1937 
1938     writer->AddObjKey("coordinate_system");
1939     formatter->setOmitTypeInImmediateChild();
1940     coordinateSystem()->_exportToJSON(formatter);
1941 
1942     ObjectUsage::baseExportToJSON(formatter);
1943 }
1944 //! @endcond
1945 
1946 // ---------------------------------------------------------------------------
1947 
1948 //! @cond Doxygen_Suppress
1949 static util::IComparable::Criterion
getStandardCriterion(util::IComparable::Criterion criterion)1950 getStandardCriterion(util::IComparable::Criterion criterion) {
1951     return criterion == util::IComparable::Criterion::
1952                             EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
1953                ? util::IComparable::Criterion::EQUIVALENT
1954                : criterion;
1955 }
1956 //! @endcond
1957 
1958 // ---------------------------------------------------------------------------
1959 
1960 //! @cond Doxygen_Suppress
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const1961 bool GeodeticCRS::_isEquivalentTo(
1962     const util::IComparable *other, util::IComparable::Criterion criterion,
1963     const io::DatabaseContextPtr &dbContext) const {
1964     if (other == nullptr || !util::isOfExactType<GeodeticCRS>(*other)) {
1965         return false;
1966     }
1967     return _isEquivalentToNoTypeCheck(other, criterion, dbContext);
1968 }
1969 
_isEquivalentToNoTypeCheck(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const1970 bool GeodeticCRS::_isEquivalentToNoTypeCheck(
1971     const util::IComparable *other, util::IComparable::Criterion criterion,
1972     const io::DatabaseContextPtr &dbContext) const {
1973     const auto standardCriterion = getStandardCriterion(criterion);
1974 
1975     // TODO test velocityModel
1976     return SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext);
1977 }
1978 //! @endcond
1979 
1980 // ---------------------------------------------------------------------------
1981 
1982 //! @cond Doxygen_Suppress
createMapNameEPSGCode(const char * name,int code)1983 static util::PropertyMap createMapNameEPSGCode(const char *name, int code) {
1984     return util::PropertyMap()
1985         .set(common::IdentifiedObject::NAME_KEY, name)
1986         .set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::EPSG)
1987         .set(metadata::Identifier::CODE_KEY, code);
1988 }
1989 //! @endcond
1990 
1991 // ---------------------------------------------------------------------------
1992 
createEPSG_4978()1993 GeodeticCRSNNPtr GeodeticCRS::createEPSG_4978() {
1994     return create(
1995         createMapNameEPSGCode("WGS 84", 4978),
1996         datum::GeodeticReferenceFrame::EPSG_6326,
1997         cs::CartesianCS::createGeocentric(common::UnitOfMeasure::METRE));
1998 }
1999 
2000 // ---------------------------------------------------------------------------
2001 
2002 //! @cond Doxygen_Suppress
2003 
hasCodeCompatibleOfAuthorityFactory(const common::IdentifiedObject * obj,const io::AuthorityFactoryPtr & authorityFactory)2004 static bool hasCodeCompatibleOfAuthorityFactory(
2005     const common::IdentifiedObject *obj,
2006     const io::AuthorityFactoryPtr &authorityFactory) {
2007     const auto &ids = obj->identifiers();
2008     if (!ids.empty() && authorityFactory->getAuthority().empty()) {
2009         return true;
2010     }
2011     for (const auto &id : ids) {
2012         if (*(id->codeSpace()) == authorityFactory->getAuthority()) {
2013             return true;
2014         }
2015     }
2016     return false;
2017 }
2018 
hasCodeCompatibleOfAuthorityFactory(const metadata::IdentifierNNPtr & id,const io::AuthorityFactoryPtr & authorityFactory)2019 static bool hasCodeCompatibleOfAuthorityFactory(
2020     const metadata::IdentifierNNPtr &id,
2021     const io::AuthorityFactoryPtr &authorityFactory) {
2022     if (authorityFactory->getAuthority().empty()) {
2023         return true;
2024     }
2025     return *(id->codeSpace()) == authorityFactory->getAuthority();
2026 }
2027 
2028 //! @endcond
2029 
2030 // ---------------------------------------------------------------------------
2031 
2032 /** \brief Identify the CRS with reference CRSs.
2033  *
2034  * The candidate CRSs are either hard-coded, or looked in the database when
2035  * authorityFactory is not null.
2036  *
2037  * Note that the implementation uses a set of heuristics to have a good
2038  * compromise of successful identifications over execution time. It might miss
2039  * legitimate matches in some circumstances.
2040  *
2041  * The method returns a list of matching reference CRS, and the percentage
2042  * (0-100) of confidence in the match:
2043  * <ul>
2044  * <li>100% means that the name of the reference entry
2045  * perfectly matches the CRS name, and both are equivalent. In which case a
2046  * single result is returned.
2047  * Note: in the case of a GeographicCRS whose axis
2048  * order is implicit in the input definition (for example ESRI WKT), then axis
2049  * order is ignored for the purpose of identification. That is the CRS built
2050  * from
2051  * GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],
2052  * PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
2053  * will be identified to EPSG:4326, but will not pass a
2054  * isEquivalentTo(EPSG_4326, util::IComparable::Criterion::EQUIVALENT) test,
2055  * but rather isEquivalentTo(EPSG_4326,
2056  * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS)
2057  * </li>
2058  * <li>90% means that CRS are equivalent, but the names are not exactly the
2059  * same.
2060  * <li>70% means that CRS are equivalent (equivalent datum and coordinate
2061  * system),
2062  * but the names are not equivalent.</li>
2063  * <li>60% means that ellipsoid, prime meridian and coordinate systems are
2064  * equivalent, but the CRS and datum names do not match.</li>
2065  * <li>25% means that the CRS are not equivalent, but there is some similarity
2066  * in
2067  * the names.</li>
2068  * </ul>
2069  *
2070  * @param authorityFactory Authority factory (or null, but degraded
2071  * functionality)
2072  * @return a list of matching reference CRS, and the percentage (0-100) of
2073  * confidence in the match.
2074  */
2075 std::list<std::pair<GeodeticCRSNNPtr, int>>
identify(const io::AuthorityFactoryPtr & authorityFactory) const2076 GeodeticCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
2077     typedef std::pair<GeodeticCRSNNPtr, int> Pair;
2078     std::list<Pair> res;
2079     const auto &thisName(nameStr());
2080 
2081     io::DatabaseContextPtr dbContext =
2082         authorityFactory ? authorityFactory->databaseContext().as_nullable()
2083                          : nullptr;
2084     const bool l_implicitCS = hasImplicitCS();
2085     const auto crsCriterion =
2086         l_implicitCS
2087             ? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
2088             : util::IComparable::Criterion::EQUIVALENT;
2089 
2090     const GeographicCRSNNPtr candidatesCRS[] = {GeographicCRS::EPSG_4326,
2091                                                 GeographicCRS::EPSG_4267,
2092                                                 GeographicCRS::EPSG_4269};
2093     for (const auto &crs : candidatesCRS) {
2094         const bool nameEquivalent = metadata::Identifier::isEquivalentName(
2095             thisName.c_str(), crs->nameStr().c_str());
2096         const bool nameEqual = thisName == crs->nameStr();
2097         const bool isEq = _isEquivalentTo(crs.get(), crsCriterion, dbContext);
2098         if (nameEquivalent && isEq && (!authorityFactory || nameEqual)) {
2099             res.emplace_back(util::nn_static_pointer_cast<GeodeticCRS>(crs),
2100                              nameEqual ? 100 : 90);
2101             return res;
2102         } else if (nameEqual && !isEq && !authorityFactory) {
2103             res.emplace_back(util::nn_static_pointer_cast<GeodeticCRS>(crs),
2104                              25);
2105             return res;
2106         } else if (isEq && !authorityFactory) {
2107             res.emplace_back(util::nn_static_pointer_cast<GeodeticCRS>(crs),
2108                              70);
2109             return res;
2110         }
2111     }
2112 
2113     std::string geodetic_crs_type;
2114     if (isGeocentric()) {
2115         geodetic_crs_type = "geocentric";
2116     } else {
2117         auto geogCRS = dynamic_cast<const GeographicCRS *>(this);
2118         if (geogCRS) {
2119             if (coordinateSystem()->axisList().size() == 2) {
2120                 geodetic_crs_type = "geographic 2D";
2121             } else {
2122                 geodetic_crs_type = "geographic 3D";
2123             }
2124         }
2125     }
2126 
2127     if (authorityFactory) {
2128 
2129         const auto thisDatum(datumNonNull(dbContext));
2130 
2131         auto searchByDatumCode = [this, &authorityFactory, &res,
2132                                   &geodetic_crs_type, crsCriterion, &dbContext](
2133             const common::IdentifiedObjectNNPtr &l_datum) {
2134             for (const auto &id : l_datum->identifiers()) {
2135                 try {
2136                     auto tempRes = authorityFactory->createGeodeticCRSFromDatum(
2137                         *id->codeSpace(), id->code(), geodetic_crs_type);
2138                     for (const auto &crs : tempRes) {
2139                         if (_isEquivalentTo(crs.get(), crsCriterion,
2140                                             dbContext)) {
2141                             res.emplace_back(crs, 70);
2142                         }
2143                     }
2144                 } catch (const std::exception &) {
2145                 }
2146             }
2147         };
2148 
2149         auto searchByEllipsoid = [this, &authorityFactory, &res, &thisDatum,
2150                                   &geodetic_crs_type, l_implicitCS,
2151                                   &dbContext]() {
2152             const auto &thisEllipsoid = thisDatum->ellipsoid();
2153             const auto ellipsoids =
2154                 thisEllipsoid->identifiers().empty()
2155                     ? authorityFactory->createEllipsoidFromExisting(
2156                           thisEllipsoid)
2157                     : std::list<datum::EllipsoidNNPtr>{thisEllipsoid};
2158             for (const auto &ellps : ellipsoids) {
2159                 for (const auto &id : ellps->identifiers()) {
2160                     try {
2161                         auto tempRes =
2162                             authorityFactory->createGeodeticCRSFromEllipsoid(
2163                                 *id->codeSpace(), id->code(),
2164                                 geodetic_crs_type);
2165                         for (const auto &crs : tempRes) {
2166                             const auto crsDatum(crs->datumNonNull(dbContext));
2167                             if (crsDatum->ellipsoid()->_isEquivalentTo(
2168                                     ellps.get(),
2169                                     util::IComparable::Criterion::EQUIVALENT,
2170                                     dbContext) &&
2171                                 crsDatum->primeMeridian()->_isEquivalentTo(
2172                                     thisDatum->primeMeridian().get(),
2173                                     util::IComparable::Criterion::EQUIVALENT,
2174                                     dbContext) &&
2175                                 (!l_implicitCS ||
2176                                  coordinateSystem()->_isEquivalentTo(
2177                                      crs->coordinateSystem().get(),
2178                                      util::IComparable::Criterion::EQUIVALENT,
2179                                      dbContext))) {
2180                                 res.emplace_back(crs, 60);
2181                             }
2182                         }
2183                     } catch (const std::exception &) {
2184                     }
2185                 }
2186             }
2187         };
2188 
2189         const auto searchByDatumOrEllipsoid = [&authorityFactory, &res,
2190                                                &thisDatum, searchByDatumCode,
2191                                                searchByEllipsoid]() {
2192             if (!thisDatum->identifiers().empty()) {
2193                 searchByDatumCode(thisDatum);
2194             } else {
2195                 auto candidateDatums = authorityFactory->createObjectsFromName(
2196                     thisDatum->nameStr(), {io::AuthorityFactory::ObjectType::
2197                                                GEODETIC_REFERENCE_FRAME},
2198                     false);
2199                 const size_t sizeBefore = res.size();
2200                 for (const auto &candidateDatum : candidateDatums) {
2201                     searchByDatumCode(candidateDatum);
2202                 }
2203                 if (sizeBefore == res.size()) {
2204                     searchByEllipsoid();
2205                 }
2206             }
2207         };
2208 
2209         const bool unsignificantName = thisName.empty() ||
2210                                        ci_equal(thisName, "unknown") ||
2211                                        ci_equal(thisName, "unnamed");
2212 
2213         if (unsignificantName) {
2214             searchByDatumOrEllipsoid();
2215         } else if (hasCodeCompatibleOfAuthorityFactory(this,
2216                                                        authorityFactory)) {
2217             // If the CRS has already an id, check in the database for the
2218             // official object, and verify that they are equivalent.
2219             for (const auto &id : identifiers()) {
2220                 if (hasCodeCompatibleOfAuthorityFactory(id, authorityFactory)) {
2221                     try {
2222                         auto crs = io::AuthorityFactory::create(
2223                                        authorityFactory->databaseContext(),
2224                                        *id->codeSpace())
2225                                        ->createGeodeticCRS(id->code());
2226                         bool match =
2227                             _isEquivalentTo(crs.get(), crsCriterion, dbContext);
2228                         res.emplace_back(crs, match ? 100 : 25);
2229                         return res;
2230                     } catch (const std::exception &) {
2231                     }
2232                 }
2233             }
2234         } else {
2235             bool gotAbove25Pct = false;
2236             for (int ipass = 0; ipass < 2; ipass++) {
2237                 const bool approximateMatch = ipass == 1;
2238                 auto objects = authorityFactory->createObjectsFromName(
2239                     thisName, {io::AuthorityFactory::ObjectType::GEODETIC_CRS},
2240                     approximateMatch);
2241                 for (const auto &obj : objects) {
2242                     auto crs = util::nn_dynamic_pointer_cast<GeodeticCRS>(obj);
2243                     assert(crs);
2244                     auto crsNN = NN_NO_CHECK(crs);
2245                     if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) {
2246                         if (crs->nameStr() == thisName) {
2247                             res.clear();
2248                             res.emplace_back(crsNN, 100);
2249                             return res;
2250                         }
2251                         const bool eqName =
2252                             metadata::Identifier::isEquivalentName(
2253                                 thisName.c_str(), crs->nameStr().c_str());
2254                         res.emplace_back(crsNN, eqName ? 90 : 70);
2255                         gotAbove25Pct = true;
2256                     } else {
2257                         res.emplace_back(crsNN, 25);
2258                     }
2259                 }
2260                 if (!res.empty()) {
2261                     break;
2262                 }
2263             }
2264             if (!gotAbove25Pct) {
2265                 searchByDatumOrEllipsoid();
2266             }
2267         }
2268 
2269         const auto &thisCS(coordinateSystem());
2270         // Sort results
2271         res.sort([&thisName, &thisDatum, &thisCS, &dbContext](const Pair &a,
2272                                                               const Pair &b) {
2273             // First consider confidence
2274             if (a.second > b.second) {
2275                 return true;
2276             }
2277             if (a.second < b.second) {
2278                 return false;
2279             }
2280 
2281             // Then consider exact name matching
2282             const auto &aName(a.first->nameStr());
2283             const auto &bName(b.first->nameStr());
2284             if (aName == thisName && bName != thisName) {
2285                 return true;
2286             }
2287             if (bName == thisName && aName != thisName) {
2288                 return false;
2289             }
2290 
2291             // Then datum matching
2292             const auto aDatum(a.first->datumNonNull(dbContext));
2293             const auto bDatum(b.first->datumNonNull(dbContext));
2294             const auto thisEquivADatum(thisDatum->_isEquivalentTo(
2295                 aDatum.get(), util::IComparable::Criterion::EQUIVALENT,
2296                 dbContext));
2297             const auto thisEquivBDatum(thisDatum->_isEquivalentTo(
2298                 bDatum.get(), util::IComparable::Criterion::EQUIVALENT,
2299                 dbContext));
2300 
2301             if (thisEquivADatum && !thisEquivBDatum) {
2302                 return true;
2303             }
2304             if (!thisEquivADatum && thisEquivBDatum) {
2305                 return false;
2306             }
2307 
2308             // Then coordinate system matching
2309             const auto &aCS(a.first->coordinateSystem());
2310             const auto &bCS(b.first->coordinateSystem());
2311             const auto thisEquivACs(thisCS->_isEquivalentTo(
2312                 aCS.get(), util::IComparable::Criterion::EQUIVALENT,
2313                 dbContext));
2314             const auto thisEquivBCs(thisCS->_isEquivalentTo(
2315                 bCS.get(), util::IComparable::Criterion::EQUIVALENT,
2316                 dbContext));
2317             if (thisEquivACs && !thisEquivBCs) {
2318                 return true;
2319             }
2320             if (!thisEquivACs && thisEquivBCs) {
2321                 return false;
2322             }
2323 
2324             // Then dimension of the coordinate system matching
2325             const auto thisCSAxisListSize = thisCS->axisList().size();
2326             const auto aCSAxistListSize = aCS->axisList().size();
2327             const auto bCSAxistListSize = bCS->axisList().size();
2328             if (thisCSAxisListSize == aCSAxistListSize &&
2329                 thisCSAxisListSize != bCSAxistListSize) {
2330                 return true;
2331             }
2332             if (thisCSAxisListSize != aCSAxistListSize &&
2333                 thisCSAxisListSize == bCSAxistListSize) {
2334                 return false;
2335             }
2336 
2337             // Favor the CRS whole ellipsoid names matches the ellipsoid
2338             // name (WGS84...)
2339             const bool aEllpsNameEqCRSName =
2340                 metadata::Identifier::isEquivalentName(
2341                     aDatum->ellipsoid()->nameStr().c_str(),
2342                     a.first->nameStr().c_str());
2343             const bool bEllpsNameEqCRSName =
2344                 metadata::Identifier::isEquivalentName(
2345                     bDatum->ellipsoid()->nameStr().c_str(),
2346                     b.first->nameStr().c_str());
2347             if (aEllpsNameEqCRSName && !bEllpsNameEqCRSName) {
2348                 return true;
2349             }
2350             if (bEllpsNameEqCRSName && !aEllpsNameEqCRSName) {
2351                 return false;
2352             }
2353 
2354             // Arbitrary final sorting criterion
2355             return aName < bName;
2356         });
2357 
2358         // If there are results with 90% confidence, only keep those
2359         if (res.size() >= 2 && res.front().second == 90) {
2360             std::list<Pair> newRes;
2361             for (const auto &pair : res) {
2362                 if (pair.second == 90) {
2363                     newRes.push_back(pair);
2364                 } else {
2365                     break;
2366                 }
2367             }
2368             return newRes;
2369         }
2370     }
2371     return res;
2372 }
2373 
2374 // ---------------------------------------------------------------------------
2375 
2376 //! @cond Doxygen_Suppress
2377 
2378 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & authorityFactory) const2379 GeodeticCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
2380     typedef std::pair<CRSNNPtr, int> Pair;
2381     std::list<Pair> res;
2382     auto resTemp = identify(authorityFactory);
2383     for (const auto &pair : resTemp) {
2384         res.emplace_back(pair.first, pair.second);
2385     }
2386     return res;
2387 }
2388 
2389 //! @endcond
2390 
2391 // ---------------------------------------------------------------------------
2392 
2393 //! @cond Doxygen_Suppress
2394 struct GeographicCRS::Private {
2395     cs::EllipsoidalCSNNPtr coordinateSystem_;
2396 
Privatecrs::GeographicCRS::Private2397     explicit Private(const cs::EllipsoidalCSNNPtr &csIn)
2398         : coordinateSystem_(csIn) {}
2399 };
2400 //! @endcond
2401 
2402 // ---------------------------------------------------------------------------
2403 
GeographicCRS(const datum::GeodeticReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::EllipsoidalCSNNPtr & csIn)2404 GeographicCRS::GeographicCRS(const datum::GeodeticReferenceFramePtr &datumIn,
2405                              const datum::DatumEnsemblePtr &datumEnsembleIn,
2406                              const cs::EllipsoidalCSNNPtr &csIn)
2407     : SingleCRS(datumIn, datumEnsembleIn, csIn),
2408       GeodeticCRS(datumIn,
2409                   checkEnsembleForGeodeticCRS(datumIn, datumEnsembleIn), csIn),
2410       d(internal::make_unique<Private>(csIn)) {}
2411 
2412 // ---------------------------------------------------------------------------
2413 
GeographicCRS(const GeographicCRS & other)2414 GeographicCRS::GeographicCRS(const GeographicCRS &other)
2415     : SingleCRS(other), GeodeticCRS(other),
2416       d(internal::make_unique<Private>(*other.d)) {}
2417 
2418 // ---------------------------------------------------------------------------
2419 
2420 //! @cond Doxygen_Suppress
2421 GeographicCRS::~GeographicCRS() = default;
2422 //! @endcond
2423 
2424 // ---------------------------------------------------------------------------
2425 
_shallowClone() const2426 CRSNNPtr GeographicCRS::_shallowClone() const {
2427     auto crs(GeographicCRS::nn_make_shared<GeographicCRS>(*this));
2428     crs->assignSelf(crs);
2429     return crs;
2430 }
2431 
2432 // ---------------------------------------------------------------------------
2433 
2434 /** \brief Return the cs::EllipsoidalCS associated with the CRS.
2435  *
2436  * @return a EllipsoidalCS.
2437  */
coordinateSystem()2438 const cs::EllipsoidalCSNNPtr &GeographicCRS::coordinateSystem() PROJ_PURE_DEFN {
2439     return d->coordinateSystem_;
2440 }
2441 
2442 // ---------------------------------------------------------------------------
2443 
2444 /** \brief Instantiate a GeographicCRS from a datum::GeodeticReferenceFrameNNPtr
2445  * and a
2446  * cs::EllipsoidalCS.
2447  *
2448  * @param properties See \ref general_properties.
2449  * At minimum the name should be defined.
2450  * @param datum The datum of the CRS.
2451  * @param cs a EllipsoidalCS.
2452  * @return new GeographicCRS.
2453  */
2454 GeographicCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFrameNNPtr & datum,const cs::EllipsoidalCSNNPtr & cs)2455 GeographicCRS::create(const util::PropertyMap &properties,
2456                       const datum::GeodeticReferenceFrameNNPtr &datum,
2457                       const cs::EllipsoidalCSNNPtr &cs) {
2458     return create(properties, datum.as_nullable(), nullptr, cs);
2459 }
2460 
2461 // ---------------------------------------------------------------------------
2462 
2463 /** \brief Instantiate a GeographicCRS from a datum::GeodeticReferenceFramePtr
2464  * or
2465  * datum::DatumEnsemble and a
2466  * cs::EllipsoidalCS.
2467  *
2468  * One and only one of datum or datumEnsemble should be set to a non-null value.
2469  *
2470  * @param properties See \ref general_properties.
2471  * At minimum the name should be defined.
2472  * @param datum The datum of the CRS, or nullptr
2473  * @param datumEnsemble The datum ensemble of the CRS, or nullptr.
2474  * @param cs a EllipsoidalCS.
2475  * @return new GeographicCRS.
2476  */
2477 GeographicCRSNNPtr
create(const util::PropertyMap & properties,const datum::GeodeticReferenceFramePtr & datum,const datum::DatumEnsemblePtr & datumEnsemble,const cs::EllipsoidalCSNNPtr & cs)2478 GeographicCRS::create(const util::PropertyMap &properties,
2479                       const datum::GeodeticReferenceFramePtr &datum,
2480                       const datum::DatumEnsemblePtr &datumEnsemble,
2481                       const cs::EllipsoidalCSNNPtr &cs) {
2482     GeographicCRSNNPtr crs(
2483         GeographicCRS::nn_make_shared<GeographicCRS>(datum, datumEnsemble, cs));
2484     crs->assignSelf(crs);
2485     crs->setProperties(properties);
2486     properties.getStringValue("EXTENSION_PROJ4",
2487                               crs->CRS::getPrivate()->extensionProj4_);
2488     crs->CRS::getPrivate()->setImplicitCS(properties);
2489     return crs;
2490 }
2491 
2492 // ---------------------------------------------------------------------------
2493 
2494 //! @cond Doxygen_Suppress
2495 
2496 /** \brief Return whether the current GeographicCRS is the 2D part of the
2497  * other 3D GeographicCRS.
2498  */
is2DPartOf3D(util::nn<const GeographicCRS * > other,const io::DatabaseContextPtr & dbContext)2499 bool GeographicCRS::is2DPartOf3D(util::nn<const GeographicCRS *> other,
2500                                  const io::DatabaseContextPtr &dbContext)
2501     PROJ_PURE_DEFN {
2502     const auto &axis = d->coordinateSystem_->axisList();
2503     const auto &otherAxis = other->d->coordinateSystem_->axisList();
2504     if (!(axis.size() == 2 && otherAxis.size() == 3)) {
2505         return false;
2506     }
2507     const auto &firstAxis = axis[0];
2508     const auto &secondAxis = axis[1];
2509     const auto &otherFirstAxis = otherAxis[0];
2510     const auto &otherSecondAxis = otherAxis[1];
2511     if (!(firstAxis->_isEquivalentTo(
2512               otherFirstAxis.get(), util::IComparable::Criterion::EQUIVALENT) &&
2513           secondAxis->_isEquivalentTo(
2514               otherSecondAxis.get(),
2515               util::IComparable::Criterion::EQUIVALENT))) {
2516         return false;
2517     }
2518     const auto thisDatum = datumNonNull(dbContext);
2519     const auto otherDatum = other->datumNonNull(dbContext);
2520     return thisDatum->_isEquivalentTo(otherDatum.get(),
2521                                       util::IComparable::Criterion::EQUIVALENT);
2522 }
2523 
2524 //! @endcond
2525 
2526 // ---------------------------------------------------------------------------
2527 
2528 //! @cond Doxygen_Suppress
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const2529 bool GeographicCRS::_isEquivalentTo(
2530     const util::IComparable *other, util::IComparable::Criterion criterion,
2531     const io::DatabaseContextPtr &dbContext) const {
2532     if (other == nullptr || !util::isOfExactType<GeographicCRS>(*other)) {
2533         return false;
2534     }
2535 
2536     const auto standardCriterion = getStandardCriterion(criterion);
2537     if (GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
2538                                                 dbContext)) {
2539         return true;
2540     }
2541     if (criterion !=
2542         util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS) {
2543         return false;
2544     }
2545     const auto axisOrder = coordinateSystem()->axisOrder();
2546     if (axisOrder == cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH ||
2547         axisOrder == cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST) {
2548         const auto &unit = coordinateSystem()->axisList()[0]->unit();
2549         return GeographicCRS::create(
2550                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
2551                                            nameStr()),
2552                    datum(), datumEnsemble(),
2553                    axisOrder ==
2554                            cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH
2555                        ? cs::EllipsoidalCS::createLatitudeLongitude(unit)
2556                        : cs::EllipsoidalCS::createLongitudeLatitude(unit))
2557             ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
2558                                                       dbContext);
2559     }
2560     if (axisOrder ==
2561             cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH_HEIGHT_UP ||
2562         axisOrder ==
2563             cs::EllipsoidalCS::AxisOrder::LAT_NORTH_LONG_EAST_HEIGHT_UP) {
2564         const auto &angularUnit = coordinateSystem()->axisList()[0]->unit();
2565         const auto &linearUnit = coordinateSystem()->axisList()[2]->unit();
2566         return GeographicCRS::create(
2567                    util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
2568                                            nameStr()),
2569                    datum(), datumEnsemble(),
2570                    axisOrder == cs::EllipsoidalCS::AxisOrder::
2571                                     LONG_EAST_LAT_NORTH_HEIGHT_UP
2572                        ? cs::EllipsoidalCS::
2573                              createLatitudeLongitudeEllipsoidalHeight(
2574                                  angularUnit, linearUnit)
2575                        : cs::EllipsoidalCS::
2576                              createLongitudeLatitudeEllipsoidalHeight(
2577                                  angularUnit, linearUnit))
2578             ->GeodeticCRS::_isEquivalentToNoTypeCheck(other, standardCriterion,
2579                                                       dbContext);
2580     }
2581     return false;
2582 }
2583 //! @endcond
2584 
2585 // ---------------------------------------------------------------------------
2586 
createEPSG_4267()2587 GeographicCRSNNPtr GeographicCRS::createEPSG_4267() {
2588     return create(createMapNameEPSGCode("NAD27", 4267),
2589                   datum::GeodeticReferenceFrame::EPSG_6267,
2590                   cs::EllipsoidalCS::createLatitudeLongitude(
2591                       common::UnitOfMeasure::DEGREE));
2592 }
2593 
2594 // ---------------------------------------------------------------------------
2595 
createEPSG_4269()2596 GeographicCRSNNPtr GeographicCRS::createEPSG_4269() {
2597     return create(createMapNameEPSGCode("NAD83", 4269),
2598                   datum::GeodeticReferenceFrame::EPSG_6269,
2599                   cs::EllipsoidalCS::createLatitudeLongitude(
2600                       common::UnitOfMeasure::DEGREE));
2601 }
2602 
2603 // ---------------------------------------------------------------------------
2604 
createEPSG_4326()2605 GeographicCRSNNPtr GeographicCRS::createEPSG_4326() {
2606     return create(createMapNameEPSGCode("WGS 84", 4326),
2607                   datum::GeodeticReferenceFrame::EPSG_6326,
2608                   cs::EllipsoidalCS::createLatitudeLongitude(
2609                       common::UnitOfMeasure::DEGREE));
2610 }
2611 
2612 // ---------------------------------------------------------------------------
2613 
createOGC_CRS84()2614 GeographicCRSNNPtr GeographicCRS::createOGC_CRS84() {
2615     util::PropertyMap propertiesCRS;
2616     propertiesCRS
2617         .set(metadata::Identifier::CODESPACE_KEY, metadata::Identifier::OGC)
2618         .set(metadata::Identifier::CODE_KEY, "CRS84")
2619         .set(common::IdentifiedObject::NAME_KEY, "WGS 84 (CRS84)");
2620     return create(propertiesCRS, datum::GeodeticReferenceFrame::EPSG_6326,
2621                   cs::EllipsoidalCS::createLongitudeLatitude( // Long Lat !
2622                       common::UnitOfMeasure::DEGREE));
2623 }
2624 
2625 // ---------------------------------------------------------------------------
2626 
createEPSG_4979()2627 GeographicCRSNNPtr GeographicCRS::createEPSG_4979() {
2628     return create(
2629         createMapNameEPSGCode("WGS 84", 4979),
2630         datum::GeodeticReferenceFrame::EPSG_6326,
2631         cs::EllipsoidalCS::createLatitudeLongitudeEllipsoidalHeight(
2632             common::UnitOfMeasure::DEGREE, common::UnitOfMeasure::METRE));
2633 }
2634 
2635 // ---------------------------------------------------------------------------
2636 
createEPSG_4807()2637 GeographicCRSNNPtr GeographicCRS::createEPSG_4807() {
2638     auto ellps(datum::Ellipsoid::createFlattenedSphere(
2639         createMapNameEPSGCode("Clarke 1880 (IGN)", 7011),
2640         common::Length(6378249.2), common::Scale(293.4660212936269)));
2641 
2642     auto cs(cs::EllipsoidalCS::createLatitudeLongitude(
2643         common::UnitOfMeasure::GRAD));
2644 
2645     auto datum(datum::GeodeticReferenceFrame::create(
2646         createMapNameEPSGCode("Nouvelle Triangulation Francaise (Paris)", 6807),
2647         ellps, util::optional<std::string>(), datum::PrimeMeridian::PARIS));
2648 
2649     return create(createMapNameEPSGCode("NTF (Paris)", 4807), datum, cs);
2650 }
2651 
2652 // ---------------------------------------------------------------------------
2653 
2654 /** \brief Return a variant of this CRS "demoted" to a 2D one, if not already
2655  * the case.
2656  *
2657  *
2658  * @param newName Name of the new CRS. If empty, nameStr() will be used.
2659  * @param dbContext Database context to look for potentially already registered
2660  *                  2D CRS. May be nullptr.
2661  * @return a new CRS demoted to 2D, or the current one if already 2D or not
2662  * applicable.
2663  * @since 6.3
2664  */
2665 GeographicCRSNNPtr
demoteTo2D(const std::string & newName,const io::DatabaseContextPtr & dbContext) const2666 GeographicCRS::demoteTo2D(const std::string &newName,
2667                           const io::DatabaseContextPtr &dbContext) const {
2668 
2669     const auto &axisList = coordinateSystem()->axisList();
2670     if (axisList.size() == 3) {
2671         const auto &l_identifiers = identifiers();
2672         // First check if there is a Geographic 2D CRS in the database
2673         // of the same name.
2674         // This is the common practice in the EPSG dataset.
2675         if (dbContext && l_identifiers.size() == 1) {
2676             auto authFactory = io::AuthorityFactory::create(
2677                 NN_NO_CHECK(dbContext), *(l_identifiers[0]->codeSpace()));
2678             auto res = authFactory->createObjectsFromName(
2679                 nameStr(),
2680                 {io::AuthorityFactory::ObjectType::GEOGRAPHIC_2D_CRS}, false);
2681             if (!res.empty()) {
2682                 const auto &firstRes = res.front();
2683                 auto firstResAsGeogCRS =
2684                     util::nn_dynamic_pointer_cast<GeographicCRS>(firstRes);
2685                 if (firstResAsGeogCRS &&
2686                     firstResAsGeogCRS->is2DPartOf3D(NN_NO_CHECK(this),
2687                                                     dbContext)) {
2688                     return NN_NO_CHECK(firstResAsGeogCRS);
2689                 }
2690             }
2691         }
2692 
2693         auto cs = cs::EllipsoidalCS::create(util::PropertyMap(), axisList[0],
2694                                             axisList[1]);
2695         return GeographicCRS::create(
2696             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
2697                                     !newName.empty() ? newName : nameStr()),
2698             datum(), datumEnsemble(), cs);
2699     }
2700 
2701     return NN_NO_CHECK(std::dynamic_pointer_cast<GeographicCRS>(
2702         shared_from_this().as_nullable()));
2703 }
2704 
2705 // ---------------------------------------------------------------------------
2706 
2707 //! @cond Doxygen_Suppress
addAngularUnitConvertAndAxisSwap(io::PROJStringFormatter * formatter) const2708 void GeographicCRS::addAngularUnitConvertAndAxisSwap(
2709     io::PROJStringFormatter *formatter) const {
2710     const auto &axisList = coordinateSystem()->axisList();
2711 
2712     formatter->addStep("unitconvert");
2713     formatter->addParam("xy_in", "rad");
2714     if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
2715         formatter->addParam("z_in", "m");
2716     }
2717     {
2718         const auto &unitHoriz = axisList[0]->unit();
2719         const auto projUnit = unitHoriz.exportToPROJString();
2720         if (projUnit.empty()) {
2721             formatter->addParam("xy_out", unitHoriz.conversionToSI());
2722         } else {
2723             formatter->addParam("xy_out", projUnit);
2724         }
2725     }
2726     if (axisList.size() == 3 && !formatter->omitZUnitConversion()) {
2727         const auto &unitZ = axisList[2]->unit();
2728         auto projVUnit = unitZ.exportToPROJString();
2729         if (projVUnit.empty()) {
2730             formatter->addParam("z_out", unitZ.conversionToSI());
2731         } else {
2732             formatter->addParam("z_out", projVUnit);
2733         }
2734     }
2735 
2736     const char *order[2] = {nullptr, nullptr};
2737     const char *one = "1";
2738     const char *two = "2";
2739     for (int i = 0; i < 2; i++) {
2740         const auto &dir = axisList[i]->direction();
2741         if (&dir == &cs::AxisDirection::WEST) {
2742             order[i] = "-1";
2743         } else if (&dir == &cs::AxisDirection::EAST) {
2744             order[i] = one;
2745         } else if (&dir == &cs::AxisDirection::SOUTH) {
2746             order[i] = "-2";
2747         } else if (&dir == &cs::AxisDirection::NORTH) {
2748             order[i] = two;
2749         }
2750     }
2751     if (order[0] && order[1] && (order[0] != one || order[1] != two)) {
2752         formatter->addStep("axisswap");
2753         char orderStr[10];
2754         sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
2755         formatter->addParam("order", orderStr);
2756     }
2757 }
2758 //! @endcond
2759 
2760 // ---------------------------------------------------------------------------
2761 
2762 //! @cond Doxygen_Suppress
_exportToPROJString(io::PROJStringFormatter * formatter) const2763 void GeographicCRS::_exportToPROJString(
2764     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
2765 {
2766     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
2767     if (!extensionProj4.empty()) {
2768         formatter->ingestPROJString(
2769             replaceAll(extensionProj4, " +type=crs", ""));
2770         formatter->addNoDefs(false);
2771         return;
2772     }
2773 
2774     if (!formatter->omitProjLongLatIfPossible() ||
2775         primeMeridian()->longitude().getSIValue() != 0.0 ||
2776         !formatter->getTOWGS84Parameters().empty() ||
2777         !formatter->getHDatumExtension().empty()) {
2778 
2779         formatter->addStep("longlat");
2780         bool done = false;
2781         if (formatter->getLegacyCRSToCRSContext() &&
2782             formatter->getHDatumExtension().empty() &&
2783             formatter->getTOWGS84Parameters().empty()) {
2784             const auto l_datum = datumNonNull(formatter->databaseContext());
2785             if (l_datum->_isEquivalentTo(
2786                     datum::GeodeticReferenceFrame::EPSG_6326.get(),
2787                     util::IComparable::Criterion::EQUIVALENT)) {
2788                 done = true;
2789                 formatter->addParam("ellps", "WGS84");
2790             } else if (l_datum->_isEquivalentTo(
2791                            datum::GeodeticReferenceFrame::EPSG_6269.get(),
2792                            util::IComparable::Criterion::EQUIVALENT)) {
2793                 done = true;
2794                 // We do not want datum=NAD83 to cause a useless towgs84=0,0,0
2795                 formatter->addParam("ellps", "GRS80");
2796             }
2797         }
2798         if (!done) {
2799             addDatumInfoToPROJString(formatter);
2800         }
2801     }
2802     if (!formatter->getCRSExport()) {
2803         addAngularUnitConvertAndAxisSwap(formatter);
2804     }
2805 }
2806 //! @endcond
2807 
2808 // ---------------------------------------------------------------------------
2809 
2810 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const2811 void GeographicCRS::_exportToJSON(
2812     io::JSONFormatter *formatter) const // throw(io::FormattingException)
2813 {
2814     auto writer = formatter->writer();
2815     auto objectContext(
2816         formatter->MakeObjectContext("GeographicCRS", !identifiers().empty()));
2817 
2818     writer->AddObjKey("name");
2819     auto l_name = nameStr();
2820     if (l_name.empty()) {
2821         writer->Add("unnamed");
2822     } else {
2823         writer->Add(l_name);
2824     }
2825 
2826     const auto &l_datum(datum());
2827     if (l_datum) {
2828         writer->AddObjKey("datum");
2829         l_datum->_exportToJSON(formatter);
2830     } else {
2831         writer->AddObjKey("datum_ensemble");
2832         formatter->setOmitTypeInImmediateChild();
2833         datumEnsemble()->_exportToJSON(formatter);
2834     }
2835 
2836     writer->AddObjKey("coordinate_system");
2837     formatter->setOmitTypeInImmediateChild();
2838     coordinateSystem()->_exportToJSON(formatter);
2839 
2840     ObjectUsage::baseExportToJSON(formatter);
2841 }
2842 //! @endcond
2843 
2844 // ---------------------------------------------------------------------------
2845 
2846 //! @cond Doxygen_Suppress
2847 struct VerticalCRS::Private {
2848     std::vector<operation::TransformationNNPtr> geoidModel{};
2849     std::vector<operation::PointMotionOperationNNPtr> velocityModel{};
2850 };
2851 //! @endcond
2852 
2853 // ---------------------------------------------------------------------------
2854 
2855 //! @cond Doxygen_Suppress
2856 static const datum::DatumEnsemblePtr &
checkEnsembleForVerticalCRS(const datum::VerticalReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & ensemble)2857 checkEnsembleForVerticalCRS(const datum::VerticalReferenceFramePtr &datumIn,
2858                             const datum::DatumEnsemblePtr &ensemble) {
2859     const char *msg = "One of Datum or DatumEnsemble should be defined";
2860     if (datumIn) {
2861         if (!ensemble) {
2862             return ensemble;
2863         }
2864         msg = "Datum and DatumEnsemble should not be defined";
2865     } else if (ensemble) {
2866         const auto &datums = ensemble->datums();
2867         assert(!datums.empty());
2868         auto grfFirst =
2869             dynamic_cast<datum::VerticalReferenceFrame *>(datums[0].get());
2870         if (grfFirst) {
2871             return ensemble;
2872         }
2873         msg = "Ensemble should contain VerticalReferenceFrame";
2874     }
2875     throw util::Exception(msg);
2876 }
2877 //! @endcond
2878 
2879 // ---------------------------------------------------------------------------
2880 
VerticalCRS(const datum::VerticalReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::VerticalCSNNPtr & csIn)2881 VerticalCRS::VerticalCRS(const datum::VerticalReferenceFramePtr &datumIn,
2882                          const datum::DatumEnsemblePtr &datumEnsembleIn,
2883                          const cs::VerticalCSNNPtr &csIn)
2884     : SingleCRS(datumIn, checkEnsembleForVerticalCRS(datumIn, datumEnsembleIn),
2885                 csIn),
2886       d(internal::make_unique<Private>()) {}
2887 
2888 // ---------------------------------------------------------------------------
2889 
VerticalCRS(const VerticalCRS & other)2890 VerticalCRS::VerticalCRS(const VerticalCRS &other)
2891     : SingleCRS(other), d(internal::make_unique<Private>(*other.d)) {}
2892 
2893 // ---------------------------------------------------------------------------
2894 
2895 //! @cond Doxygen_Suppress
2896 VerticalCRS::~VerticalCRS() = default;
2897 //! @endcond
2898 
2899 // ---------------------------------------------------------------------------
2900 
_shallowClone() const2901 CRSNNPtr VerticalCRS::_shallowClone() const {
2902     auto crs(VerticalCRS::nn_make_shared<VerticalCRS>(*this));
2903     crs->assignSelf(crs);
2904     return crs;
2905 }
2906 
2907 // ---------------------------------------------------------------------------
2908 
2909 /** \brief Return the datum::VerticalReferenceFrame associated with the CRS.
2910  *
2911  * @return a VerticalReferenceFrame.
2912  */
datum() const2913 const datum::VerticalReferenceFramePtr VerticalCRS::datum() const {
2914     return std::static_pointer_cast<datum::VerticalReferenceFrame>(
2915         SingleCRS::getPrivate()->datum);
2916 }
2917 
2918 // ---------------------------------------------------------------------------
2919 
2920 /** \brief Return the geoid model associated with the CRS.
2921  *
2922  * Geoid height model or height correction model linked to a geoid-based
2923  * vertical CRS.
2924  *
2925  * @return a geoid model. might be null
2926  */
2927 const std::vector<operation::TransformationNNPtr> &
geoidModel()2928 VerticalCRS::geoidModel() PROJ_PURE_DEFN {
2929     return d->geoidModel;
2930 }
2931 
2932 // ---------------------------------------------------------------------------
2933 
2934 /** \brief Return the velocity model associated with the CRS.
2935  *
2936  * @return a velocity model. might be null.
2937  */
2938 const std::vector<operation::PointMotionOperationNNPtr> &
velocityModel()2939 VerticalCRS::velocityModel() PROJ_PURE_DEFN {
2940     return d->velocityModel;
2941 }
2942 
2943 // ---------------------------------------------------------------------------
2944 
2945 /** \brief Return the cs::VerticalCS associated with the CRS.
2946  *
2947  * @return a VerticalCS.
2948  */
coordinateSystem() const2949 const cs::VerticalCSNNPtr VerticalCRS::coordinateSystem() const {
2950     return util::nn_static_pointer_cast<cs::VerticalCS>(
2951         SingleCRS::getPrivate()->coordinateSystem);
2952 }
2953 
2954 // ---------------------------------------------------------------------------
2955 
2956 //! @cond Doxygen_Suppress
2957 /** \brief Return the real datum or a synthetized one if a datumEnsemble.
2958  */
2959 const datum::VerticalReferenceFrameNNPtr
datumNonNull(const io::DatabaseContextPtr & dbContext) const2960 VerticalCRS::datumNonNull(const io::DatabaseContextPtr &dbContext) const {
2961     return NN_NO_CHECK(
2962         util::nn_dynamic_pointer_cast<datum::VerticalReferenceFrame>(
2963             SingleCRS::datumNonNull(dbContext)));
2964 }
2965 //! @endcond
2966 
2967 // ---------------------------------------------------------------------------
2968 
2969 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const2970 void VerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
2971     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
2972     formatter->startNode(isWKT2 ? io::WKTConstants::VERTCRS
2973                                 : formatter->useESRIDialect()
2974                                       ? io::WKTConstants::VERTCS
2975                                       : io::WKTConstants::VERT_CS,
2976                          !identifiers().empty());
2977 
2978     auto l_name = nameStr();
2979     const auto &dbContext = formatter->databaseContext();
2980     if (formatter->useESRIDialect()) {
2981         bool aliasFound = false;
2982         if (dbContext) {
2983             auto l_alias = dbContext->getAliasFromOfficialName(
2984                 l_name, "vertical_crs", "ESRI");
2985             if (!l_alias.empty()) {
2986                 l_name = l_alias;
2987                 aliasFound = true;
2988             }
2989         }
2990         if (!aliasFound) {
2991             l_name = io::WKTFormatter::morphNameToESRI(l_name);
2992         }
2993     }
2994 
2995     formatter->addQuotedString(l_name);
2996 
2997     const auto l_datum = datum();
2998     if (formatter->useESRIDialect() && l_datum &&
2999         l_datum->getWKT1DatumType() == "2002") {
3000         bool foundMatch = false;
3001         if (dbContext) {
3002             auto authFactory = io::AuthorityFactory::create(
3003                 NN_NO_CHECK(dbContext), std::string());
3004             auto list = authFactory->createObjectsFromName(
3005                 l_datum->nameStr(),
3006                 {io::AuthorityFactory::ObjectType::GEODETIC_REFERENCE_FRAME},
3007                 false /* approximate=false*/);
3008             if (!list.empty()) {
3009                 auto gdatum =
3010                     util::nn_dynamic_pointer_cast<datum::Datum>(list.front());
3011                 if (gdatum) {
3012                     gdatum->_exportToWKT(formatter);
3013                     foundMatch = true;
3014                 }
3015             }
3016         }
3017         if (!foundMatch) {
3018             // We should export a geodetic datum, but we cannot really do better
3019             l_datum->_exportToWKT(formatter);
3020         }
3021     } else {
3022         exportDatumOrDatumEnsembleToWkt(formatter);
3023     }
3024     const auto &cs = SingleCRS::getPrivate()->coordinateSystem;
3025     const auto &axisList = cs->axisList();
3026 
3027     if (formatter->useESRIDialect()) {
3028         // Seems to be a constant value...
3029         formatter->startNode(io::WKTConstants::PARAMETER, false);
3030         formatter->addQuotedString("Vertical_Shift");
3031         formatter->add(0.0);
3032         formatter->endNode();
3033 
3034         formatter->startNode(io::WKTConstants::PARAMETER, false);
3035         formatter->addQuotedString("Direction");
3036         formatter->add(
3037             axisList[0]->direction() == cs::AxisDirection::UP ? 1.0 : -1.0);
3038         formatter->endNode();
3039     }
3040 
3041     if (!isWKT2) {
3042         axisList[0]->unit()._exportToWKT(formatter);
3043     }
3044 
3045     const auto oldAxisOutputRule = formatter->outputAxis();
3046     if (oldAxisOutputRule ==
3047         io::WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE) {
3048         formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES);
3049     }
3050     cs->_exportToWKT(formatter);
3051     formatter->setOutputAxis(oldAxisOutputRule);
3052 
3053     if (isWKT2 && formatter->use2019Keywords() && !d->geoidModel.empty()) {
3054         const auto &model = d->geoidModel[0];
3055         formatter->startNode(io::WKTConstants::GEOIDMODEL, false);
3056         formatter->addQuotedString(model->nameStr());
3057         model->formatID(formatter);
3058         formatter->endNode();
3059     }
3060 
3061     ObjectUsage::baseExportToWKT(formatter);
3062     formatter->endNode();
3063 }
3064 //! @endcond
3065 
3066 // ---------------------------------------------------------------------------
3067 
3068 //! @cond Doxygen_Suppress
_exportToPROJString(io::PROJStringFormatter * formatter) const3069 void VerticalCRS::_exportToPROJString(
3070     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
3071 {
3072     auto geoidgrids = formatter->getVDatumExtension();
3073     if (!geoidgrids.empty()) {
3074         formatter->addParam("geoidgrids", geoidgrids);
3075     }
3076 
3077     auto &axisList = coordinateSystem()->axisList();
3078     if (!axisList.empty()) {
3079         auto projUnit = axisList[0]->unit().exportToPROJString();
3080         if (projUnit.empty()) {
3081             formatter->addParam("vto_meter",
3082                                 axisList[0]->unit().conversionToSI());
3083         } else {
3084             formatter->addParam("vunits", projUnit);
3085         }
3086     }
3087 }
3088 //! @endcond
3089 
3090 // ---------------------------------------------------------------------------
3091 
3092 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const3093 void VerticalCRS::_exportToJSON(
3094     io::JSONFormatter *formatter) const // throw(io::FormattingException)
3095 {
3096     auto writer = formatter->writer();
3097     auto objectContext(
3098         formatter->MakeObjectContext("VerticalCRS", !identifiers().empty()));
3099 
3100     writer->AddObjKey("name");
3101     auto l_name = nameStr();
3102     if (l_name.empty()) {
3103         writer->Add("unnamed");
3104     } else {
3105         writer->Add(l_name);
3106     }
3107 
3108     const auto &l_datum(datum());
3109     if (l_datum) {
3110         writer->AddObjKey("datum");
3111         l_datum->_exportToJSON(formatter);
3112     } else {
3113         writer->AddObjKey("datum_ensemble");
3114         formatter->setOmitTypeInImmediateChild();
3115         datumEnsemble()->_exportToJSON(formatter);
3116     }
3117 
3118     writer->AddObjKey("coordinate_system");
3119     formatter->setOmitTypeInImmediateChild();
3120     coordinateSystem()->_exportToJSON(formatter);
3121 
3122     if (!d->geoidModel.empty()) {
3123         const auto &model = d->geoidModel[0];
3124         writer->AddObjKey("geoid_model");
3125         auto objectContext2(formatter->MakeObjectContext(nullptr, false));
3126         writer->AddObjKey("name");
3127         writer->Add(model->nameStr());
3128 
3129         if (model->identifiers().empty()) {
3130             const auto &interpCRS = model->interpolationCRS();
3131             if (interpCRS) {
3132                 writer->AddObjKey("interpolation_crs");
3133                 interpCRS->_exportToJSON(formatter);
3134             }
3135         }
3136 
3137         model->formatID(formatter);
3138     }
3139 
3140     ObjectUsage::baseExportToJSON(formatter);
3141 }
3142 //! @endcond
3143 
3144 // ---------------------------------------------------------------------------
3145 
3146 //! @cond Doxygen_Suppress
addLinearUnitConvert(io::PROJStringFormatter * formatter) const3147 void VerticalCRS::addLinearUnitConvert(
3148     io::PROJStringFormatter *formatter) const {
3149     auto &axisList = coordinateSystem()->axisList();
3150 
3151     if (!axisList.empty()) {
3152         if (axisList[0]->unit().conversionToSI() != 1.0) {
3153             formatter->addStep("unitconvert");
3154             formatter->addParam("z_in", "m");
3155             auto projVUnit = axisList[0]->unit().exportToPROJString();
3156             if (projVUnit.empty()) {
3157                 formatter->addParam("z_out",
3158                                     axisList[0]->unit().conversionToSI());
3159             } else {
3160                 formatter->addParam("z_out", projVUnit);
3161             }
3162         }
3163     }
3164 }
3165 //! @endcond
3166 
3167 // ---------------------------------------------------------------------------
3168 
3169 /** \brief Instantiate a VerticalCRS from a datum::VerticalReferenceFrame and a
3170  * cs::VerticalCS.
3171  *
3172  * @param properties See \ref general_properties.
3173  * At minimum the name should be defined. The GEOID_MODEL property can be set
3174  * to a TransformationNNPtr object.
3175  * @param datumIn The datum of the CRS.
3176  * @param csIn a VerticalCS.
3177  * @return new VerticalCRS.
3178  */
3179 VerticalCRSNNPtr
create(const util::PropertyMap & properties,const datum::VerticalReferenceFrameNNPtr & datumIn,const cs::VerticalCSNNPtr & csIn)3180 VerticalCRS::create(const util::PropertyMap &properties,
3181                     const datum::VerticalReferenceFrameNNPtr &datumIn,
3182                     const cs::VerticalCSNNPtr &csIn) {
3183     return create(properties, datumIn.as_nullable(), nullptr, csIn);
3184 }
3185 
3186 // ---------------------------------------------------------------------------
3187 
3188 /** \brief Instantiate a VerticalCRS from a datum::VerticalReferenceFrame or
3189  * datum::DatumEnsemble and a cs::VerticalCS.
3190  *
3191  * One and only one of datum or datumEnsemble should be set to a non-null value.
3192  *
3193  * @param properties See \ref general_properties.
3194  * At minimum the name should be defined. The GEOID_MODEL property can be set
3195  * to a TransformationNNPtr object.
3196  * @param datumIn The datum of the CRS, or nullptr
3197  * @param datumEnsembleIn The datum ensemble of the CRS, or nullptr.
3198  * @param csIn a VerticalCS.
3199  * @return new VerticalCRS.
3200  */
3201 VerticalCRSNNPtr
create(const util::PropertyMap & properties,const datum::VerticalReferenceFramePtr & datumIn,const datum::DatumEnsemblePtr & datumEnsembleIn,const cs::VerticalCSNNPtr & csIn)3202 VerticalCRS::create(const util::PropertyMap &properties,
3203                     const datum::VerticalReferenceFramePtr &datumIn,
3204                     const datum::DatumEnsemblePtr &datumEnsembleIn,
3205                     const cs::VerticalCSNNPtr &csIn) {
3206     auto crs(VerticalCRS::nn_make_shared<VerticalCRS>(datumIn, datumEnsembleIn,
3207                                                       csIn));
3208     crs->assignSelf(crs);
3209     crs->setProperties(properties);
3210     const auto geoidModelPtr = properties.get("GEOID_MODEL");
3211     if (geoidModelPtr) {
3212         auto transf = util::nn_dynamic_pointer_cast<operation::Transformation>(
3213             *geoidModelPtr);
3214         if (transf) {
3215             crs->d->geoidModel.emplace_back(NN_NO_CHECK(transf));
3216         }
3217     }
3218     return crs;
3219 }
3220 
3221 // ---------------------------------------------------------------------------
3222 
3223 //! @cond Doxygen_Suppress
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const3224 bool VerticalCRS::_isEquivalentTo(
3225     const util::IComparable *other, util::IComparable::Criterion criterion,
3226     const io::DatabaseContextPtr &dbContext) const {
3227     auto otherVertCRS = dynamic_cast<const VerticalCRS *>(other);
3228     // TODO test geoidModel and velocityModel
3229     return otherVertCRS != nullptr &&
3230            SingleCRS::baseIsEquivalentTo(other, criterion, dbContext);
3231 }
3232 //! @endcond
3233 
3234 // ---------------------------------------------------------------------------
3235 
3236 /** \brief Identify the CRS with reference CRSs.
3237  *
3238  * The candidate CRSs are looked in the database when
3239  * authorityFactory is not null.
3240  *
3241  * Note that the implementation uses a set of heuristics to have a good
3242  * compromise of successful identifications over execution time. It might miss
3243  * legitimate matches in some circumstances.
3244  *
3245  * The method returns a list of matching reference CRS, and the percentage
3246  * (0-100) of confidence in the match.
3247  * 100% means that the name of the reference entry
3248  * perfectly matches the CRS name, and both are equivalent. In which case a
3249  * single result is returned.
3250  * 90% means that CRS are equivalent, but the names are not exactly the same.
3251  * 70% means that CRS are equivalent (equivalent datum and coordinate system),
3252  * but the names are not equivalent.
3253  * 25% means that the CRS are not equivalent, but there is some similarity in
3254  * the names.
3255  *
3256  * @param authorityFactory Authority factory (if null, will return an empty
3257  * list)
3258  * @return a list of matching reference CRS, and the percentage (0-100) of
3259  * confidence in the match.
3260  */
3261 std::list<std::pair<VerticalCRSNNPtr, int>>
identify(const io::AuthorityFactoryPtr & authorityFactory) const3262 VerticalCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
3263     typedef std::pair<VerticalCRSNNPtr, int> Pair;
3264     std::list<Pair> res;
3265 
3266     const auto &thisName(nameStr());
3267 
3268     if (authorityFactory) {
3269         const io::DatabaseContextNNPtr &dbContext =
3270             authorityFactory->databaseContext();
3271 
3272         const bool unsignificantName = thisName.empty() ||
3273                                        ci_equal(thisName, "unknown") ||
3274                                        ci_equal(thisName, "unnamed");
3275         if (hasCodeCompatibleOfAuthorityFactory(this, authorityFactory)) {
3276             // If the CRS has already an id, check in the database for the
3277             // official object, and verify that they are equivalent.
3278             for (const auto &id : identifiers()) {
3279                 if (hasCodeCompatibleOfAuthorityFactory(id, authorityFactory)) {
3280                     try {
3281                         auto crs = io::AuthorityFactory::create(
3282                                        dbContext, *id->codeSpace())
3283                                        ->createVerticalCRS(id->code());
3284                         bool match = _isEquivalentTo(
3285                             crs.get(), util::IComparable::Criterion::EQUIVALENT,
3286                             dbContext);
3287                         res.emplace_back(crs, match ? 100 : 25);
3288                         return res;
3289                     } catch (const std::exception &) {
3290                     }
3291                 }
3292             }
3293         } else if (!unsignificantName) {
3294             for (int ipass = 0; ipass < 2; ipass++) {
3295                 const bool approximateMatch = ipass == 1;
3296                 auto objects = authorityFactory->createObjectsFromName(
3297                     thisName, {io::AuthorityFactory::ObjectType::VERTICAL_CRS},
3298                     approximateMatch);
3299                 for (const auto &obj : objects) {
3300                     auto crs = util::nn_dynamic_pointer_cast<VerticalCRS>(obj);
3301                     assert(crs);
3302                     auto crsNN = NN_NO_CHECK(crs);
3303                     if (_isEquivalentTo(
3304                             crs.get(), util::IComparable::Criterion::EQUIVALENT,
3305                             dbContext)) {
3306                         if (crs->nameStr() == thisName) {
3307                             res.clear();
3308                             res.emplace_back(crsNN, 100);
3309                             return res;
3310                         }
3311                         res.emplace_back(crsNN, 90);
3312                     } else {
3313                         res.emplace_back(crsNN, 25);
3314                     }
3315                 }
3316                 if (!res.empty()) {
3317                     break;
3318                 }
3319             }
3320         }
3321 
3322         // Sort results
3323         res.sort([&thisName](const Pair &a, const Pair &b) {
3324             // First consider confidence
3325             if (a.second > b.second) {
3326                 return true;
3327             }
3328             if (a.second < b.second) {
3329                 return false;
3330             }
3331 
3332             // Then consider exact name matching
3333             const auto &aName(a.first->nameStr());
3334             const auto &bName(b.first->nameStr());
3335             if (aName == thisName && bName != thisName) {
3336                 return true;
3337             }
3338             if (bName == thisName && aName != thisName) {
3339                 return false;
3340             }
3341 
3342             // Arbitrary final sorting criterion
3343             return aName < bName;
3344         });
3345 
3346         // Keep only results of the highest confidence
3347         if (res.size() >= 2) {
3348             const auto highestConfidence = res.front().second;
3349             std::list<Pair> newRes;
3350             for (const auto &pair : res) {
3351                 if (pair.second == highestConfidence) {
3352                     newRes.push_back(pair);
3353                 } else {
3354                     break;
3355                 }
3356             }
3357             return newRes;
3358         }
3359     }
3360 
3361     return res;
3362 }
3363 
3364 // ---------------------------------------------------------------------------
3365 
3366 //! @cond Doxygen_Suppress
3367 
3368 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & authorityFactory) const3369 VerticalCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
3370     typedef std::pair<CRSNNPtr, int> Pair;
3371     std::list<Pair> res;
3372     auto resTemp = identify(authorityFactory);
3373     for (const auto &pair : resTemp) {
3374         res.emplace_back(pair.first, pair.second);
3375     }
3376     return res;
3377 }
3378 
3379 //! @endcond
3380 
3381 // ---------------------------------------------------------------------------
3382 
3383 //! @cond Doxygen_Suppress
3384 struct DerivedCRS::Private {
3385     SingleCRSNNPtr baseCRS_;
3386     operation::ConversionNNPtr derivingConversion_;
3387 
Privatecrs::DerivedCRS::Private3388     Private(const SingleCRSNNPtr &baseCRSIn,
3389             const operation::ConversionNNPtr &derivingConversionIn)
3390         : baseCRS_(baseCRSIn), derivingConversion_(derivingConversionIn) {}
3391 
3392     // For the conversion make a _shallowClone(), so that we can later set
3393     // its targetCRS to this.
Privatecrs::DerivedCRS::Private3394     Private(const Private &other)
3395         : baseCRS_(other.baseCRS_),
3396           derivingConversion_(other.derivingConversion_->shallowClone()) {}
3397 };
3398 
3399 //! @endcond
3400 
3401 // ---------------------------------------------------------------------------
3402 
3403 // DerivedCRS is an abstract class, that virtually inherits from SingleCRS
3404 // Consequently the base constructor in SingleCRS will never be called by
3405 // that constructor. clang -Wabstract-vbase-init and VC++ underline this, but
3406 // other
3407 // compilers will complain if we don't call the base constructor.
3408 
DerivedCRS(const SingleCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CoordinateSystemNNPtr & cs)3409 DerivedCRS::DerivedCRS(const SingleCRSNNPtr &baseCRSIn,
3410                        const operation::ConversionNNPtr &derivingConversionIn,
3411                        const cs::CoordinateSystemNNPtr &
3412 #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
3413                            cs
3414 #endif
3415                        )
3416     :
3417 #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
3418       SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), cs),
3419 #endif
3420       d(internal::make_unique<Private>(baseCRSIn, derivingConversionIn)) {
3421 }
3422 
3423 // ---------------------------------------------------------------------------
3424 
DerivedCRS(const DerivedCRS & other)3425 DerivedCRS::DerivedCRS(const DerivedCRS &other)
3426     :
3427 #if !defined(COMPILER_WARNS_ABOUT_ABSTRACT_VBASE_INIT)
3428       SingleCRS(other),
3429 #endif
3430       d(internal::make_unique<Private>(*other.d)) {
3431 }
3432 
3433 // ---------------------------------------------------------------------------
3434 
3435 //! @cond Doxygen_Suppress
3436 DerivedCRS::~DerivedCRS() = default;
3437 //! @endcond
3438 
3439 // ---------------------------------------------------------------------------
3440 
3441 /** \brief Return the base CRS of a DerivedCRS.
3442  *
3443  * @return the base CRS.
3444  */
baseCRS()3445 const SingleCRSNNPtr &DerivedCRS::baseCRS() PROJ_PURE_DEFN {
3446     return d->baseCRS_;
3447 }
3448 
3449 // ---------------------------------------------------------------------------
3450 
3451 /** \brief Return the deriving conversion from the base CRS to this CRS.
3452  *
3453  * @return the deriving conversion.
3454  */
derivingConversion() const3455 const operation::ConversionNNPtr DerivedCRS::derivingConversion() const {
3456     return d->derivingConversion_->shallowClone();
3457 }
3458 
3459 // ---------------------------------------------------------------------------
3460 
3461 //! @cond Doxygen_Suppress
3462 const operation::ConversionNNPtr &
derivingConversionRef()3463 DerivedCRS::derivingConversionRef() PROJ_PURE_DEFN {
3464     return d->derivingConversion_;
3465 }
3466 //! @endcond
3467 
3468 // ---------------------------------------------------------------------------
3469 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const3470 bool DerivedCRS::_isEquivalentTo(
3471     const util::IComparable *other, util::IComparable::Criterion criterion,
3472     const io::DatabaseContextPtr &dbContext) const {
3473     auto otherDerivedCRS = dynamic_cast<const DerivedCRS *>(other);
3474     const auto standardCriterion = getStandardCriterion(criterion);
3475     if (otherDerivedCRS == nullptr ||
3476         !SingleCRS::baseIsEquivalentTo(other, standardCriterion, dbContext)) {
3477         return false;
3478     }
3479     return d->baseCRS_->_isEquivalentTo(otherDerivedCRS->d->baseCRS_.get(),
3480                                         criterion, dbContext) &&
3481            d->derivingConversion_->_isEquivalentTo(
3482                otherDerivedCRS->d->derivingConversion_.get(), standardCriterion,
3483                dbContext);
3484 }
3485 
3486 // ---------------------------------------------------------------------------
3487 
setDerivingConversionCRS()3488 void DerivedCRS::setDerivingConversionCRS() {
3489     derivingConversionRef()->setWeakSourceTargetCRS(
3490         baseCRS().as_nullable(),
3491         std::static_pointer_cast<CRS>(shared_from_this().as_nullable()));
3492 }
3493 
3494 // ---------------------------------------------------------------------------
3495 
baseExportToWKT(io::WKTFormatter * formatter,const std::string & keyword,const std::string & baseKeyword) const3496 void DerivedCRS::baseExportToWKT(io::WKTFormatter *formatter,
3497                                  const std::string &keyword,
3498                                  const std::string &baseKeyword) const {
3499     formatter->startNode(keyword, !identifiers().empty());
3500     formatter->addQuotedString(nameStr());
3501 
3502     const auto &l_baseCRS = d->baseCRS_;
3503     formatter->startNode(baseKeyword, formatter->use2019Keywords() &&
3504                                           !l_baseCRS->identifiers().empty());
3505     formatter->addQuotedString(l_baseCRS->nameStr());
3506     l_baseCRS->exportDatumOrDatumEnsembleToWkt(formatter);
3507     if (formatter->use2019Keywords() &&
3508         !(formatter->idOnTopLevelOnly() && formatter->topLevelHasId())) {
3509         l_baseCRS->formatID(formatter);
3510     }
3511     formatter->endNode();
3512 
3513     formatter->setUseDerivingConversion(true);
3514     derivingConversionRef()->_exportToWKT(formatter);
3515     formatter->setUseDerivingConversion(false);
3516 
3517     coordinateSystem()->_exportToWKT(formatter);
3518     ObjectUsage::baseExportToWKT(formatter);
3519     formatter->endNode();
3520 }
3521 
3522 // ---------------------------------------------------------------------------
3523 
3524 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const3525 void DerivedCRS::_exportToJSON(
3526     io::JSONFormatter *formatter) const // throw(io::FormattingException)
3527 {
3528     auto writer = formatter->writer();
3529     auto objectContext(
3530         formatter->MakeObjectContext(className(), !identifiers().empty()));
3531 
3532     writer->AddObjKey("name");
3533     auto l_name = nameStr();
3534     if (l_name.empty()) {
3535         writer->Add("unnamed");
3536     } else {
3537         writer->Add(l_name);
3538     }
3539 
3540     writer->AddObjKey("base_crs");
3541     baseCRS()->_exportToJSON(formatter);
3542 
3543     writer->AddObjKey("conversion");
3544     formatter->setOmitTypeInImmediateChild();
3545     derivingConversionRef()->_exportToJSON(formatter);
3546 
3547     writer->AddObjKey("coordinate_system");
3548     formatter->setOmitTypeInImmediateChild();
3549     coordinateSystem()->_exportToJSON(formatter);
3550 
3551     ObjectUsage::baseExportToJSON(formatter);
3552 }
3553 //! @endcond
3554 
3555 // ---------------------------------------------------------------------------
3556 
3557 //! @cond Doxygen_Suppress
3558 struct ProjectedCRS::Private {
3559     GeodeticCRSNNPtr baseCRS_;
3560     cs::CartesianCSNNPtr cs_;
Privatecrs::ProjectedCRS::Private3561     Private(const GeodeticCRSNNPtr &baseCRSIn, const cs::CartesianCSNNPtr &csIn)
3562         : baseCRS_(baseCRSIn), cs_(csIn) {}
3563 
baseCRScrs::ProjectedCRS::Private3564     inline const GeodeticCRSNNPtr &baseCRS() const { return baseCRS_; }
3565 
coordinateSystemcrs::ProjectedCRS::Private3566     inline const cs::CartesianCSNNPtr &coordinateSystem() const { return cs_; }
3567 };
3568 //! @endcond
3569 
3570 // ---------------------------------------------------------------------------
3571 
ProjectedCRS(const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CartesianCSNNPtr & csIn)3572 ProjectedCRS::ProjectedCRS(
3573     const GeodeticCRSNNPtr &baseCRSIn,
3574     const operation::ConversionNNPtr &derivingConversionIn,
3575     const cs::CartesianCSNNPtr &csIn)
3576     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
3577       DerivedCRS(baseCRSIn, derivingConversionIn, csIn),
3578       d(internal::make_unique<Private>(baseCRSIn, csIn)) {}
3579 
3580 // ---------------------------------------------------------------------------
3581 
ProjectedCRS(const ProjectedCRS & other)3582 ProjectedCRS::ProjectedCRS(const ProjectedCRS &other)
3583     : SingleCRS(other), DerivedCRS(other),
3584       d(internal::make_unique<Private>(other.baseCRS(),
3585                                        other.coordinateSystem())) {}
3586 
3587 // ---------------------------------------------------------------------------
3588 
3589 //! @cond Doxygen_Suppress
3590 ProjectedCRS::~ProjectedCRS() = default;
3591 //! @endcond
3592 
3593 // ---------------------------------------------------------------------------
3594 
_shallowClone() const3595 CRSNNPtr ProjectedCRS::_shallowClone() const {
3596     auto crs(ProjectedCRS::nn_make_shared<ProjectedCRS>(*this));
3597     crs->assignSelf(crs);
3598     crs->setDerivingConversionCRS();
3599     return crs;
3600 }
3601 
3602 // ---------------------------------------------------------------------------
3603 
3604 /** \brief Return the base CRS (a GeodeticCRS, which is generally a
3605  * GeographicCRS) of the ProjectedCRS.
3606  *
3607  * @return the base CRS.
3608  */
baseCRS()3609 const GeodeticCRSNNPtr &ProjectedCRS::baseCRS() PROJ_PURE_DEFN {
3610     return d->baseCRS();
3611 }
3612 
3613 // ---------------------------------------------------------------------------
3614 
3615 /** \brief Return the cs::CartesianCS associated with the CRS.
3616  *
3617  * @return a CartesianCS
3618  */
coordinateSystem()3619 const cs::CartesianCSNNPtr &ProjectedCRS::coordinateSystem() PROJ_PURE_DEFN {
3620     return d->coordinateSystem();
3621 }
3622 
3623 // ---------------------------------------------------------------------------
3624 
3625 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const3626 void ProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
3627     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
3628 
3629     const auto &l_identifiers = identifiers();
3630     // Try to perfectly round-trip ESRI projectedCRS if the current object
3631     // perfectly matches the database definition
3632     const auto &dbContext = formatter->databaseContext();
3633 
3634     auto l_name = nameStr();
3635     const auto &l_coordinateSystem = d->coordinateSystem();
3636     const auto &axisList = l_coordinateSystem->axisList();
3637     if (axisList.size() == 3 && !(isWKT2 && formatter->use2019Keywords())) {
3638         auto projCRS2D = demoteTo2D(std::string(), dbContext);
3639         if (dbContext) {
3640             const auto res = projCRS2D->identify(
3641                 io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "EPSG"));
3642             if (res.size() == 1) {
3643                 const auto &front = res.front();
3644                 if (front.second == 100) {
3645                     projCRS2D = front.first;
3646                 }
3647             }
3648         }
3649 
3650         if (formatter->useESRIDialect() && dbContext) {
3651             // Try to format the ProjecteD 3D CRS as a
3652             // PROJCS[],VERTCS[...,DATUM[]]
3653             // if we find corresponding objects
3654             if (exportAsESRIWktCompoundCRSWithEllipsoidalHeight(
3655                     this, baseCRS().as_nullable().get(), formatter)) {
3656                 return;
3657             }
3658         }
3659 
3660         if (!formatter->useESRIDialect() &&
3661             CRS::getPrivate()->allowNonConformantWKT1Export_) {
3662             formatter->startNode(io::WKTConstants::COMPD_CS, false);
3663             formatter->addQuotedString(l_name + " + " + baseCRS()->nameStr());
3664             projCRS2D->_exportToWKT(formatter);
3665             baseCRS()
3666                 ->demoteTo2D(std::string(), dbContext)
3667                 ->_exportToWKT(formatter);
3668             formatter->endNode();
3669             return;
3670         }
3671 
3672         auto &originalCompoundCRS = CRS::getPrivate()->originalCompoundCRS_;
3673         if (!formatter->useESRIDialect() && originalCompoundCRS) {
3674             originalCompoundCRS->_exportToWKT(formatter);
3675             return;
3676         }
3677 
3678         if (!formatter->useESRIDialect() &&
3679             formatter->isAllowedEllipsoidalHeightAsVerticalCRS()) {
3680             if (exportAsWKT1CompoundCRSWithEllipsoidalHeight(
3681                     projCRS2D, axisList[2], formatter)) {
3682                 return;
3683             }
3684         }
3685 
3686         io::FormattingException::Throw(
3687             "Projected 3D CRS can only be exported since WKT2:2019");
3688     }
3689 
3690     std::string l_alias;
3691     if (formatter->useESRIDialect() && dbContext) {
3692         l_alias = dbContext->getAliasFromOfficialName(l_name, "projected_crs",
3693                                                       "ESRI");
3694     }
3695 
3696     if (!isWKT2 && formatter->useESRIDialect() && !l_identifiers.empty() &&
3697         *(l_identifiers[0]->codeSpace()) == "ESRI" && dbContext) {
3698         try {
3699             const auto definition = dbContext->getTextDefinition(
3700                 "projected_crs", "ESRI", l_identifiers[0]->code());
3701             if (starts_with(definition, "PROJCS")) {
3702                 auto crsFromFromDef = io::WKTParser()
3703                                           .attachDatabaseContext(dbContext)
3704                                           .createFromWKT(definition);
3705                 if (_isEquivalentTo(
3706                         dynamic_cast<IComparable *>(crsFromFromDef.get()),
3707                         util::IComparable::Criterion::EQUIVALENT)) {
3708                     formatter->ingestWKTNode(
3709                         io::WKTNode::createFrom(definition));
3710                     return;
3711                 }
3712             }
3713         } catch (const std::exception &) {
3714         }
3715     } else if (!isWKT2 && formatter->useESRIDialect() && !l_alias.empty()) {
3716         try {
3717             auto res =
3718                 io::AuthorityFactory::create(NN_NO_CHECK(dbContext), "ESRI")
3719                     ->createObjectsFromName(
3720                         l_alias,
3721                         {io::AuthorityFactory::ObjectType::PROJECTED_CRS},
3722                         false);
3723             if (res.size() == 1) {
3724                 const auto definition = dbContext->getTextDefinition(
3725                     "projected_crs", "ESRI",
3726                     res.front()->identifiers()[0]->code());
3727                 if (starts_with(definition, "PROJCS")) {
3728                     if (_isEquivalentTo(
3729                             dynamic_cast<IComparable *>(res.front().get()),
3730                             util::IComparable::Criterion::EQUIVALENT)) {
3731                         formatter->ingestWKTNode(
3732                             io::WKTNode::createFrom(definition));
3733                         return;
3734                     }
3735                 }
3736             }
3737         } catch (const std::exception &) {
3738         }
3739     }
3740 
3741     const auto exportAxis = [&l_coordinateSystem, &axisList, &formatter]() {
3742         const auto oldAxisOutputRule = formatter->outputAxis();
3743         if (oldAxisOutputRule ==
3744             io::WKTFormatter::OutputAxisRule::WKT1_GDAL_EPSG_STYLE) {
3745             if (&axisList[0]->direction() == &cs::AxisDirection::EAST &&
3746                 &axisList[1]->direction() == &cs::AxisDirection::NORTH) {
3747                 formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES);
3748             }
3749         }
3750         l_coordinateSystem->_exportToWKT(formatter);
3751         formatter->setOutputAxis(oldAxisOutputRule);
3752     };
3753 
3754     if (!isWKT2 && !formatter->useESRIDialect() &&
3755         starts_with(nameStr(), "Popular Visualisation CRS / Mercator")) {
3756         formatter->startNode(io::WKTConstants::PROJCS, !l_identifiers.empty());
3757         formatter->addQuotedString(nameStr());
3758         formatter->setTOWGS84Parameters({0, 0, 0, 0, 0, 0, 0});
3759         baseCRS()->_exportToWKT(formatter);
3760         formatter->setTOWGS84Parameters({});
3761 
3762         formatter->startNode(io::WKTConstants::PROJECTION, false);
3763         formatter->addQuotedString("Mercator_1SP");
3764         formatter->endNode();
3765 
3766         formatter->startNode(io::WKTConstants::PARAMETER, false);
3767         formatter->addQuotedString("central_meridian");
3768         formatter->add(0.0);
3769         formatter->endNode();
3770 
3771         formatter->startNode(io::WKTConstants::PARAMETER, false);
3772         formatter->addQuotedString("scale_factor");
3773         formatter->add(1.0);
3774         formatter->endNode();
3775 
3776         formatter->startNode(io::WKTConstants::PARAMETER, false);
3777         formatter->addQuotedString("false_easting");
3778         formatter->add(0.0);
3779         formatter->endNode();
3780 
3781         formatter->startNode(io::WKTConstants::PARAMETER, false);
3782         formatter->addQuotedString("false_northing");
3783         formatter->add(0.0);
3784         formatter->endNode();
3785 
3786         axisList[0]->unit()._exportToWKT(formatter);
3787         exportAxis();
3788         derivingConversionRef()->addWKTExtensionNode(formatter);
3789         ObjectUsage::baseExportToWKT(formatter);
3790         formatter->endNode();
3791         return;
3792     }
3793 
3794     formatter->startNode(isWKT2 ? io::WKTConstants::PROJCRS
3795                                 : io::WKTConstants::PROJCS,
3796                          !l_identifiers.empty());
3797 
3798     if (formatter->useESRIDialect()) {
3799         if (l_alias.empty()) {
3800             l_name = io::WKTFormatter::morphNameToESRI(l_name);
3801         } else {
3802             l_name = l_alias;
3803         }
3804     }
3805     if (!isWKT2 && !formatter->useESRIDialect() && isDeprecated()) {
3806         l_name += " (deprecated)";
3807     }
3808     formatter->addQuotedString(l_name);
3809 
3810     const auto &l_baseCRS = d->baseCRS();
3811     const auto &geodeticCRSAxisList = l_baseCRS->coordinateSystem()->axisList();
3812 
3813     if (isWKT2) {
3814         formatter->startNode(
3815             (formatter->use2019Keywords() &&
3816              dynamic_cast<const GeographicCRS *>(l_baseCRS.get()))
3817                 ? io::WKTConstants::BASEGEOGCRS
3818                 : io::WKTConstants::BASEGEODCRS,
3819             formatter->use2019Keywords() && !l_baseCRS->identifiers().empty());
3820         formatter->addQuotedString(l_baseCRS->nameStr());
3821         l_baseCRS->exportDatumOrDatumEnsembleToWkt(formatter);
3822         // insert ellipsoidal cs unit when the units of the map
3823         // projection angular parameters are not explicitly given within those
3824         // parameters. See
3825         // http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#61
3826         if (formatter->primeMeridianOrParameterUnitOmittedIfSameAsAxis()) {
3827             geodeticCRSAxisList[0]->unit()._exportToWKT(formatter);
3828         }
3829         l_baseCRS->primeMeridian()->_exportToWKT(formatter);
3830         if (formatter->use2019Keywords() &&
3831             !(formatter->idOnTopLevelOnly() && formatter->topLevelHasId())) {
3832             l_baseCRS->formatID(formatter);
3833         }
3834         formatter->endNode();
3835     } else {
3836         const auto oldAxisOutputRule = formatter->outputAxis();
3837         formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::NO);
3838         l_baseCRS->_exportToWKT(formatter);
3839         formatter->setOutputAxis(oldAxisOutputRule);
3840     }
3841 
3842     formatter->pushAxisLinearUnit(
3843         common::UnitOfMeasure::create(axisList[0]->unit()));
3844 
3845     formatter->pushAxisAngularUnit(
3846         common::UnitOfMeasure::create(geodeticCRSAxisList[0]->unit()));
3847 
3848     derivingConversionRef()->_exportToWKT(formatter);
3849 
3850     formatter->popAxisAngularUnit();
3851 
3852     formatter->popAxisLinearUnit();
3853 
3854     if (!isWKT2) {
3855         axisList[0]->unit()._exportToWKT(formatter);
3856     }
3857 
3858     exportAxis();
3859 
3860     if (!isWKT2 && !formatter->useESRIDialect()) {
3861         const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
3862         if (!extensionProj4.empty()) {
3863             formatter->startNode(io::WKTConstants::EXTENSION, false);
3864             formatter->addQuotedString("PROJ4");
3865             formatter->addQuotedString(extensionProj4);
3866             formatter->endNode();
3867         } else {
3868             derivingConversionRef()->addWKTExtensionNode(formatter);
3869         }
3870     }
3871 
3872     ObjectUsage::baseExportToWKT(formatter);
3873     formatter->endNode();
3874     return;
3875 }
3876 //! @endcond
3877 
3878 // ---------------------------------------------------------------------------
3879 
3880 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const3881 void ProjectedCRS::_exportToJSON(
3882     io::JSONFormatter *formatter) const // throw(io::FormattingException)
3883 {
3884     auto writer = formatter->writer();
3885     auto objectContext(
3886         formatter->MakeObjectContext("ProjectedCRS", !identifiers().empty()));
3887 
3888     writer->AddObjKey("name");
3889     auto l_name = nameStr();
3890     if (l_name.empty()) {
3891         writer->Add("unnamed");
3892     } else {
3893         writer->Add(l_name);
3894     }
3895 
3896     writer->AddObjKey("base_crs");
3897     formatter->setAllowIDInImmediateChild();
3898     formatter->setOmitTypeInImmediateChild();
3899     baseCRS()->_exportToJSON(formatter);
3900 
3901     writer->AddObjKey("conversion");
3902     formatter->setOmitTypeInImmediateChild();
3903     derivingConversionRef()->_exportToJSON(formatter);
3904 
3905     writer->AddObjKey("coordinate_system");
3906     formatter->setOmitTypeInImmediateChild();
3907     coordinateSystem()->_exportToJSON(formatter);
3908 
3909     ObjectUsage::baseExportToJSON(formatter);
3910 }
3911 //! @endcond
3912 
3913 // ---------------------------------------------------------------------------
3914 
_exportToPROJString(io::PROJStringFormatter * formatter) const3915 void ProjectedCRS::_exportToPROJString(
3916     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
3917 {
3918     const auto &extensionProj4 = CRS::getPrivate()->extensionProj4_;
3919     if (!extensionProj4.empty()) {
3920         formatter->ingestPROJString(
3921             replaceAll(extensionProj4, " +type=crs", ""));
3922         formatter->addNoDefs(false);
3923         return;
3924     }
3925 
3926     derivingConversionRef()->_exportToPROJString(formatter);
3927 }
3928 
3929 // ---------------------------------------------------------------------------
3930 
3931 /** \brief Instantiate a ProjectedCRS from a base CRS, a deriving
3932  * operation::Conversion
3933  * and a coordinate system.
3934  *
3935  * @param properties See \ref general_properties.
3936  * At minimum the name should be defined.
3937  * @param baseCRSIn The base CRS, a GeodeticCRS that is generally a
3938  * GeographicCRS.
3939  * @param derivingConversionIn The deriving operation::Conversion (typically
3940  * using a map
3941  * projection method)
3942  * @param csIn The coordniate system.
3943  * @return new ProjectedCRS.
3944  */
3945 ProjectedCRSNNPtr
create(const util::PropertyMap & properties,const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CartesianCSNNPtr & csIn)3946 ProjectedCRS::create(const util::PropertyMap &properties,
3947                      const GeodeticCRSNNPtr &baseCRSIn,
3948                      const operation::ConversionNNPtr &derivingConversionIn,
3949                      const cs::CartesianCSNNPtr &csIn) {
3950     auto crs = ProjectedCRS::nn_make_shared<ProjectedCRS>(
3951         baseCRSIn, derivingConversionIn, csIn);
3952     crs->assignSelf(crs);
3953     crs->setProperties(properties);
3954     crs->setDerivingConversionCRS();
3955     properties.getStringValue("EXTENSION_PROJ4",
3956                               crs->CRS::getPrivate()->extensionProj4_);
3957     crs->CRS::getPrivate()->setImplicitCS(properties);
3958     return crs;
3959 }
3960 
3961 // ---------------------------------------------------------------------------
3962 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const3963 bool ProjectedCRS::_isEquivalentTo(
3964     const util::IComparable *other, util::IComparable::Criterion criterion,
3965     const io::DatabaseContextPtr &dbContext) const {
3966     return other != nullptr && util::isOfExactType<ProjectedCRS>(*other) &&
3967            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
3968 }
3969 
3970 // ---------------------------------------------------------------------------
3971 
3972 //! @cond Doxygen_Suppress
3973 ProjectedCRSNNPtr
alterParametersLinearUnit(const common::UnitOfMeasure & unit,bool convertToNewUnit) const3974 ProjectedCRS::alterParametersLinearUnit(const common::UnitOfMeasure &unit,
3975                                         bool convertToNewUnit) const {
3976     return create(
3977         createPropertyMap(this), baseCRS(),
3978         derivingConversion()->alterParametersLinearUnit(unit, convertToNewUnit),
3979         coordinateSystem());
3980 }
3981 //! @endcond
3982 
3983 // ---------------------------------------------------------------------------
3984 
3985 //! @cond Doxygen_Suppress
addUnitConvertAndAxisSwap(io::PROJStringFormatter * formatter,bool axisSpecFound) const3986 void ProjectedCRS::addUnitConvertAndAxisSwap(io::PROJStringFormatter *formatter,
3987                                              bool axisSpecFound) const {
3988     const auto &axisList = d->coordinateSystem()->axisList();
3989     const auto &unit = axisList[0]->unit();
3990     const auto *zUnit = axisList.size() == 3 ? &(axisList[2]->unit()) : nullptr;
3991     if (!unit._isEquivalentTo(common::UnitOfMeasure::METRE,
3992                               util::IComparable::Criterion::EQUIVALENT) ||
3993         (zUnit &&
3994          !zUnit->_isEquivalentTo(common::UnitOfMeasure::METRE,
3995                                  util::IComparable::Criterion::EQUIVALENT))) {
3996         auto projUnit = unit.exportToPROJString();
3997         const double toSI = unit.conversionToSI();
3998         if (!formatter->getCRSExport()) {
3999             formatter->addStep("unitconvert");
4000             formatter->addParam("xy_in", "m");
4001             if (zUnit)
4002                 formatter->addParam("z_in", "m");
4003 
4004             if (projUnit.empty()) {
4005                 formatter->addParam("xy_out", toSI);
4006             } else {
4007                 formatter->addParam("xy_out", projUnit);
4008             }
4009             if (zUnit) {
4010                 auto projZUnit = zUnit->exportToPROJString();
4011                 const double zToSI = zUnit->conversionToSI();
4012                 if (projZUnit.empty()) {
4013                     formatter->addParam("z_out", zToSI);
4014                 } else {
4015                     formatter->addParam("z_out", projZUnit);
4016                 }
4017             }
4018         } else {
4019             if (projUnit.empty()) {
4020                 formatter->addParam("to_meter", toSI);
4021             } else {
4022                 formatter->addParam("units", projUnit);
4023             }
4024         }
4025     } else if (formatter->getCRSExport() &&
4026                !formatter->getLegacyCRSToCRSContext()) {
4027         formatter->addParam("units", "m");
4028     }
4029 
4030     if (!axisSpecFound && !formatter->getCRSExport()) {
4031         const auto &dir0 = axisList[0]->direction();
4032         const auto &dir1 = axisList[1]->direction();
4033         if (!(&dir0 == &cs::AxisDirection::EAST &&
4034               &dir1 == &cs::AxisDirection::NORTH) &&
4035             // For polar projections, that have south+south direction,
4036             // we don't want to mess with axes.
4037             dir0 != dir1) {
4038 
4039             const char *order[2] = {nullptr, nullptr};
4040             for (int i = 0; i < 2; i++) {
4041                 const auto &dir = axisList[i]->direction();
4042                 if (&dir == &cs::AxisDirection::WEST)
4043                     order[i] = "-1";
4044                 else if (&dir == &cs::AxisDirection::EAST)
4045                     order[i] = "1";
4046                 else if (&dir == &cs::AxisDirection::SOUTH)
4047                     order[i] = "-2";
4048                 else if (&dir == &cs::AxisDirection::NORTH)
4049                     order[i] = "2";
4050             }
4051 
4052             if (order[0] && order[1]) {
4053                 formatter->addStep("axisswap");
4054                 char orderStr[10];
4055                 sprintf(orderStr, "%.2s,%.2s", order[0], order[1]);
4056                 formatter->addParam("order", orderStr);
4057             }
4058         } else {
4059             const auto &name0 = axisList[0]->nameStr();
4060             const auto &name1 = axisList[1]->nameStr();
4061             const bool northingEasting = ci_starts_with(name0, "northing") &&
4062                                          ci_starts_with(name1, "easting");
4063             // case of EPSG:32661 ["WGS 84 / UPS North (N,E)]"
4064             // case of EPSG:32761 ["WGS 84 / UPS South (N,E)]"
4065             if (((&dir0 == &cs::AxisDirection::SOUTH &&
4066                   &dir1 == &cs::AxisDirection::SOUTH) ||
4067                  (&dir0 == &cs::AxisDirection::NORTH &&
4068                   &dir1 == &cs::AxisDirection::NORTH)) &&
4069                 northingEasting) {
4070                 formatter->addStep("axisswap");
4071                 formatter->addParam("order", "2,1");
4072             }
4073         }
4074     }
4075 }
4076 //! @endcond
4077 
4078 // ---------------------------------------------------------------------------
4079 
4080 /** \brief Identify the CRS with reference CRSs.
4081  *
4082  * The candidate CRSs are either hard-coded, or looked in the database when
4083  * authorityFactory is not null.
4084  *
4085  * Note that the implementation uses a set of heuristics to have a good
4086  * compromise of successful identifications over execution time. It might miss
4087  * legitimate matches in some circumstances.
4088  *
4089  * The method returns a list of matching reference CRS, and the percentage
4090  * (0-100) of confidence in the match. The list is sorted by decreasing
4091  * confidence.
4092  *
4093  * 100% means that the name of the reference entry
4094  * perfectly matches the CRS name, and both are equivalent. In which case a
4095  * single result is returned.
4096  * 90% means that CRS are equivalent, but the names are not exactly the same.
4097  * 70% means that CRS are equivalent (equivalent base CRS, conversion and
4098  * coordinate system), but the names are not equivalent.
4099  * 60% means that CRS have strong similarity (equivalent base datum, conversion
4100  * and coordinate system), but the names are not equivalent.
4101  * 50% means that CRS have similarity (equivalent base ellipsoid and
4102  * conversion),
4103  * but the coordinate system do not match (e.g. different axis ordering or
4104  * axis unit).
4105  * 25% means that the CRS are not equivalent, but there is some similarity in
4106  * the names.
4107  *
4108  * For the purpose of this function, equivalence is tested with the
4109  * util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS, that is
4110  * to say that the axis order of the base GeographicCRS is ignored.
4111  *
4112  * @param authorityFactory Authority factory (or null, but degraded
4113  * functionality)
4114  * @return a list of matching reference CRS, and the percentage (0-100) of
4115  * confidence in the match.
4116  */
4117 std::list<std::pair<ProjectedCRSNNPtr, int>>
identify(const io::AuthorityFactoryPtr & authorityFactory) const4118 ProjectedCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
4119     typedef std::pair<ProjectedCRSNNPtr, int> Pair;
4120     std::list<Pair> res;
4121 
4122     const auto &thisName(nameStr());
4123     io::DatabaseContextPtr dbContext =
4124         authorityFactory ? authorityFactory->databaseContext().as_nullable()
4125                          : nullptr;
4126 
4127     std::list<std::pair<GeodeticCRSNNPtr, int>> baseRes;
4128     const auto &l_baseCRS(baseCRS());
4129     const auto l_datum = l_baseCRS->datumNonNull(dbContext);
4130     const bool significantNameForDatum =
4131         !ci_starts_with(l_datum->nameStr(), "unknown") &&
4132         l_datum->nameStr() != "unnamed";
4133     const auto &ellipsoid = l_baseCRS->ellipsoid();
4134     auto geogCRS = dynamic_cast<const GeographicCRS *>(l_baseCRS.get());
4135     if (geogCRS &&
4136         geogCRS->coordinateSystem()->axisOrder() ==
4137             cs::EllipsoidalCS::AxisOrder::LONG_EAST_LAT_NORTH) {
4138         baseRes =
4139             GeographicCRS::create(
4140                 util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
4141                                         geogCRS->nameStr()),
4142                 geogCRS->datum(), geogCRS->datumEnsemble(),
4143                 cs::EllipsoidalCS::createLatitudeLongitude(
4144                     geogCRS->coordinateSystem()->axisList()[0]->unit()))
4145                 ->identify(authorityFactory);
4146     } else {
4147         baseRes = l_baseCRS->identify(authorityFactory);
4148     }
4149 
4150     int zone = 0;
4151     bool north = false;
4152 
4153     auto computeConfidence = [&thisName](const std::string &crsName) {
4154         return crsName == thisName ? 100
4155                                    : metadata::Identifier::isEquivalentName(
4156                                          crsName.c_str(), thisName.c_str())
4157                                          ? 90
4158                                          : 70;
4159     };
4160     auto computeUTMCRSName = [](const char *base, int l_zone, bool l_north) {
4161         return base + toString(l_zone) + (l_north ? "N" : "S");
4162     };
4163 
4164     const auto &conv = derivingConversionRef();
4165     const auto &cs = coordinateSystem();
4166 
4167     if (baseRes.size() == 1 && baseRes.front().second >= 70 &&
4168         conv->isUTM(zone, north) &&
4169         cs->_isEquivalentTo(
4170             cs::CartesianCS::createEastingNorthing(common::UnitOfMeasure::METRE)
4171                 .get(),
4172             util::IComparable::Criterion::EQUIVALENT, dbContext)) {
4173         if (baseRes.front().first->_isEquivalentTo(
4174                 GeographicCRS::EPSG_4326.get(),
4175                 util::IComparable::Criterion::EQUIVALENT, dbContext)) {
4176             std::string crsName(
4177                 computeUTMCRSName("WGS 84 / UTM zone ", zone, north));
4178             res.emplace_back(
4179                 ProjectedCRS::create(
4180                     createMapNameEPSGCode(crsName.c_str(),
4181                                           (north ? 32600 : 32700) + zone),
4182                     GeographicCRS::EPSG_4326, conv->identify(), cs),
4183                 computeConfidence(crsName));
4184             return res;
4185         } else if (((zone >= 1 && zone <= 22) || zone == 59 || zone == 60) &&
4186                    north &&
4187                    baseRes.front().first->_isEquivalentTo(
4188                        GeographicCRS::EPSG_4267.get(),
4189                        util::IComparable::Criterion::EQUIVALENT, dbContext)) {
4190             std::string crsName(
4191                 computeUTMCRSName("NAD27 / UTM zone ", zone, north));
4192             res.emplace_back(
4193                 ProjectedCRS::create(
4194                     createMapNameEPSGCode(crsName.c_str(),
4195                                           (zone >= 59) ? 3370 + zone - 59
4196                                                        : 26700 + zone),
4197                     GeographicCRS::EPSG_4267, conv->identify(), cs),
4198                 computeConfidence(crsName));
4199             return res;
4200         } else if (((zone >= 1 && zone <= 23) || zone == 59 || zone == 60) &&
4201                    north &&
4202                    baseRes.front().first->_isEquivalentTo(
4203                        GeographicCRS::EPSG_4269.get(),
4204                        util::IComparable::Criterion::EQUIVALENT, dbContext)) {
4205             std::string crsName(
4206                 computeUTMCRSName("NAD83 / UTM zone ", zone, north));
4207             res.emplace_back(
4208                 ProjectedCRS::create(
4209                     createMapNameEPSGCode(crsName.c_str(),
4210                                           (zone >= 59) ? 3372 + zone - 59
4211                                                        : 26900 + zone),
4212                     GeographicCRS::EPSG_4269, conv->identify(), cs),
4213                 computeConfidence(crsName));
4214             return res;
4215         }
4216     }
4217 
4218     const bool l_implicitCS = hasImplicitCS();
4219     const auto addCRS = [&](const ProjectedCRSNNPtr &crs, const bool eqName) {
4220         const auto &l_unit = cs->axisList()[0]->unit();
4221         if (_isEquivalentTo(crs.get(), util::IComparable::Criterion::
4222                                            EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
4223                             dbContext) ||
4224             (l_implicitCS &&
4225              l_unit._isEquivalentTo(
4226                  crs->coordinateSystem()->axisList()[0]->unit(),
4227                  util::IComparable::Criterion::EQUIVALENT) &&
4228              l_baseCRS->_isEquivalentTo(
4229                  crs->baseCRS().get(), util::IComparable::Criterion::
4230                                            EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
4231                  dbContext) &&
4232              derivingConversionRef()->_isEquivalentTo(
4233                  crs->derivingConversionRef().get(),
4234                  util::IComparable::Criterion::EQUIVALENT, dbContext))) {
4235             if (crs->nameStr() == thisName) {
4236                 res.clear();
4237                 res.emplace_back(crs, 100);
4238             } else {
4239                 res.emplace_back(crs, eqName ? 90 : 70);
4240             }
4241         } else if (ellipsoid->_isEquivalentTo(
4242                        crs->baseCRS()->ellipsoid().get(),
4243                        util::IComparable::Criterion::EQUIVALENT, dbContext) &&
4244                    derivingConversionRef()->_isEquivalentTo(
4245                        crs->derivingConversionRef().get(),
4246                        util::IComparable::Criterion::EQUIVALENT, dbContext)) {
4247             if ((l_implicitCS &&
4248                  l_unit._isEquivalentTo(
4249                      crs->coordinateSystem()->axisList()[0]->unit(),
4250                      util::IComparable::Criterion::EQUIVALENT)) ||
4251                 cs->_isEquivalentTo(crs->coordinateSystem().get(),
4252                                     util::IComparable::Criterion::EQUIVALENT,
4253                                     dbContext)) {
4254                 if (!significantNameForDatum ||
4255                     l_datum->_isEquivalentTo(
4256                         crs->baseCRS()->datumNonNull(dbContext).get(),
4257                         util::IComparable::Criterion::EQUIVALENT)) {
4258                     res.emplace_back(crs, 70);
4259                 } else {
4260                     res.emplace_back(crs, 60);
4261                 }
4262             } else {
4263                 res.emplace_back(crs, 50);
4264             }
4265         } else {
4266             res.emplace_back(crs, 25);
4267         }
4268     };
4269 
4270     if (authorityFactory) {
4271 
4272         const bool unsignificantName = thisName.empty() ||
4273                                        ci_equal(thisName, "unknown") ||
4274                                        ci_equal(thisName, "unnamed");
4275         bool foundEquivalentName = false;
4276 
4277         if (hasCodeCompatibleOfAuthorityFactory(this, authorityFactory)) {
4278             // If the CRS has already an id, check in the database for the
4279             // official object, and verify that they are equivalent.
4280             for (const auto &id : identifiers()) {
4281                 if (hasCodeCompatibleOfAuthorityFactory(id, authorityFactory)) {
4282                     try {
4283                         auto crs = io::AuthorityFactory::create(
4284                                        authorityFactory->databaseContext(),
4285                                        *id->codeSpace())
4286                                        ->createProjectedCRS(id->code());
4287                         bool match = _isEquivalentTo(
4288                             crs.get(), util::IComparable::Criterion::
4289                                            EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS,
4290                             dbContext);
4291                         res.emplace_back(crs, match ? 100 : 25);
4292                         return res;
4293                     } catch (const std::exception &) {
4294                     }
4295                 }
4296             }
4297         } else if (!unsignificantName) {
4298             for (int ipass = 0; ipass < 2; ipass++) {
4299                 const bool approximateMatch = ipass == 1;
4300                 auto objects = authorityFactory->createObjectsFromNameEx(
4301                     thisName, {io::AuthorityFactory::ObjectType::PROJECTED_CRS},
4302                     approximateMatch);
4303                 for (const auto &pairObjName : objects) {
4304                     auto crs = util::nn_dynamic_pointer_cast<ProjectedCRS>(
4305                         pairObjName.first);
4306                     assert(crs);
4307                     auto crsNN = NN_NO_CHECK(crs);
4308                     const bool eqName = metadata::Identifier::isEquivalentName(
4309                         thisName.c_str(), pairObjName.second.c_str());
4310                     foundEquivalentName |= eqName;
4311 
4312                     addCRS(crsNN, eqName);
4313                     if (res.back().second == 100) {
4314                         return res;
4315                     }
4316                 }
4317                 if (!res.empty()) {
4318                     break;
4319                 }
4320             }
4321         }
4322 
4323         const auto lambdaSort = [&thisName](const Pair &a, const Pair &b) {
4324             // First consider confidence
4325             if (a.second > b.second) {
4326                 return true;
4327             }
4328             if (a.second < b.second) {
4329                 return false;
4330             }
4331 
4332             // Then consider exact name matching
4333             const auto &aName(a.first->nameStr());
4334             const auto &bName(b.first->nameStr());
4335             if (aName == thisName && bName != thisName) {
4336                 return true;
4337             }
4338             if (bName == thisName && aName != thisName) {
4339                 return false;
4340             }
4341 
4342             // Arbitrary final sorting criterion
4343             return aName < bName;
4344         };
4345 
4346         // Sort results
4347         res.sort(lambdaSort);
4348 
4349         if (!hasCodeCompatibleOfAuthorityFactory(this, authorityFactory) &&
4350             !foundEquivalentName && (res.empty() || res.front().second < 50)) {
4351             std::set<std::pair<std::string, std::string>> alreadyKnown;
4352             for (const auto &pair : res) {
4353                 const auto &ids = pair.first->identifiers();
4354                 assert(!ids.empty());
4355                 alreadyKnown.insert(std::pair<std::string, std::string>(
4356                     *(ids[0]->codeSpace()), ids[0]->code()));
4357             }
4358 
4359             auto self = NN_NO_CHECK(std::dynamic_pointer_cast<ProjectedCRS>(
4360                 shared_from_this().as_nullable()));
4361             auto candidates =
4362                 authorityFactory->createProjectedCRSFromExisting(self);
4363             for (const auto &crs : candidates) {
4364                 const auto &ids = crs->identifiers();
4365                 assert(!ids.empty());
4366                 if (alreadyKnown.find(std::pair<std::string, std::string>(
4367                         *(ids[0]->codeSpace()), ids[0]->code())) !=
4368                     alreadyKnown.end()) {
4369                     continue;
4370                 }
4371 
4372                 addCRS(crs, unsignificantName);
4373             }
4374 
4375             res.sort(lambdaSort);
4376         }
4377 
4378         // Keep only results of the highest confidence
4379         if (res.size() >= 2) {
4380             const auto highestConfidence = res.front().second;
4381             std::list<Pair> newRes;
4382             for (const auto &pair : res) {
4383                 if (pair.second == highestConfidence) {
4384                     newRes.push_back(pair);
4385                 } else {
4386                     break;
4387                 }
4388             }
4389             return newRes;
4390         }
4391     }
4392 
4393     return res;
4394 }
4395 
4396 // ---------------------------------------------------------------------------
4397 
4398 /** \brief Return a variant of this CRS "demoted" to a 2D one, if not already
4399  * the case.
4400  *
4401  *
4402  * @param newName Name of the new CRS. If empty, nameStr() will be used.
4403  * @param dbContext Database context to look for potentially already registered
4404  *                  2D CRS. May be nullptr.
4405  * @return a new CRS demoted to 2D, or the current one if already 2D or not
4406  * applicable.
4407  * @since 6.3
4408  */
4409 ProjectedCRSNNPtr
demoteTo2D(const std::string & newName,const io::DatabaseContextPtr & dbContext) const4410 ProjectedCRS::demoteTo2D(const std::string &newName,
4411                          const io::DatabaseContextPtr &dbContext) const {
4412 
4413     const auto &axisList = coordinateSystem()->axisList();
4414     if (axisList.size() == 3) {
4415         auto cs = cs::CartesianCS::create(util::PropertyMap(), axisList[0],
4416                                           axisList[1]);
4417         const auto &l_baseCRS = baseCRS();
4418         const auto geogCRS =
4419             dynamic_cast<const GeographicCRS *>(l_baseCRS.get());
4420         const auto newBaseCRS =
4421             geogCRS ? util::nn_static_pointer_cast<GeodeticCRS>(
4422                           geogCRS->demoteTo2D(std::string(), dbContext))
4423                     : l_baseCRS;
4424         return ProjectedCRS::create(
4425             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
4426                                     !newName.empty() ? newName : nameStr()),
4427             newBaseCRS, derivingConversion(), cs);
4428     }
4429 
4430     return NN_NO_CHECK(std::dynamic_pointer_cast<ProjectedCRS>(
4431         shared_from_this().as_nullable()));
4432 }
4433 
4434 // ---------------------------------------------------------------------------
4435 
4436 //! @cond Doxygen_Suppress
4437 
4438 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & authorityFactory) const4439 ProjectedCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
4440     typedef std::pair<CRSNNPtr, int> Pair;
4441     std::list<Pair> res;
4442     auto resTemp = identify(authorityFactory);
4443     for (const auto &pair : resTemp) {
4444         res.emplace_back(pair.first, pair.second);
4445     }
4446     return res;
4447 }
4448 
4449 //! @endcond
4450 
4451 // ---------------------------------------------------------------------------
4452 
4453 //! @cond Doxygen_Suppress
InvalidCompoundCRSException(const char * message)4454 InvalidCompoundCRSException::InvalidCompoundCRSException(const char *message)
4455     : Exception(message) {}
4456 
4457 // ---------------------------------------------------------------------------
4458 
InvalidCompoundCRSException(const std::string & message)4459 InvalidCompoundCRSException::InvalidCompoundCRSException(
4460     const std::string &message)
4461     : Exception(message) {}
4462 
4463 // ---------------------------------------------------------------------------
4464 
4465 InvalidCompoundCRSException::~InvalidCompoundCRSException() = default;
4466 
4467 // ---------------------------------------------------------------------------
4468 
4469 InvalidCompoundCRSException::InvalidCompoundCRSException(
4470     const InvalidCompoundCRSException &) = default;
4471 //! @endcond
4472 
4473 // ---------------------------------------------------------------------------
4474 
4475 //! @cond Doxygen_Suppress
4476 struct CompoundCRS::Private {
4477     std::vector<CRSNNPtr> components_{};
4478 };
4479 //! @endcond
4480 
4481 // ---------------------------------------------------------------------------
4482 
CompoundCRS(const std::vector<CRSNNPtr> & components)4483 CompoundCRS::CompoundCRS(const std::vector<CRSNNPtr> &components)
4484     : CRS(), d(internal::make_unique<Private>()) {
4485     d->components_ = components;
4486 }
4487 
4488 // ---------------------------------------------------------------------------
4489 
CompoundCRS(const CompoundCRS & other)4490 CompoundCRS::CompoundCRS(const CompoundCRS &other)
4491     : CRS(other), d(internal::make_unique<Private>(*other.d)) {}
4492 
4493 // ---------------------------------------------------------------------------
4494 
4495 //! @cond Doxygen_Suppress
4496 CompoundCRS::~CompoundCRS() = default;
4497 //! @endcond
4498 
4499 // ---------------------------------------------------------------------------
4500 
_shallowClone() const4501 CRSNNPtr CompoundCRS::_shallowClone() const {
4502     auto crs(CompoundCRS::nn_make_shared<CompoundCRS>(*this));
4503     crs->assignSelf(crs);
4504     return crs;
4505 }
4506 
4507 // ---------------------------------------------------------------------------
4508 
4509 /** \brief Return the components of a CompoundCRS.
4510  *
4511  * @return the components.
4512  */
4513 const std::vector<CRSNNPtr> &
componentReferenceSystems()4514 CompoundCRS::componentReferenceSystems() PROJ_PURE_DEFN {
4515     return d->components_;
4516 }
4517 
4518 // ---------------------------------------------------------------------------
4519 
4520 /** \brief Instantiate a CompoundCRS from a vector of CRS.
4521  *
4522  * @param properties See \ref general_properties.
4523  * At minimum the name should be defined.
4524  * @param components the component CRS of the CompoundCRS.
4525  * @return new CompoundCRS.
4526  * @throw InvalidCompoundCRSException
4527  */
create(const util::PropertyMap & properties,const std::vector<CRSNNPtr> & components)4528 CompoundCRSNNPtr CompoundCRS::create(const util::PropertyMap &properties,
4529                                      const std::vector<CRSNNPtr> &components) {
4530 
4531     if (components.size() < 2) {
4532         throw InvalidCompoundCRSException(
4533             "compound CRS should have at least 2 components");
4534     }
4535 
4536     auto comp0 = components[0].get();
4537     auto comp0Bound = dynamic_cast<const BoundCRS *>(comp0);
4538     if (comp0Bound) {
4539         comp0 = comp0Bound->baseCRS().get();
4540     }
4541     auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0);
4542     auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0);
4543     auto comp0Eng = dynamic_cast<const EngineeringCRS *>(comp0);
4544 
4545     auto comp1 = components[1].get();
4546     auto comp1Bound = dynamic_cast<const BoundCRS *>(comp1);
4547     if (comp1Bound) {
4548         comp1 = comp1Bound->baseCRS().get();
4549     }
4550     auto comp1Vert = dynamic_cast<const VerticalCRS *>(comp1);
4551     auto comp1Eng = dynamic_cast<const EngineeringCRS *>(comp1);
4552     // Loose validation based on
4553     // http://docs.opengeospatial.org/as/18-005r4/18-005r4.html#34
4554     bool ok = false;
4555     if ((comp0Geog && comp0Geog->coordinateSystem()->axisList().size() == 2 &&
4556          (comp1Vert ||
4557           (comp1Eng &&
4558            comp1Eng->coordinateSystem()->axisList().size() == 1))) ||
4559         (comp0Proj && comp0Proj->coordinateSystem()->axisList().size() == 2 &&
4560          (comp1Vert ||
4561           (comp1Eng &&
4562            comp1Eng->coordinateSystem()->axisList().size() == 1))) ||
4563         (comp0Eng && comp0Eng->coordinateSystem()->axisList().size() <= 2 &&
4564          comp1Vert)) {
4565         // Spatial compound coordinate reference system
4566         ok = true;
4567     } else {
4568         bool isComp0Spatial = comp0Geog || comp0Proj || comp0Eng ||
4569                               dynamic_cast<const GeodeticCRS *>(comp0) ||
4570                               dynamic_cast<const VerticalCRS *>(comp0);
4571         if (isComp0Spatial && dynamic_cast<const TemporalCRS *>(comp1)) {
4572             // Spatio-temporal compound coordinate reference system
4573             ok = true;
4574         } else if (isComp0Spatial &&
4575                    dynamic_cast<const ParametricCRS *>(comp1)) {
4576             // Spatio-parametric compound coordinate reference system
4577             ok = true;
4578         }
4579     }
4580     if (!ok) {
4581         throw InvalidCompoundCRSException(
4582             "components of the compound CRS do not belong to one of the "
4583             "allowed combinations of "
4584             "http://docs.opengeospatial.org/as/18-005r4/18-005r4.html#34");
4585     }
4586 
4587     auto compoundCRS(CompoundCRS::nn_make_shared<CompoundCRS>(components));
4588     compoundCRS->assignSelf(compoundCRS);
4589     compoundCRS->setProperties(properties);
4590     if (!properties.get(common::IdentifiedObject::NAME_KEY)) {
4591         std::string name;
4592         for (const auto &crs : components) {
4593             if (!name.empty()) {
4594                 name += " + ";
4595             }
4596             const auto &l_name = crs->nameStr();
4597             if (!l_name.empty()) {
4598                 name += l_name;
4599             } else {
4600                 name += "unnamed";
4601             }
4602         }
4603         util::PropertyMap propertyName;
4604         propertyName.set(common::IdentifiedObject::NAME_KEY, name);
4605         compoundCRS->setProperties(propertyName);
4606     }
4607 
4608     return compoundCRS;
4609 }
4610 
4611 // ---------------------------------------------------------------------------
4612 
4613 //! @cond Doxygen_Suppress
4614 
4615 /** \brief Instantiate a CompoundCRS, a Geographic 3D CRS or a Projected CRS
4616  * from a vector of CRS.
4617  *
4618  * Be a bit "lax", in allowing formulations like EPSG:4326+4326 or
4619  * EPSG:32631+4326 to express Geographic 3D CRS / Projected3D CRS.
4620  *
4621  * @param properties See \ref general_properties.
4622  * At minimum the name should be defined.
4623  * @param components the component CRS of the CompoundCRS.
4624  * @return new CRS.
4625  * @throw InvalidCompoundCRSException
4626  */
createLax(const util::PropertyMap & properties,const std::vector<CRSNNPtr> & components,const io::DatabaseContextPtr & dbContext)4627 CRSNNPtr CompoundCRS::createLax(const util::PropertyMap &properties,
4628                                 const std::vector<CRSNNPtr> &components,
4629                                 const io::DatabaseContextPtr &dbContext) {
4630 
4631     if (components.size() == 2) {
4632         auto comp0 = components[0].get();
4633         auto comp1 = components[1].get();
4634         auto comp0Geog = dynamic_cast<const GeographicCRS *>(comp0);
4635         auto comp0Proj = dynamic_cast<const ProjectedCRS *>(comp0);
4636         auto comp0Bound = dynamic_cast<const BoundCRS *>(comp0);
4637         if (comp0Geog == nullptr && comp0Proj == nullptr) {
4638             if (comp0Bound) {
4639                 const auto *baseCRS = comp0Bound->baseCRS().get();
4640                 comp0Geog = dynamic_cast<const GeographicCRS *>(baseCRS);
4641                 comp0Proj = dynamic_cast<const ProjectedCRS *>(baseCRS);
4642             }
4643         }
4644         auto comp1Geog = dynamic_cast<const GeographicCRS *>(comp1);
4645         if ((comp0Geog != nullptr || comp0Proj != nullptr) &&
4646             comp1Geog != nullptr) {
4647             const auto horizGeog =
4648                 (comp0Proj != nullptr)
4649                     ? comp0Proj->baseCRS().as_nullable().get()
4650                     : comp0Geog;
4651             if (horizGeog->_isEquivalentTo(
4652                     comp1Geog->demoteTo2D(std::string(), dbContext).get())) {
4653                 return components[0]
4654                     ->promoteTo3D(std::string(), dbContext)
4655                     ->allowNonConformantWKT1Export();
4656             }
4657             throw InvalidCompoundCRSException(
4658                 "The 'vertical' geographic CRS is not equivalent to the "
4659                 "geographic CRS of the horizontal part");
4660         }
4661 
4662         // Detect a COMPD_CS whose VERT_CS is for ellipoidal heights
4663         auto comp1Vert =
4664             util::nn_dynamic_pointer_cast<VerticalCRS>(components[1]);
4665         if (comp1Vert != nullptr && comp1Vert->datum() &&
4666             comp1Vert->datum()->getWKT1DatumType() == "2002") {
4667             const auto &axis = comp1Vert->coordinateSystem()->axisList()[0];
4668             std::string name(components[0]->nameStr());
4669             if (!(axis->unit()._isEquivalentTo(
4670                       common::UnitOfMeasure::METRE,
4671                       util::IComparable::Criterion::EQUIVALENT) &&
4672                   &(axis->direction()) == &(cs::AxisDirection::UP))) {
4673                 name += " (" + comp1Vert->nameStr() + ')';
4674             }
4675             auto newVertAxis = cs::CoordinateSystemAxis::create(
4676                 util::PropertyMap().set(IdentifiedObject::NAME_KEY,
4677                                         cs::AxisName::Ellipsoidal_height),
4678                 cs::AxisAbbreviation::h, axis->direction(), axis->unit());
4679             return components[0]
4680                 ->promoteTo3D(name, dbContext, newVertAxis)
4681                 ->attachOriginalCompoundCRS(create(
4682                     properties,
4683                     comp0Bound ? std::vector<CRSNNPtr>{comp0Bound->baseCRS(),
4684                                                        components[1]}
4685                                : components));
4686         }
4687     }
4688 
4689     return create(properties, components);
4690 }
4691 //! @endcond
4692 
4693 // ---------------------------------------------------------------------------
4694 
4695 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const4696 void CompoundCRS::_exportToWKT(io::WKTFormatter *formatter) const {
4697     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
4698     const auto &l_components = componentReferenceSystems();
4699     if (!isWKT2 && formatter->useESRIDialect() && l_components.size() == 2) {
4700         l_components[0]->_exportToWKT(formatter);
4701         l_components[1]->_exportToWKT(formatter);
4702     } else {
4703         formatter->startNode(isWKT2 ? io::WKTConstants::COMPOUNDCRS
4704                                     : io::WKTConstants::COMPD_CS,
4705                              !identifiers().empty());
4706         formatter->addQuotedString(nameStr());
4707         for (const auto &crs : l_components) {
4708             crs->_exportToWKT(formatter);
4709         }
4710         ObjectUsage::baseExportToWKT(formatter);
4711         formatter->endNode();
4712     }
4713 }
4714 //! @endcond
4715 
4716 // ---------------------------------------------------------------------------
4717 
4718 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const4719 void CompoundCRS::_exportToJSON(
4720     io::JSONFormatter *formatter) const // throw(io::FormattingException)
4721 {
4722     auto writer = formatter->writer();
4723     auto objectContext(
4724         formatter->MakeObjectContext("CompoundCRS", !identifiers().empty()));
4725 
4726     writer->AddObjKey("name");
4727     auto l_name = nameStr();
4728     if (l_name.empty()) {
4729         writer->Add("unnamed");
4730     } else {
4731         writer->Add(l_name);
4732     }
4733 
4734     writer->AddObjKey("components");
4735     {
4736         auto componentsContext(writer->MakeArrayContext(false));
4737         for (const auto &crs : componentReferenceSystems()) {
4738             crs->_exportToJSON(formatter);
4739         }
4740     }
4741 
4742     ObjectUsage::baseExportToJSON(formatter);
4743 }
4744 //! @endcond
4745 
4746 // ---------------------------------------------------------------------------
4747 
_exportToPROJString(io::PROJStringFormatter * formatter) const4748 void CompoundCRS::_exportToPROJString(
4749     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
4750 {
4751     for (const auto &crs : componentReferenceSystems()) {
4752         auto crs_exportable =
4753             dynamic_cast<const IPROJStringExportable *>(crs.get());
4754         if (crs_exportable) {
4755             crs_exportable->_exportToPROJString(formatter);
4756         }
4757     }
4758 }
4759 
4760 // ---------------------------------------------------------------------------
4761 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const4762 bool CompoundCRS::_isEquivalentTo(
4763     const util::IComparable *other, util::IComparable::Criterion criterion,
4764     const io::DatabaseContextPtr &dbContext) const {
4765     auto otherCompoundCRS = dynamic_cast<const CompoundCRS *>(other);
4766     if (otherCompoundCRS == nullptr ||
4767         (criterion == util::IComparable::Criterion::STRICT &&
4768          !ObjectUsage::_isEquivalentTo(other, criterion, dbContext))) {
4769         return false;
4770     }
4771     const auto &components = componentReferenceSystems();
4772     const auto &otherComponents = otherCompoundCRS->componentReferenceSystems();
4773     if (components.size() != otherComponents.size()) {
4774         return false;
4775     }
4776     for (size_t i = 0; i < components.size(); i++) {
4777         if (!components[i]->_isEquivalentTo(otherComponents[i].get(), criterion,
4778                                             dbContext)) {
4779             return false;
4780         }
4781     }
4782     return true;
4783 }
4784 
4785 // ---------------------------------------------------------------------------
4786 
4787 /** \brief Identify the CRS with reference CRSs.
4788  *
4789  * The candidate CRSs are looked in the database when
4790  * authorityFactory is not null.
4791  *
4792  * Note that the implementation uses a set of heuristics to have a good
4793  * compromise of successful identifications over execution time. It might miss
4794  * legitimate matches in some circumstances.
4795  *
4796  * The method returns a list of matching reference CRS, and the percentage
4797  * (0-100) of confidence in the match. The list is sorted by decreasing
4798  * confidence.
4799  *
4800  * 100% means that the name of the reference entry
4801  * perfectly matches the CRS name, and both are equivalent. In which case a
4802  * single result is returned.
4803  * 90% means that CRS are equivalent, but the names are not exactly the same.
4804  * 70% means that CRS are equivalent (equivalent horizontal and vertical CRS),
4805  * but the names are not equivalent.
4806  * 25% means that the CRS are not equivalent, but there is some similarity in
4807  * the names.
4808  *
4809  * @param authorityFactory Authority factory (if null, will return an empty
4810  * list)
4811  * @return a list of matching reference CRS, and the percentage (0-100) of
4812  * confidence in the match.
4813  */
4814 std::list<std::pair<CompoundCRSNNPtr, int>>
identify(const io::AuthorityFactoryPtr & authorityFactory) const4815 CompoundCRS::identify(const io::AuthorityFactoryPtr &authorityFactory) const {
4816     typedef std::pair<CompoundCRSNNPtr, int> Pair;
4817     std::list<Pair> res;
4818 
4819     const auto &thisName(nameStr());
4820 
4821     const auto &components = componentReferenceSystems();
4822     const bool l_implicitCS = components[0]->hasImplicitCS();
4823     const auto crsCriterion =
4824         l_implicitCS
4825             ? util::IComparable::Criterion::EQUIVALENT_EXCEPT_AXIS_ORDER_GEOGCRS
4826             : util::IComparable::Criterion::EQUIVALENT;
4827 
4828     if (authorityFactory) {
4829         const io::DatabaseContextNNPtr &dbContext =
4830             authorityFactory->databaseContext();
4831 
4832         const bool unsignificantName = thisName.empty() ||
4833                                        ci_equal(thisName, "unknown") ||
4834                                        ci_equal(thisName, "unnamed");
4835         bool foundEquivalentName = false;
4836 
4837         if (hasCodeCompatibleOfAuthorityFactory(this, authorityFactory)) {
4838             // If the CRS has already an id, check in the database for the
4839             // official object, and verify that they are equivalent.
4840             for (const auto &id : identifiers()) {
4841                 if (hasCodeCompatibleOfAuthorityFactory(id, authorityFactory)) {
4842                     try {
4843                         auto crs = io::AuthorityFactory::create(
4844                                        dbContext, *id->codeSpace())
4845                                        ->createCompoundCRS(id->code());
4846                         bool match =
4847                             _isEquivalentTo(crs.get(), crsCriterion, dbContext);
4848                         res.emplace_back(crs, match ? 100 : 25);
4849                         return res;
4850                     } catch (const std::exception &) {
4851                     }
4852                 }
4853             }
4854         } else if (!unsignificantName) {
4855             for (int ipass = 0; ipass < 2; ipass++) {
4856                 const bool approximateMatch = ipass == 1;
4857                 auto objects = authorityFactory->createObjectsFromName(
4858                     thisName, {io::AuthorityFactory::ObjectType::COMPOUND_CRS},
4859                     approximateMatch);
4860                 for (const auto &obj : objects) {
4861                     auto crs = util::nn_dynamic_pointer_cast<CompoundCRS>(obj);
4862                     assert(crs);
4863                     auto crsNN = NN_NO_CHECK(crs);
4864                     const bool eqName = metadata::Identifier::isEquivalentName(
4865                         thisName.c_str(), crs->nameStr().c_str());
4866                     foundEquivalentName |= eqName;
4867                     if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) {
4868                         if (crs->nameStr() == thisName) {
4869                             res.clear();
4870                             res.emplace_back(crsNN, 100);
4871                             return res;
4872                         }
4873                         res.emplace_back(crsNN, eqName ? 90 : 70);
4874                     } else {
4875 
4876                         res.emplace_back(crsNN, 25);
4877                     }
4878                 }
4879                 if (!res.empty()) {
4880                     break;
4881                 }
4882             }
4883         }
4884 
4885         const auto lambdaSort = [&thisName](const Pair &a, const Pair &b) {
4886             // First consider confidence
4887             if (a.second > b.second) {
4888                 return true;
4889             }
4890             if (a.second < b.second) {
4891                 return false;
4892             }
4893 
4894             // Then consider exact name matching
4895             const auto &aName(a.first->nameStr());
4896             const auto &bName(b.first->nameStr());
4897             if (aName == thisName && bName != thisName) {
4898                 return true;
4899             }
4900             if (bName == thisName && aName != thisName) {
4901                 return false;
4902             }
4903 
4904             // Arbitrary final sorting criterion
4905             return aName < bName;
4906         };
4907 
4908         // Sort results
4909         res.sort(lambdaSort);
4910 
4911         if (identifiers().empty() && !foundEquivalentName &&
4912             (res.empty() || res.front().second < 50)) {
4913             std::set<std::pair<std::string, std::string>> alreadyKnown;
4914             for (const auto &pair : res) {
4915                 const auto &ids = pair.first->identifiers();
4916                 assert(!ids.empty());
4917                 alreadyKnown.insert(std::pair<std::string, std::string>(
4918                     *(ids[0]->codeSpace()), ids[0]->code()));
4919             }
4920 
4921             auto self = NN_NO_CHECK(std::dynamic_pointer_cast<CompoundCRS>(
4922                 shared_from_this().as_nullable()));
4923             auto candidates =
4924                 authorityFactory->createCompoundCRSFromExisting(self);
4925             for (const auto &crs : candidates) {
4926                 const auto &ids = crs->identifiers();
4927                 assert(!ids.empty());
4928                 if (alreadyKnown.find(std::pair<std::string, std::string>(
4929                         *(ids[0]->codeSpace()), ids[0]->code())) !=
4930                     alreadyKnown.end()) {
4931                     continue;
4932                 }
4933 
4934                 if (_isEquivalentTo(crs.get(), crsCriterion, dbContext)) {
4935                     res.emplace_back(crs, unsignificantName ? 90 : 70);
4936                 } else {
4937                     res.emplace_back(crs, 25);
4938                 }
4939             }
4940 
4941             res.sort(lambdaSort);
4942         }
4943 
4944         // If we didn't find a match for the CompoundCRS, check if the
4945         // horizontal and vertical parts are not themselves well known.
4946         if (identifiers().empty() && res.empty() && components.size() == 2) {
4947             auto candidatesHorizCRS = components[0]->identify(authorityFactory);
4948             auto candidatesVertCRS = components[1]->identify(authorityFactory);
4949             if (candidatesHorizCRS.size() == 1 &&
4950                 candidatesVertCRS.size() == 1 &&
4951                 candidatesHorizCRS.front().second >= 70 &&
4952                 candidatesVertCRS.front().second >= 70) {
4953                 auto newCRS = CompoundCRS::create(
4954                     util::PropertyMap().set(
4955                         common::IdentifiedObject::NAME_KEY,
4956                         candidatesHorizCRS.front().first->nameStr() + " + " +
4957                             candidatesVertCRS.front().first->nameStr()),
4958                     {candidatesHorizCRS.front().first,
4959                      candidatesVertCRS.front().first});
4960                 const bool eqName = metadata::Identifier::isEquivalentName(
4961                     thisName.c_str(), newCRS->nameStr().c_str());
4962                 res.emplace_back(
4963                     newCRS,
4964                     std::min(thisName == newCRS->nameStr() ? 100 : eqName ? 90
4965                                                                           : 70,
4966                              std::min(candidatesHorizCRS.front().second,
4967                                       candidatesVertCRS.front().second)));
4968             }
4969         }
4970 
4971         // Keep only results of the highest confidence
4972         if (res.size() >= 2) {
4973             const auto highestConfidence = res.front().second;
4974             std::list<Pair> newRes;
4975             for (const auto &pair : res) {
4976                 if (pair.second == highestConfidence) {
4977                     newRes.push_back(pair);
4978                 } else {
4979                     break;
4980                 }
4981             }
4982             return newRes;
4983         }
4984     }
4985 
4986     return res;
4987 }
4988 
4989 // ---------------------------------------------------------------------------
4990 
4991 //! @cond Doxygen_Suppress
4992 
4993 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & authorityFactory) const4994 CompoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
4995     typedef std::pair<CRSNNPtr, int> Pair;
4996     std::list<Pair> res;
4997     auto resTemp = identify(authorityFactory);
4998     for (const auto &pair : resTemp) {
4999         res.emplace_back(pair.first, pair.second);
5000     }
5001     return res;
5002 }
5003 
5004 //! @endcond
5005 
5006 // ---------------------------------------------------------------------------
5007 
5008 //! @cond Doxygen_Suppress
5009 struct PROJ_INTERNAL BoundCRS::Private {
5010     CRSNNPtr baseCRS_;
5011     CRSNNPtr hubCRS_;
5012     operation::TransformationNNPtr transformation_;
5013 
5014     Private(const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn,
5015             const operation::TransformationNNPtr &transformationIn);
5016 
baseCRScrs::BoundCRS::Private5017     inline const CRSNNPtr &baseCRS() const { return baseCRS_; }
hubCRScrs::BoundCRS::Private5018     inline const CRSNNPtr &hubCRS() const { return hubCRS_; }
transformationcrs::BoundCRS::Private5019     inline const operation::TransformationNNPtr &transformation() const {
5020         return transformation_;
5021     }
5022 };
5023 
Private(const CRSNNPtr & baseCRSIn,const CRSNNPtr & hubCRSIn,const operation::TransformationNNPtr & transformationIn)5024 BoundCRS::Private::Private(
5025     const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn,
5026     const operation::TransformationNNPtr &transformationIn)
5027     : baseCRS_(baseCRSIn), hubCRS_(hubCRSIn),
5028       transformation_(transformationIn) {}
5029 
5030 //! @endcond
5031 
5032 // ---------------------------------------------------------------------------
5033 
BoundCRS(const CRSNNPtr & baseCRSIn,const CRSNNPtr & hubCRSIn,const operation::TransformationNNPtr & transformationIn)5034 BoundCRS::BoundCRS(const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn,
5035                    const operation::TransformationNNPtr &transformationIn)
5036     : d(internal::make_unique<Private>(baseCRSIn, hubCRSIn, transformationIn)) {
5037 }
5038 
5039 // ---------------------------------------------------------------------------
5040 
BoundCRS(const BoundCRS & other)5041 BoundCRS::BoundCRS(const BoundCRS &other)
5042     : CRS(other),
5043       d(internal::make_unique<Private>(other.d->baseCRS(), other.d->hubCRS(),
5044                                        other.d->transformation())) {}
5045 
5046 // ---------------------------------------------------------------------------
5047 
5048 //! @cond Doxygen_Suppress
5049 BoundCRS::~BoundCRS() = default;
5050 //! @endcond
5051 
5052 // ---------------------------------------------------------------------------
5053 
shallowCloneAsBoundCRS() const5054 BoundCRSNNPtr BoundCRS::shallowCloneAsBoundCRS() const {
5055     auto crs(BoundCRS::nn_make_shared<BoundCRS>(*this));
5056     crs->assignSelf(crs);
5057     return crs;
5058 }
5059 
5060 // ---------------------------------------------------------------------------
5061 
_shallowClone() const5062 CRSNNPtr BoundCRS::_shallowClone() const { return shallowCloneAsBoundCRS(); }
5063 
5064 // ---------------------------------------------------------------------------
5065 
5066 /** \brief Return the base CRS.
5067  *
5068  * This is the CRS into which coordinates of the BoundCRS are expressed.
5069  *
5070  * @return the base CRS.
5071  */
baseCRS()5072 const CRSNNPtr &BoundCRS::baseCRS() PROJ_PURE_DEFN { return d->baseCRS_; }
5073 
5074 // ---------------------------------------------------------------------------
5075 
5076 // The only legit caller is BoundCRS::baseCRSWithCanonicalBoundCRS()
setCanonicalBoundCRS(const BoundCRSNNPtr & boundCRS)5077 void CRS::setCanonicalBoundCRS(const BoundCRSNNPtr &boundCRS) {
5078 
5079     d->canonicalBoundCRS_ = boundCRS;
5080 }
5081 
5082 // ---------------------------------------------------------------------------
5083 
5084 /** \brief Return a shallow clone of the base CRS that points to a
5085  * shallow clone of this BoundCRS.
5086  *
5087  * The base CRS is the CRS into which coordinates of the BoundCRS are expressed.
5088  *
5089  * The returned CRS will actually be a shallow clone of the actual base CRS,
5090  * with the extra property that CRS::canonicalBoundCRS() will point to a
5091  * shallow clone of this BoundCRS. Use this only if you want to work with
5092  * the base CRS object rather than the BoundCRS, but wanting to be able to
5093  * retrieve the BoundCRS later.
5094  *
5095  * @return the base CRS.
5096  */
baseCRSWithCanonicalBoundCRS() const5097 CRSNNPtr BoundCRS::baseCRSWithCanonicalBoundCRS() const {
5098     auto baseCRSClone = baseCRS()->_shallowClone();
5099     baseCRSClone->setCanonicalBoundCRS(shallowCloneAsBoundCRS());
5100     return baseCRSClone;
5101 }
5102 
5103 // ---------------------------------------------------------------------------
5104 
5105 /** \brief Return the target / hub CRS.
5106  *
5107  * @return the hub CRS.
5108  */
hubCRS()5109 const CRSNNPtr &BoundCRS::hubCRS() PROJ_PURE_DEFN { return d->hubCRS_; }
5110 
5111 // ---------------------------------------------------------------------------
5112 
5113 /** \brief Return the transformation to the hub RS.
5114  *
5115  * @return transformation.
5116  */
5117 const operation::TransformationNNPtr &
transformation()5118 BoundCRS::transformation() PROJ_PURE_DEFN {
5119     return d->transformation_;
5120 }
5121 
5122 // ---------------------------------------------------------------------------
5123 
5124 /** \brief Instantiate a BoundCRS from a base CRS, a hub CRS and a
5125  * transformation.
5126  *
5127  * @param baseCRSIn base CRS.
5128  * @param hubCRSIn hub CRS.
5129  * @param transformationIn transformation from base CRS to hub CRS.
5130  * @return new BoundCRS.
5131  */
5132 BoundCRSNNPtr
create(const CRSNNPtr & baseCRSIn,const CRSNNPtr & hubCRSIn,const operation::TransformationNNPtr & transformationIn)5133 BoundCRS::create(const CRSNNPtr &baseCRSIn, const CRSNNPtr &hubCRSIn,
5134                  const operation::TransformationNNPtr &transformationIn) {
5135     auto crs = BoundCRS::nn_make_shared<BoundCRS>(baseCRSIn, hubCRSIn,
5136                                                   transformationIn);
5137     crs->assignSelf(crs);
5138     const auto &l_name = baseCRSIn->nameStr();
5139     if (!l_name.empty()) {
5140         crs->setProperties(util::PropertyMap().set(
5141             common::IdentifiedObject::NAME_KEY, l_name));
5142     }
5143     return crs;
5144 }
5145 
5146 // ---------------------------------------------------------------------------
5147 
5148 /** \brief Instantiate a BoundCRS from a base CRS and TOWGS84 parameters
5149  *
5150  * @param baseCRSIn base CRS.
5151  * @param TOWGS84Parameters a vector of 3 or 7 double values representing WKT1
5152  * TOWGS84 parameter.
5153  * @return new BoundCRS.
5154  */
5155 BoundCRSNNPtr
createFromTOWGS84(const CRSNNPtr & baseCRSIn,const std::vector<double> & TOWGS84Parameters)5156 BoundCRS::createFromTOWGS84(const CRSNNPtr &baseCRSIn,
5157                             const std::vector<double> &TOWGS84Parameters) {
5158     auto transf =
5159         operation::Transformation::createTOWGS84(baseCRSIn, TOWGS84Parameters);
5160     return create(baseCRSIn, transf->targetCRS(), transf);
5161 }
5162 
5163 // ---------------------------------------------------------------------------
5164 
5165 /** \brief Instantiate a BoundCRS from a base CRS and nadgrids parameters
5166  *
5167  * @param baseCRSIn base CRS.
5168  * @param filename Horizontal grid filename
5169  * @return new BoundCRS.
5170  */
createFromNadgrids(const CRSNNPtr & baseCRSIn,const std::string & filename)5171 BoundCRSNNPtr BoundCRS::createFromNadgrids(const CRSNNPtr &baseCRSIn,
5172                                            const std::string &filename) {
5173     const auto sourceGeographicCRS = baseCRSIn->extractGeographicCRS();
5174     auto transformationSourceCRS =
5175         sourceGeographicCRS
5176             ? NN_NO_CHECK(std::static_pointer_cast<CRS>(sourceGeographicCRS))
5177             : baseCRSIn;
5178     if (sourceGeographicCRS != nullptr &&
5179         sourceGeographicCRS->primeMeridian()->longitude().getSIValue() != 0.0) {
5180         transformationSourceCRS = GeographicCRS::create(
5181             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
5182                                     sourceGeographicCRS->nameStr() +
5183                                         " (with Greenwich prime meridian)"),
5184             datum::GeodeticReferenceFrame::create(
5185                 util::PropertyMap().set(
5186                     common::IdentifiedObject::NAME_KEY,
5187                     sourceGeographicCRS->datumNonNull(nullptr)->nameStr() +
5188                         " (with Greenwich prime meridian)"),
5189                 sourceGeographicCRS->datumNonNull(nullptr)->ellipsoid(),
5190                 util::optional<std::string>(), datum::PrimeMeridian::GREENWICH),
5191             cs::EllipsoidalCS::createLatitudeLongitude(
5192                 common::UnitOfMeasure::DEGREE));
5193     }
5194     std::string transformationName = transformationSourceCRS->nameStr();
5195     transformationName += " to WGS84";
5196 
5197     return create(
5198         baseCRSIn, GeographicCRS::EPSG_4326,
5199         operation::Transformation::createNTv2(
5200             util::PropertyMap().set(common::IdentifiedObject::NAME_KEY,
5201                                     transformationName),
5202             transformationSourceCRS, GeographicCRS::EPSG_4326, filename,
5203             std::vector<metadata::PositionalAccuracyNNPtr>()));
5204 }
5205 
5206 // ---------------------------------------------------------------------------
5207 
isTOWGS84Compatible() const5208 bool BoundCRS::isTOWGS84Compatible() const {
5209     return dynamic_cast<GeodeticCRS *>(d->hubCRS().get()) != nullptr &&
5210            ci_equal(d->hubCRS()->nameStr(), "WGS 84");
5211 }
5212 
5213 // ---------------------------------------------------------------------------
5214 
getHDatumPROJ4GRIDS() const5215 std::string BoundCRS::getHDatumPROJ4GRIDS() const {
5216     if (ci_equal(d->hubCRS()->nameStr(), "WGS 84")) {
5217         return d->transformation()->getNTv2Filename();
5218     }
5219     return std::string();
5220 }
5221 
5222 // ---------------------------------------------------------------------------
5223 
getVDatumPROJ4GRIDS() const5224 std::string BoundCRS::getVDatumPROJ4GRIDS() const {
5225     if (dynamic_cast<VerticalCRS *>(d->baseCRS().get()) &&
5226         ci_equal(d->hubCRS()->nameStr(), "WGS 84")) {
5227         return d->transformation()->getHeightToGeographic3DFilename();
5228     }
5229     return std::string();
5230 }
5231 
5232 // ---------------------------------------------------------------------------
5233 
5234 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const5235 void BoundCRS::_exportToWKT(io::WKTFormatter *formatter) const {
5236     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
5237     if (isWKT2) {
5238         formatter->startNode(io::WKTConstants::BOUNDCRS, false);
5239         formatter->startNode(io::WKTConstants::SOURCECRS, false);
5240         d->baseCRS()->_exportToWKT(formatter);
5241         formatter->endNode();
5242         formatter->startNode(io::WKTConstants::TARGETCRS, false);
5243         d->hubCRS()->_exportToWKT(formatter);
5244         formatter->endNode();
5245         formatter->setAbridgedTransformation(true);
5246         d->transformation()->_exportToWKT(formatter);
5247         formatter->setAbridgedTransformation(false);
5248         formatter->endNode();
5249     } else {
5250 
5251         auto vdatumProj4GridName = getVDatumPROJ4GRIDS();
5252         if (!vdatumProj4GridName.empty()) {
5253             formatter->setVDatumExtension(vdatumProj4GridName);
5254             d->baseCRS()->_exportToWKT(formatter);
5255             formatter->setVDatumExtension(std::string());
5256             return;
5257         }
5258 
5259         auto hdatumProj4GridName = getHDatumPROJ4GRIDS();
5260         if (!hdatumProj4GridName.empty()) {
5261             formatter->setHDatumExtension(hdatumProj4GridName);
5262             d->baseCRS()->_exportToWKT(formatter);
5263             formatter->setHDatumExtension(std::string());
5264             return;
5265         }
5266 
5267         if (!isTOWGS84Compatible()) {
5268             io::FormattingException::Throw(
5269                 "Cannot export BoundCRS with non-WGS 84 hub CRS in WKT1");
5270         }
5271         auto params = d->transformation()->getTOWGS84Parameters();
5272         if (!formatter->useESRIDialect()) {
5273             formatter->setTOWGS84Parameters(params);
5274         }
5275         d->baseCRS()->_exportToWKT(formatter);
5276         formatter->setTOWGS84Parameters(std::vector<double>());
5277     }
5278 }
5279 //! @endcond
5280 
5281 // ---------------------------------------------------------------------------
5282 
5283 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const5284 void BoundCRS::_exportToJSON(
5285     io::JSONFormatter *formatter) const // throw(io::FormattingException)
5286 {
5287     auto writer = formatter->writer();
5288     auto objectContext(
5289         formatter->MakeObjectContext("BoundCRS", !identifiers().empty()));
5290 
5291     writer->AddObjKey("source_crs");
5292     d->baseCRS()->_exportToJSON(formatter);
5293 
5294     writer->AddObjKey("target_crs");
5295     d->hubCRS()->_exportToJSON(formatter);
5296 
5297     writer->AddObjKey("transformation");
5298     formatter->setOmitTypeInImmediateChild();
5299     formatter->setAbridgedTransformation(true);
5300     d->transformation()->_exportToJSON(formatter);
5301     formatter->setAbridgedTransformation(false);
5302 }
5303 //! @endcond
5304 
5305 // ---------------------------------------------------------------------------
5306 
_exportToPROJString(io::PROJStringFormatter * formatter) const5307 void BoundCRS::_exportToPROJString(
5308     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
5309 {
5310     auto crs_exportable =
5311         dynamic_cast<const io::IPROJStringExportable *>(d->baseCRS_.get());
5312     if (!crs_exportable) {
5313         io::FormattingException::Throw(
5314             "baseCRS of BoundCRS cannot be exported as a PROJ string");
5315     }
5316 
5317     auto vdatumProj4GridName = getVDatumPROJ4GRIDS();
5318     if (!vdatumProj4GridName.empty()) {
5319         formatter->setVDatumExtension(vdatumProj4GridName);
5320         crs_exportable->_exportToPROJString(formatter);
5321         formatter->setVDatumExtension(std::string());
5322     } else {
5323         auto hdatumProj4GridName = getHDatumPROJ4GRIDS();
5324         if (!hdatumProj4GridName.empty()) {
5325             formatter->setHDatumExtension(hdatumProj4GridName);
5326             crs_exportable->_exportToPROJString(formatter);
5327             formatter->setHDatumExtension(std::string());
5328         } else {
5329             if (isTOWGS84Compatible()) {
5330                 auto params = transformation()->getTOWGS84Parameters();
5331                 formatter->setTOWGS84Parameters(params);
5332             }
5333             crs_exportable->_exportToPROJString(formatter);
5334             formatter->setTOWGS84Parameters(std::vector<double>());
5335         }
5336     }
5337 }
5338 
5339 // ---------------------------------------------------------------------------
5340 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const5341 bool BoundCRS::_isEquivalentTo(const util::IComparable *other,
5342                                util::IComparable::Criterion criterion,
5343                                const io::DatabaseContextPtr &dbContext) const {
5344     auto otherBoundCRS = dynamic_cast<const BoundCRS *>(other);
5345     if (otherBoundCRS == nullptr ||
5346         (criterion == util::IComparable::Criterion::STRICT &&
5347          !ObjectUsage::_isEquivalentTo(other, criterion, dbContext))) {
5348         return false;
5349     }
5350     const auto standardCriterion = getStandardCriterion(criterion);
5351     return d->baseCRS_->_isEquivalentTo(otherBoundCRS->d->baseCRS_.get(),
5352                                         criterion, dbContext) &&
5353            d->hubCRS_->_isEquivalentTo(otherBoundCRS->d->hubCRS_.get(),
5354                                        criterion, dbContext) &&
5355            d->transformation_->_isEquivalentTo(
5356                otherBoundCRS->d->transformation_.get(), standardCriterion,
5357                dbContext);
5358 }
5359 
5360 // ---------------------------------------------------------------------------
5361 
5362 //! @cond Doxygen_Suppress
5363 
5364 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & authorityFactory) const5365 BoundCRS::_identify(const io::AuthorityFactoryPtr &authorityFactory) const {
5366     typedef std::pair<CRSNNPtr, int> Pair;
5367     std::list<Pair> res;
5368     if (!authorityFactory)
5369         return res;
5370     std::list<Pair> resMatchOfTransfToWGS84;
5371     const io::DatabaseContextNNPtr &dbContext =
5372         authorityFactory->databaseContext();
5373     if (d->hubCRS_->_isEquivalentTo(GeographicCRS::EPSG_4326.get(),
5374                                     util::IComparable::Criterion::EQUIVALENT,
5375                                     dbContext)) {
5376         auto resTemp = d->baseCRS_->identify(authorityFactory);
5377 
5378         std::string refTransfPROJString;
5379         bool refTransfPROJStringValid = false;
5380         auto refTransf = d->transformation_->normalizeForVisualization();
5381         try {
5382             refTransfPROJString = refTransf->exportToPROJString(
5383                 io::PROJStringFormatter::create().get());
5384             refTransfPROJString = replaceAll(
5385                 refTransfPROJString,
5386                 " +rx=0 +ry=0 +rz=0 +s=0 +convention=position_vector", "");
5387             refTransfPROJStringValid = true;
5388         } catch (const std::exception &) {
5389         }
5390         bool refIsNullTransform = false;
5391         if (isTOWGS84Compatible()) {
5392             auto params = transformation()->getTOWGS84Parameters();
5393             if (params == std::vector<double>{0, 0, 0, 0, 0, 0, 0}) {
5394                 refIsNullTransform = true;
5395             }
5396         }
5397 
5398         for (const auto &pair : resTemp) {
5399             const auto &candidateBaseCRS = pair.first;
5400             auto projCRS =
5401                 dynamic_cast<const ProjectedCRS *>(candidateBaseCRS.get());
5402             auto geodCRS = projCRS ? projCRS->baseCRS().as_nullable()
5403                                    : util::nn_dynamic_pointer_cast<GeodeticCRS>(
5404                                          candidateBaseCRS);
5405             if (geodCRS) {
5406                 auto context = operation::CoordinateOperationContext::create(
5407                     authorityFactory, nullptr, 0.0);
5408                 context->setSpatialCriterion(
5409                     operation::CoordinateOperationContext::SpatialCriterion::
5410                         PARTIAL_INTERSECTION);
5411                 auto ops =
5412                     operation::CoordinateOperationFactory::create()
5413                         ->createOperations(NN_NO_CHECK(geodCRS),
5414                                            GeographicCRS::EPSG_4326, context);
5415 
5416                 bool foundOp = false;
5417                 for (const auto &op : ops) {
5418                     auto opNormalized = op->normalizeForVisualization();
5419                     std::string opTransfPROJString;
5420                     bool opTransfPROJStringValid = false;
5421                     if (op->nameStr().find("Ballpark geographic") == 0) {
5422                         if (refIsNullTransform) {
5423                             res.emplace_back(create(candidateBaseCRS,
5424                                                     d->hubCRS_,
5425                                                     transformation()),
5426                                              pair.second);
5427                             foundOp = true;
5428                             break;
5429                         }
5430                         continue;
5431                     }
5432                     try {
5433                         opTransfPROJString = opNormalized->exportToPROJString(
5434                             io::PROJStringFormatter::create().get());
5435                         opTransfPROJStringValid = true;
5436                         opTransfPROJString = replaceAll(
5437                             opTransfPROJString, " +rx=0 +ry=0 +rz=0 +s=0 "
5438                                                 "+convention=position_vector",
5439                             "");
5440                     } catch (const std::exception &) {
5441                     }
5442                     if ((refTransfPROJStringValid && opTransfPROJStringValid &&
5443                          refTransfPROJString == opTransfPROJString) ||
5444                         opNormalized->_isEquivalentTo(
5445                             refTransf.get(),
5446                             util::IComparable::Criterion::EQUIVALENT,
5447                             dbContext)) {
5448                         resMatchOfTransfToWGS84.emplace_back(
5449                             create(candidateBaseCRS, d->hubCRS_,
5450                                    NN_NO_CHECK(util::nn_dynamic_pointer_cast<
5451                                                operation::Transformation>(op))),
5452                             pair.second);
5453                         foundOp = true;
5454                         break;
5455                     }
5456                 }
5457                 if (!foundOp) {
5458                     res.emplace_back(
5459                         create(candidateBaseCRS, d->hubCRS_, transformation()),
5460                         std::min(70, pair.second));
5461                 }
5462             }
5463         }
5464     }
5465     return !resMatchOfTransfToWGS84.empty() ? resMatchOfTransfToWGS84 : res;
5466 }
5467 
5468 // ---------------------------------------------------------------------------
5469 
5470 //! @cond Doxygen_Suppress
5471 struct DerivedGeodeticCRS::Private {};
5472 //! @endcond
5473 
5474 // ---------------------------------------------------------------------------
5475 
5476 //! @cond Doxygen_Suppress
5477 DerivedGeodeticCRS::~DerivedGeodeticCRS() = default;
5478 //! @endcond
5479 
5480 // ---------------------------------------------------------------------------
5481 
DerivedGeodeticCRS(const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CartesianCSNNPtr & csIn)5482 DerivedGeodeticCRS::DerivedGeodeticCRS(
5483     const GeodeticCRSNNPtr &baseCRSIn,
5484     const operation::ConversionNNPtr &derivingConversionIn,
5485     const cs::CartesianCSNNPtr &csIn)
5486     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5487       GeodeticCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5488       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
5489 
5490 // ---------------------------------------------------------------------------
5491 
DerivedGeodeticCRS(const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::SphericalCSNNPtr & csIn)5492 DerivedGeodeticCRS::DerivedGeodeticCRS(
5493     const GeodeticCRSNNPtr &baseCRSIn,
5494     const operation::ConversionNNPtr &derivingConversionIn,
5495     const cs::SphericalCSNNPtr &csIn)
5496     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5497       GeodeticCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5498       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
5499 
5500 // ---------------------------------------------------------------------------
5501 
DerivedGeodeticCRS(const DerivedGeodeticCRS & other)5502 DerivedGeodeticCRS::DerivedGeodeticCRS(const DerivedGeodeticCRS &other)
5503     : SingleCRS(other), GeodeticCRS(other), DerivedCRS(other), d(nullptr) {}
5504 
5505 // ---------------------------------------------------------------------------
5506 
_shallowClone() const5507 CRSNNPtr DerivedGeodeticCRS::_shallowClone() const {
5508     auto crs(DerivedGeodeticCRS::nn_make_shared<DerivedGeodeticCRS>(*this));
5509     crs->assignSelf(crs);
5510     crs->setDerivingConversionCRS();
5511     return crs;
5512 }
5513 
5514 // ---------------------------------------------------------------------------
5515 
5516 /** \brief Return the base CRS (a GeodeticCRS) of a DerivedGeodeticCRS.
5517  *
5518  * @return the base CRS.
5519  */
baseCRS() const5520 const GeodeticCRSNNPtr DerivedGeodeticCRS::baseCRS() const {
5521     return NN_NO_CHECK(util::nn_dynamic_pointer_cast<GeodeticCRS>(
5522         DerivedCRS::getPrivate()->baseCRS_));
5523 }
5524 
5525 // ---------------------------------------------------------------------------
5526 
5527 /** \brief Instantiate a DerivedGeodeticCRS from a base CRS, a deriving
5528  * conversion and a cs::CartesianCS.
5529  *
5530  * @param properties See \ref general_properties.
5531  * At minimum the name should be defined.
5532  * @param baseCRSIn base CRS.
5533  * @param derivingConversionIn the deriving conversion from the base CRS to this
5534  * CRS.
5535  * @param csIn the coordinate system.
5536  * @return new DerivedGeodeticCRS.
5537  */
create(const util::PropertyMap & properties,const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CartesianCSNNPtr & csIn)5538 DerivedGeodeticCRSNNPtr DerivedGeodeticCRS::create(
5539     const util::PropertyMap &properties, const GeodeticCRSNNPtr &baseCRSIn,
5540     const operation::ConversionNNPtr &derivingConversionIn,
5541     const cs::CartesianCSNNPtr &csIn) {
5542     auto crs(DerivedGeodeticCRS::nn_make_shared<DerivedGeodeticCRS>(
5543         baseCRSIn, derivingConversionIn, csIn));
5544     crs->assignSelf(crs);
5545     crs->setProperties(properties);
5546     crs->setDerivingConversionCRS();
5547     return crs;
5548 }
5549 
5550 // ---------------------------------------------------------------------------
5551 
5552 /** \brief Instantiate a DerivedGeodeticCRS from a base CRS, a deriving
5553  * conversion and a cs::SphericalCS.
5554  *
5555  * @param properties See \ref general_properties.
5556  * At minimum the name should be defined.
5557  * @param baseCRSIn base CRS.
5558  * @param derivingConversionIn the deriving conversion from the base CRS to this
5559  * CRS.
5560  * @param csIn the coordinate system.
5561  * @return new DerivedGeodeticCRS.
5562  */
create(const util::PropertyMap & properties,const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::SphericalCSNNPtr & csIn)5563 DerivedGeodeticCRSNNPtr DerivedGeodeticCRS::create(
5564     const util::PropertyMap &properties, const GeodeticCRSNNPtr &baseCRSIn,
5565     const operation::ConversionNNPtr &derivingConversionIn,
5566     const cs::SphericalCSNNPtr &csIn) {
5567     auto crs(DerivedGeodeticCRS::nn_make_shared<DerivedGeodeticCRS>(
5568         baseCRSIn, derivingConversionIn, csIn));
5569     crs->assignSelf(crs);
5570     crs->setProperties(properties);
5571     crs->setDerivingConversionCRS();
5572     return crs;
5573 }
5574 
5575 // ---------------------------------------------------------------------------
5576 
5577 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const5578 void DerivedGeodeticCRS::_exportToWKT(io::WKTFormatter *formatter) const {
5579     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
5580     if (!isWKT2) {
5581         io::FormattingException::Throw(
5582             "DerivedGeodeticCRS can only be exported to WKT2");
5583     }
5584     formatter->startNode(io::WKTConstants::GEODCRS, !identifiers().empty());
5585     formatter->addQuotedString(nameStr());
5586 
5587     auto l_baseCRS = baseCRS();
5588     formatter->startNode((formatter->use2019Keywords() &&
5589                           dynamic_cast<const GeographicCRS *>(l_baseCRS.get()))
5590                              ? io::WKTConstants::BASEGEOGCRS
5591                              : io::WKTConstants::BASEGEODCRS,
5592                          !baseCRS()->identifiers().empty());
5593     formatter->addQuotedString(l_baseCRS->nameStr());
5594     auto l_datum = l_baseCRS->datum();
5595     if (l_datum) {
5596         l_datum->_exportToWKT(formatter);
5597     } else {
5598         auto l_datumEnsemble = datumEnsemble();
5599         assert(l_datumEnsemble);
5600         l_datumEnsemble->_exportToWKT(formatter);
5601     }
5602     l_baseCRS->primeMeridian()->_exportToWKT(formatter);
5603     formatter->endNode();
5604 
5605     formatter->setUseDerivingConversion(true);
5606     derivingConversionRef()->_exportToWKT(formatter);
5607     formatter->setUseDerivingConversion(false);
5608 
5609     coordinateSystem()->_exportToWKT(formatter);
5610     ObjectUsage::baseExportToWKT(formatter);
5611     formatter->endNode();
5612 }
5613 //! @endcond
5614 
5615 // ---------------------------------------------------------------------------
5616 
_exportToPROJString(io::PROJStringFormatter *) const5617 void DerivedGeodeticCRS::_exportToPROJString(
5618     io::PROJStringFormatter *) const // throw(io::FormattingException)
5619 {
5620     throw io::FormattingException(
5621         "DerivedGeodeticCRS cannot be exported to PROJ string");
5622 }
5623 
5624 // ---------------------------------------------------------------------------
5625 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const5626 bool DerivedGeodeticCRS::_isEquivalentTo(
5627     const util::IComparable *other, util::IComparable::Criterion criterion,
5628     const io::DatabaseContextPtr &dbContext) const {
5629     auto otherDerivedCRS = dynamic_cast<const DerivedGeodeticCRS *>(other);
5630     return otherDerivedCRS != nullptr &&
5631            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
5632 }
5633 
5634 // ---------------------------------------------------------------------------
5635 
5636 //! @cond Doxygen_Suppress
5637 
5638 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & factory) const5639 DerivedGeodeticCRS::_identify(const io::AuthorityFactoryPtr &factory) const {
5640     return CRS::_identify(factory);
5641 }
5642 
5643 //! @endcond
5644 
5645 // ---------------------------------------------------------------------------
5646 
5647 //! @cond Doxygen_Suppress
5648 struct DerivedGeographicCRS::Private {};
5649 //! @endcond
5650 
5651 // ---------------------------------------------------------------------------
5652 
5653 //! @cond Doxygen_Suppress
5654 DerivedGeographicCRS::~DerivedGeographicCRS() = default;
5655 //! @endcond
5656 
5657 // ---------------------------------------------------------------------------
5658 
DerivedGeographicCRS(const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::EllipsoidalCSNNPtr & csIn)5659 DerivedGeographicCRS::DerivedGeographicCRS(
5660     const GeodeticCRSNNPtr &baseCRSIn,
5661     const operation::ConversionNNPtr &derivingConversionIn,
5662     const cs::EllipsoidalCSNNPtr &csIn)
5663     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5664       GeographicCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5665       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
5666 
5667 // ---------------------------------------------------------------------------
5668 
DerivedGeographicCRS(const DerivedGeographicCRS & other)5669 DerivedGeographicCRS::DerivedGeographicCRS(const DerivedGeographicCRS &other)
5670     : SingleCRS(other), GeographicCRS(other), DerivedCRS(other), d(nullptr) {}
5671 
5672 // ---------------------------------------------------------------------------
5673 
_shallowClone() const5674 CRSNNPtr DerivedGeographicCRS::_shallowClone() const {
5675     auto crs(DerivedGeographicCRS::nn_make_shared<DerivedGeographicCRS>(*this));
5676     crs->assignSelf(crs);
5677     crs->setDerivingConversionCRS();
5678     return crs;
5679 }
5680 
5681 // ---------------------------------------------------------------------------
5682 
5683 /** \brief Return the base CRS (a GeodeticCRS) of a DerivedGeographicCRS.
5684  *
5685  * @return the base CRS.
5686  */
baseCRS() const5687 const GeodeticCRSNNPtr DerivedGeographicCRS::baseCRS() const {
5688     return NN_NO_CHECK(util::nn_dynamic_pointer_cast<GeodeticCRS>(
5689         DerivedCRS::getPrivate()->baseCRS_));
5690 }
5691 
5692 // ---------------------------------------------------------------------------
5693 
5694 /** \brief Instantiate a DerivedGeographicCRS from a base CRS, a deriving
5695  * conversion and a cs::EllipsoidalCS.
5696  *
5697  * @param properties See \ref general_properties.
5698  * At minimum the name should be defined.
5699  * @param baseCRSIn base CRS.
5700  * @param derivingConversionIn the deriving conversion from the base CRS to this
5701  * CRS.
5702  * @param csIn the coordinate system.
5703  * @return new DerivedGeographicCRS.
5704  */
create(const util::PropertyMap & properties,const GeodeticCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::EllipsoidalCSNNPtr & csIn)5705 DerivedGeographicCRSNNPtr DerivedGeographicCRS::create(
5706     const util::PropertyMap &properties, const GeodeticCRSNNPtr &baseCRSIn,
5707     const operation::ConversionNNPtr &derivingConversionIn,
5708     const cs::EllipsoidalCSNNPtr &csIn) {
5709     auto crs(DerivedGeographicCRS::nn_make_shared<DerivedGeographicCRS>(
5710         baseCRSIn, derivingConversionIn, csIn));
5711     crs->assignSelf(crs);
5712     crs->setProperties(properties);
5713     crs->setDerivingConversionCRS();
5714     return crs;
5715 }
5716 
5717 // ---------------------------------------------------------------------------
5718 
5719 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const5720 void DerivedGeographicCRS::_exportToWKT(io::WKTFormatter *formatter) const {
5721     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
5722     if (!isWKT2) {
5723         io::FormattingException::Throw(
5724             "DerivedGeographicCRS can only be exported to WKT2");
5725     }
5726     formatter->startNode(formatter->use2019Keywords()
5727                              ? io::WKTConstants::GEOGCRS
5728                              : io::WKTConstants::GEODCRS,
5729                          !identifiers().empty());
5730     formatter->addQuotedString(nameStr());
5731 
5732     auto l_baseCRS = baseCRS();
5733     formatter->startNode((formatter->use2019Keywords() &&
5734                           dynamic_cast<const GeographicCRS *>(l_baseCRS.get()))
5735                              ? io::WKTConstants::BASEGEOGCRS
5736                              : io::WKTConstants::BASEGEODCRS,
5737                          !l_baseCRS->identifiers().empty());
5738     formatter->addQuotedString(l_baseCRS->nameStr());
5739     l_baseCRS->exportDatumOrDatumEnsembleToWkt(formatter);
5740     l_baseCRS->primeMeridian()->_exportToWKT(formatter);
5741     formatter->endNode();
5742 
5743     formatter->setUseDerivingConversion(true);
5744     derivingConversionRef()->_exportToWKT(formatter);
5745     formatter->setUseDerivingConversion(false);
5746 
5747     coordinateSystem()->_exportToWKT(formatter);
5748     ObjectUsage::baseExportToWKT(formatter);
5749     formatter->endNode();
5750 }
5751 //! @endcond
5752 
5753 // ---------------------------------------------------------------------------
5754 
_exportToPROJString(io::PROJStringFormatter * formatter) const5755 void DerivedGeographicCRS::_exportToPROJString(
5756     io::PROJStringFormatter *formatter) const // throw(io::FormattingException)
5757 {
5758     const auto &l_conv = derivingConversionRef();
5759     const auto &methodName = l_conv->method()->nameStr();
5760     if (methodName == "PROJ ob_tran o_proj=longlat" ||
5761         methodName == "PROJ ob_tran o_proj=lonlat" ||
5762         methodName == "PROJ ob_tran o_proj=latlong" ||
5763         methodName == "PROJ ob_tran o_proj=latlon" ||
5764         ci_equal(methodName,
5765                  PROJ_WKT2_NAME_METHOD_POLE_ROTATION_GRIB_CONVENTION)) {
5766         l_conv->_exportToPROJString(formatter);
5767         return;
5768     }
5769 
5770     throw io::FormattingException(
5771         "DerivedGeographicCRS cannot be exported to PROJ string");
5772 }
5773 
5774 // ---------------------------------------------------------------------------
5775 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const5776 bool DerivedGeographicCRS::_isEquivalentTo(
5777     const util::IComparable *other, util::IComparable::Criterion criterion,
5778     const io::DatabaseContextPtr &dbContext) const {
5779     auto otherDerivedCRS = dynamic_cast<const DerivedGeographicCRS *>(other);
5780     return otherDerivedCRS != nullptr &&
5781            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
5782 }
5783 
5784 // ---------------------------------------------------------------------------
5785 
5786 //! @cond Doxygen_Suppress
5787 
5788 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & factory) const5789 DerivedGeographicCRS::_identify(const io::AuthorityFactoryPtr &factory) const {
5790     return CRS::_identify(factory);
5791 }
5792 
5793 //! @endcond
5794 
5795 // ---------------------------------------------------------------------------
5796 
5797 //! @cond Doxygen_Suppress
5798 struct DerivedProjectedCRS::Private {};
5799 //! @endcond
5800 
5801 // ---------------------------------------------------------------------------
5802 
5803 //! @cond Doxygen_Suppress
5804 DerivedProjectedCRS::~DerivedProjectedCRS() = default;
5805 //! @endcond
5806 
5807 // ---------------------------------------------------------------------------
5808 
DerivedProjectedCRS(const ProjectedCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CoordinateSystemNNPtr & csIn)5809 DerivedProjectedCRS::DerivedProjectedCRS(
5810     const ProjectedCRSNNPtr &baseCRSIn,
5811     const operation::ConversionNNPtr &derivingConversionIn,
5812     const cs::CoordinateSystemNNPtr &csIn)
5813     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
5814       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
5815 
5816 // ---------------------------------------------------------------------------
5817 
DerivedProjectedCRS(const DerivedProjectedCRS & other)5818 DerivedProjectedCRS::DerivedProjectedCRS(const DerivedProjectedCRS &other)
5819     : SingleCRS(other), DerivedCRS(other), d(nullptr) {}
5820 
5821 // ---------------------------------------------------------------------------
5822 
_shallowClone() const5823 CRSNNPtr DerivedProjectedCRS::_shallowClone() const {
5824     auto crs(DerivedProjectedCRS::nn_make_shared<DerivedProjectedCRS>(*this));
5825     crs->assignSelf(crs);
5826     crs->setDerivingConversionCRS();
5827     return crs;
5828 }
5829 
5830 // ---------------------------------------------------------------------------
5831 
5832 /** \brief Return the base CRS (a ProjectedCRS) of a DerivedProjectedCRS.
5833  *
5834  * @return the base CRS.
5835  */
baseCRS() const5836 const ProjectedCRSNNPtr DerivedProjectedCRS::baseCRS() const {
5837     return NN_NO_CHECK(util::nn_dynamic_pointer_cast<ProjectedCRS>(
5838         DerivedCRS::getPrivate()->baseCRS_));
5839 }
5840 
5841 // ---------------------------------------------------------------------------
5842 
5843 /** \brief Instantiate a DerivedProjectedCRS from a base CRS, a deriving
5844  * conversion and a cs::CS.
5845  *
5846  * @param properties See \ref general_properties.
5847  * At minimum the name should be defined.
5848  * @param baseCRSIn base CRS.
5849  * @param derivingConversionIn the deriving conversion from the base CRS to this
5850  * CRS.
5851  * @param csIn the coordinate system.
5852  * @return new DerivedProjectedCRS.
5853  */
create(const util::PropertyMap & properties,const ProjectedCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::CoordinateSystemNNPtr & csIn)5854 DerivedProjectedCRSNNPtr DerivedProjectedCRS::create(
5855     const util::PropertyMap &properties, const ProjectedCRSNNPtr &baseCRSIn,
5856     const operation::ConversionNNPtr &derivingConversionIn,
5857     const cs::CoordinateSystemNNPtr &csIn) {
5858     auto crs(DerivedProjectedCRS::nn_make_shared<DerivedProjectedCRS>(
5859         baseCRSIn, derivingConversionIn, csIn));
5860     crs->assignSelf(crs);
5861     crs->setProperties(properties);
5862     crs->setDerivingConversionCRS();
5863     return crs;
5864 }
5865 
5866 // ---------------------------------------------------------------------------
5867 
5868 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const5869 void DerivedProjectedCRS::_exportToWKT(io::WKTFormatter *formatter) const {
5870     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
5871     if (!isWKT2 || !formatter->use2019Keywords()) {
5872         io::FormattingException::Throw(
5873             "DerivedProjectedCRS can only be exported to WKT2:2019");
5874     }
5875     formatter->startNode(io::WKTConstants::DERIVEDPROJCRS,
5876                          !identifiers().empty());
5877     formatter->addQuotedString(nameStr());
5878 
5879     {
5880         auto l_baseProjCRS = baseCRS();
5881         formatter->startNode(io::WKTConstants::BASEPROJCRS,
5882                              !l_baseProjCRS->identifiers().empty());
5883         formatter->addQuotedString(l_baseProjCRS->nameStr());
5884 
5885         auto l_baseGeodCRS = l_baseProjCRS->baseCRS();
5886         auto &geodeticCRSAxisList =
5887             l_baseGeodCRS->coordinateSystem()->axisList();
5888 
5889         formatter->startNode(
5890             dynamic_cast<const GeographicCRS *>(l_baseGeodCRS.get())
5891                 ? io::WKTConstants::BASEGEOGCRS
5892                 : io::WKTConstants::BASEGEODCRS,
5893             !l_baseGeodCRS->identifiers().empty());
5894         formatter->addQuotedString(l_baseGeodCRS->nameStr());
5895         l_baseGeodCRS->exportDatumOrDatumEnsembleToWkt(formatter);
5896         // insert ellipsoidal cs unit when the units of the map
5897         // projection angular parameters are not explicitly given within those
5898         // parameters. See
5899         // http://docs.opengeospatial.org/is/12-063r5/12-063r5.html#61
5900         if (formatter->primeMeridianOrParameterUnitOmittedIfSameAsAxis() &&
5901             !geodeticCRSAxisList.empty()) {
5902             geodeticCRSAxisList[0]->unit()._exportToWKT(formatter);
5903         }
5904         l_baseGeodCRS->primeMeridian()->_exportToWKT(formatter);
5905         formatter->endNode();
5906 
5907         l_baseProjCRS->derivingConversionRef()->_exportToWKT(formatter);
5908         formatter->endNode();
5909     }
5910 
5911     formatter->setUseDerivingConversion(true);
5912     derivingConversionRef()->_exportToWKT(formatter);
5913     formatter->setUseDerivingConversion(false);
5914 
5915     coordinateSystem()->_exportToWKT(formatter);
5916     ObjectUsage::baseExportToWKT(formatter);
5917     formatter->endNode();
5918 }
5919 //! @endcond
5920 
5921 // ---------------------------------------------------------------------------
5922 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const5923 bool DerivedProjectedCRS::_isEquivalentTo(
5924     const util::IComparable *other, util::IComparable::Criterion criterion,
5925     const io::DatabaseContextPtr &dbContext) const {
5926     auto otherDerivedCRS = dynamic_cast<const DerivedProjectedCRS *>(other);
5927     return otherDerivedCRS != nullptr &&
5928            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
5929 }
5930 
5931 // ---------------------------------------------------------------------------
5932 
5933 //! @cond Doxygen_Suppress
5934 struct TemporalCRS::Private {};
5935 //! @endcond
5936 
5937 // ---------------------------------------------------------------------------
5938 
5939 //! @cond Doxygen_Suppress
5940 TemporalCRS::~TemporalCRS() = default;
5941 //! @endcond
5942 
5943 // ---------------------------------------------------------------------------
5944 
TemporalCRS(const datum::TemporalDatumNNPtr & datumIn,const cs::TemporalCSNNPtr & csIn)5945 TemporalCRS::TemporalCRS(const datum::TemporalDatumNNPtr &datumIn,
5946                          const cs::TemporalCSNNPtr &csIn)
5947     : SingleCRS(datumIn.as_nullable(), nullptr, csIn), d(nullptr) {}
5948 
5949 // ---------------------------------------------------------------------------
5950 
TemporalCRS(const TemporalCRS & other)5951 TemporalCRS::TemporalCRS(const TemporalCRS &other)
5952     : SingleCRS(other), d(nullptr) {}
5953 
5954 // ---------------------------------------------------------------------------
5955 
_shallowClone() const5956 CRSNNPtr TemporalCRS::_shallowClone() const {
5957     auto crs(TemporalCRS::nn_make_shared<TemporalCRS>(*this));
5958     crs->assignSelf(crs);
5959     return crs;
5960 }
5961 
5962 // ---------------------------------------------------------------------------
5963 
5964 /** \brief Return the datum::TemporalDatum associated with the CRS.
5965  *
5966  * @return a TemporalDatum
5967  */
datum() const5968 const datum::TemporalDatumNNPtr TemporalCRS::datum() const {
5969     return NN_NO_CHECK(std::static_pointer_cast<datum::TemporalDatum>(
5970         SingleCRS::getPrivate()->datum));
5971 }
5972 
5973 // ---------------------------------------------------------------------------
5974 
5975 /** \brief Return the cs::TemporalCS associated with the CRS.
5976  *
5977  * @return a TemporalCS
5978  */
coordinateSystem() const5979 const cs::TemporalCSNNPtr TemporalCRS::coordinateSystem() const {
5980     return util::nn_static_pointer_cast<cs::TemporalCS>(
5981         SingleCRS::getPrivate()->coordinateSystem);
5982 }
5983 
5984 // ---------------------------------------------------------------------------
5985 
5986 /** \brief Instantiate a TemporalCRS from a datum and a coordinate system.
5987  *
5988  * @param properties See \ref general_properties.
5989  * At minimum the name should be defined.
5990  * @param datumIn the datum.
5991  * @param csIn the coordinate system.
5992  * @return new TemporalCRS.
5993  */
create(const util::PropertyMap & properties,const datum::TemporalDatumNNPtr & datumIn,const cs::TemporalCSNNPtr & csIn)5994 TemporalCRSNNPtr TemporalCRS::create(const util::PropertyMap &properties,
5995                                      const datum::TemporalDatumNNPtr &datumIn,
5996                                      const cs::TemporalCSNNPtr &csIn) {
5997     auto crs(TemporalCRS::nn_make_shared<TemporalCRS>(datumIn, csIn));
5998     crs->assignSelf(crs);
5999     crs->setProperties(properties);
6000     return crs;
6001 }
6002 
6003 // ---------------------------------------------------------------------------
6004 
6005 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const6006 void TemporalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
6007     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
6008     if (!isWKT2) {
6009         io::FormattingException::Throw(
6010             "TemporalCRS can only be exported to WKT2");
6011     }
6012     formatter->startNode(io::WKTConstants::TIMECRS, !identifiers().empty());
6013     formatter->addQuotedString(nameStr());
6014     datum()->_exportToWKT(formatter);
6015     coordinateSystem()->_exportToWKT(formatter);
6016     ObjectUsage::baseExportToWKT(formatter);
6017     formatter->endNode();
6018 }
6019 //! @endcond
6020 
6021 // ---------------------------------------------------------------------------
6022 
6023 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const6024 void TemporalCRS::_exportToJSON(
6025     io::JSONFormatter *formatter) const // throw(io::FormattingException)
6026 {
6027     auto writer = formatter->writer();
6028     auto objectContext(
6029         formatter->MakeObjectContext("TemporalCRS", !identifiers().empty()));
6030 
6031     writer->AddObjKey("name");
6032     auto l_name = nameStr();
6033     if (l_name.empty()) {
6034         writer->Add("unnamed");
6035     } else {
6036         writer->Add(l_name);
6037     }
6038 
6039     writer->AddObjKey("datum");
6040     formatter->setOmitTypeInImmediateChild();
6041     datum()->_exportToJSON(formatter);
6042 
6043     writer->AddObjKey("coordinate_system");
6044     formatter->setOmitTypeInImmediateChild();
6045     coordinateSystem()->_exportToJSON(formatter);
6046 
6047     ObjectUsage::baseExportToJSON(formatter);
6048 }
6049 //! @endcond
6050 
6051 // ---------------------------------------------------------------------------
6052 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const6053 bool TemporalCRS::_isEquivalentTo(
6054     const util::IComparable *other, util::IComparable::Criterion criterion,
6055     const io::DatabaseContextPtr &dbContext) const {
6056     auto otherTemporalCRS = dynamic_cast<const TemporalCRS *>(other);
6057     return otherTemporalCRS != nullptr &&
6058            SingleCRS::baseIsEquivalentTo(other, criterion, dbContext);
6059 }
6060 
6061 // ---------------------------------------------------------------------------
6062 
6063 //! @cond Doxygen_Suppress
6064 struct EngineeringCRS::Private {};
6065 //! @endcond
6066 
6067 // ---------------------------------------------------------------------------
6068 
6069 //! @cond Doxygen_Suppress
6070 EngineeringCRS::~EngineeringCRS() = default;
6071 //! @endcond
6072 
6073 // ---------------------------------------------------------------------------
6074 
EngineeringCRS(const datum::EngineeringDatumNNPtr & datumIn,const cs::CoordinateSystemNNPtr & csIn)6075 EngineeringCRS::EngineeringCRS(const datum::EngineeringDatumNNPtr &datumIn,
6076                                const cs::CoordinateSystemNNPtr &csIn)
6077     : SingleCRS(datumIn.as_nullable(), nullptr, csIn),
6078       d(internal::make_unique<Private>()) {}
6079 
6080 // ---------------------------------------------------------------------------
6081 
EngineeringCRS(const EngineeringCRS & other)6082 EngineeringCRS::EngineeringCRS(const EngineeringCRS &other)
6083     : SingleCRS(other), d(internal::make_unique<Private>(*(other.d))) {}
6084 
6085 // ---------------------------------------------------------------------------
6086 
_shallowClone() const6087 CRSNNPtr EngineeringCRS::_shallowClone() const {
6088     auto crs(EngineeringCRS::nn_make_shared<EngineeringCRS>(*this));
6089     crs->assignSelf(crs);
6090     return crs;
6091 }
6092 
6093 // ---------------------------------------------------------------------------
6094 
6095 /** \brief Return the datum::EngineeringDatum associated with the CRS.
6096  *
6097  * @return a EngineeringDatum
6098  */
datum() const6099 const datum::EngineeringDatumNNPtr EngineeringCRS::datum() const {
6100     return NN_NO_CHECK(std::static_pointer_cast<datum::EngineeringDatum>(
6101         SingleCRS::getPrivate()->datum));
6102 }
6103 
6104 // ---------------------------------------------------------------------------
6105 
6106 /** \brief Instantiate a EngineeringCRS from a datum and a coordinate system.
6107  *
6108  * @param properties See \ref general_properties.
6109  * At minimum the name should be defined.
6110  * @param datumIn the datum.
6111  * @param csIn the coordinate system.
6112  * @return new EngineeringCRS.
6113  */
6114 EngineeringCRSNNPtr
create(const util::PropertyMap & properties,const datum::EngineeringDatumNNPtr & datumIn,const cs::CoordinateSystemNNPtr & csIn)6115 EngineeringCRS::create(const util::PropertyMap &properties,
6116                        const datum::EngineeringDatumNNPtr &datumIn,
6117                        const cs::CoordinateSystemNNPtr &csIn) {
6118     auto crs(EngineeringCRS::nn_make_shared<EngineeringCRS>(datumIn, csIn));
6119     crs->assignSelf(crs);
6120     crs->setProperties(properties);
6121 
6122     return crs;
6123 }
6124 
6125 // ---------------------------------------------------------------------------
6126 
6127 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const6128 void EngineeringCRS::_exportToWKT(io::WKTFormatter *formatter) const {
6129     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
6130     formatter->startNode(isWKT2 ? io::WKTConstants::ENGCRS
6131                                 : io::WKTConstants::LOCAL_CS,
6132                          !identifiers().empty());
6133     formatter->addQuotedString(nameStr());
6134     if (isWKT2 || !datum()->nameStr().empty()) {
6135         datum()->_exportToWKT(formatter);
6136     }
6137     if (!isWKT2) {
6138         coordinateSystem()->axisList()[0]->unit()._exportToWKT(formatter);
6139     }
6140 
6141     const auto oldAxisOutputRule = formatter->outputAxis();
6142     formatter->setOutputAxis(io::WKTFormatter::OutputAxisRule::YES);
6143     coordinateSystem()->_exportToWKT(formatter);
6144     formatter->setOutputAxis(oldAxisOutputRule);
6145 
6146     ObjectUsage::baseExportToWKT(formatter);
6147     formatter->endNode();
6148 }
6149 //! @endcond
6150 
6151 // ---------------------------------------------------------------------------
6152 
6153 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const6154 void EngineeringCRS::_exportToJSON(
6155     io::JSONFormatter *formatter) const // throw(io::FormattingException)
6156 {
6157     auto writer = formatter->writer();
6158     auto objectContext(
6159         formatter->MakeObjectContext("EngineeringCRS", !identifiers().empty()));
6160 
6161     writer->AddObjKey("name");
6162     auto l_name = nameStr();
6163     if (l_name.empty()) {
6164         writer->Add("unnamed");
6165     } else {
6166         writer->Add(l_name);
6167     }
6168 
6169     writer->AddObjKey("datum");
6170     formatter->setOmitTypeInImmediateChild();
6171     datum()->_exportToJSON(formatter);
6172 
6173     writer->AddObjKey("coordinate_system");
6174     formatter->setOmitTypeInImmediateChild();
6175     coordinateSystem()->_exportToJSON(formatter);
6176 
6177     ObjectUsage::baseExportToJSON(formatter);
6178 }
6179 //! @endcond
6180 
6181 // ---------------------------------------------------------------------------
6182 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const6183 bool EngineeringCRS::_isEquivalentTo(
6184     const util::IComparable *other, util::IComparable::Criterion criterion,
6185     const io::DatabaseContextPtr &dbContext) const {
6186     auto otherEngineeringCRS = dynamic_cast<const EngineeringCRS *>(other);
6187     return otherEngineeringCRS != nullptr &&
6188            SingleCRS::baseIsEquivalentTo(other, criterion, dbContext);
6189 }
6190 
6191 // ---------------------------------------------------------------------------
6192 
6193 //! @cond Doxygen_Suppress
6194 struct ParametricCRS::Private {};
6195 //! @endcond
6196 
6197 // ---------------------------------------------------------------------------
6198 
6199 //! @cond Doxygen_Suppress
6200 ParametricCRS::~ParametricCRS() = default;
6201 //! @endcond
6202 
6203 // ---------------------------------------------------------------------------
6204 
ParametricCRS(const datum::ParametricDatumNNPtr & datumIn,const cs::ParametricCSNNPtr & csIn)6205 ParametricCRS::ParametricCRS(const datum::ParametricDatumNNPtr &datumIn,
6206                              const cs::ParametricCSNNPtr &csIn)
6207     : SingleCRS(datumIn.as_nullable(), nullptr, csIn), d(nullptr) {}
6208 
6209 // ---------------------------------------------------------------------------
6210 
ParametricCRS(const ParametricCRS & other)6211 ParametricCRS::ParametricCRS(const ParametricCRS &other)
6212     : SingleCRS(other), d(nullptr) {}
6213 
6214 // ---------------------------------------------------------------------------
6215 
_shallowClone() const6216 CRSNNPtr ParametricCRS::_shallowClone() const {
6217     auto crs(ParametricCRS::nn_make_shared<ParametricCRS>(*this));
6218     crs->assignSelf(crs);
6219     return crs;
6220 }
6221 
6222 // ---------------------------------------------------------------------------
6223 
6224 /** \brief Return the datum::ParametricDatum associated with the CRS.
6225  *
6226  * @return a ParametricDatum
6227  */
datum() const6228 const datum::ParametricDatumNNPtr ParametricCRS::datum() const {
6229     return NN_NO_CHECK(std::static_pointer_cast<datum::ParametricDatum>(
6230         SingleCRS::getPrivate()->datum));
6231 }
6232 
6233 // ---------------------------------------------------------------------------
6234 
6235 /** \brief Return the cs::TemporalCS associated with the CRS.
6236  *
6237  * @return a TemporalCS
6238  */
coordinateSystem() const6239 const cs::ParametricCSNNPtr ParametricCRS::coordinateSystem() const {
6240     return util::nn_static_pointer_cast<cs::ParametricCS>(
6241         SingleCRS::getPrivate()->coordinateSystem);
6242 }
6243 
6244 // ---------------------------------------------------------------------------
6245 
6246 /** \brief Instantiate a ParametricCRS from a datum and a coordinate system.
6247  *
6248  * @param properties See \ref general_properties.
6249  * At minimum the name should be defined.
6250  * @param datumIn the datum.
6251  * @param csIn the coordinate system.
6252  * @return new ParametricCRS.
6253  */
6254 ParametricCRSNNPtr
create(const util::PropertyMap & properties,const datum::ParametricDatumNNPtr & datumIn,const cs::ParametricCSNNPtr & csIn)6255 ParametricCRS::create(const util::PropertyMap &properties,
6256                       const datum::ParametricDatumNNPtr &datumIn,
6257                       const cs::ParametricCSNNPtr &csIn) {
6258     auto crs(ParametricCRS::nn_make_shared<ParametricCRS>(datumIn, csIn));
6259     crs->assignSelf(crs);
6260     crs->setProperties(properties);
6261     return crs;
6262 }
6263 
6264 // ---------------------------------------------------------------------------
6265 
6266 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const6267 void ParametricCRS::_exportToWKT(io::WKTFormatter *formatter) const {
6268     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
6269     if (!isWKT2) {
6270         io::FormattingException::Throw(
6271             "ParametricCRS can only be exported to WKT2");
6272     }
6273     formatter->startNode(io::WKTConstants::PARAMETRICCRS,
6274                          !identifiers().empty());
6275     formatter->addQuotedString(nameStr());
6276     datum()->_exportToWKT(formatter);
6277     coordinateSystem()->_exportToWKT(formatter);
6278     ObjectUsage::baseExportToWKT(formatter);
6279     formatter->endNode();
6280 }
6281 //! @endcond
6282 
6283 // ---------------------------------------------------------------------------
6284 
6285 //! @cond Doxygen_Suppress
_exportToJSON(io::JSONFormatter * formatter) const6286 void ParametricCRS::_exportToJSON(
6287     io::JSONFormatter *formatter) const // throw(io::FormattingException)
6288 {
6289     auto writer = formatter->writer();
6290     auto objectContext(
6291         formatter->MakeObjectContext("ParametricCRS", !identifiers().empty()));
6292 
6293     writer->AddObjKey("name");
6294     auto l_name = nameStr();
6295     if (l_name.empty()) {
6296         writer->Add("unnamed");
6297     } else {
6298         writer->Add(l_name);
6299     }
6300 
6301     writer->AddObjKey("datum");
6302     formatter->setOmitTypeInImmediateChild();
6303     datum()->_exportToJSON(formatter);
6304 
6305     writer->AddObjKey("coordinate_system");
6306     formatter->setOmitTypeInImmediateChild();
6307     coordinateSystem()->_exportToJSON(formatter);
6308 
6309     ObjectUsage::baseExportToJSON(formatter);
6310 }
6311 //! @endcond
6312 
6313 // ---------------------------------------------------------------------------
6314 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const6315 bool ParametricCRS::_isEquivalentTo(
6316     const util::IComparable *other, util::IComparable::Criterion criterion,
6317     const io::DatabaseContextPtr &dbContext) const {
6318     auto otherParametricCRS = dynamic_cast<const ParametricCRS *>(other);
6319     return otherParametricCRS != nullptr &&
6320            SingleCRS::baseIsEquivalentTo(other, criterion, dbContext);
6321 }
6322 
6323 // ---------------------------------------------------------------------------
6324 
6325 //! @cond Doxygen_Suppress
6326 struct DerivedVerticalCRS::Private {};
6327 //! @endcond
6328 
6329 // ---------------------------------------------------------------------------
6330 
6331 //! @cond Doxygen_Suppress
6332 DerivedVerticalCRS::~DerivedVerticalCRS() = default;
6333 //! @endcond
6334 
6335 // ---------------------------------------------------------------------------
6336 
DerivedVerticalCRS(const VerticalCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::VerticalCSNNPtr & csIn)6337 DerivedVerticalCRS::DerivedVerticalCRS(
6338     const VerticalCRSNNPtr &baseCRSIn,
6339     const operation::ConversionNNPtr &derivingConversionIn,
6340     const cs::VerticalCSNNPtr &csIn)
6341     : SingleCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
6342       VerticalCRS(baseCRSIn->datum(), baseCRSIn->datumEnsemble(), csIn),
6343       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
6344 
6345 // ---------------------------------------------------------------------------
6346 
DerivedVerticalCRS(const DerivedVerticalCRS & other)6347 DerivedVerticalCRS::DerivedVerticalCRS(const DerivedVerticalCRS &other)
6348     : SingleCRS(other), VerticalCRS(other), DerivedCRS(other), d(nullptr) {}
6349 
6350 // ---------------------------------------------------------------------------
6351 
_shallowClone() const6352 CRSNNPtr DerivedVerticalCRS::_shallowClone() const {
6353     auto crs(DerivedVerticalCRS::nn_make_shared<DerivedVerticalCRS>(*this));
6354     crs->assignSelf(crs);
6355     crs->setDerivingConversionCRS();
6356     return crs;
6357 }
6358 
6359 // ---------------------------------------------------------------------------
6360 
6361 /** \brief Return the base CRS (a VerticalCRS) of a DerivedVerticalCRS.
6362  *
6363  * @return the base CRS.
6364  */
baseCRS() const6365 const VerticalCRSNNPtr DerivedVerticalCRS::baseCRS() const {
6366     return NN_NO_CHECK(util::nn_dynamic_pointer_cast<VerticalCRS>(
6367         DerivedCRS::getPrivate()->baseCRS_));
6368 }
6369 
6370 // ---------------------------------------------------------------------------
6371 
6372 /** \brief Instantiate a DerivedVerticalCRS from a base CRS, a deriving
6373  * conversion and a cs::VerticalCS.
6374  *
6375  * @param properties See \ref general_properties.
6376  * At minimum the name should be defined.
6377  * @param baseCRSIn base CRS.
6378  * @param derivingConversionIn the deriving conversion from the base CRS to this
6379  * CRS.
6380  * @param csIn the coordinate system.
6381  * @return new DerivedVerticalCRS.
6382  */
create(const util::PropertyMap & properties,const VerticalCRSNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const cs::VerticalCSNNPtr & csIn)6383 DerivedVerticalCRSNNPtr DerivedVerticalCRS::create(
6384     const util::PropertyMap &properties, const VerticalCRSNNPtr &baseCRSIn,
6385     const operation::ConversionNNPtr &derivingConversionIn,
6386     const cs::VerticalCSNNPtr &csIn) {
6387     auto crs(DerivedVerticalCRS::nn_make_shared<DerivedVerticalCRS>(
6388         baseCRSIn, derivingConversionIn, csIn));
6389     crs->assignSelf(crs);
6390     crs->setProperties(properties);
6391     crs->setDerivingConversionCRS();
6392     return crs;
6393 }
6394 
6395 // ---------------------------------------------------------------------------
6396 
6397 //! @cond Doxygen_Suppress
_exportToWKT(io::WKTFormatter * formatter) const6398 void DerivedVerticalCRS::_exportToWKT(io::WKTFormatter *formatter) const {
6399     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
6400     if (!isWKT2) {
6401         io::FormattingException::Throw(
6402             "DerivedVerticalCRS can only be exported to WKT2");
6403     }
6404     baseExportToWKT(formatter, io::WKTConstants::VERTCRS,
6405                     io::WKTConstants::BASEVERTCRS);
6406 }
6407 //! @endcond
6408 
6409 // ---------------------------------------------------------------------------
6410 
_exportToPROJString(io::PROJStringFormatter *) const6411 void DerivedVerticalCRS::_exportToPROJString(
6412     io::PROJStringFormatter *) const // throw(io::FormattingException)
6413 {
6414     throw io::FormattingException(
6415         "DerivedVerticalCRS cannot be exported to PROJ string");
6416 }
6417 
6418 // ---------------------------------------------------------------------------
6419 
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const6420 bool DerivedVerticalCRS::_isEquivalentTo(
6421     const util::IComparable *other, util::IComparable::Criterion criterion,
6422     const io::DatabaseContextPtr &dbContext) const {
6423     auto otherDerivedCRS = dynamic_cast<const DerivedVerticalCRS *>(other);
6424     return otherDerivedCRS != nullptr &&
6425            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
6426 }
6427 
6428 // ---------------------------------------------------------------------------
6429 
6430 //! @cond Doxygen_Suppress
6431 
6432 std::list<std::pair<CRSNNPtr, int>>
_identify(const io::AuthorityFactoryPtr & factory) const6433 DerivedVerticalCRS::_identify(const io::AuthorityFactoryPtr &factory) const {
6434     return CRS::_identify(factory);
6435 }
6436 
6437 //! @endcond
6438 
6439 // ---------------------------------------------------------------------------
6440 
6441 //! @cond Doxygen_Suppress
6442 template <class DerivedCRSTraits>
6443 struct DerivedCRSTemplate<DerivedCRSTraits>::Private {};
6444 //! @endcond
6445 
6446 // ---------------------------------------------------------------------------
6447 
6448 //! @cond Doxygen_Suppress
6449 template <class DerivedCRSTraits>
6450 DerivedCRSTemplate<DerivedCRSTraits>::~DerivedCRSTemplate() = default;
6451 //! @endcond
6452 
6453 // ---------------------------------------------------------------------------
6454 
6455 template <class DerivedCRSTraits>
DerivedCRSTemplate(const BaseNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const CSNNPtr & csIn)6456 DerivedCRSTemplate<DerivedCRSTraits>::DerivedCRSTemplate(
6457     const BaseNNPtr &baseCRSIn,
6458     const operation::ConversionNNPtr &derivingConversionIn, const CSNNPtr &csIn)
6459     : SingleCRS(baseCRSIn->datum().as_nullable(), nullptr, csIn),
6460       BaseType(baseCRSIn->datum(), csIn),
6461       DerivedCRS(baseCRSIn, derivingConversionIn, csIn), d(nullptr) {}
6462 
6463 // ---------------------------------------------------------------------------
6464 
6465 template <class DerivedCRSTraits>
DerivedCRSTemplate(const DerivedCRSTemplate & other)6466 DerivedCRSTemplate<DerivedCRSTraits>::DerivedCRSTemplate(
6467     const DerivedCRSTemplate &other)
6468     : SingleCRS(other), BaseType(other), DerivedCRS(other), d(nullptr) {}
6469 
6470 // ---------------------------------------------------------------------------
6471 
6472 template <class DerivedCRSTraits>
6473 const typename DerivedCRSTemplate<DerivedCRSTraits>::BaseNNPtr
baseCRS() const6474 DerivedCRSTemplate<DerivedCRSTraits>::baseCRS() const {
6475     auto l_baseCRS = DerivedCRS::getPrivate()->baseCRS_;
6476     return NN_NO_CHECK(util::nn_dynamic_pointer_cast<BaseType>(l_baseCRS));
6477 }
6478 
6479 // ---------------------------------------------------------------------------
6480 
6481 //! @cond Doxygen_Suppress
6482 
6483 template <class DerivedCRSTraits>
_shallowClone() const6484 CRSNNPtr DerivedCRSTemplate<DerivedCRSTraits>::_shallowClone() const {
6485     auto crs(DerivedCRSTemplate::nn_make_shared<DerivedCRSTemplate>(*this));
6486     crs->assignSelf(crs);
6487     crs->setDerivingConversionCRS();
6488     return crs;
6489 }
6490 
6491 // ---------------------------------------------------------------------------
6492 
6493 template <class DerivedCRSTraits>
6494 typename DerivedCRSTemplate<DerivedCRSTraits>::NNPtr
create(const util::PropertyMap & properties,const BaseNNPtr & baseCRSIn,const operation::ConversionNNPtr & derivingConversionIn,const CSNNPtr & csIn)6495 DerivedCRSTemplate<DerivedCRSTraits>::create(
6496     const util::PropertyMap &properties, const BaseNNPtr &baseCRSIn,
6497     const operation::ConversionNNPtr &derivingConversionIn,
6498     const CSNNPtr &csIn) {
6499     auto crs(DerivedCRSTemplate::nn_make_shared<DerivedCRSTemplate>(
6500         baseCRSIn, derivingConversionIn, csIn));
6501     crs->assignSelf(crs);
6502     crs->setProperties(properties);
6503     crs->setDerivingConversionCRS();
6504     return crs;
6505 }
6506 
6507 // ---------------------------------------------------------------------------
6508 
6509 template <class DerivedCRSTraits>
className() const6510 const char *DerivedCRSTemplate<DerivedCRSTraits>::className() const {
6511     return DerivedCRSTraits::CRSName().c_str();
6512 }
6513 
6514 // ---------------------------------------------------------------------------
6515 
DerivedCRSTemplateCheckExportToWKT(io::WKTFormatter * formatter,const std::string & crsName,bool wkt2_2019_only)6516 static void DerivedCRSTemplateCheckExportToWKT(io::WKTFormatter *formatter,
6517                                                const std::string &crsName,
6518                                                bool wkt2_2019_only) {
6519     const bool isWKT2 = formatter->version() == io::WKTFormatter::Version::WKT2;
6520     if (!isWKT2 || (wkt2_2019_only && !formatter->use2019Keywords())) {
6521         io::FormattingException::Throw(crsName +
6522                                        " can only be exported to WKT2" +
6523                                        (wkt2_2019_only ? ":2019" : ""));
6524     }
6525 }
6526 
6527 // ---------------------------------------------------------------------------
6528 
6529 template <class DerivedCRSTraits>
_exportToWKT(io::WKTFormatter * formatter) const6530 void DerivedCRSTemplate<DerivedCRSTraits>::_exportToWKT(
6531     io::WKTFormatter *formatter) const {
6532     DerivedCRSTemplateCheckExportToWKT(formatter, DerivedCRSTraits::CRSName(),
6533                                        DerivedCRSTraits::wkt2_2019_only);
6534     baseExportToWKT(formatter, DerivedCRSTraits::WKTKeyword(),
6535                     DerivedCRSTraits::WKTBaseKeyword());
6536 }
6537 
6538 // ---------------------------------------------------------------------------
6539 
6540 template <class DerivedCRSTraits>
_isEquivalentTo(const util::IComparable * other,util::IComparable::Criterion criterion,const io::DatabaseContextPtr & dbContext) const6541 bool DerivedCRSTemplate<DerivedCRSTraits>::_isEquivalentTo(
6542     const util::IComparable *other, util::IComparable::Criterion criterion,
6543     const io::DatabaseContextPtr &dbContext) const {
6544     auto otherDerivedCRS = dynamic_cast<const DerivedCRSTemplate *>(other);
6545     return otherDerivedCRS != nullptr &&
6546            DerivedCRS::_isEquivalentTo(other, criterion, dbContext);
6547 }
6548 
6549 //! @endcond
6550 
6551 // ---------------------------------------------------------------------------
6552 
6553 //! @cond Doxygen_Suppress
6554 static const std::string STRING_DerivedEngineeringCRS("DerivedEngineeringCRS");
CRSName()6555 const std::string &DerivedEngineeringCRSTraits::CRSName() {
6556     return STRING_DerivedEngineeringCRS;
6557 }
WKTKeyword()6558 const std::string &DerivedEngineeringCRSTraits::WKTKeyword() {
6559     return io::WKTConstants::ENGCRS;
6560 }
WKTBaseKeyword()6561 const std::string &DerivedEngineeringCRSTraits::WKTBaseKeyword() {
6562     return io::WKTConstants::BASEENGCRS;
6563 }
6564 
6565 template class DerivedCRSTemplate<DerivedEngineeringCRSTraits>;
6566 //! @endcond
6567 
6568 // ---------------------------------------------------------------------------
6569 
6570 //! @cond Doxygen_Suppress
6571 static const std::string STRING_DerivedParametricCRS("DerivedParametricCRS");
CRSName()6572 const std::string &DerivedParametricCRSTraits::CRSName() {
6573     return STRING_DerivedParametricCRS;
6574 }
WKTKeyword()6575 const std::string &DerivedParametricCRSTraits::WKTKeyword() {
6576     return io::WKTConstants::PARAMETRICCRS;
6577 }
WKTBaseKeyword()6578 const std::string &DerivedParametricCRSTraits::WKTBaseKeyword() {
6579     return io::WKTConstants::BASEPARAMCRS;
6580 }
6581 
6582 template class DerivedCRSTemplate<DerivedParametricCRSTraits>;
6583 //! @endcond
6584 
6585 // ---------------------------------------------------------------------------
6586 
6587 //! @cond Doxygen_Suppress
6588 static const std::string STRING_DerivedTemporalCRS("DerivedTemporalCRS");
CRSName()6589 const std::string &DerivedTemporalCRSTraits::CRSName() {
6590     return STRING_DerivedTemporalCRS;
6591 }
WKTKeyword()6592 const std::string &DerivedTemporalCRSTraits::WKTKeyword() {
6593     return io::WKTConstants::TIMECRS;
6594 }
WKTBaseKeyword()6595 const std::string &DerivedTemporalCRSTraits::WKTBaseKeyword() {
6596     return io::WKTConstants::BASETIMECRS;
6597 }
6598 
6599 template class DerivedCRSTemplate<DerivedTemporalCRSTraits>;
6600 //! @endcond
6601 
6602 // ---------------------------------------------------------------------------
6603 
6604 } // namespace crs
6605 NS_PROJ_END
6606