/****************************************************************************** * Project: PROJ * Purpose: Functionality related to deformation model * Author: Even Rouault, * ****************************************************************************** * Copyright (c) 2020, Even Rouault, * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. *****************************************************************************/ /** This file implements the gridded deformation model proposol of * https://docs.google.com/document/d/1wiyrAmzqh8MZlzHSp3wf594Ob_M1LeFtDA5swuzvLZY * It is written in a generic way, independent of the rest of PROJ * infrastructure. * * Verbose debugging info can be turned on by setting the DEBUG_DEFMODEL macro */ #ifndef DEFMODEL_HPP #define DEFMODEL_HPP #ifdef PROJ_COMPILATION #include "proj/internal/include_nlohmann_json.hpp" #else #include "nlohmann/json.hpp" #endif #include #include #include #include #include #include #include #ifndef DEFORMATON_MODEL_NAMESPACE #define DEFORMATON_MODEL_NAMESPACE DeformationModel #endif #include "defmodel_exceptions.hpp" namespace DEFORMATON_MODEL_NAMESPACE { using json = nlohmann::json; // --------------------------------------------------------------------------- /** Spatial extent as a bounding box. */ class SpatialExtent { public: /** Parse the provided object as an extent. * * @throws ParsingException in case of error. */ static SpatialExtent parse(const json &j); double minx() const { return mMinx; } double miny() const { return mMiny; } double maxx() const { return mMaxx; } double maxy() const { return mMaxy; } double minxNormalized(bool bIsGeographic) const { return bIsGeographic ? mMinxRad : mMinx; } double minyNormalized(bool bIsGeographic) const { return bIsGeographic ? mMinyRad : mMiny; } double maxxNormalized(bool bIsGeographic) const { return bIsGeographic ? mMaxxRad : mMaxx; } double maxyNormalized(bool bIsGeographic) const { return bIsGeographic ? mMaxyRad : mMaxy; } protected: friend class MasterFile; friend class Component; SpatialExtent() = default; private: double mMinx = std::numeric_limits::quiet_NaN(); double mMiny = std::numeric_limits::quiet_NaN(); double mMaxx = std::numeric_limits::quiet_NaN(); double mMaxy = std::numeric_limits::quiet_NaN(); double mMinxRad = std::numeric_limits::quiet_NaN(); double mMinyRad = std::numeric_limits::quiet_NaN(); double mMaxxRad = std::numeric_limits::quiet_NaN(); double mMaxyRad = std::numeric_limits::quiet_NaN(); }; // --------------------------------------------------------------------------- /** Epoch */ class Epoch { public: /** Constructor from a ISO 8601 date-time. May throw ParsingException */ explicit Epoch(const std::string &dt = std::string()); /** Return ISO 8601 date-time */ const std::string &toString() const { return mDt; } /** Return decimal year */ double toDecimalYear() const; private: std::string mDt{}; double mDecimalYear = 0; }; // --------------------------------------------------------------------------- /** Component of a deformation model. */ class Component { public: /** Parse the provided object as a component. * * @throws ParsingException in case of error. */ static Component parse(const json &j); /** Get a text description of the component. */ const std::string &description() const { return mDescription; } /** Get the region within which the component is defined. Outside this * region the component evaluates to 0. */ const SpatialExtent &extent() const { return mSpatialExtent; } /** Get the displacement parameters defined by the component, one of * "none", "horizontal", "vertical", and "3d". The "none" option allows * for a component which defines uncertainty with different grids to those * defining displacement. */ const std::string &displacementType() const { return mDisplacementType; } /** Get the uncertainty parameters defined by the component, * one of "none", "horizontal", "vertical", "3d". */ const std::string &uncertaintyType() const { return mUncertaintyType; } /** Get the horizontal uncertainty to use if it is not defined explicitly * in the spatial model. */ double horizontalUncertainty() const { return mHorizontalUncertainty; } /** Get the vertical uncertainty to use if it is not defined explicitly in * the spatial model. */ double verticalUncertainty() const { return mVerticalUncertainty; } struct SpatialModel { /** Specifies the type of the spatial model data file. Initially it * is proposed that only "GeoTIFF" is supported. */ std::string type{}; /** How values in model should be interpolated. This proposal will * support "bilinear" and "geocentric_bilinear". */ std::string interpolationMethod{}; /** Specifies location of the spatial model GeoTIFF file relative to * the master JSON file. */ std::string filename{}; /** A hex encoded MD5 checksum of the grid file that can be used to * validate that it is the correct version of the file. */ std::string md5Checksum{}; }; /** Get the spatial model. */ const SpatialModel &spatialModel() const { return mSpatialModel; } /** Generic type for a type function */ struct TimeFunction { std::string type{}; virtual ~TimeFunction(); virtual double evaluateAt(double dt) const = 0; protected: TimeFunction() = default; }; struct ConstantTimeFunction : public TimeFunction { virtual double evaluateAt(double dt) const override; }; struct VelocityTimeFunction : public TimeFunction { /** Date/time at which the velocity function is zero. */ Epoch referenceEpoch{}; virtual double evaluateAt(double dt) const override; }; struct StepTimeFunction : public TimeFunction { /** Epoch at which the step function transitions from 0 to 1. */ Epoch stepEpoch{}; virtual double evaluateAt(double dt) const override; }; struct ReverseStepTimeFunction : public TimeFunction { /** Epoch at which the reverse step function transitions from 1. to 0 */ Epoch stepEpoch{}; virtual double evaluateAt(double dt) const override; }; struct PiecewiseTimeFunction : public TimeFunction { /** One of "zero", "constant", and "linear", defines the behavior of * the function before the first defined epoch */ std::string beforeFirst{}; /** One of "zero", "constant", and "linear", defines the behavior of * the function after the last defined epoch */ std::string afterLast{}; struct EpochScaleFactorTuple { /** Defines the date/time of the data point. */ Epoch epoch{}; /** Function value at the above epoch */ double scaleFactor = std::numeric_limits::quiet_NaN(); }; /** A sorted array data points each defined by two elements. * The array is sorted in order of increasing epoch. * Note: where the time function includes a step it is represented by * two consecutive data points with the same epoch. The first defines * the scale factor that applies before the epoch and the second the * scale factor that applies after the epoch. */ std::vector model{}; virtual double evaluateAt(double dt) const override; }; struct ExponentialTimeFunction : public TimeFunction { /** The date/time at which the exponential decay starts. */ Epoch referenceEpoch{}; /** The date/time at which the exponential decay ends. */ Epoch endEpoch{}; /** The relaxation constant in years. */ double relaxationConstant = std::numeric_limits::quiet_NaN(); /** The scale factor that applies before the reference epoch. */ double beforeScaleFactor = std::numeric_limits::quiet_NaN(); /** Initial scale factor. */ double initialScaleFactor = std::numeric_limits::quiet_NaN(); /** The scale factor the exponential function approaches. */ double finalScaleFactor = std::numeric_limits::quiet_NaN(); virtual double evaluateAt(double dt) const override; }; /** Get the time function. */ const TimeFunction *timeFunction() const { return mTimeFunction.get(); } private: Component() = default; std::string mDescription{}; SpatialExtent mSpatialExtent{}; std::string mDisplacementType{}; std::string mUncertaintyType{}; double mHorizontalUncertainty = std::numeric_limits::quiet_NaN(); double mVerticalUncertainty = std::numeric_limits::quiet_NaN(); SpatialModel mSpatialModel{}; std::unique_ptr mTimeFunction{}; }; Component::TimeFunction::~TimeFunction() = default; // --------------------------------------------------------------------------- /** Master file of a deformation model. */ class MasterFile { public: /** Parse the provided serialized JSON content and return an object. * * @throws ParsingException in case of error. */ static std::unique_ptr parse(const std::string &text); /** Get file type. Should always be "deformation_model_master_file" */ const std::string &fileType() const { return mFileType; } /** Get the version of the format. At time of writing, only "1.0" is known */ const std::string &formatVersion() const { return mFormatVersion; } /** Get brief descriptive name of the deformation model. */ const std::string &name() const { return mName; } /** Get a string identifying the version of the deformation model. * The format for specifying version is defined by the agency * responsible for the deformation model. */ const std::string &version() const { return mVersion; } /** Get a string identifying the license of the file. * e.g "Create Commons Attribution 4.0 International" */ const std::string &license() const { return mLicense; } /** Get a text description of the model. Intended to be longer than name() */ const std::string &description() const { return mDescription; } /** Get a text description of the model. Intended to be longer than name() */ const std::string &publicationDate() const { return mPublicationDate; } /** Basic information on the agency responsible for the model. */ struct Authority { std::string name{}; std::string url{}; std::string address{}; std::string email{}; }; /** Get basic information on the agency responsible for the model. */ const Authority &authority() const { return mAuthority; } /** Hyperlink related to the model. */ struct Link { /** URL holding the information */ std::string href{}; /** Relationship to the dataset. e.g. "about", "source", "license", * "metadata" */ std::string rel{}; /** Mime type */ std::string type{}; /** Description of the link */ std::string title{}; }; /** Get links to related information. */ const std::vector links() const { return mLinks; } /** Get a string identifying the source CRS. That is the coordinate * reference system to which the deformation model applies. Typically * "EPSG:XXXX" */ const std::string &sourceCRS() const { return mSourceCRS; } /** Get a string identifying the target CRS. That is, for a time * dependent coordinate transformation, the coordinate reference * system resulting from applying the deformation. * Typically "EPSG:XXXX" */ const std::string &targetCRS() const { return mTargetCRS; } /** Get a string identifying the definition CRS. That is, the * coordinate reference system used to define the component spatial * models. Typically "EPSG:XXXX" */ const std::string &definitionCRS() const { return mDefinitionCRS; } /** Get the nominal reference epoch of the deformation model. Formatted * as a ISO-8601 date-time. This is not necessarily used to calculate * the deformation model - each component defines its own time function. */ const std::string &referenceEpoch() const { return mReferenceEpoch; } /** Get the epoch at which the uncertainties of the deformation model * are calculated. Formatted as a ISO-8601 date-time. */ const std::string &uncertaintyReferenceEpoch() const { return mUncertaintyReferenceEpoch; } /** Unit of horizontal offsets. Only "metre" and "degree" are supported. */ const std::string &horizontalOffsetUnit() const { return mHorizontalOffsetUnit; } /** Unit of vertical offsets. Only "metre" is supported. */ const std::string &verticalOffsetUnit() const { return mVerticalOffsetUnit; } /** Type of horizontal uncertainty. e.g "circular 95% confidence limit" */ const std::string &horizontalUncertaintyType() const { return mHorizontalUncertaintyType; } /** Unit of horizontal uncertainty. Only "metre" is supported. */ const std::string &horizontalUncertaintyUnit() const { return mHorizontalUncertaintyUnit; } /** Type of vertical uncertainty. e.g "circular 95% confidence limit" */ const std::string &verticalUncertaintyType() const { return mVerticalUncertaintyType; } /** Unit of vertical uncertainty. Only "metre" is supported. */ const std::string &verticalUncertaintyUnit() const { return mVerticalUncertaintyUnit; } /** Defines how the horizontal offsets are applied to geographic * coordinates. Only "addition" and "geocentric" are supported */ const std::string &horizontalOffsetMethod() const { return mHorizontalOffsetMethod; } /** Get the region within which the deformation model is defined. * It cannot be calculated outside this region */ const SpatialExtent &extent() const { return mSpatialExtent; } /** Defines the range of times for which the model is valid, specified * by a first and a last value. The deformation model is undefined for * dates outside this range. */ struct TimeExtent { Epoch first{}; Epoch last{}; }; /** Get the range of times for which the model is valid. */ const TimeExtent &timeExtent() const { return mTimeExtent; } /** Get an array of the components comprising the deformation model. */ const std::vector &components() const { return mComponents; } private: MasterFile() = default; std::string mFileType{}; std::string mFormatVersion{}; std::string mName{}; std::string mVersion{}; std::string mLicense{}; std::string mDescription{}; std::string mPublicationDate{}; Authority mAuthority{}; std::vector mLinks{}; std::string mSourceCRS{}; std::string mTargetCRS{}; std::string mDefinitionCRS{}; std::string mReferenceEpoch{}; std::string mUncertaintyReferenceEpoch{}; std::string mHorizontalOffsetUnit{}; std::string mVerticalOffsetUnit{}; std::string mHorizontalUncertaintyType{}; std::string mHorizontalUncertaintyUnit{}; std::string mVerticalUncertaintyType{}; std::string mVerticalUncertaintyUnit{}; std::string mHorizontalOffsetMethod{}; SpatialExtent mSpatialExtent{}; TimeExtent mTimeExtent{}; std::vector mComponents{}; }; // --------------------------------------------------------------------------- /** Prototype for a Grid used by GridSet. Intended to be implemented * by user code */ struct GridPrototype { double minx = 0; double miny = 0; double resx = 0; double resy = 0; int width = 0; int height = 0; // cppcheck-suppress functionStatic bool getLongLatOffset(int /*ix*/, int /*iy*/, double & /*longOffsetRadian*/, double & /*latOffsetRadian*/) const { throw UnimplementedException("getLongLatOffset unimplemented"); } // cppcheck-suppress functionStatic bool getZOffset(int /*ix*/, int /*iy*/, double & /*zOffset*/) const { throw UnimplementedException("getZOffset unimplemented"); } // cppcheck-suppress functionStatic bool getEastingNorthingOffset(int /*ix*/, int /*iy*/, double & /*eastingOffset*/, double & /*northingOffset*/) const { throw UnimplementedException("getEastingNorthingOffset unimplemented"); } // cppcheck-suppress functionStatic bool getLongLatZOffset(int /*ix*/, int /*iy*/, double & /*longOffsetRadian*/, double & /*latOffsetRadian*/, double & /*zOffset*/) const { throw UnimplementedException("getLongLatZOffset unimplemented"); #if 0 return getLongLatOffset(ix, iy, longOffsetRadian, latOffsetRadian) && getZOffset(ix, iy, zOffset); #endif } // cppcheck-suppress functionStatic bool getEastingNorthingZOffset(int /*ix*/, int /*iy*/, double & /*eastingOffset*/, double & /*northingOffset*/, double & /*zOffset*/) const { throw UnimplementedException("getEastingNorthingOffset unimplemented"); #if 0 return getEastingNorthingOffset(ix, iy, eastingOffset, northingOffset) && getZOffset(ix, iy, zOffset); #endif } #ifdef DEBUG_DEFMODEL std::string name() const { throw UnimplementedException("name() unimplemented"); } #endif }; // --------------------------------------------------------------------------- /** Prototype for a GridSet used by EvaluatorIface. Intended to be implemented * by user code */ template struct GridSetPrototype { // The return pointer should remain "stable" over time for a given grid // of a GridSet. // cppcheck-suppress functionStatic const Grid *gridAt(double /*x */, double /* y */) { throw UnimplementedException("gridAt unimplemented"); } }; // --------------------------------------------------------------------------- /** Prototype for a EvaluatorIface used by Evaluator. Intended to be implemented * by user code */ template > struct EvaluatorIfacePrototype { std::unique_ptr open(const std::string & /* filename*/) { throw UnimplementedException("open unimplemented"); } // cppcheck-suppress functionStatic void geographicToGeocentric(double /* lam */, double /* phi */, double /* height*/, double /* a */, double /* b */, double /*es*/, double & /* X */, double & /* Y */, double & /* Z */) { throw UnimplementedException("geographicToGeocentric unimplemented"); } // cppcheck-suppress functionStatic void geocentricToGeographic(double /* X */, double /* Y */, double /* Z */, double /* a */, double /* b */, double /*es*/, double & /* lam */, double & /* phi */, double & /* height*/) { throw UnimplementedException("geocentricToGeographic unimplemented"); } // cppcheck-suppress functionStatic bool isGeographicCRS(const std::string & /* crsDef */) { throw UnimplementedException("isGeographicCRS unimplemented"); } #ifdef DEBUG_DEFMODEL void log(const std::string & /* msg */) { throw UnimplementedException("log unimplemented"); } #endif }; // --------------------------------------------------------------------------- /** Internal class to offer caching services over a Component */ template struct ComponentEx; // --------------------------------------------------------------------------- /** Class to evaluate the transformation of a coordinate */ template , class EvaluatorIface = EvaluatorIfacePrototype<>> class Evaluator { public: /** Constructor. May throw EvaluatorException */ explicit Evaluator(std::unique_ptr &&model, EvaluatorIface &iface, double a, double b); /** Evaluate displacement of a position given by (x,y,z,t) and * return it in (x_out,y_out_,z_out). * For geographic CRS (only supported at that time), x must be a * longitude, and y a latitude. */ bool forward(EvaluatorIface &iface, double x, double y, double z, double t, double &x_out, double &y_out, double &z_out) { return forward(iface, x, y, z, t, false, x_out, y_out, z_out); } /** Apply inverse transformation. */ bool inverse(EvaluatorIface &iface, double x, double y, double z, double t, double &x_out, double &y_out, double &z_out); /** Clear grid cache */ void clearGridCache(); /** Return whether the definition CRS is a geographic CRS */ bool isGeographicCRS() const { return mIsGeographicCRS; } private: std::unique_ptr mModel; const double mA; const double mB; const double mEs; const bool mIsHorizontalUnitDegree; /* degree vs metre */ const bool mIsAddition; /* addition vs geocentric */ const bool mIsGeographicCRS; bool forward(EvaluatorIface &iface, double x, double y, double z, double t, bool forInverseComputation, double &x_out, double &y_out, double &z_out); std::vector>> mComponents{}; }; // --------------------------------------------------------------------------- } // namespace DEFORMATON_MODEL_NAMESPACE // --------------------------------------------------------------------------- #include "defmodel_impl.hpp" #endif // DEFMODEL_HPP