diff options
Diffstat (limited to 'src/transformations/defmodel_impl.hpp')
| -rw-r--r-- | src/transformations/defmodel_impl.hpp | 1265 |
1 files changed, 1265 insertions, 0 deletions
diff --git a/src/transformations/defmodel_impl.hpp b/src/transformations/defmodel_impl.hpp new file mode 100644 index 00000000..a15137d7 --- /dev/null +++ b/src/transformations/defmodel_impl.hpp @@ -0,0 +1,1265 @@ +/****************************************************************************** + * 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. + *****************************************************************************/ + +#ifndef DEFORMATON_MODEL_NAMESPACE +#error "Should be included only by defmodel.hpp" +#endif + +namespace DEFORMATON_MODEL_NAMESPACE { + +// --------------------------------------------------------------------------- + +static const std::string STR_DEGREE("degree"); +static const std::string STR_METRE("metre"); + +static const std::string STR_ADDITION("addition"); +static const std::string STR_GEOCENTRIC("geocentric"); + +static const std::string STR_BILINEAR("bilinear"); +static const std::string STR_GEOCENTRIC_BILINEAR("geocentric_bilinear"); + +static const std::string STR_NONE("none"); +static const std::string STR_HORIZONTAL("horizontal"); +static const std::string STR_VERTICAL("vertical"); +static const std::string STR_3D("3d"); + +constexpr double DEFMODEL_PI = 3.14159265358979323846; +constexpr double DEG_TO_RAD_CONSTANT = 3.14159265358979323846 / 180.; +inline constexpr double DegToRad(double d) { return d * DEG_TO_RAD_CONSTANT; } + +// --------------------------------------------------------------------------- + +enum class DisplacementType { NONE, HORIZONTAL, VERTICAL, THREE_D }; + +// --------------------------------------------------------------------------- + +/** Internal class to offer caching services over a Grid */ +template <class Grid> struct GridEx { + const Grid *grid; + + bool smallResx; // lesser than one degree + + double sinhalfresx; + double coshalfresx; + + double sinresy; + double cosresy; + + int last_ix0 = -1; + int last_iy0 = -1; + double dX00 = 0; + double dY00 = 0; + double dZ00 = 0; + double dX01 = 0; + double dY01 = 0; + double dZ01 = 0; + double dX10 = 0; + double dY10 = 0; + double dZ10 = 0; + double dX11 = 0; + double dY11 = 0; + double dZ11 = 0; + double sinphi0 = 0; + double cosphi0 = 0; + double sinphi1 = 0; + double cosphi1 = 0; + + explicit GridEx(const Grid *gridIn) + : grid(gridIn), smallResx(grid->resx < DegToRad(1)), + sinhalfresx(sin(grid->resx / 2)), coshalfresx(cos(grid->resx / 2)), + sinresy(sin(grid->resy)), cosresy(cos(grid->resy)) {} + + // Return geocentric offset (dX, dY, dZ) relative to a point + // where x0 = -resx / 2 + inline void getBilinearGeocentric(int ix0, int iy0, double de00, + double dn00, double de01, double dn01, + double de10, double dn10, double de11, + double dn11, double m00, double m01, + double m10, double m11, double &dX, + double &dY, double &dZ) { + // If interpolating in the same cell as before, then we can skip + // the recomputation of dXij, dYij and dZij + if (ix0 != last_ix0 || iy0 != last_iy0) { + + last_ix0 = ix0; + if (iy0 != last_iy0) { + const double y0 = grid->miny + iy0 * grid->resy; + sinphi0 = sin(y0); + cosphi0 = cos(y0); + + // Use trigonometric formulas to avoid new calls to sin/cos + sinphi1 = + /* sin(y0+grid->resyRad) = */ sinphi0 * cosresy + + cosphi0 * sinresy; + cosphi1 = + /* cos(y0+grid->resyRad) = */ cosphi0 * cosresy - + sinphi0 * sinresy; + + last_iy0 = iy0; + } + + // Using "traditional" formulas to convert from easting, northing + // offsets to geocentric offsets + const double sinlam00 = -sinhalfresx; + const double coslam00 = coshalfresx; + const double sinphi00 = sinphi0; + const double cosphi00 = cosphi0; + const double dn00sinphi00 = dn00 * sinphi00; + dX00 = -de00 * sinlam00 - dn00sinphi00 * coslam00; + dY00 = de00 * coslam00 - dn00sinphi00 * sinlam00; + dZ00 = dn00 * cosphi00; + + const double sinlam01 = -sinhalfresx; + const double coslam01 = coshalfresx; + const double sinphi01 = sinphi1; + const double cosphi01 = cosphi1; + const double dn01sinphi01 = dn01 * sinphi01; + dX01 = -de01 * sinlam01 - dn01sinphi01 * coslam01; + dY01 = de01 * coslam01 - dn01sinphi01 * sinlam01; + dZ01 = dn01 * cosphi01; + + const double sinlam10 = sinhalfresx; + const double coslam10 = coshalfresx; + const double sinphi10 = sinphi0; + const double cosphi10 = cosphi0; + const double dn10sinphi10 = dn10 * sinphi10; + dX10 = -de10 * sinlam10 - dn10sinphi10 * coslam10; + dY10 = de10 * coslam10 - dn10sinphi10 * sinlam10; + dZ10 = dn10 * cosphi10; + + const double sinlam11 = sinhalfresx; + const double coslam11 = coshalfresx; + const double sinphi11 = sinphi1; + const double cosphi11 = cosphi1; + const double dn11sinphi11 = dn11 * sinphi11; + dX11 = -de11 * sinlam11 - dn11sinphi11 * coslam11; + dY11 = de11 * coslam11 - dn11sinphi11 * sinlam11; + dZ11 = dn11 * cosphi11; + } + + dX = m00 * dX00 + m01 * dX01 + m10 * dX10 + m11 * dX11; + dY = m00 * dY00 + m01 * dY01 + m10 * dY10 + m11 * dY11; + dZ = m00 * dZ00 + m01 * dZ01 + m10 * dZ10 + m11 * dZ11; + } +}; + +// --------------------------------------------------------------------------- + +/** Internal class to offer caching services over a Component */ +template <class Grid, class GridSet> struct ComponentEx { + const Component &component; + + const bool isBilinearInterpolation; /* bilinear vs geocentric_bilinear */ + + const DisplacementType displacementType; + + // Cache + std::unique_ptr<GridSet> gridSet{}; + std::map<const Grid *, GridEx<Grid>> mapGrids{}; + + private: + mutable double mCachedDt = 0; + mutable double mCachedValue = 0; + + static DisplacementType getDisplacementType(const std::string &s) { + if (s == STR_HORIZONTAL) + return DisplacementType::HORIZONTAL; + if (s == STR_VERTICAL) + return DisplacementType::VERTICAL; + if (s == STR_3D) + return DisplacementType::THREE_D; + return DisplacementType::NONE; + } + + public: + explicit ComponentEx(const Component &componentIn) + : component(componentIn), + isBilinearInterpolation( + componentIn.spatialModel().interpolationMethod == STR_BILINEAR), + displacementType(getDisplacementType(component.displacementType())) {} + + double evaluateAt(double dt) const { + if (dt == mCachedDt) + return mCachedValue; + mCachedDt = dt; + mCachedValue = component.timeFunction()->evaluateAt(dt); + return mCachedValue; + } + + void clearGridCache() { + gridSet.reset(); + mapGrids.clear(); + } +}; + +// --------------------------------------------------------------------------- + +/** Converts a ISO8601 date-time string, formatted as "YYYY-MM-DDTHH:MM:SSZ", + * into a decimal year. + * Leap years are taken into account, but not leap seconds. + */ +static double ISO8601ToDecimalYear(const std::string &dt) { + int year, month, day, hour, min, sec; + if (sscanf(dt.c_str(), "%04d-%02d-%02dT%02d:%02d:%02dZ", &year, &month, + &day, &hour, &min, &sec) != 6 || + year < 1582 || // Start of Gregorian calendar + month < 1 || month > 12 || day < 1 || day > 31 || hour < 0 || + hour >= 24 || min < 0 || min >= 60 || sec < 0 || sec >= 61) { + throw ParsingException("Wrong formatting / invalid date-time for " + + dt); + } + const bool isLeapYear = + (((year % 4) == 0 && (year % 100) != 0) || (year % 400) == 0); + // Given the intended use, we omit leap seconds... + const int month_table[2][12] = { + {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}, + {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}}; + int dayInYear = day - 1; + for (int m = 1; m < month; m++) { + dayInYear += month_table[isLeapYear ? 1 : 0][m - 1]; + } + if (day > month_table[isLeapYear ? 1 : 0][month - 1]) { + throw ParsingException("Wrong formatting / invalid date-time for " + + dt); + } + return year + + (dayInYear * 86400 + hour * 3600 + min * 60 + sec) / + (isLeapYear ? 86400. * 366 : 86400. * 365); +} + +// --------------------------------------------------------------------------- + +Epoch::Epoch(const std::string &dt) : mDt(dt) { + if (!dt.empty()) { + mDecimalYear = ISO8601ToDecimalYear(dt); + } +} + +// --------------------------------------------------------------------------- + +double Epoch::toDecimalYear() const { return mDecimalYear; } + +// --------------------------------------------------------------------------- + +static std::string getString(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::string(); + } + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json v = j[key]; + if (!v.is_string()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a string"); + } + return v.get<std::string>(); +} + +static std::string getReqString(const json &j, const char *key) { + return getString(j, key, false); +} + +static std::string getOptString(const json &j, const char *key) { + return getString(j, key, true); +} + +// --------------------------------------------------------------------------- + +static double getDouble(const json &j, const char *key, bool optional) { + if (!j.contains(key)) { + if (optional) { + return std::numeric_limits<double>::quiet_NaN(); + } + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json v = j[key]; + if (!v.is_number()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a number"); + } + return v.get<double>(); +} + +static double getReqDouble(const json &j, const char *key) { + return getDouble(j, key, false); +} +static double getOptDouble(const json &j, const char *key) { + return getDouble(j, key, true); +} + +// --------------------------------------------------------------------------- + +static json getObjectMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json obj = j[key]; + if (!obj.is_object()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a object"); + } + return obj; +} + +// --------------------------------------------------------------------------- + +static json getArrayMember(const json &j, const char *key) { + if (!j.contains(key)) { + throw ParsingException(std::string("Missing \"") + key + "\" key"); + } + const json obj = j[key]; + if (!obj.is_array()) { + throw ParsingException(std::string("The value of \"") + key + + "\" should be a array"); + } + return obj; +} + +// --------------------------------------------------------------------------- + +std::unique_ptr<MasterFile> MasterFile::parse(const std::string &text) { + std::unique_ptr<MasterFile> dmmf(new MasterFile()); + json j; + try { + j = json::parse(text); + } catch (const std::exception &e) { + throw ParsingException(e.what()); + } + if (!j.is_object()) { + throw ParsingException("Not an object"); + } + dmmf->mFileType = getReqString(j, "file_type"); + dmmf->mFormatVersion = getReqString(j, "format_version"); + dmmf->mName = getOptString(j, "name"); + dmmf->mVersion = getOptString(j, "version"); + dmmf->mLicense = getOptString(j, "license"); + dmmf->mDescription = getOptString(j, "description"); + dmmf->mPublicationDate = getOptString(j, "publication_date"); + + if (j.contains("authority")) { + const json jAuthority = j["authority"]; + if (!jAuthority.is_object()) { + throw ParsingException("authority is not a object"); + } + dmmf->mAuthority.name = getOptString(jAuthority, "name"); + dmmf->mAuthority.url = getOptString(jAuthority, "url"); + dmmf->mAuthority.address = getOptString(jAuthority, "address"); + dmmf->mAuthority.email = getOptString(jAuthority, "email"); + } + + if (j.contains("links")) { + const json jLinks = j["links"]; + if (!jLinks.is_array()) { + throw ParsingException("links is not an array"); + } + for (const json &jLink : jLinks) { + if (!jLink.is_object()) { + throw ParsingException("links[] item is not an object"); + } + Link link; + link.href = getOptString(jLink, "href"); + link.rel = getOptString(jLink, "rel"); + link.type = getOptString(jLink, "type"); + link.title = getOptString(jLink, "title"); + dmmf->mLinks.emplace_back(std::move(link)); + } + } + dmmf->mSourceCRS = getReqString(j, "source_crs"); + dmmf->mTargetCRS = getReqString(j, "target_crs"); + dmmf->mDefinitionCRS = getReqString(j, "definition_crs"); + if (dmmf->mSourceCRS != dmmf->mDefinitionCRS) { + throw ParsingException( + "source_crs != definition_crs not currently supported"); + } + dmmf->mReferenceEpoch = getOptString(j, "reference_epoch"); + dmmf->mUncertaintyReferenceEpoch = + getOptString(j, "uncertainty_reference_epoch"); + dmmf->mHorizontalOffsetUnit = getOptString(j, "horizontal_offset_unit"); + if (!dmmf->mHorizontalOffsetUnit.empty() && + dmmf->mHorizontalOffsetUnit != STR_METRE && + dmmf->mHorizontalOffsetUnit != STR_DEGREE) { + throw ParsingException("Unsupported value for horizontal_offset_unit"); + } + dmmf->mVerticalOffsetUnit = getOptString(j, "vertical_offset_unit"); + if (!dmmf->mVerticalOffsetUnit.empty() && + dmmf->mVerticalOffsetUnit != STR_METRE) { + throw ParsingException("Unsupported value for vertical_offset_unit"); + } + dmmf->mHorizontalUncertaintyType = + getOptString(j, "horizontal_uncertainty_type"); + dmmf->mHorizontalUncertaintyUnit = + getOptString(j, "horizontal_uncertainty_unit"); + dmmf->mVerticalUncertaintyType = + getOptString(j, "vertical_uncertainty_type"); + dmmf->mVerticalUncertaintyUnit = + getOptString(j, "vertical_uncertainty_unit"); + dmmf->mHorizontalOffsetMethod = getOptString(j, "horizontal_offset_method"); + if (!dmmf->mHorizontalOffsetMethod.empty() && + dmmf->mHorizontalOffsetMethod != STR_ADDITION && + dmmf->mHorizontalOffsetMethod != STR_GEOCENTRIC) { + throw ParsingException( + "Unsupported value for horizontal_offset_method"); + } + dmmf->mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); + + const json jTimeExtent = getObjectMember(j, "time_extent"); + dmmf->mTimeExtent.first = Epoch(getReqString(jTimeExtent, "first")); + dmmf->mTimeExtent.last = Epoch(getReqString(jTimeExtent, "last")); + + const json jComponents = getArrayMember(j, "components"); + for (const json &jComponent : jComponents) { + dmmf->mComponents.emplace_back(Component::parse(jComponent)); + const auto &comp = dmmf->mComponents.back(); + if (comp.displacementType() == STR_HORIZONTAL || + comp.displacementType() == STR_3D) { + if (dmmf->mHorizontalOffsetUnit.empty()) { + throw ParsingException("horizontal_offset_unit should be " + "defined as there is a component with " + "displacement_type = horizontal/3d"); + } + if (dmmf->mHorizontalOffsetMethod.empty()) { + throw ParsingException("horizontal_offset_method should be " + "defined as there is a component with " + "displacement_type = horizontal/3d"); + } + } + if (comp.displacementType() == STR_VERTICAL || + comp.displacementType() == STR_3D) { + if (dmmf->mVerticalOffsetUnit.empty()) { + throw ParsingException("vertical_offset_unit should be defined " + "as there is a component with " + "displacement_type = vertical/3d"); + } + } + if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && + comp.spatialModel().interpolationMethod != STR_BILINEAR) { + throw ParsingException("horizontal_offset_unit = degree can only " + "be used with interpolation_method = " + "bilinear"); + } + } + + if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && + dmmf->mHorizontalOffsetMethod != STR_ADDITION) { + throw ParsingException("horizontal_offset_unit = degree can only be " + "used with horizontal_offset_method = addition"); + } + + return dmmf; +} + +// --------------------------------------------------------------------------- + +SpatialExtent SpatialExtent::parse(const json &j) { + SpatialExtent ex; + + const std::string type = getReqString(j, "type"); + if (type != "bbox") { + throw ParsingException("unsupported type of extent"); + } + + const json jParameter = getObjectMember(j, "parameters"); + const json jBbox = getArrayMember(jParameter, "bbox"); + if (jBbox.size() != 4) { + throw ParsingException("bbox is not an array of 4 numeric elements"); + } + for (int i = 0; i < 4; i++) { + if (!jBbox[i].is_number()) { + throw ParsingException( + "bbox is not an array of 4 numeric elements"); + } + } + ex.mMinx = jBbox[0].get<double>(); + ex.mMiny = jBbox[1].get<double>(); + ex.mMaxx = jBbox[2].get<double>(); + ex.mMaxy = jBbox[3].get<double>(); + + ex.mMinxRad = DegToRad(ex.mMinx); + ex.mMinyRad = DegToRad(ex.mMiny); + ex.mMaxxRad = DegToRad(ex.mMaxx); + ex.mMaxyRad = DegToRad(ex.mMaxy); + + return ex; +} + +// --------------------------------------------------------------------------- + +Component Component::parse(const json &j) { + Component comp; + if (!j.is_object()) { + throw ParsingException("component is not an object"); + } + comp.mDescription = getOptString(j, "description"); + comp.mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); + comp.mDisplacementType = getReqString(j, "displacement_type"); + if (comp.mDisplacementType != STR_NONE && + comp.mDisplacementType != STR_HORIZONTAL && + comp.mDisplacementType != STR_VERTICAL && + comp.mDisplacementType != STR_3D) { + throw ParsingException("Unsupported value for displacement_type"); + } + comp.mUncertaintyType = getReqString(j, "uncertainty_type"); + comp.mHorizontalUncertainty = getOptDouble(j, "horizontal_uncertainty"); + comp.mVerticalUncertainty = getOptDouble(j, "vertical_uncertainty"); + + const json jSpatialModel = getObjectMember(j, "spatial_model"); + comp.mSpatialModel.type = getReqString(jSpatialModel, "type"); + comp.mSpatialModel.interpolationMethod = + getReqString(jSpatialModel, "interpolation_method"); + if (comp.mSpatialModel.interpolationMethod != STR_BILINEAR && + comp.mSpatialModel.interpolationMethod != STR_GEOCENTRIC_BILINEAR) { + throw ParsingException("Unsupported value for interpolation_method"); + } + comp.mSpatialModel.filename = getReqString(jSpatialModel, "filename"); + comp.mSpatialModel.md5Checksum = + getOptString(jSpatialModel, "md5_checksum"); + + const json jTimeFunction = getObjectMember(j, "time_function"); + const std::string timeFunctionType = getReqString(jTimeFunction, "type"); + const json jParameters = timeFunctionType == "constant" + ? json() + : getObjectMember(jTimeFunction, "parameters"); + + if (timeFunctionType == "constant") { + std::unique_ptr<ConstantTimeFunction> tf(new ConstantTimeFunction()); + tf->type = timeFunctionType; + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "velocity") { + std::unique_ptr<VelocityTimeFunction> tf(new VelocityTimeFunction()); + tf->type = timeFunctionType; + tf->referenceEpoch = + Epoch(getReqString(jParameters, "reference_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "step") { + std::unique_ptr<StepTimeFunction> tf(new StepTimeFunction()); + tf->type = timeFunctionType; + tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "reverse_step") { + std::unique_ptr<ReverseStepTimeFunction> tf( + new ReverseStepTimeFunction()); + tf->type = timeFunctionType; + tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "piecewise") { + std::unique_ptr<PiecewiseTimeFunction> tf(new PiecewiseTimeFunction()); + tf->type = timeFunctionType; + tf->beforeFirst = getReqString(jParameters, "before_first"); + if (tf->beforeFirst != "zero" && tf->beforeFirst != "constant" && + tf->beforeFirst != "linear") { + throw ParsingException("Unsupported value for before_first"); + } + tf->afterLast = getReqString(jParameters, "after_last"); + if (tf->afterLast != "zero" && tf->afterLast != "constant" && + tf->afterLast != "linear") { + throw ParsingException("Unsupported value for afterLast"); + } + const json jModel = getArrayMember(jParameters, "model"); + for (const json &jModelElt : jModel) { + if (!jModelElt.is_object()) { + throw ParsingException("model[] element is not an object"); + } + PiecewiseTimeFunction::EpochScaleFactorTuple tuple; + tuple.epoch = Epoch(getReqString(jModelElt, "epoch")); + tuple.scaleFactor = getReqDouble(jModelElt, "scale_factor"); + tf->model.emplace_back(std::move(tuple)); + } + + comp.mTimeFunction = std::move(tf); + } else if (timeFunctionType == "exponential") { + std::unique_ptr<ExponentialTimeFunction> tf( + new ExponentialTimeFunction()); + tf->type = timeFunctionType; + tf->referenceEpoch = + Epoch(getReqString(jParameters, "reference_epoch")); + tf->endEpoch = Epoch(getOptString(jParameters, "end_epoch")); + tf->relaxationConstant = + getReqDouble(jParameters, "relaxation_constant"); + if (tf->relaxationConstant <= 0.0) { + throw ParsingException("Invalid value for relaxation_constant"); + } + tf->beforeScaleFactor = + getReqDouble(jParameters, "before_scale_factor"); + tf->initialScaleFactor = + getReqDouble(jParameters, "initial_scale_factor"); + tf->finalScaleFactor = getReqDouble(jParameters, "final_scale_factor"); + comp.mTimeFunction = std::move(tf); + } else { + throw ParsingException("Unsupported type of time function: " + + timeFunctionType); + } + + return comp; +} + +// --------------------------------------------------------------------------- + +double Component::ConstantTimeFunction::evaluateAt(double) const { return 1.0; } + +// --------------------------------------------------------------------------- + +double Component::VelocityTimeFunction::evaluateAt(double dt) const { + return dt - referenceEpoch.toDecimalYear(); +} + +// --------------------------------------------------------------------------- + +double Component::StepTimeFunction::evaluateAt(double dt) const { + if (dt < stepEpoch.toDecimalYear()) + return 0.0; + return 1.0; +} + +// --------------------------------------------------------------------------- + +double Component::ReverseStepTimeFunction::evaluateAt(double dt) const { + if (dt < stepEpoch.toDecimalYear()) + return -1.0; + return 0.0; +} + +// --------------------------------------------------------------------------- + +double Component::PiecewiseTimeFunction::evaluateAt(double dt) const { + if (model.empty()) { + return 0.0; + } + + const double dt1 = model[0].epoch.toDecimalYear(); + if (dt < dt1) { + if (beforeFirst == "zero") + return 0.0; + if (beforeFirst == "constant" || model.size() == 1) + return model[0].scaleFactor; + + // linear + const double f1 = model[0].scaleFactor; + const double dt2 = model[1].epoch.toDecimalYear(); + const double f2 = model[1].scaleFactor; + if (dt1 == dt2) + return f1; + return (f1 * (dt2 - dt) + f2 * (dt - dt1)) / (dt2 - dt1); + } + for (size_t i = 1; i < model.size(); i++) { + const double dtip1 = model[i].epoch.toDecimalYear(); + if (dt < dtip1) { + const double dti = model[i - 1].epoch.toDecimalYear(); + const double fip1 = model[i].scaleFactor; + const double fi = model[i - 1].scaleFactor; + return (fi * (dtip1 - dt) + fip1 * (dt - dti)) / (dtip1 - dti); + } + } + if (afterLast == "zero") { + return 0.0; + } + if (afterLast == "constant" || model.size() == 1) + return model.back().scaleFactor; + + // linear + const double dtnm1 = model[model.size() - 2].epoch.toDecimalYear(); + const double fnm1 = model[model.size() - 2].scaleFactor; + const double dtn = model.back().epoch.toDecimalYear(); + const double fn = model.back().scaleFactor; + if (dtnm1 == dtn) + return fn; + return (fnm1 * (dtn - dt) + fn * (dt - dtnm1)) / (dtn - dtnm1); +} + +// --------------------------------------------------------------------------- + +double Component::ExponentialTimeFunction::evaluateAt(double dt) const { + const double t0 = referenceEpoch.toDecimalYear(); + if (dt < t0) + return beforeScaleFactor; + if (!endEpoch.toString().empty()) { + dt = std::min(dt, endEpoch.toDecimalYear()); + } + return initialScaleFactor + + (finalScaleFactor - initialScaleFactor) * + (1.0 - std::exp(-(dt - t0) / relaxationConstant)); +} + +// --------------------------------------------------------------------------- + +inline void DeltaEastingNorthingToLongLat(double cosphi, double de, double dn, + double a, double b, double es, + double &dlam, double &dphi) { + const double oneMinuX = es * (1 - cosphi * cosphi); + const double X = 1 - oneMinuX; + const double sqrtX = sqrt(X); +#if 0 + // With es of Earth, absolute/relative error is at most 2e-8 + // const double sqrtX = 1 - 0.5 * oneMinuX * (1 + 0.25 * oneMinuX); +#endif + dlam = de * sqrtX / (a * cosphi); + dphi = dn * a * sqrtX * X / (b * b); +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +Evaluator<Grid, GridSet, EvaluatorIface>::Evaluator( + std::unique_ptr<MasterFile> &&model, EvaluatorIface &iface, double a, + double b) + : mModel(std::move(model)), mA(a), mB(b), mEs(1 - (b * b) / (a * a)), + mIsHorizontalUnitDegree(mModel->horizontalOffsetUnit() == STR_DEGREE), + mIsAddition(mModel->horizontalOffsetMethod() == STR_ADDITION), + mIsGeographicCRS(iface.isGeographicCRS(mModel->definitionCRS())) { + if (!mIsGeographicCRS && mIsHorizontalUnitDegree) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "horizontal_offset_unit = degree are incompatible"); + } + if (!mIsGeographicCRS && !mIsAddition) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "horizontal_offset_method = geocentric are incompatible"); + } + mComponents.reserve(mModel->components().size()); + for (const auto &comp : mModel->components()) { + mComponents.emplace_back(std::unique_ptr<ComponentEx<Grid, GridSet>>( + new ComponentEx<Grid, GridSet>(comp))); + if (!mIsGeographicCRS && !mComponents.back()->isBilinearInterpolation) { + throw EvaluatorException( + "definition_crs = projected CRS and " + "interpolation_method = geocentric_bilinear are incompatible"); + } + } +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +void Evaluator<Grid, GridSet, EvaluatorIface>::clearGridCache() { + for (auto &comp : mComponents) { + comp->clearGridCache(); + } +} + +// --------------------------------------------------------------------------- + +#ifdef DEBUG_DEFMODEL + +static std::string shortName(const Component &comp) { + const auto &desc = comp.description(); + return desc.substr(0, desc.find('\n')) + " (" + + comp.spatialModel().filename + ")"; +} + +static std::string toString(double val) { + char buffer[32]; + snprintf(buffer, sizeof(buffer), "%.9g", val); + return buffer; +} + +#endif + +// --------------------------------------------------------------------------- + +static bool bboxCheck(double &x, double &y, bool forInverseComputation, + const double minx, const double miny, const double maxx, + const double maxy, const double EPS, + const double extraMarginForInverse) { + if (x < minx - EPS || x > maxx + EPS || y < miny - EPS || y > maxy + EPS) { + if (!forInverseComputation) { + return false; + } + // In case of iterative computation for inverse, allow to be a + // slightly bit outside of the grid and clamp to the edges + bool xOk = false; + if (x >= minx - EPS && x <= maxx + EPS) { + xOk = true; + } else if (x > minx - extraMarginForInverse && x < minx) { + x = minx; + xOk = true; + } else if (x < maxx + extraMarginForInverse && x > maxx) { + x = maxx; + xOk = true; + } + + bool yOk = false; + if (y >= miny - EPS && y <= maxy + EPS) { + yOk = true; + } else if (y > miny - extraMarginForInverse && y < miny) { + y = miny; + yOk = true; + } else if (y < maxy + extraMarginForInverse && y > maxy) { + y = maxy; + yOk = true; + } + + return xOk && yOk; + } + return true; +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +bool Evaluator<Grid, GridSet, EvaluatorIface>::forward( + EvaluatorIface &iface, double x, double y, double z, double t, + bool forInverseComputation, double &x_out, double &y_out, double &z_out) + +{ + x_out = x; + y_out = y; + z_out = z; + + const double EPS = mIsGeographicCRS ? 1e-10 : 1e-5; + + // Check against global model spatial extent, potentially wrapping + // longitude to match + { + const auto &extent = mModel->extent(); + const double minx = extent.minxNormalized(mIsGeographicCRS); + const double maxx = extent.maxxNormalized(mIsGeographicCRS); + if (mIsGeographicCRS) { + while (x < minx - EPS) { + x += 2.0 * DEFMODEL_PI; + } + while (x > maxx + EPS) { + x -= 2.0 * DEFMODEL_PI; + } + } + const double miny = extent.minyNormalized(mIsGeographicCRS); + const double maxy = extent.maxyNormalized(mIsGeographicCRS); + const double extraMarginForInverse = + mIsGeographicCRS ? DegToRad(0.1) : 10000; + if (!bboxCheck(x, y, forInverseComputation, minx, miny, maxx, maxy, EPS, + extraMarginForInverse)) { +#ifdef DEBUG_DEFMODEL + iface.log("Calculation point " + toString(x) + "," + toString(y) + + " is outside the extents of the deformation model"); +#endif + return false; + } + } + + // Check against global model temporal extent + { + const auto &timeExtent = mModel->timeExtent(); + if (t < timeExtent.first.toDecimalYear() || + t > timeExtent.last.toDecimalYear()) { +#ifdef DEBUG_DEFMODEL + iface.log("Calculation epoch " + toString(t) + + " is not valid for the deformation model"); +#endif + return false; + } + } + + // For mIsHorizontalUnitDegree + double dlam = 0; + double dphi = 0; + + // For !mIsHorizontalUnitDegree + double de = 0; + double dn = 0; + + double dz = 0; + + bool sincosphiInitialized = false; + double sinphi = 0; + double cosphi = 0; + + for (auto &compEx : mComponents) { + const auto &comp = compEx->component; + if (compEx->displacementType == DisplacementType::NONE) { + continue; + } + const auto &extent = comp.extent(); + double xForGrid = x; + double yForGrid = y; + const double minx = extent.minxNormalized(mIsGeographicCRS); + const double maxx = extent.maxxNormalized(mIsGeographicCRS); + const double miny = extent.minyNormalized(mIsGeographicCRS); + const double maxy = extent.maxyNormalized(mIsGeographicCRS); + const double extraMarginForInverse = 0; + if (!bboxCheck(xForGrid, yForGrid, forInverseComputation, minx, miny, + maxx, maxy, EPS, extraMarginForInverse)) { +#ifdef DEBUG_DEFMODEL + iface.log( + "Skipping component " + shortName(comp) + + " due to point being outside of its declared spatial extent."); +#endif + continue; + } + xForGrid = std::max(xForGrid, minx); + yForGrid = std::max(yForGrid, miny); + xForGrid = std::min(xForGrid, maxx); + yForGrid = std::min(yForGrid, maxy); + const auto tfactor = compEx->evaluateAt(t); + if (tfactor == 0.0) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to time function evaluating to 0."); +#endif + continue; + } + +#ifdef DEBUG_DEFMODEL + iface.log("Entering component " + shortName(comp) + + " with time function evaluating to " + toString(tfactor) + + "."); +#endif + + if (compEx->gridSet == nullptr) { + compEx->gridSet = iface.open(comp.spatialModel().filename); + if (compEx->gridSet == nullptr) { + return false; + } + } + const Grid *grid = compEx->gridSet->gridAt(xForGrid, yForGrid); + if (grid == nullptr) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to no grid found for this point in the grid set."); +#endif + continue; + } + if (grid->width < 2 || grid->height < 2) { + return false; + } + const double ix_d = (xForGrid - grid->minx) / grid->resx; + const double iy_d = (yForGrid - grid->miny) / grid->resy; + if (ix_d < -EPS || iy_d < -EPS || ix_d + 1 >= grid->width + EPS || + iy_d + 1 >= grid->height + EPS) { +#ifdef DEBUG_DEFMODEL + iface.log("Skipping component " + shortName(comp) + + " due to point being outside of actual spatial extent of " + "grid " + + grid->name() + "."); +#endif + continue; + } + const int ix0 = std::min(static_cast<int>(ix_d), grid->width - 2); + const int iy0 = std::min(static_cast<int>(iy_d), grid->height - 2); + const int ix1 = ix0 + 1; + const int iy1 = iy0 + 1; + const double frct_x = ix_d - ix0; + const double frct_y = iy_d - iy0; + const double one_minus_frct_x = 1. - frct_x; + const double one_minus_frct_y = 1. - frct_y; + const double m00 = one_minus_frct_x * one_minus_frct_y; + const double m10 = frct_x * one_minus_frct_y; + const double m01 = one_minus_frct_x * frct_y; + const double m11 = frct_x * frct_y; + + if (compEx->displacementType == DisplacementType::VERTICAL) { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getZOffset(ix0, iy0, dz00) || + !grid->getZOffset(ix1, iy0, dz10) || + !grid->getZOffset(ix0, iy1, dz01) || + !grid->getZOffset(ix1, iy1, dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } else if (mIsHorizontalUnitDegree) { + double dx00 = 0; + double dy00 = 0; + double dx01 = 0; + double dy01 = 0; + double dx10 = 0; + double dy10 = 0; + double dx11 = 0; + double dy11 = 0; + if (compEx->displacementType == DisplacementType::HORIZONTAL) { + if (!grid->getLonLatOffset(ix0, iy0, dx00, dy00) || + !grid->getLonLatOffset(ix1, iy0, dx10, dy10) || + !grid->getLonLatOffset(ix0, iy1, dx01, dy01) || + !grid->getLonLatOffset(ix1, iy1, dx11, dy11)) { + return false; + } + } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getLonLatZOffset(ix0, iy0, dx00, dy00, dz00) || + !grid->getLonLatZOffset(ix1, iy0, dx10, dy10, dz10) || + !grid->getLonLatZOffset(ix0, iy1, dx01, dy01, dz01) || + !grid->getLonLatZOffset(ix1, iy1, dx11, dy11, dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } + const double dlamInterp = + dx00 * m00 + dx01 * m01 + dx10 * m10 + dx11 * m11; + const double dphiInterp = + dy00 * m00 + dy01 * m01 + dy10 * m10 + dy11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dlamInterp = " + toString(tfactor) + " * " + + toString(dlamInterp) + "."); + iface.log("tfactor * dphiInterp = " + toString(tfactor) + " * " + + toString(dphiInterp) + "."); +#endif + dlam += tfactor * dlamInterp; + dphi += tfactor * dphiInterp; + } else /* horizontal unit is metre */ { + double de00 = 0; + double dn00 = 0; + double de01 = 0; + double dn01 = 0; + double de10 = 0; + double dn10 = 0; + double de11 = 0; + double dn11 = 0; + if (compEx->displacementType == DisplacementType::HORIZONTAL) { + if (!grid->getEastingNorthingOffset(ix0, iy0, de00, dn00) || + !grid->getEastingNorthingOffset(ix1, iy0, de10, dn10) || + !grid->getEastingNorthingOffset(ix0, iy1, de01, dn01) || + !grid->getEastingNorthingOffset(ix1, iy1, de11, dn11)) { + return false; + } + } else /* if (compEx->displacementType == DisplacementType::THREE_D) */ { + double dz00 = 0; + double dz01 = 0; + double dz10 = 0; + double dz11 = 0; + if (!grid->getEastingNorthingZOffset(ix0, iy0, de00, dn00, + dz00) || + !grid->getEastingNorthingZOffset(ix1, iy0, de10, dn10, + dz10) || + !grid->getEastingNorthingZOffset(ix0, iy1, de01, dn01, + dz01) || + !grid->getEastingNorthingZOffset(ix1, iy1, de11, dn11, + dz11)) { + return false; + } + const double dzInterp = + dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + + toString(dzInterp) + "."); +#endif + dz += tfactor * dzInterp; + } + if (compEx->isBilinearInterpolation) { + const double deInterp = + de00 * m00 + de01 * m01 + de10 * m10 + de11 * m11; + const double dnInterp = + dn00 * m00 + dn01 * m01 + dn10 * m10 + dn11 * m11; +#ifdef DEBUG_DEFMODEL + iface.log("tfactor * deInterp = " + toString(tfactor) + " * " + + toString(deInterp) + "."); + iface.log("tfactor * dnInterp = " + toString(tfactor) + " * " + + toString(dnInterp) + "."); +#endif + de += tfactor * deInterp; + dn += tfactor * dnInterp; + } else /* geocentric_bilinear */ { + double dX; + double dY; + double dZ; + + auto iter = compEx->mapGrids.find(grid); + if (iter == compEx->mapGrids.end()) { + GridEx<Grid> gridWithCache(grid); + iter = compEx->mapGrids + .insert(std::pair<const Grid *, GridEx<Grid>>( + grid, std::move(gridWithCache))) + .first; + } + GridEx<Grid> &gridwithCacheRef = iter->second; + + gridwithCacheRef.getBilinearGeocentric( + ix0, iy0, de00, dn00, de01, dn01, de10, dn10, de11, dn11, + m00, m01, m10, m11, dX, dY, dZ); + if (!sincosphiInitialized) { + sincosphiInitialized = true; + sinphi = sin(y); + cosphi = cos(y); + } + const double lam_rel_to_cell_center = + (frct_x - 0.5) * grid->resx; + // Use small-angle approximation of sin/cos when reasonable + // Max abs/rel error on cos is 3.9e-9 and on sin 1.3e-11 + const double sinlam = + gridwithCacheRef.smallResx + ? lam_rel_to_cell_center * + (1 - + (1. / 6) * (lam_rel_to_cell_center * + lam_rel_to_cell_center)) + : sin(lam_rel_to_cell_center); + const double coslam = gridwithCacheRef.smallResx + ? (1 - + 0.5 * (lam_rel_to_cell_center * + lam_rel_to_cell_center)) + : cos(lam_rel_to_cell_center); + + // Convert back from geocentric deltas to easting, northing + // deltas + const double deInterp = -dX * sinlam + dY * coslam; + const double dnInterp = + (-dX * coslam - dY * sinlam) * sinphi + dZ * cosphi; +#ifdef DEBUG_DEFMODEL + iface.log("After geocentric_bilinear interpolation: tfactor * " + "deInterp = " + + toString(tfactor) + " * " + toString(deInterp) + "."); + iface.log("After geocentric_bilinear interpolation: tfactor * " + "dnInterp = " + + toString(tfactor) + " * " + toString(dnInterp) + "."); +#endif + de += tfactor * deInterp; + dn += tfactor * dnInterp; + } + } + } + + // Apply shifts depending on horizontal_offset_unit and + // horizontal_offset_method + if (mIsHorizontalUnitDegree) { + x_out += dlam; + y_out += dphi; + } else { +#ifdef DEBUG_DEFMODEL + iface.log("Total sum of de: " + toString(de)); + iface.log("Total sum of dn: " + toString(dn)); +#endif + if (mIsAddition && !mIsGeographicCRS) { + x_out += de; + y_out += dn; + } else if (mIsAddition) { + // Simple way of adding the offset + if (!sincosphiInitialized) { + cosphi = cos(y); + } + DeltaEastingNorthingToLongLat(cosphi, de, dn, mA, mB, mEs, dlam, + dphi); +#ifdef DEBUG_DEFMODEL + iface.log("Result dlam: " + toString(dlam)); + iface.log("Result dphi: " + toString(dphi)); +#endif + x_out += dlam; + y_out += dphi; + } else { + // Geocentric way of adding the offset + if (!sincosphiInitialized) { + sinphi = sin(y); + cosphi = cos(y); + } + const double sinlam = sin(x); + const double coslam = cos(x); + const double dnsinphi = dn * sinphi; + const double dX = -de * sinlam - dnsinphi * coslam; + const double dY = de * coslam - dnsinphi * sinlam; + const double dZ = dn * cosphi; + + double X; + double Y; + double Z; + iface.geographicToGeocentric(x, y, 0, mA, mB, mEs, X, Y, Z); +#ifdef DEBUG_DEFMODEL + iface.log("Geocentric coordinate before: " + toString(X) + "," + + toString(Y) + "," + toString(Z)); + iface.log("Geocentric shift: " + toString(dX) + "," + toString(dY) + + "," + toString(dZ)); +#endif + X += dX; + Y += dY; + Z += dZ; +#ifdef DEBUG_DEFMODEL + iface.log("Geocentric coordinate after: " + toString(X) + "," + + toString(Y) + "," + toString(Z)); +#endif + + double h_out_ignored; + iface.geocentricToGeographic(X, Y, Z, mA, mB, mEs, x_out, y_out, + h_out_ignored); + } + } +#ifdef DEBUG_DEFMODEL + iface.log("Total sum of dz: " + toString(dz)); +#endif + z_out += dz; + + return true; +} + +// --------------------------------------------------------------------------- + +template <class Grid, class GridSet, class EvaluatorIface> +bool Evaluator<Grid, GridSet, EvaluatorIface>::inverse( + EvaluatorIface &iface, double x, double y, double z, double t, + double &x_out, double &y_out, double &z_out) + +{ + x_out = x; + y_out = y; + z_out = z; + constexpr double EPS_HORIZ = 1e-12; + constexpr double EPS_VERT = 1e-3; + constexpr bool forInverseComputation = true; + for (int i = 0; i < 10; i++) { +#ifdef DEBUG_DEFMODEL + iface.log("Iteration " + std::to_string(i) + ": before foward: x=" + + toString(x_out) + ", y=" + toString(y_out)); +#endif + double x_new; + double y_new; + double z_new; + if (!forward(iface, x_out, y_out, z_out, t, forInverseComputation, + x_new, y_new, z_new)) { + return false; + } +#ifdef DEBUG_DEFMODEL + iface.log("After foward: x=" + toString(x_new) + ", y=" + + toString(y_new)); +#endif + const double dx = x_new - x; + const double dy = y_new - y; + const double dz = z_new - z; + x_out -= dx; + y_out -= dy; + z_out -= dz; + if (std::max(std::fabs(dx), std::fabs(dy)) < EPS_HORIZ && + std::fabs(dz) < EPS_VERT) { + return true; + } + } + return false; +} + +// --------------------------------------------------------------------------- + +} // namespace DEFORMATON_MODEL_NAMESPACE |