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