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