/src/PROJ/src/transformations/defmodel_impl.hpp
Line | Count | Source (jump to first uncovered line) |
1 | | /****************************************************************************** |
2 | | * Project: PROJ |
3 | | * Purpose: Functionality related to deformation model |
4 | | * Author: Even Rouault, <even.rouault at spatialys.com> |
5 | | * |
6 | | ****************************************************************************** |
7 | | * Copyright (c) 2020, Even Rouault, <even.rouault at spatialys.com> |
8 | | * |
9 | | * Permission is hereby granted, free of charge, to any person obtaining a |
10 | | * copy of this software and associated documentation files (the "Software"), |
11 | | * to deal in the Software without restriction, including without limitation |
12 | | * the rights to use, copy, modify, merge, publish, distribute, sublicense, |
13 | | * and/or sell copies of the Software, and to permit persons to whom the |
14 | | * Software is furnished to do so, subject to the following conditions: |
15 | | * |
16 | | * The above copyright notice and this permission notice shall be included |
17 | | * in all copies or substantial portions of the Software. |
18 | | * |
19 | | * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS |
20 | | * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
21 | | * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL |
22 | | * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
23 | | * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
24 | | * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER |
25 | | * DEALINGS IN THE SOFTWARE. |
26 | | *****************************************************************************/ |
27 | | |
28 | | #ifndef DEFORMATON_MODEL_NAMESPACE |
29 | | #error "Should be included only by defmodel.hpp" |
30 | | #endif |
31 | | |
32 | | namespace DEFORMATON_MODEL_NAMESPACE { |
33 | | |
34 | | // --------------------------------------------------------------------------- |
35 | | |
36 | | static const std::string STR_DEGREE("degree"); |
37 | | static const std::string STR_METRE("metre"); |
38 | | |
39 | | static const std::string STR_ADDITION("addition"); |
40 | | static const std::string STR_GEOCENTRIC("geocentric"); |
41 | | |
42 | | static const std::string STR_BILINEAR("bilinear"); |
43 | | static const std::string STR_GEOCENTRIC_BILINEAR("geocentric_bilinear"); |
44 | | |
45 | | static const std::string STR_NONE("none"); |
46 | | static const std::string STR_HORIZONTAL("horizontal"); |
47 | | static const std::string STR_VERTICAL("vertical"); |
48 | | static const std::string STR_3D("3d"); |
49 | | |
50 | | constexpr double DEFMODEL_PI = 3.14159265358979323846; |
51 | | constexpr double DEG_TO_RAD_CONSTANT = 3.14159265358979323846 / 180.; |
52 | 0 | inline constexpr double DegToRad(double d) { return d * DEG_TO_RAD_CONSTANT; } |
53 | | |
54 | | // --------------------------------------------------------------------------- |
55 | | |
56 | | enum class DisplacementType { NONE, HORIZONTAL, VERTICAL, THREE_D }; |
57 | | |
58 | | // --------------------------------------------------------------------------- |
59 | | |
60 | | /** Internal class to offer caching services over a Grid */ |
61 | | template <class Grid> struct GridEx { |
62 | | const Grid *grid; |
63 | | |
64 | | bool smallResx; // lesser than one degree |
65 | | |
66 | | double sinhalfresx; |
67 | | double coshalfresx; |
68 | | |
69 | | double sinresy; |
70 | | double cosresy; |
71 | | |
72 | | int last_ix0 = -1; |
73 | | int last_iy0 = -1; |
74 | | double dX00 = 0; |
75 | | double dY00 = 0; |
76 | | double dZ00 = 0; |
77 | | double dX01 = 0; |
78 | | double dY01 = 0; |
79 | | double dZ01 = 0; |
80 | | double dX10 = 0; |
81 | | double dY10 = 0; |
82 | | double dZ10 = 0; |
83 | | double dX11 = 0; |
84 | | double dY11 = 0; |
85 | | double dZ11 = 0; |
86 | | double sinphi0 = 0; |
87 | | double cosphi0 = 0; |
88 | | double sinphi1 = 0; |
89 | | double cosphi1 = 0; |
90 | | |
91 | | explicit GridEx(const Grid *gridIn) |
92 | | : grid(gridIn), smallResx(grid->resx < DegToRad(1)), |
93 | | sinhalfresx(sin(grid->resx / 2)), coshalfresx(cos(grid->resx / 2)), |
94 | 0 | sinresy(sin(grid->resy)), cosresy(cos(grid->resy)) {} |
95 | | |
96 | | // Return geocentric offset (dX, dY, dZ) relative to a point |
97 | | // where x0 = -resx / 2 |
98 | | inline void getBilinearGeocentric(int ix0, int iy0, double de00, |
99 | | double dn00, double de01, double dn01, |
100 | | double de10, double dn10, double de11, |
101 | | double dn11, double m00, double m01, |
102 | | double m10, double m11, double &dX, |
103 | 0 | double &dY, double &dZ) { |
104 | | // If interpolating in the same cell as before, then we can skip |
105 | | // the recomputation of dXij, dYij and dZij |
106 | 0 | if (ix0 != last_ix0 || iy0 != last_iy0) { |
107 | |
|
108 | 0 | last_ix0 = ix0; |
109 | 0 | if (iy0 != last_iy0) { |
110 | 0 | const double y0 = grid->miny + iy0 * grid->resy; |
111 | 0 | sinphi0 = sin(y0); |
112 | 0 | cosphi0 = cos(y0); |
113 | | |
114 | | // Use trigonometric formulas to avoid new calls to sin/cos |
115 | 0 | sinphi1 = |
116 | 0 | /* sin(y0+grid->resyRad) = */ sinphi0 * cosresy + |
117 | 0 | cosphi0 * sinresy; |
118 | 0 | cosphi1 = |
119 | 0 | /* cos(y0+grid->resyRad) = */ cosphi0 * cosresy - |
120 | 0 | sinphi0 * sinresy; |
121 | |
|
122 | 0 | last_iy0 = iy0; |
123 | 0 | } |
124 | | |
125 | | // Using "traditional" formulas to convert from easting, northing |
126 | | // offsets to geocentric offsets |
127 | 0 | const double sinlam00 = -sinhalfresx; |
128 | 0 | const double coslam00 = coshalfresx; |
129 | 0 | const double sinphi00 = sinphi0; |
130 | 0 | const double cosphi00 = cosphi0; |
131 | 0 | const double dn00sinphi00 = dn00 * sinphi00; |
132 | 0 | dX00 = -de00 * sinlam00 - dn00sinphi00 * coslam00; |
133 | 0 | dY00 = de00 * coslam00 - dn00sinphi00 * sinlam00; |
134 | 0 | dZ00 = dn00 * cosphi00; |
135 | |
|
136 | 0 | const double sinlam01 = -sinhalfresx; |
137 | 0 | const double coslam01 = coshalfresx; |
138 | 0 | const double sinphi01 = sinphi1; |
139 | 0 | const double cosphi01 = cosphi1; |
140 | 0 | const double dn01sinphi01 = dn01 * sinphi01; |
141 | 0 | dX01 = -de01 * sinlam01 - dn01sinphi01 * coslam01; |
142 | 0 | dY01 = de01 * coslam01 - dn01sinphi01 * sinlam01; |
143 | 0 | dZ01 = dn01 * cosphi01; |
144 | |
|
145 | 0 | const double sinlam10 = sinhalfresx; |
146 | 0 | const double coslam10 = coshalfresx; |
147 | 0 | const double sinphi10 = sinphi0; |
148 | 0 | const double cosphi10 = cosphi0; |
149 | 0 | const double dn10sinphi10 = dn10 * sinphi10; |
150 | 0 | dX10 = -de10 * sinlam10 - dn10sinphi10 * coslam10; |
151 | 0 | dY10 = de10 * coslam10 - dn10sinphi10 * sinlam10; |
152 | 0 | dZ10 = dn10 * cosphi10; |
153 | |
|
154 | 0 | const double sinlam11 = sinhalfresx; |
155 | 0 | const double coslam11 = coshalfresx; |
156 | 0 | const double sinphi11 = sinphi1; |
157 | 0 | const double cosphi11 = cosphi1; |
158 | 0 | const double dn11sinphi11 = dn11 * sinphi11; |
159 | 0 | dX11 = -de11 * sinlam11 - dn11sinphi11 * coslam11; |
160 | 0 | dY11 = de11 * coslam11 - dn11sinphi11 * sinlam11; |
161 | 0 | dZ11 = dn11 * cosphi11; |
162 | 0 | } |
163 | |
|
164 | 0 | dX = m00 * dX00 + m01 * dX01 + m10 * dX10 + m11 * dX11; |
165 | 0 | dY = m00 * dY00 + m01 * dY01 + m10 * dY10 + m11 * dY11; |
166 | 0 | dZ = m00 * dZ00 + m01 * dZ01 + m10 * dZ10 + m11 * dZ11; |
167 | 0 | } |
168 | | }; |
169 | | |
170 | | // --------------------------------------------------------------------------- |
171 | | |
172 | | /** Internal class to offer caching services over a Component */ |
173 | | template <class Grid, class GridSet> struct ComponentEx { |
174 | | const Component &component; |
175 | | |
176 | | const bool isBilinearInterpolation; /* bilinear vs geocentric_bilinear */ |
177 | | |
178 | | const DisplacementType displacementType; |
179 | | |
180 | | // Cache |
181 | | std::unique_ptr<GridSet> gridSet{}; |
182 | | std::map<const Grid *, GridEx<Grid>> mapGrids{}; |
183 | | |
184 | | private: |
185 | | mutable double mCachedDt = 0; |
186 | | mutable double mCachedValue = 0; |
187 | | |
188 | 0 | static DisplacementType getDisplacementType(const std::string &s) { |
189 | 0 | if (s == STR_HORIZONTAL) |
190 | 0 | return DisplacementType::HORIZONTAL; |
191 | 0 | if (s == STR_VERTICAL) |
192 | 0 | return DisplacementType::VERTICAL; |
193 | 0 | if (s == STR_3D) |
194 | 0 | return DisplacementType::THREE_D; |
195 | 0 | return DisplacementType::NONE; |
196 | 0 | } |
197 | | |
198 | | public: |
199 | | explicit ComponentEx(const Component &componentIn) |
200 | | : component(componentIn), |
201 | | isBilinearInterpolation( |
202 | | componentIn.spatialModel().interpolationMethod == STR_BILINEAR), |
203 | 0 | displacementType(getDisplacementType(component.displacementType())) {} |
204 | | |
205 | 0 | double evaluateAt(double dt) const { |
206 | 0 | if (dt == mCachedDt) |
207 | 0 | return mCachedValue; |
208 | 0 | mCachedDt = dt; |
209 | 0 | mCachedValue = component.timeFunction()->evaluateAt(dt); |
210 | 0 | return mCachedValue; |
211 | 0 | } |
212 | | |
213 | 0 | void clearGridCache() { |
214 | 0 | gridSet.reset(); |
215 | 0 | mapGrids.clear(); |
216 | 0 | } |
217 | | }; |
218 | | |
219 | | // --------------------------------------------------------------------------- |
220 | | |
221 | | /** Converts a ISO8601 date-time string, formatted as "YYYY-MM-DDTHH:MM:SSZ", |
222 | | * into a decimal year. |
223 | | * Leap years are taken into account, but not leap seconds. |
224 | | */ |
225 | 0 | static double ISO8601ToDecimalYear(const std::string &dt) { |
226 | 0 | int year, month, day, hour, min, sec; |
227 | 0 | if (sscanf(dt.c_str(), "%04d-%02d-%02dT%02d:%02d:%02dZ", &year, &month, |
228 | 0 | &day, &hour, &min, &sec) != 6 || |
229 | 0 | year < 1582 || // Start of Gregorian calendar |
230 | 0 | month < 1 || month > 12 || day < 1 || day > 31 || hour < 0 || |
231 | 0 | hour >= 24 || min < 0 || min >= 60 || sec < 0 || sec >= 61) { |
232 | 0 | throw ParsingException("Wrong formatting / invalid date-time for " + |
233 | 0 | dt); |
234 | 0 | } |
235 | 0 | const bool isLeapYear = |
236 | 0 | (((year % 4) == 0 && (year % 100) != 0) || (year % 400) == 0); |
237 | | // Given the intended use, we omit leap seconds... |
238 | 0 | const int month_table[2][12] = { |
239 | 0 | {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}, |
240 | 0 | {31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}}; |
241 | 0 | int dayInYear = day - 1; |
242 | 0 | for (int m = 1; m < month; m++) { |
243 | 0 | dayInYear += month_table[isLeapYear ? 1 : 0][m - 1]; |
244 | 0 | } |
245 | 0 | if (day > month_table[isLeapYear ? 1 : 0][month - 1]) { |
246 | 0 | throw ParsingException("Wrong formatting / invalid date-time for " + |
247 | 0 | dt); |
248 | 0 | } |
249 | 0 | return year + (dayInYear * 86400 + hour * 3600 + min * 60 + sec) / |
250 | 0 | (isLeapYear ? 86400. * 366 : 86400. * 365); |
251 | 0 | } |
252 | | |
253 | | // --------------------------------------------------------------------------- |
254 | | |
255 | 0 | Epoch::Epoch(const std::string &dt) : mDt(dt) { |
256 | 0 | if (!dt.empty()) { |
257 | 0 | mDecimalYear = ISO8601ToDecimalYear(dt); |
258 | 0 | } |
259 | 0 | } |
260 | | |
261 | | // --------------------------------------------------------------------------- |
262 | | |
263 | 0 | double Epoch::toDecimalYear() const { return mDecimalYear; } |
264 | | |
265 | | // --------------------------------------------------------------------------- |
266 | | |
267 | 0 | static std::string getString(const json &j, const char *key, bool optional) { |
268 | 0 | if (!j.contains(key)) { |
269 | 0 | if (optional) { |
270 | 0 | return std::string(); |
271 | 0 | } |
272 | 0 | throw ParsingException(std::string("Missing \"") + key + "\" key"); |
273 | 0 | } |
274 | 0 | const json v = j[key]; |
275 | 0 | if (!v.is_string()) { |
276 | 0 | throw ParsingException(std::string("The value of \"") + key + |
277 | 0 | "\" should be a string"); |
278 | 0 | } |
279 | 0 | return v.get<std::string>(); |
280 | 0 | } |
281 | | |
282 | 0 | static std::string getReqString(const json &j, const char *key) { |
283 | 0 | return getString(j, key, false); |
284 | 0 | } |
285 | | |
286 | 0 | static std::string getOptString(const json &j, const char *key) { |
287 | 0 | return getString(j, key, true); |
288 | 0 | } |
289 | | |
290 | | // --------------------------------------------------------------------------- |
291 | | |
292 | 0 | static double getDouble(const json &j, const char *key, bool optional) { |
293 | 0 | if (!j.contains(key)) { |
294 | 0 | if (optional) { |
295 | 0 | return std::numeric_limits<double>::quiet_NaN(); |
296 | 0 | } |
297 | 0 | throw ParsingException(std::string("Missing \"") + key + "\" key"); |
298 | 0 | } |
299 | 0 | const json v = j[key]; |
300 | 0 | if (!v.is_number()) { |
301 | 0 | throw ParsingException(std::string("The value of \"") + key + |
302 | 0 | "\" should be a number"); |
303 | 0 | } |
304 | 0 | return v.get<double>(); |
305 | 0 | } |
306 | | |
307 | 0 | static double getReqDouble(const json &j, const char *key) { |
308 | 0 | return getDouble(j, key, false); |
309 | 0 | } |
310 | 0 | static double getOptDouble(const json &j, const char *key) { |
311 | 0 | return getDouble(j, key, true); |
312 | 0 | } |
313 | | |
314 | | // --------------------------------------------------------------------------- |
315 | | |
316 | 0 | static json getObjectMember(const json &j, const char *key) { |
317 | 0 | if (!j.contains(key)) { |
318 | 0 | throw ParsingException(std::string("Missing \"") + key + "\" key"); |
319 | 0 | } |
320 | 0 | const json obj = j[key]; |
321 | 0 | if (!obj.is_object()) { |
322 | 0 | throw ParsingException(std::string("The value of \"") + key + |
323 | 0 | "\" should be a object"); |
324 | 0 | } |
325 | 0 | return obj; |
326 | 0 | } |
327 | | |
328 | | // --------------------------------------------------------------------------- |
329 | | |
330 | 0 | static json getArrayMember(const json &j, const char *key) { |
331 | 0 | if (!j.contains(key)) { |
332 | 0 | throw ParsingException(std::string("Missing \"") + key + "\" key"); |
333 | 0 | } |
334 | 0 | const json obj = j[key]; |
335 | 0 | if (!obj.is_array()) { |
336 | 0 | throw ParsingException(std::string("The value of \"") + key + |
337 | 0 | "\" should be a array"); |
338 | 0 | } |
339 | 0 | return obj; |
340 | 0 | } |
341 | | |
342 | | // --------------------------------------------------------------------------- |
343 | | |
344 | 0 | std::unique_ptr<MasterFile> MasterFile::parse(const std::string &text) { |
345 | 0 | std::unique_ptr<MasterFile> dmmf(new MasterFile()); |
346 | 0 | json j; |
347 | 0 | try { |
348 | 0 | j = json::parse(text); |
349 | 0 | } catch (const std::exception &e) { |
350 | 0 | throw ParsingException(e.what()); |
351 | 0 | } |
352 | 0 | if (!j.is_object()) { |
353 | 0 | throw ParsingException("Not an object"); |
354 | 0 | } |
355 | 0 | dmmf->mFileType = getReqString(j, "file_type"); |
356 | 0 | dmmf->mFormatVersion = getReqString(j, "format_version"); |
357 | 0 | dmmf->mName = getOptString(j, "name"); |
358 | 0 | dmmf->mVersion = getOptString(j, "version"); |
359 | 0 | dmmf->mLicense = getOptString(j, "license"); |
360 | 0 | dmmf->mDescription = getOptString(j, "description"); |
361 | 0 | dmmf->mPublicationDate = getOptString(j, "publication_date"); |
362 | |
|
363 | 0 | if (j.contains("authority")) { |
364 | 0 | const json jAuthority = j["authority"]; |
365 | 0 | if (!jAuthority.is_object()) { |
366 | 0 | throw ParsingException("authority is not a object"); |
367 | 0 | } |
368 | 0 | dmmf->mAuthority.name = getOptString(jAuthority, "name"); |
369 | 0 | dmmf->mAuthority.url = getOptString(jAuthority, "url"); |
370 | 0 | dmmf->mAuthority.address = getOptString(jAuthority, "address"); |
371 | 0 | dmmf->mAuthority.email = getOptString(jAuthority, "email"); |
372 | 0 | } |
373 | | |
374 | 0 | if (j.contains("links")) { |
375 | 0 | const json jLinks = j["links"]; |
376 | 0 | if (!jLinks.is_array()) { |
377 | 0 | throw ParsingException("links is not an array"); |
378 | 0 | } |
379 | 0 | for (const json &jLink : jLinks) { |
380 | 0 | if (!jLink.is_object()) { |
381 | 0 | throw ParsingException("links[] item is not an object"); |
382 | 0 | } |
383 | 0 | Link link; |
384 | 0 | link.href = getOptString(jLink, "href"); |
385 | 0 | link.rel = getOptString(jLink, "rel"); |
386 | 0 | link.type = getOptString(jLink, "type"); |
387 | 0 | link.title = getOptString(jLink, "title"); |
388 | 0 | dmmf->mLinks.emplace_back(std::move(link)); |
389 | 0 | } |
390 | 0 | } |
391 | 0 | dmmf->mSourceCRS = getReqString(j, "source_crs"); |
392 | 0 | dmmf->mTargetCRS = getReqString(j, "target_crs"); |
393 | 0 | dmmf->mDefinitionCRS = getReqString(j, "definition_crs"); |
394 | 0 | if (dmmf->mSourceCRS != dmmf->mDefinitionCRS) { |
395 | 0 | throw ParsingException( |
396 | 0 | "source_crs != definition_crs not currently supported"); |
397 | 0 | } |
398 | 0 | dmmf->mReferenceEpoch = getOptString(j, "reference_epoch"); |
399 | 0 | dmmf->mUncertaintyReferenceEpoch = |
400 | 0 | getOptString(j, "uncertainty_reference_epoch"); |
401 | 0 | dmmf->mHorizontalOffsetUnit = getOptString(j, "horizontal_offset_unit"); |
402 | 0 | if (!dmmf->mHorizontalOffsetUnit.empty() && |
403 | 0 | dmmf->mHorizontalOffsetUnit != STR_METRE && |
404 | 0 | dmmf->mHorizontalOffsetUnit != STR_DEGREE) { |
405 | 0 | throw ParsingException("Unsupported value for horizontal_offset_unit"); |
406 | 0 | } |
407 | 0 | dmmf->mVerticalOffsetUnit = getOptString(j, "vertical_offset_unit"); |
408 | 0 | if (!dmmf->mVerticalOffsetUnit.empty() && |
409 | 0 | dmmf->mVerticalOffsetUnit != STR_METRE) { |
410 | 0 | throw ParsingException("Unsupported value for vertical_offset_unit"); |
411 | 0 | } |
412 | 0 | dmmf->mHorizontalUncertaintyType = |
413 | 0 | getOptString(j, "horizontal_uncertainty_type"); |
414 | 0 | dmmf->mHorizontalUncertaintyUnit = |
415 | 0 | getOptString(j, "horizontal_uncertainty_unit"); |
416 | 0 | dmmf->mVerticalUncertaintyType = |
417 | 0 | getOptString(j, "vertical_uncertainty_type"); |
418 | 0 | dmmf->mVerticalUncertaintyUnit = |
419 | 0 | getOptString(j, "vertical_uncertainty_unit"); |
420 | 0 | dmmf->mHorizontalOffsetMethod = getOptString(j, "horizontal_offset_method"); |
421 | 0 | if (!dmmf->mHorizontalOffsetMethod.empty() && |
422 | 0 | dmmf->mHorizontalOffsetMethod != STR_ADDITION && |
423 | 0 | dmmf->mHorizontalOffsetMethod != STR_GEOCENTRIC) { |
424 | 0 | throw ParsingException( |
425 | 0 | "Unsupported value for horizontal_offset_method"); |
426 | 0 | } |
427 | 0 | dmmf->mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); |
428 | |
|
429 | 0 | const json jTimeExtent = getObjectMember(j, "time_extent"); |
430 | 0 | dmmf->mTimeExtent.first = Epoch(getReqString(jTimeExtent, "first")); |
431 | 0 | dmmf->mTimeExtent.last = Epoch(getReqString(jTimeExtent, "last")); |
432 | |
|
433 | 0 | const json jComponents = getArrayMember(j, "components"); |
434 | 0 | for (const json &jComponent : jComponents) { |
435 | 0 | dmmf->mComponents.emplace_back(Component::parse(jComponent)); |
436 | 0 | const auto &comp = dmmf->mComponents.back(); |
437 | 0 | if (comp.displacementType() == STR_HORIZONTAL || |
438 | 0 | comp.displacementType() == STR_3D) { |
439 | 0 | if (dmmf->mHorizontalOffsetUnit.empty()) { |
440 | 0 | throw ParsingException("horizontal_offset_unit should be " |
441 | 0 | "defined as there is a component with " |
442 | 0 | "displacement_type = horizontal/3d"); |
443 | 0 | } |
444 | 0 | if (dmmf->mHorizontalOffsetMethod.empty()) { |
445 | 0 | throw ParsingException("horizontal_offset_method should be " |
446 | 0 | "defined as there is a component with " |
447 | 0 | "displacement_type = horizontal/3d"); |
448 | 0 | } |
449 | 0 | } |
450 | 0 | if (comp.displacementType() == STR_VERTICAL || |
451 | 0 | comp.displacementType() == STR_3D) { |
452 | 0 | if (dmmf->mVerticalOffsetUnit.empty()) { |
453 | 0 | throw ParsingException("vertical_offset_unit should be defined " |
454 | 0 | "as there is a component with " |
455 | 0 | "displacement_type = vertical/3d"); |
456 | 0 | } |
457 | 0 | } |
458 | 0 | if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && |
459 | 0 | comp.spatialModel().interpolationMethod != STR_BILINEAR) { |
460 | 0 | throw ParsingException("horizontal_offset_unit = degree can only " |
461 | 0 | "be used with interpolation_method = " |
462 | 0 | "bilinear"); |
463 | 0 | } |
464 | 0 | } |
465 | | |
466 | 0 | if (dmmf->mHorizontalOffsetUnit == STR_DEGREE && |
467 | 0 | dmmf->mHorizontalOffsetMethod != STR_ADDITION) { |
468 | 0 | throw ParsingException("horizontal_offset_unit = degree can only be " |
469 | 0 | "used with horizontal_offset_method = addition"); |
470 | 0 | } |
471 | | |
472 | 0 | return dmmf; |
473 | 0 | } |
474 | | |
475 | | // --------------------------------------------------------------------------- |
476 | | |
477 | 0 | SpatialExtent SpatialExtent::parse(const json &j) { |
478 | 0 | SpatialExtent ex; |
479 | |
|
480 | 0 | const std::string type = getReqString(j, "type"); |
481 | 0 | if (type != "bbox") { |
482 | 0 | throw ParsingException("unsupported type of extent"); |
483 | 0 | } |
484 | | |
485 | 0 | const json jParameter = getObjectMember(j, "parameters"); |
486 | 0 | const json jBbox = getArrayMember(jParameter, "bbox"); |
487 | 0 | if (jBbox.size() != 4) { |
488 | 0 | throw ParsingException("bbox is not an array of 4 numeric elements"); |
489 | 0 | } |
490 | 0 | for (int i = 0; i < 4; i++) { |
491 | 0 | if (!jBbox[i].is_number()) { |
492 | 0 | throw ParsingException( |
493 | 0 | "bbox is not an array of 4 numeric elements"); |
494 | 0 | } |
495 | 0 | } |
496 | 0 | ex.mMinx = jBbox[0].get<double>(); |
497 | 0 | ex.mMiny = jBbox[1].get<double>(); |
498 | 0 | ex.mMaxx = jBbox[2].get<double>(); |
499 | 0 | ex.mMaxy = jBbox[3].get<double>(); |
500 | |
|
501 | 0 | ex.mMinxRad = DegToRad(ex.mMinx); |
502 | 0 | ex.mMinyRad = DegToRad(ex.mMiny); |
503 | 0 | ex.mMaxxRad = DegToRad(ex.mMaxx); |
504 | 0 | ex.mMaxyRad = DegToRad(ex.mMaxy); |
505 | |
|
506 | 0 | return ex; |
507 | 0 | } |
508 | | |
509 | | // --------------------------------------------------------------------------- |
510 | | |
511 | 0 | Component Component::parse(const json &j) { |
512 | 0 | Component comp; |
513 | 0 | if (!j.is_object()) { |
514 | 0 | throw ParsingException("component is not an object"); |
515 | 0 | } |
516 | 0 | comp.mDescription = getOptString(j, "description"); |
517 | 0 | comp.mSpatialExtent = SpatialExtent::parse(getObjectMember(j, "extent")); |
518 | 0 | comp.mDisplacementType = getReqString(j, "displacement_type"); |
519 | 0 | if (comp.mDisplacementType != STR_NONE && |
520 | 0 | comp.mDisplacementType != STR_HORIZONTAL && |
521 | 0 | comp.mDisplacementType != STR_VERTICAL && |
522 | 0 | comp.mDisplacementType != STR_3D) { |
523 | 0 | throw ParsingException("Unsupported value for displacement_type"); |
524 | 0 | } |
525 | 0 | comp.mUncertaintyType = getReqString(j, "uncertainty_type"); |
526 | 0 | comp.mHorizontalUncertainty = getOptDouble(j, "horizontal_uncertainty"); |
527 | 0 | comp.mVerticalUncertainty = getOptDouble(j, "vertical_uncertainty"); |
528 | |
|
529 | 0 | const json jSpatialModel = getObjectMember(j, "spatial_model"); |
530 | 0 | comp.mSpatialModel.type = getReqString(jSpatialModel, "type"); |
531 | 0 | comp.mSpatialModel.interpolationMethod = |
532 | 0 | getReqString(jSpatialModel, "interpolation_method"); |
533 | 0 | if (comp.mSpatialModel.interpolationMethod != STR_BILINEAR && |
534 | 0 | comp.mSpatialModel.interpolationMethod != STR_GEOCENTRIC_BILINEAR) { |
535 | 0 | throw ParsingException("Unsupported value for interpolation_method"); |
536 | 0 | } |
537 | 0 | comp.mSpatialModel.filename = getReqString(jSpatialModel, "filename"); |
538 | 0 | comp.mSpatialModel.md5Checksum = |
539 | 0 | getOptString(jSpatialModel, "md5_checksum"); |
540 | |
|
541 | 0 | const json jTimeFunction = getObjectMember(j, "time_function"); |
542 | 0 | const std::string timeFunctionType = getReqString(jTimeFunction, "type"); |
543 | 0 | const json jParameters = timeFunctionType == "constant" |
544 | 0 | ? json() |
545 | 0 | : getObjectMember(jTimeFunction, "parameters"); |
546 | |
|
547 | 0 | if (timeFunctionType == "constant") { |
548 | 0 | std::unique_ptr<ConstantTimeFunction> tf(new ConstantTimeFunction()); |
549 | 0 | tf->type = timeFunctionType; |
550 | 0 | comp.mTimeFunction = std::move(tf); |
551 | 0 | } else if (timeFunctionType == "velocity") { |
552 | 0 | std::unique_ptr<VelocityTimeFunction> tf(new VelocityTimeFunction()); |
553 | 0 | tf->type = timeFunctionType; |
554 | 0 | tf->referenceEpoch = |
555 | 0 | Epoch(getReqString(jParameters, "reference_epoch")); |
556 | 0 | comp.mTimeFunction = std::move(tf); |
557 | 0 | } else if (timeFunctionType == "step") { |
558 | 0 | std::unique_ptr<StepTimeFunction> tf(new StepTimeFunction()); |
559 | 0 | tf->type = timeFunctionType; |
560 | 0 | tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); |
561 | 0 | comp.mTimeFunction = std::move(tf); |
562 | 0 | } else if (timeFunctionType == "reverse_step") { |
563 | 0 | std::unique_ptr<ReverseStepTimeFunction> tf( |
564 | 0 | new ReverseStepTimeFunction()); |
565 | 0 | tf->type = timeFunctionType; |
566 | 0 | tf->stepEpoch = Epoch(getReqString(jParameters, "step_epoch")); |
567 | 0 | comp.mTimeFunction = std::move(tf); |
568 | 0 | } else if (timeFunctionType == "piecewise") { |
569 | 0 | std::unique_ptr<PiecewiseTimeFunction> tf(new PiecewiseTimeFunction()); |
570 | 0 | tf->type = timeFunctionType; |
571 | 0 | tf->beforeFirst = getReqString(jParameters, "before_first"); |
572 | 0 | if (tf->beforeFirst != "zero" && tf->beforeFirst != "constant" && |
573 | 0 | tf->beforeFirst != "linear") { |
574 | 0 | throw ParsingException("Unsupported value for before_first"); |
575 | 0 | } |
576 | 0 | tf->afterLast = getReqString(jParameters, "after_last"); |
577 | 0 | if (tf->afterLast != "zero" && tf->afterLast != "constant" && |
578 | 0 | tf->afterLast != "linear") { |
579 | 0 | throw ParsingException("Unsupported value for afterLast"); |
580 | 0 | } |
581 | 0 | const json jModel = getArrayMember(jParameters, "model"); |
582 | 0 | for (const json &jModelElt : jModel) { |
583 | 0 | if (!jModelElt.is_object()) { |
584 | 0 | throw ParsingException("model[] element is not an object"); |
585 | 0 | } |
586 | 0 | PiecewiseTimeFunction::EpochScaleFactorTuple tuple; |
587 | 0 | tuple.epoch = Epoch(getReqString(jModelElt, "epoch")); |
588 | 0 | tuple.scaleFactor = getReqDouble(jModelElt, "scale_factor"); |
589 | 0 | tf->model.emplace_back(std::move(tuple)); |
590 | 0 | } |
591 | | |
592 | 0 | comp.mTimeFunction = std::move(tf); |
593 | 0 | } else if (timeFunctionType == "exponential") { |
594 | 0 | std::unique_ptr<ExponentialTimeFunction> tf( |
595 | 0 | new ExponentialTimeFunction()); |
596 | 0 | tf->type = timeFunctionType; |
597 | 0 | tf->referenceEpoch = |
598 | 0 | Epoch(getReqString(jParameters, "reference_epoch")); |
599 | 0 | tf->endEpoch = Epoch(getOptString(jParameters, "end_epoch")); |
600 | 0 | tf->relaxationConstant = |
601 | 0 | getReqDouble(jParameters, "relaxation_constant"); |
602 | 0 | if (tf->relaxationConstant <= 0.0) { |
603 | 0 | throw ParsingException("Invalid value for relaxation_constant"); |
604 | 0 | } |
605 | 0 | tf->beforeScaleFactor = |
606 | 0 | getReqDouble(jParameters, "before_scale_factor"); |
607 | 0 | tf->initialScaleFactor = |
608 | 0 | getReqDouble(jParameters, "initial_scale_factor"); |
609 | 0 | tf->finalScaleFactor = getReqDouble(jParameters, "final_scale_factor"); |
610 | 0 | comp.mTimeFunction = std::move(tf); |
611 | 0 | } else { |
612 | 0 | throw ParsingException("Unsupported type of time function: " + |
613 | 0 | timeFunctionType); |
614 | 0 | } |
615 | | |
616 | 0 | return comp; |
617 | 0 | } |
618 | | |
619 | | // --------------------------------------------------------------------------- |
620 | | |
621 | 0 | double Component::ConstantTimeFunction::evaluateAt(double) const { return 1.0; } |
622 | | |
623 | | // --------------------------------------------------------------------------- |
624 | | |
625 | 0 | double Component::VelocityTimeFunction::evaluateAt(double dt) const { |
626 | 0 | return dt - referenceEpoch.toDecimalYear(); |
627 | 0 | } |
628 | | |
629 | | // --------------------------------------------------------------------------- |
630 | | |
631 | 0 | double Component::StepTimeFunction::evaluateAt(double dt) const { |
632 | 0 | if (dt < stepEpoch.toDecimalYear()) |
633 | 0 | return 0.0; |
634 | 0 | return 1.0; |
635 | 0 | } |
636 | | |
637 | | // --------------------------------------------------------------------------- |
638 | | |
639 | 0 | double Component::ReverseStepTimeFunction::evaluateAt(double dt) const { |
640 | 0 | if (dt < stepEpoch.toDecimalYear()) |
641 | 0 | return -1.0; |
642 | 0 | return 0.0; |
643 | 0 | } |
644 | | |
645 | | // --------------------------------------------------------------------------- |
646 | | |
647 | 0 | double Component::PiecewiseTimeFunction::evaluateAt(double dt) const { |
648 | 0 | if (model.empty()) { |
649 | 0 | return 0.0; |
650 | 0 | } |
651 | | |
652 | 0 | const double dt1 = model[0].epoch.toDecimalYear(); |
653 | 0 | if (dt < dt1) { |
654 | 0 | if (beforeFirst == "zero") |
655 | 0 | return 0.0; |
656 | 0 | if (beforeFirst == "constant" || model.size() == 1) |
657 | 0 | return model[0].scaleFactor; |
658 | | |
659 | | // linear |
660 | 0 | const double f1 = model[0].scaleFactor; |
661 | 0 | const double dt2 = model[1].epoch.toDecimalYear(); |
662 | 0 | const double f2 = model[1].scaleFactor; |
663 | 0 | if (dt1 == dt2) |
664 | 0 | return f1; |
665 | 0 | return (f1 * (dt2 - dt) + f2 * (dt - dt1)) / (dt2 - dt1); |
666 | 0 | } |
667 | 0 | for (size_t i = 1; i < model.size(); i++) { |
668 | 0 | const double dtip1 = model[i].epoch.toDecimalYear(); |
669 | 0 | if (dt < dtip1) { |
670 | 0 | const double dti = model[i - 1].epoch.toDecimalYear(); |
671 | 0 | const double fip1 = model[i].scaleFactor; |
672 | 0 | const double fi = model[i - 1].scaleFactor; |
673 | 0 | return (fi * (dtip1 - dt) + fip1 * (dt - dti)) / (dtip1 - dti); |
674 | 0 | } |
675 | 0 | } |
676 | 0 | if (afterLast == "zero") { |
677 | 0 | return 0.0; |
678 | 0 | } |
679 | 0 | if (afterLast == "constant" || model.size() == 1) |
680 | 0 | return model.back().scaleFactor; |
681 | | |
682 | | // linear |
683 | 0 | const double dtnm1 = model[model.size() - 2].epoch.toDecimalYear(); |
684 | 0 | const double fnm1 = model[model.size() - 2].scaleFactor; |
685 | 0 | const double dtn = model.back().epoch.toDecimalYear(); |
686 | 0 | const double fn = model.back().scaleFactor; |
687 | 0 | if (dtnm1 == dtn) |
688 | 0 | return fn; |
689 | 0 | return (fnm1 * (dtn - dt) + fn * (dt - dtnm1)) / (dtn - dtnm1); |
690 | 0 | } |
691 | | |
692 | | // --------------------------------------------------------------------------- |
693 | | |
694 | 0 | double Component::ExponentialTimeFunction::evaluateAt(double dt) const { |
695 | 0 | const double t0 = referenceEpoch.toDecimalYear(); |
696 | 0 | if (dt < t0) |
697 | 0 | return beforeScaleFactor; |
698 | 0 | if (!endEpoch.toString().empty()) { |
699 | 0 | dt = std::min(dt, endEpoch.toDecimalYear()); |
700 | 0 | } |
701 | 0 | return initialScaleFactor + |
702 | 0 | (finalScaleFactor - initialScaleFactor) * |
703 | 0 | (1.0 - std::exp(-(dt - t0) / relaxationConstant)); |
704 | 0 | } |
705 | | |
706 | | // --------------------------------------------------------------------------- |
707 | | |
708 | | inline void DeltaEastingNorthingToLongLat(double cosphi, double de, double dn, |
709 | | double a, double b, double es, |
710 | 0 | double &dlam, double &dphi) { |
711 | 0 | const double oneMinuX = es * (1 - cosphi * cosphi); |
712 | 0 | const double X = 1 - oneMinuX; |
713 | 0 | const double sqrtX = sqrt(X); |
714 | | #if 0 |
715 | | // With es of Earth, absolute/relative error is at most 2e-8 |
716 | | // const double sqrtX = 1 - 0.5 * oneMinuX * (1 + 0.25 * oneMinuX); |
717 | | #endif |
718 | 0 | dlam = de * sqrtX / (a * cosphi); |
719 | 0 | dphi = dn * a * sqrtX * X / (b * b); |
720 | 0 | } |
721 | | |
722 | | // --------------------------------------------------------------------------- |
723 | | |
724 | | template <class Grid, class GridSet, class EvaluatorIface> |
725 | | Evaluator<Grid, GridSet, EvaluatorIface>::Evaluator( |
726 | | std::unique_ptr<MasterFile> &&model, EvaluatorIface &iface, double a, |
727 | | double b) |
728 | | : mModel(std::move(model)), mA(a), mB(b), mEs(1 - (b * b) / (a * a)), |
729 | | mIsHorizontalUnitDegree(mModel->horizontalOffsetUnit() == STR_DEGREE), |
730 | | mIsAddition(mModel->horizontalOffsetMethod() == STR_ADDITION), |
731 | 0 | mIsGeographicCRS(iface.isGeographicCRS(mModel->definitionCRS())) { |
732 | 0 | if (!mIsGeographicCRS && mIsHorizontalUnitDegree) { |
733 | 0 | throw EvaluatorException( |
734 | 0 | "definition_crs = projected CRS and " |
735 | 0 | "horizontal_offset_unit = degree are incompatible"); |
736 | 0 | } |
737 | 0 | if (!mIsGeographicCRS && !mIsAddition) { |
738 | 0 | throw EvaluatorException( |
739 | 0 | "definition_crs = projected CRS and " |
740 | 0 | "horizontal_offset_method = geocentric are incompatible"); |
741 | 0 | } |
742 | 0 | mComponents.reserve(mModel->components().size()); |
743 | 0 | for (const auto &comp : mModel->components()) { |
744 | 0 | mComponents.emplace_back(std::unique_ptr<ComponentEx<Grid, GridSet>>( |
745 | 0 | new ComponentEx<Grid, GridSet>(comp))); |
746 | 0 | if (!mIsGeographicCRS && !mComponents.back()->isBilinearInterpolation) { |
747 | 0 | throw EvaluatorException( |
748 | 0 | "definition_crs = projected CRS and " |
749 | 0 | "interpolation_method = geocentric_bilinear are incompatible"); |
750 | 0 | } |
751 | 0 | } |
752 | 0 | } |
753 | | |
754 | | // --------------------------------------------------------------------------- |
755 | | |
756 | | template <class Grid, class GridSet, class EvaluatorIface> |
757 | 0 | void Evaluator<Grid, GridSet, EvaluatorIface>::clearGridCache() { |
758 | 0 | for (auto &comp : mComponents) { |
759 | 0 | comp->clearGridCache(); |
760 | 0 | } |
761 | 0 | } |
762 | | |
763 | | // --------------------------------------------------------------------------- |
764 | | |
765 | | #ifdef DEBUG_DEFMODEL |
766 | | |
767 | | static std::string shortName(const Component &comp) { |
768 | | const auto &desc = comp.description(); |
769 | | return desc.substr(0, desc.find('\n')) + " (" + |
770 | | comp.spatialModel().filename + ")"; |
771 | | } |
772 | | |
773 | | static std::string toString(double val) { |
774 | | char buffer[32]; |
775 | | snprintf(buffer, sizeof(buffer), "%.9g", val); |
776 | | return buffer; |
777 | | } |
778 | | |
779 | | #endif |
780 | | |
781 | | // --------------------------------------------------------------------------- |
782 | | |
783 | | static bool bboxCheck(double &x, double &y, bool forInverseComputation, |
784 | | const double minx, const double miny, const double maxx, |
785 | | const double maxy, const double EPS, |
786 | 0 | const double extraMarginForInverse) { |
787 | 0 | if (x < minx - EPS || x > maxx + EPS || y < miny - EPS || y > maxy + EPS) { |
788 | 0 | if (!forInverseComputation) { |
789 | 0 | return false; |
790 | 0 | } |
791 | | // In case of iterative computation for inverse, allow to be a |
792 | | // slightly bit outside of the grid and clamp to the edges |
793 | 0 | bool xOk = false; |
794 | 0 | if (x >= minx - EPS && x <= maxx + EPS) { |
795 | 0 | xOk = true; |
796 | 0 | } else if (x > minx - extraMarginForInverse && x < minx) { |
797 | 0 | x = minx; |
798 | 0 | xOk = true; |
799 | 0 | } else if (x < maxx + extraMarginForInverse && x > maxx) { |
800 | 0 | x = maxx; |
801 | 0 | xOk = true; |
802 | 0 | } |
803 | |
|
804 | 0 | bool yOk = false; |
805 | 0 | if (y >= miny - EPS && y <= maxy + EPS) { |
806 | 0 | yOk = true; |
807 | 0 | } else if (y > miny - extraMarginForInverse && y < miny) { |
808 | 0 | y = miny; |
809 | 0 | yOk = true; |
810 | 0 | } else if (y < maxy + extraMarginForInverse && y > maxy) { |
811 | 0 | y = maxy; |
812 | 0 | yOk = true; |
813 | 0 | } |
814 | |
|
815 | 0 | return xOk && yOk; |
816 | 0 | } |
817 | 0 | return true; |
818 | 0 | } |
819 | | |
820 | | // --------------------------------------------------------------------------- |
821 | | |
822 | | template <class Grid, class GridSet, class EvaluatorIface> |
823 | | bool Evaluator<Grid, GridSet, EvaluatorIface>::forward( |
824 | | EvaluatorIface &iface, double x, double y, double z, double t, |
825 | | bool forInverseComputation, double &x_out, double &y_out, double &z_out) |
826 | | |
827 | 0 | { |
828 | 0 | x_out = x; |
829 | 0 | y_out = y; |
830 | 0 | z_out = z; |
831 | |
|
832 | 0 | const double EPS = mIsGeographicCRS ? 1e-10 : 1e-5; |
833 | | |
834 | | // Check against global model spatial extent, potentially wrapping |
835 | | // longitude to match |
836 | 0 | { |
837 | 0 | const auto &extent = mModel->extent(); |
838 | 0 | const double minx = extent.minxNormalized(mIsGeographicCRS); |
839 | 0 | const double maxx = extent.maxxNormalized(mIsGeographicCRS); |
840 | 0 | if (mIsGeographicCRS) { |
841 | 0 | while (x < minx - EPS) { |
842 | 0 | x += 2.0 * DEFMODEL_PI; |
843 | 0 | } |
844 | 0 | while (x > maxx + EPS) { |
845 | 0 | x -= 2.0 * DEFMODEL_PI; |
846 | 0 | } |
847 | 0 | } |
848 | 0 | const double miny = extent.minyNormalized(mIsGeographicCRS); |
849 | 0 | const double maxy = extent.maxyNormalized(mIsGeographicCRS); |
850 | 0 | const double extraMarginForInverse = |
851 | 0 | mIsGeographicCRS ? DegToRad(0.1) : 10000; |
852 | 0 | if (!bboxCheck(x, y, forInverseComputation, minx, miny, maxx, maxy, EPS, |
853 | 0 | extraMarginForInverse)) { |
854 | | #ifdef DEBUG_DEFMODEL |
855 | | iface.log("Calculation point " + toString(x) + "," + toString(y) + |
856 | | " is outside the extents of the deformation model"); |
857 | | #endif |
858 | 0 | return false; |
859 | 0 | } |
860 | 0 | } |
861 | | |
862 | | // Check against global model temporal extent |
863 | 0 | { |
864 | 0 | const auto &timeExtent = mModel->timeExtent(); |
865 | 0 | if (t < timeExtent.first.toDecimalYear() || |
866 | 0 | t > timeExtent.last.toDecimalYear()) { |
867 | | #ifdef DEBUG_DEFMODEL |
868 | | iface.log("Calculation epoch " + toString(t) + |
869 | | " is not valid for the deformation model"); |
870 | | #endif |
871 | 0 | return false; |
872 | 0 | } |
873 | 0 | } |
874 | | |
875 | | // For mIsHorizontalUnitDegree |
876 | 0 | double dlam = 0; |
877 | 0 | double dphi = 0; |
878 | | |
879 | | // For !mIsHorizontalUnitDegree |
880 | 0 | double de = 0; |
881 | 0 | double dn = 0; |
882 | |
|
883 | 0 | double dz = 0; |
884 | |
|
885 | 0 | bool sincosphiInitialized = false; |
886 | 0 | double sinphi = 0; |
887 | 0 | double cosphi = 0; |
888 | |
|
889 | 0 | for (auto &compEx : mComponents) { |
890 | 0 | const auto &comp = compEx->component; |
891 | 0 | if (compEx->displacementType == DisplacementType::NONE) { |
892 | 0 | continue; |
893 | 0 | } |
894 | 0 | const auto &extent = comp.extent(); |
895 | 0 | double xForGrid = x; |
896 | 0 | double yForGrid = y; |
897 | 0 | const double minx = extent.minxNormalized(mIsGeographicCRS); |
898 | 0 | const double maxx = extent.maxxNormalized(mIsGeographicCRS); |
899 | 0 | const double miny = extent.minyNormalized(mIsGeographicCRS); |
900 | 0 | const double maxy = extent.maxyNormalized(mIsGeographicCRS); |
901 | 0 | const double extraMarginForInverse = 0; |
902 | 0 | if (!bboxCheck(xForGrid, yForGrid, forInverseComputation, minx, miny, |
903 | 0 | maxx, maxy, EPS, extraMarginForInverse)) { |
904 | | #ifdef DEBUG_DEFMODEL |
905 | | iface.log( |
906 | | "Skipping component " + shortName(comp) + |
907 | | " due to point being outside of its declared spatial extent."); |
908 | | #endif |
909 | 0 | continue; |
910 | 0 | } |
911 | 0 | xForGrid = std::max(xForGrid, minx); |
912 | 0 | yForGrid = std::max(yForGrid, miny); |
913 | 0 | xForGrid = std::min(xForGrid, maxx); |
914 | 0 | yForGrid = std::min(yForGrid, maxy); |
915 | 0 | const auto tfactor = compEx->evaluateAt(t); |
916 | 0 | if (tfactor == 0.0) { |
917 | | #ifdef DEBUG_DEFMODEL |
918 | | iface.log("Skipping component " + shortName(comp) + |
919 | | " due to time function evaluating to 0."); |
920 | | #endif |
921 | 0 | continue; |
922 | 0 | } |
923 | | |
924 | | #ifdef DEBUG_DEFMODEL |
925 | | iface.log("Entering component " + shortName(comp) + |
926 | | " with time function evaluating to " + toString(tfactor) + |
927 | | "."); |
928 | | #endif |
929 | | |
930 | 0 | if (compEx->gridSet == nullptr) { |
931 | 0 | compEx->gridSet = iface.open(comp.spatialModel().filename); |
932 | 0 | if (compEx->gridSet == nullptr) { |
933 | 0 | return false; |
934 | 0 | } |
935 | 0 | } |
936 | 0 | const Grid *grid = compEx->gridSet->gridAt(xForGrid, yForGrid); |
937 | 0 | if (grid == nullptr) { |
938 | | #ifdef DEBUG_DEFMODEL |
939 | | iface.log("Skipping component " + shortName(comp) + |
940 | | " due to no grid found for this point in the grid set."); |
941 | | #endif |
942 | 0 | continue; |
943 | 0 | } |
944 | 0 | if (grid->width < 2 || grid->height < 2) { |
945 | 0 | return false; |
946 | 0 | } |
947 | 0 | const double ix_d = (xForGrid - grid->minx) / grid->resx; |
948 | 0 | const double iy_d = (yForGrid - grid->miny) / grid->resy; |
949 | 0 | if (ix_d < -EPS || iy_d < -EPS || ix_d + 1 >= grid->width + EPS || |
950 | 0 | iy_d + 1 >= grid->height + EPS) { |
951 | | #ifdef DEBUG_DEFMODEL |
952 | | iface.log("Skipping component " + shortName(comp) + |
953 | | " due to point being outside of actual spatial extent of " |
954 | | "grid " + |
955 | | grid->name() + "."); |
956 | | #endif |
957 | 0 | continue; |
958 | 0 | } |
959 | 0 | const int ix0 = std::min(static_cast<int>(ix_d), grid->width - 2); |
960 | 0 | const int iy0 = std::min(static_cast<int>(iy_d), grid->height - 2); |
961 | 0 | const int ix1 = ix0 + 1; |
962 | 0 | const int iy1 = iy0 + 1; |
963 | 0 | const double frct_x = ix_d - ix0; |
964 | 0 | const double frct_y = iy_d - iy0; |
965 | 0 | const double one_minus_frct_x = 1. - frct_x; |
966 | 0 | const double one_minus_frct_y = 1. - frct_y; |
967 | 0 | const double m00 = one_minus_frct_x * one_minus_frct_y; |
968 | 0 | const double m10 = frct_x * one_minus_frct_y; |
969 | 0 | const double m01 = one_minus_frct_x * frct_y; |
970 | 0 | const double m11 = frct_x * frct_y; |
971 | |
|
972 | 0 | if (compEx->displacementType == DisplacementType::VERTICAL) { |
973 | 0 | double dz00 = 0; |
974 | 0 | double dz01 = 0; |
975 | 0 | double dz10 = 0; |
976 | 0 | double dz11 = 0; |
977 | 0 | if (!grid->getZOffset(ix0, iy0, dz00) || |
978 | 0 | !grid->getZOffset(ix1, iy0, dz10) || |
979 | 0 | !grid->getZOffset(ix0, iy1, dz01) || |
980 | 0 | !grid->getZOffset(ix1, iy1, dz11)) { |
981 | 0 | return false; |
982 | 0 | } |
983 | 0 | const double dzInterp = |
984 | 0 | dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; |
985 | | #ifdef DEBUG_DEFMODEL |
986 | | iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + |
987 | | toString(dzInterp) + "."); |
988 | | #endif |
989 | 0 | dz += tfactor * dzInterp; |
990 | 0 | } else if (mIsHorizontalUnitDegree) { |
991 | 0 | double dx00 = 0; |
992 | 0 | double dy00 = 0; |
993 | 0 | double dx01 = 0; |
994 | 0 | double dy01 = 0; |
995 | 0 | double dx10 = 0; |
996 | 0 | double dy10 = 0; |
997 | 0 | double dx11 = 0; |
998 | 0 | double dy11 = 0; |
999 | 0 | if (compEx->displacementType == DisplacementType::HORIZONTAL) { |
1000 | 0 | if (!grid->getLongLatOffset(ix0, iy0, dx00, dy00) || |
1001 | 0 | !grid->getLongLatOffset(ix1, iy0, dx10, dy10) || |
1002 | 0 | !grid->getLongLatOffset(ix0, iy1, dx01, dy01) || |
1003 | 0 | !grid->getLongLatOffset(ix1, iy1, dx11, dy11)) { |
1004 | 0 | return false; |
1005 | 0 | } |
1006 | 0 | } else /* if (compEx->displacementType == DisplacementType::THREE_D) |
1007 | | */ |
1008 | 0 | { |
1009 | 0 | double dz00 = 0; |
1010 | 0 | double dz01 = 0; |
1011 | 0 | double dz10 = 0; |
1012 | 0 | double dz11 = 0; |
1013 | 0 | if (!grid->getLongLatZOffset(ix0, iy0, dx00, dy00, dz00) || |
1014 | 0 | !grid->getLongLatZOffset(ix1, iy0, dx10, dy10, dz10) || |
1015 | 0 | !grid->getLongLatZOffset(ix0, iy1, dx01, dy01, dz01) || |
1016 | 0 | !grid->getLongLatZOffset(ix1, iy1, dx11, dy11, dz11)) { |
1017 | 0 | return false; |
1018 | 0 | } |
1019 | 0 | const double dzInterp = |
1020 | 0 | dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; |
1021 | | #ifdef DEBUG_DEFMODEL |
1022 | | iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + |
1023 | | toString(dzInterp) + "."); |
1024 | | #endif |
1025 | 0 | dz += tfactor * dzInterp; |
1026 | 0 | } |
1027 | 0 | const double dlamInterp = |
1028 | 0 | dx00 * m00 + dx01 * m01 + dx10 * m10 + dx11 * m11; |
1029 | 0 | const double dphiInterp = |
1030 | 0 | dy00 * m00 + dy01 * m01 + dy10 * m10 + dy11 * m11; |
1031 | | #ifdef DEBUG_DEFMODEL |
1032 | | iface.log("tfactor * dlamInterp = " + toString(tfactor) + " * " + |
1033 | | toString(dlamInterp) + "."); |
1034 | | iface.log("tfactor * dphiInterp = " + toString(tfactor) + " * " + |
1035 | | toString(dphiInterp) + "."); |
1036 | | #endif |
1037 | 0 | dlam += tfactor * dlamInterp; |
1038 | 0 | dphi += tfactor * dphiInterp; |
1039 | 0 | } else /* horizontal unit is metre */ { |
1040 | 0 | double de00 = 0; |
1041 | 0 | double dn00 = 0; |
1042 | 0 | double de01 = 0; |
1043 | 0 | double dn01 = 0; |
1044 | 0 | double de10 = 0; |
1045 | 0 | double dn10 = 0; |
1046 | 0 | double de11 = 0; |
1047 | 0 | double dn11 = 0; |
1048 | 0 | if (compEx->displacementType == DisplacementType::HORIZONTAL) { |
1049 | 0 | if (!grid->getEastingNorthingOffset(ix0, iy0, de00, dn00) || |
1050 | 0 | !grid->getEastingNorthingOffset(ix1, iy0, de10, dn10) || |
1051 | 0 | !grid->getEastingNorthingOffset(ix0, iy1, de01, dn01) || |
1052 | 0 | !grid->getEastingNorthingOffset(ix1, iy1, de11, dn11)) { |
1053 | 0 | return false; |
1054 | 0 | } |
1055 | 0 | } else /* if (compEx->displacementType == DisplacementType::THREE_D) |
1056 | | */ |
1057 | 0 | { |
1058 | 0 | double dz00 = 0; |
1059 | 0 | double dz01 = 0; |
1060 | 0 | double dz10 = 0; |
1061 | 0 | double dz11 = 0; |
1062 | 0 | if (!grid->getEastingNorthingZOffset(ix0, iy0, de00, dn00, |
1063 | 0 | dz00) || |
1064 | 0 | !grid->getEastingNorthingZOffset(ix1, iy0, de10, dn10, |
1065 | 0 | dz10) || |
1066 | 0 | !grid->getEastingNorthingZOffset(ix0, iy1, de01, dn01, |
1067 | 0 | dz01) || |
1068 | 0 | !grid->getEastingNorthingZOffset(ix1, iy1, de11, dn11, |
1069 | 0 | dz11)) { |
1070 | 0 | return false; |
1071 | 0 | } |
1072 | 0 | const double dzInterp = |
1073 | 0 | dz00 * m00 + dz01 * m01 + dz10 * m10 + dz11 * m11; |
1074 | | #ifdef DEBUG_DEFMODEL |
1075 | | iface.log("tfactor * dzInterp = " + toString(tfactor) + " * " + |
1076 | | toString(dzInterp) + "."); |
1077 | | #endif |
1078 | 0 | dz += tfactor * dzInterp; |
1079 | 0 | } |
1080 | 0 | if (compEx->isBilinearInterpolation) { |
1081 | 0 | const double deInterp = |
1082 | 0 | de00 * m00 + de01 * m01 + de10 * m10 + de11 * m11; |
1083 | 0 | const double dnInterp = |
1084 | 0 | dn00 * m00 + dn01 * m01 + dn10 * m10 + dn11 * m11; |
1085 | | #ifdef DEBUG_DEFMODEL |
1086 | | iface.log("tfactor * deInterp = " + toString(tfactor) + " * " + |
1087 | | toString(deInterp) + "."); |
1088 | | iface.log("tfactor * dnInterp = " + toString(tfactor) + " * " + |
1089 | | toString(dnInterp) + "."); |
1090 | | #endif |
1091 | 0 | de += tfactor * deInterp; |
1092 | 0 | dn += tfactor * dnInterp; |
1093 | 0 | } else /* geocentric_bilinear */ { |
1094 | 0 | double dX; |
1095 | 0 | double dY; |
1096 | 0 | double dZ; |
1097 | |
|
1098 | 0 | auto iter = compEx->mapGrids.find(grid); |
1099 | 0 | if (iter == compEx->mapGrids.end()) { |
1100 | 0 | GridEx<Grid> gridWithCache(grid); |
1101 | 0 | iter = compEx->mapGrids |
1102 | 0 | .insert(std::pair<const Grid *, GridEx<Grid>>( |
1103 | 0 | grid, std::move(gridWithCache))) |
1104 | 0 | .first; |
1105 | 0 | } |
1106 | 0 | GridEx<Grid> &gridwithCacheRef = iter->second; |
1107 | |
|
1108 | 0 | gridwithCacheRef.getBilinearGeocentric( |
1109 | 0 | ix0, iy0, de00, dn00, de01, dn01, de10, dn10, de11, dn11, |
1110 | 0 | m00, m01, m10, m11, dX, dY, dZ); |
1111 | 0 | if (!sincosphiInitialized) { |
1112 | 0 | sincosphiInitialized = true; |
1113 | 0 | sinphi = sin(y); |
1114 | 0 | cosphi = cos(y); |
1115 | 0 | } |
1116 | 0 | const double lam_rel_to_cell_center = |
1117 | 0 | (frct_x - 0.5) * grid->resx; |
1118 | | // Use small-angle approximation of sin/cos when reasonable |
1119 | | // Max abs/rel error on cos is 3.9e-9 and on sin 1.3e-11 |
1120 | 0 | const double sinlam = |
1121 | 0 | gridwithCacheRef.smallResx |
1122 | 0 | ? lam_rel_to_cell_center * |
1123 | 0 | (1 - (1. / 6) * (lam_rel_to_cell_center * |
1124 | 0 | lam_rel_to_cell_center)) |
1125 | 0 | : sin(lam_rel_to_cell_center); |
1126 | 0 | const double coslam = gridwithCacheRef.smallResx |
1127 | 0 | ? (1 - 0.5 * (lam_rel_to_cell_center * |
1128 | 0 | lam_rel_to_cell_center)) |
1129 | 0 | : cos(lam_rel_to_cell_center); |
1130 | | |
1131 | | // Convert back from geocentric deltas to easting, northing |
1132 | | // deltas |
1133 | 0 | const double deInterp = -dX * sinlam + dY * coslam; |
1134 | 0 | const double dnInterp = |
1135 | 0 | (-dX * coslam - dY * sinlam) * sinphi + dZ * cosphi; |
1136 | | #ifdef DEBUG_DEFMODEL |
1137 | | iface.log("After geocentric_bilinear interpolation: tfactor * " |
1138 | | "deInterp = " + |
1139 | | toString(tfactor) + " * " + toString(deInterp) + "."); |
1140 | | iface.log("After geocentric_bilinear interpolation: tfactor * " |
1141 | | "dnInterp = " + |
1142 | | toString(tfactor) + " * " + toString(dnInterp) + "."); |
1143 | | #endif |
1144 | 0 | de += tfactor * deInterp; |
1145 | 0 | dn += tfactor * dnInterp; |
1146 | 0 | } |
1147 | 0 | } |
1148 | 0 | } |
1149 | | |
1150 | | // Apply shifts depending on horizontal_offset_unit and |
1151 | | // horizontal_offset_method |
1152 | 0 | if (mIsHorizontalUnitDegree) { |
1153 | 0 | x_out += dlam; |
1154 | 0 | y_out += dphi; |
1155 | 0 | } else { |
1156 | | #ifdef DEBUG_DEFMODEL |
1157 | | iface.log("Total sum of de: " + toString(de)); |
1158 | | iface.log("Total sum of dn: " + toString(dn)); |
1159 | | #endif |
1160 | 0 | if (mIsAddition && !mIsGeographicCRS) { |
1161 | 0 | x_out += de; |
1162 | 0 | y_out += dn; |
1163 | 0 | } else if (mIsAddition) { |
1164 | | // Simple way of adding the offset |
1165 | 0 | if (!sincosphiInitialized) { |
1166 | 0 | cosphi = cos(y); |
1167 | 0 | } |
1168 | 0 | DeltaEastingNorthingToLongLat(cosphi, de, dn, mA, mB, mEs, dlam, |
1169 | 0 | dphi); |
1170 | | #ifdef DEBUG_DEFMODEL |
1171 | | iface.log("Result dlam: " + toString(dlam)); |
1172 | | iface.log("Result dphi: " + toString(dphi)); |
1173 | | #endif |
1174 | 0 | x_out += dlam; |
1175 | 0 | y_out += dphi; |
1176 | 0 | } else { |
1177 | | // Geocentric way of adding the offset |
1178 | 0 | if (!sincosphiInitialized) { |
1179 | 0 | sinphi = sin(y); |
1180 | 0 | cosphi = cos(y); |
1181 | 0 | } |
1182 | 0 | const double sinlam = sin(x); |
1183 | 0 | const double coslam = cos(x); |
1184 | 0 | const double dnsinphi = dn * sinphi; |
1185 | 0 | const double dX = -de * sinlam - dnsinphi * coslam; |
1186 | 0 | const double dY = de * coslam - dnsinphi * sinlam; |
1187 | 0 | const double dZ = dn * cosphi; |
1188 | |
|
1189 | 0 | double X; |
1190 | 0 | double Y; |
1191 | 0 | double Z; |
1192 | 0 | iface.geographicToGeocentric(x, y, 0, mA, mB, mEs, X, Y, Z); |
1193 | | #ifdef DEBUG_DEFMODEL |
1194 | | iface.log("Geocentric coordinate before: " + toString(X) + "," + |
1195 | | toString(Y) + "," + toString(Z)); |
1196 | | iface.log("Geocentric shift: " + toString(dX) + "," + toString(dY) + |
1197 | | "," + toString(dZ)); |
1198 | | #endif |
1199 | 0 | X += dX; |
1200 | 0 | Y += dY; |
1201 | 0 | Z += dZ; |
1202 | | #ifdef DEBUG_DEFMODEL |
1203 | | iface.log("Geocentric coordinate after: " + toString(X) + "," + |
1204 | | toString(Y) + "," + toString(Z)); |
1205 | | #endif |
1206 | |
|
1207 | 0 | double h_out_ignored; |
1208 | 0 | iface.geocentricToGeographic(X, Y, Z, mA, mB, mEs, x_out, y_out, |
1209 | 0 | h_out_ignored); |
1210 | 0 | } |
1211 | 0 | } |
1212 | | #ifdef DEBUG_DEFMODEL |
1213 | | iface.log("Total sum of dz: " + toString(dz)); |
1214 | | #endif |
1215 | 0 | z_out += dz; |
1216 | |
|
1217 | 0 | return true; |
1218 | 0 | } |
1219 | | |
1220 | | // --------------------------------------------------------------------------- |
1221 | | |
1222 | | template <class Grid, class GridSet, class EvaluatorIface> |
1223 | | bool Evaluator<Grid, GridSet, EvaluatorIface>::inverse( |
1224 | | EvaluatorIface &iface, double x, double y, double z, double t, |
1225 | | double &x_out, double &y_out, double &z_out) |
1226 | | |
1227 | 0 | { |
1228 | 0 | x_out = x; |
1229 | 0 | y_out = y; |
1230 | 0 | z_out = z; |
1231 | 0 | constexpr double EPS_HORIZ = 1e-12; |
1232 | 0 | constexpr double EPS_VERT = 1e-3; |
1233 | 0 | constexpr bool forInverseComputation = true; |
1234 | 0 | for (int i = 0; i < 10; i++) { |
1235 | | #ifdef DEBUG_DEFMODEL |
1236 | | iface.log("Iteration " + std::to_string(i) + ": before forward: x=" + |
1237 | | toString(x_out) + ", y=" + toString(y_out)); |
1238 | | #endif |
1239 | 0 | double x_new; |
1240 | 0 | double y_new; |
1241 | 0 | double z_new; |
1242 | 0 | if (!forward(iface, x_out, y_out, z_out, t, forInverseComputation, |
1243 | 0 | x_new, y_new, z_new)) { |
1244 | 0 | return false; |
1245 | 0 | } |
1246 | | #ifdef DEBUG_DEFMODEL |
1247 | | iface.log("After forward: x=" + toString(x_new) + |
1248 | | ", y=" + toString(y_new)); |
1249 | | #endif |
1250 | 0 | const double dx = x_new - x; |
1251 | 0 | const double dy = y_new - y; |
1252 | 0 | const double dz = z_new - z; |
1253 | 0 | x_out -= dx; |
1254 | 0 | y_out -= dy; |
1255 | 0 | z_out -= dz; |
1256 | 0 | if (std::max(std::fabs(dx), std::fabs(dy)) < EPS_HORIZ && |
1257 | 0 | std::fabs(dz) < EPS_VERT) { |
1258 | 0 | return true; |
1259 | 0 | } |
1260 | 0 | } |
1261 | 0 | return false; |
1262 | 0 | } |
1263 | | |
1264 | | // --------------------------------------------------------------------------- |
1265 | | |
1266 | | } // namespace DEFORMATON_MODEL_NAMESPACE |