path: root/src/transformations/defmodel.hpp

/******************************************************************************
 * Project:  PROJ
 * Purpose:  Functionality related to deformation model
 * Author:   Even Rouault, <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com>
 *
 * 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 <algorithm>
#include <cmath>
#include <exception>
#include <limits>
#include <memory>
#include <string>
#include <vector>

#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
     */
    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<double>::quiet_NaN();
    double mMiny = std::numeric_limits<double>::quiet_NaN();
    double mMaxx = std::numeric_limits<double>::quiet_NaN();
    double mMaxy = std::numeric_limits<double>::quiet_NaN();
    double mMinxRad = std::numeric_limits<double>::quiet_NaN();
    double mMinyRad = std::numeric_limits<double>::quiet_NaN();
    double mMaxxRad = std::numeric_limits<double>::quiet_NaN();
    double mMaxyRad = std::numeric_limits<double>::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
     */
    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<double>::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<EpochScaleFactorTuple> 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<double>::quiet_NaN();

        /** The scale factor that applies before the reference epoch. */
        double beforeScaleFactor = std::numeric_limits<double>::quiet_NaN();

        /** Initial scale factor. */
        double initialScaleFactor = std::numeric_limits<double>::quiet_NaN();

        /** The scale factor the exponential function approaches. */
        double finalScaleFactor = std::numeric_limits<double>::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<double>::quiet_NaN();
    double mVerticalUncertainty = std::numeric_limits<double>::quiet_NaN();
    SpatialModel mSpatialModel{};
    std::unique_ptr<TimeFunction> 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
     */
    static std::unique_ptr<MasterFile> 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<Link> 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<Component> &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<Link> 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<Component> 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 getLonLatOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/,
                         double & /*latOffsetRadian*/) const {
        throw UnimplementedException("getLonLatOffset 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 getLonLatZOffset(int /*ix*/, int /*iy*/, double & /*lonOffsetRadian*/,
                          double & /*latOffsetRadian*/,
                          double & /*zOffset*/) const {
        throw UnimplementedException("getLonLatZOffset unimplemented");
#if 0
        return getLonLatOffset(ix, iy, lonOffsetRadian, 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 <class Grid = GridPrototype> 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 <class Grid = GridPrototype, class GridSet = GridSetPrototype<>>
struct EvaluatorIfacePrototype {

    std::unique_ptr<GridSet> 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 <class Grid, class GridSet> struct ComponentEx;

// ---------------------------------------------------------------------------

/** Class to evaluate the transformation of a coordinate */
template <class Grid = GridPrototype, class GridSet = GridSetPrototype<>,
          class EvaluatorIface = EvaluatorIfacePrototype<>>
class Evaluator {
  public:
    /** Constructor. May throw EvaluatorException */
    explicit Evaluator(std::unique_ptr<MasterFile> &&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<MasterFile> 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<std::unique_ptr<ComponentEx<Grid, GridSet>>> mComponents{};
};

// ---------------------------------------------------------------------------

} // namespace DEFORMATON_MODEL_NAMESPACE

// ---------------------------------------------------------------------------

#include "defmodel_impl.hpp"

#endif // DEFMODEL_HPP