1 #define PJ_LIB__
2
3 #include <errno.h>
4 #include <stddef.h>
5 #include <string.h>
6 #include <time.h>
7
8 #include "proj_internal.h"
9 #include "grids.hpp"
10
11 PROJ_HEAD(vgridshift, "Vertical grid shift");
12
13 using namespace NS_PROJ;
14
15
16 #ifdef __MINGW32__
17 // mingw32-win32 doesn't implement std::mutex
18 namespace {
19 class MyMutex {
20 public:
21 // cppcheck-suppress functionStatic
lock()22 void lock() { pj_acquire_lock(); }
23 // cppcheck-suppress functionStatic
unlock()24 void unlock() { pj_release_lock(); }
25 };
26 }
27 #else
28 #include <mutex>
29 #define MyMutex std::mutex
30 #endif
31
32 static MyMutex gMutex{};
33 static std::set<std::string> gKnownGrids{};
34
35
36 namespace { // anonymous namespace
37 struct vgridshiftData {
38 double t_final = 0;
39 double t_epoch = 0;
40 double forward_multiplier = 0;
41 ListOfVGrids grids{};
42 bool defer_grid_opening = false;
43 };
44 } // anonymous namespace
45
deal_with_vertcon_gtx_hack(PJ * P)46 static void deal_with_vertcon_gtx_hack(PJ *P)
47 {
48 struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
49 // The .gtx VERTCON files stored millimeters, but the .tif files
50 // are in metres.
51 if( Q->forward_multiplier != 0.001 ) {
52 return;
53 }
54 const char* gridname = pj_param(P->ctx, P->params, "sgrids").s;
55 if( !gridname ) {
56 return;
57 }
58 if( strcmp(gridname, "vertconw.gtx") != 0 &&
59 strcmp(gridname, "vertconc.gtx") != 0 &&
60 strcmp(gridname, "vertcone.gtx") != 0 ) {
61 return;
62 }
63 if( Q->grids.empty() ) {
64 return;
65 }
66 const auto& grids = Q->grids[0]->grids();
67 if( !grids.empty() &&
68 grids[0]->name().find(".tif") != std::string::npos ) {
69 Q->forward_multiplier = 1.0;
70 }
71 }
72
forward_3d(PJ_LPZ lpz,PJ * P)73 static PJ_XYZ forward_3d(PJ_LPZ lpz, PJ *P) {
74 struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
75 PJ_COORD point = {{0,0,0,0}};
76 point.lpz = lpz;
77
78 if ( Q->defer_grid_opening ) {
79 Q->defer_grid_opening = false;
80 Q->grids = pj_vgrid_init(P, "grids");
81 deal_with_vertcon_gtx_hack(P);
82 if ( proj_errno(P) ) {
83 return proj_coord_error().xyz;
84 }
85 }
86
87 if (!Q->grids.empty()) {
88 /* Only try the gridshift if at least one grid is loaded,
89 * otherwise just pass the coordinate through unchanged. */
90 point.xyz.z += pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
91 }
92
93 return point.xyz;
94 }
95
96
reverse_3d(PJ_XYZ xyz,PJ * P)97 static PJ_LPZ reverse_3d(PJ_XYZ xyz, PJ *P) {
98 struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
99 PJ_COORD point = {{0,0,0,0}};
100 point.xyz = xyz;
101
102 if ( Q->defer_grid_opening ) {
103 Q->defer_grid_opening = false;
104 Q->grids = pj_vgrid_init(P, "grids");
105 deal_with_vertcon_gtx_hack(P);
106 if ( proj_errno(P) ) {
107 return proj_coord_error().lpz;
108 }
109 }
110
111 if (!Q->grids.empty()) {
112 /* Only try the gridshift if at least one grid is loaded,
113 * otherwise just pass the coordinate through unchanged. */
114 point.xyz.z -= pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
115 }
116
117 return point.lpz;
118 }
119
120
forward_4d(PJ_COORD obs,PJ * P)121 static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
122 struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
123 PJ_COORD point = obs;
124
125 /* If transformation is not time restricted, we always call it */
126 if (Q->t_final==0 || Q->t_epoch==0) {
127 point.xyz = forward_3d (obs.lpz, P);
128 return point;
129 }
130
131 /* Time restricted - only apply transform if within time bracket */
132 if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
133 point.xyz = forward_3d (obs.lpz, P);
134
135
136 return point;
137 }
138
reverse_4d(PJ_COORD obs,PJ * P)139 static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
140 struct vgridshiftData *Q = (struct vgridshiftData *) P->opaque;
141 PJ_COORD point = obs;
142
143 /* If transformation is not time restricted, we always call it */
144 if (Q->t_final==0 || Q->t_epoch==0) {
145 point.lpz = reverse_3d (obs.xyz, P);
146 return point;
147 }
148
149 /* Time restricted - only apply transform if within time bracket */
150 if (obs.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch)
151 point.lpz = reverse_3d (obs.xyz, P);
152
153 return point;
154 }
155
destructor(PJ * P,int errlev)156 static PJ *destructor (PJ *P, int errlev) {
157 if (nullptr==P)
158 return nullptr;
159
160 delete static_cast<struct vgridshiftData*>(P->opaque);
161 P->opaque = nullptr;
162
163 return pj_default_destructor(P, errlev);
164 }
165
reassign_context(PJ * P,PJ_CONTEXT * ctx)166 static void reassign_context( PJ* P, PJ_CONTEXT* ctx )
167 {
168 auto Q = (struct vgridshiftData *) P->opaque;
169 for( auto& grid: Q->grids ) {
170 grid->reassign_context(ctx);
171 }
172 }
173
174
175 PJ *TRANSFORMATION(vgridshift,0) {
176 auto Q = new vgridshiftData;
177 P->opaque = (void *) Q;
178 P->destructor = destructor;
179 P->reassign_context = reassign_context;
180
181 if (!pj_param(P->ctx, P->params, "tgrids").i) {
182 proj_log_error(P, "vgridshift: +grids parameter missing.");
183 return destructor(P, PJD_ERR_NO_ARGS);
184 }
185
186 /* TODO: Refactor into shared function that can be used */
187 /* by both vgridshift and hgridshift */
188 if (pj_param(P->ctx, P->params, "tt_final").i) {
189 Q->t_final = pj_param (P->ctx, P->params, "dt_final").f;
190 if (Q->t_final == 0) {
191 /* a number wasn't passed to +t_final, let's see if it was "now" */
192 /* and set the time accordingly. */
193 if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) {
194 time_t now;
195 struct tm *date;
196 time(&now);
197 date = localtime(&now);
198 Q->t_final = 1900.0 + date->tm_year + date->tm_yday/365.0;
199 }
200 }
201 }
202
203 if (pj_param(P->ctx, P->params, "tt_epoch").i)
204 Q->t_epoch = pj_param (P->ctx, P->params, "dt_epoch").f;
205
206 /* historical: the forward direction subtracts the grid offset. */
207 Q->forward_multiplier = -1.0;
208 if (pj_param(P->ctx, P->params, "tmultiplier").i) {
209 Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
210 }
211
212 if( P->ctx->defer_grid_opening ) {
213 Q->defer_grid_opening = true;
214 }
215 else {
216 const char *gridnames = pj_param(P->ctx, P->params, "sgrids").s;
217 gMutex.lock();
218 const bool isKnownGrid = gKnownGrids.find(gridnames) != gKnownGrids.end();
219 gMutex.unlock();
220
221 if( isKnownGrid ) {
222 Q->defer_grid_opening = true;
223 }
224 else {
225 /* Build gridlist. P->vgridlist_geoid can be empty if +grids only ask for optional grids. */
226 Q->grids = pj_vgrid_init(P, "grids");
227
228 /* Was gridlist compiled properly? */
229 if ( proj_errno(P) ) {
230 proj_log_error(P, "vgridshift: could not find required grid(s).");
231 return destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
232 }
233
234 gMutex.lock();
235 gKnownGrids.insert(gridnames);
236 gMutex.unlock();
237 }
238 }
239
240 P->fwd4d = forward_4d;
241 P->inv4d = reverse_4d;
242 P->fwd3d = forward_3d;
243 P->inv3d = reverse_3d;
244 P->fwd = nullptr;
245 P->inv = nullptr;
246
247 P->left = PJ_IO_UNITS_RADIANS;
248 P->right = PJ_IO_UNITS_RADIANS;
249
250 return P;
251 }
252
pj_clear_vgridshift_knowngrids_cache()253 void pj_clear_vgridshift_knowngrids_cache() {
254 gMutex.lock();
255 gKnownGrids.clear();
256 gMutex.unlock();
257 }
258