/src/proj/src/transformations/vgridshift.cpp

Source (jump to first uncovered line)


#include <errno.h>
#include <mutex>
#include <stddef.h>
#include <string.h>
#include <time.h>

#include "grids.hpp"
#include "proj_internal.h"

PROJ_HEAD(vgridshift, "Vertical grid shift");

static std::mutex gMutexVGridShift{};
static std::set<std::string> gKnownGridsVGridShift{};

using namespace NS_PROJ;

namespace { // anonymous namespace
struct vgridshiftData {
    double t_final = 0;
    double t_epoch = 0;
    double forward_multiplier = 0;
    ListOfVGrids grids{};
    bool defer_grid_opening = false;
    int error_code_in_defer_grid_opening = 0;
};
} // anonymous namespace

static void deal_with_vertcon_gtx_hack(PJ *P) {
    struct vgridshiftData *Q = (struct vgridshiftData *)P->opaque;
    // The .gtx VERTCON files stored millimeters, but the .tif files
    // are in metres.
    if (Q->forward_multiplier != 0.001) {
        return;
    }
    const char *gridname = pj_param(P->ctx, P->params, "sgrids").s;
    if (!gridname) {
        return;
    }
    if (strcmp(gridname, "vertconw.gtx") != 0 &&
        strcmp(gridname, "vertconc.gtx") != 0 &&
        strcmp(gridname, "vertcone.gtx") != 0) {
        return;
    }
    if (Q->grids.empty()) {
        return;
    }
    const auto &grids = Q->grids[0]->grids();
    if (!grids.empty() && grids[0]->name().find(".tif") != std::string::npos) {
        Q->forward_multiplier = 1.0;
    }
}

static PJ_XYZ pj_vgridshift_forward_3d(PJ_LPZ lpz, PJ *P) {
    struct vgridshiftData *Q = (struct vgridshiftData *)P->opaque;
    PJ_COORD point = {{0, 0, 0, 0}};
    point.lpz = lpz;

    if (Q->defer_grid_opening) {
        Q->defer_grid_opening = false;
        Q->grids = pj_vgrid_init(P, "grids");
        deal_with_vertcon_gtx_hack(P);
        Q->error_code_in_defer_grid_opening = proj_errno(P);
    }
    if (Q->error_code_in_defer_grid_opening) {
        proj_errno_set(P, Q->error_code_in_defer_grid_opening);
        return proj_coord_error().xyz;
    }

    if (!Q->grids.empty()) {
        /* Only try the gridshift if at least one grid is loaded,
         * otherwise just pass the coordinate through unchanged. */
        point.xyz.z +=
            pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
    }

    return point.xyz;
}

static PJ_LPZ pj_vgridshift_reverse_3d(PJ_XYZ xyz, PJ *P) {
    struct vgridshiftData *Q = (struct vgridshiftData *)P->opaque;
    PJ_COORD point = {{0, 0, 0, 0}};
    point.xyz = xyz;

    if (Q->defer_grid_opening) {
        Q->defer_grid_opening = false;
        Q->grids = pj_vgrid_init(P, "grids");
        deal_with_vertcon_gtx_hack(P);
        Q->error_code_in_defer_grid_opening = proj_errno(P);
    }
    if (Q->error_code_in_defer_grid_opening) {
        proj_errno_set(P, Q->error_code_in_defer_grid_opening);
        return proj_coord_error().lpz;
    }

    if (!Q->grids.empty()) {
        /* Only try the gridshift if at least one grid is loaded,
         * otherwise just pass the coordinate through unchanged. */
        point.xyz.z -=
            pj_vgrid_value(P, Q->grids, point.lp, Q->forward_multiplier);
    }

    return point.lpz;
}

static void pj_vgridshift_forward_4d(PJ_COORD &coo, PJ *P) {
    struct vgridshiftData *Q = (struct vgridshiftData *)P->opaque;

    /* If transformation is not time restricted, we always call it */
    if (Q->t_final == 0 || Q->t_epoch == 0) {
        // Assigning in 2 steps avoids cppcheck warning
        // "Overlapping read/write of union is undefined behavior"
        // Cf
        // https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
        const auto xyz = pj_vgridshift_forward_3d(coo.lpz, P);
        coo.xyz = xyz;
        return;
    }

    /* Time restricted - only apply transform if within time bracket */
    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) {
        // Assigning in 2 steps avoids cppcheck warning
        // "Overlapping read/write of union is undefined behavior"
        // Cf
        // https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
        const auto xyz = pj_vgridshift_forward_3d(coo.lpz, P);
        coo.xyz = xyz;
    }
}

static void pj_vgridshift_reverse_4d(PJ_COORD &coo, PJ *P) {
    struct vgridshiftData *Q = (struct vgridshiftData *)P->opaque;

    /* If transformation is not time restricted, we always call it */
    if (Q->t_final == 0 || Q->t_epoch == 0) {
        // Assigning in 2 steps avoids cppcheck warning
        // "Overlapping read/write of union is undefined behavior"
        // Cf
        // https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
        const auto lpz = pj_vgridshift_reverse_3d(coo.xyz, P);
        coo.lpz = lpz;
        return;
    }

    /* Time restricted - only apply transform if within time bracket */
    if (coo.lpzt.t < Q->t_epoch && Q->t_final > Q->t_epoch) {
        // Assigning in 2 steps avoids cppcheck warning
        // "Overlapping read/write of union is undefined behavior"
        // Cf
        // https://github.com/OSGeo/PROJ/pull/3527#pullrequestreview-1233332710
        const auto lpz = pj_vgridshift_reverse_3d(coo.xyz, P);
        coo.lpz = lpz;
    }
}

static PJ *pj_vgridshift_destructor(PJ *P, int errlev) {
    if (nullptr == P)
        return nullptr;

    delete static_cast<struct vgridshiftData *>(P->opaque);
    P->opaque = nullptr;

    return pj_default_destructor(P, errlev);
}

static void pj_vgridshift_reassign_context(PJ *P, PJ_CONTEXT *ctx) {
    auto Q = (struct vgridshiftData *)P->opaque;
    for (auto &grid : Q->grids) {
        grid->reassign_context(ctx);
    }
}

PJ *PJ_TRANSFORMATION(vgridshift, 0) {
    auto Q = new vgridshiftData;
    P->opaque = (void *)Q;
    P->destructor = pj_vgridshift_destructor;
    P->reassign_context = pj_vgridshift_reassign_context;

    if (!pj_param(P->ctx, P->params, "tgrids").i) {
        proj_log_error(P, _("+grids parameter missing."));
        return pj_vgridshift_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
    }

    /* TODO: Refactor into shared function that can be used  */
    /* by both vgridshift and hgridshift                     */
    if (pj_param(P->ctx, P->params, "tt_final").i) {
        Q->t_final = pj_param(P->ctx, P->params, "dt_final").f;
        if (Q->t_final == 0) {
            /* a number wasn't passed to +t_final, let's see if it was "now" */
            /* and set the time accordingly.                                 */
            if (!strcmp("now", pj_param(P->ctx, P->params, "st_final").s)) {
                time_t now;
                struct tm *date;
                time(&now);
                date = localtime(&now);
                Q->t_final = 1900.0 + date->tm_year + date->tm_yday / 365.0;
            }
        }
    }

    if (pj_param(P->ctx, P->params, "tt_epoch").i)
        Q->t_epoch = pj_param(P->ctx, P->params, "dt_epoch").f;

    /* historical: the forward direction subtracts the grid offset. */
    Q->forward_multiplier = -1.0;
    if (pj_param(P->ctx, P->params, "tmultiplier").i) {
        Q->forward_multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
    }

    if (P->ctx->defer_grid_opening) {
        Q->defer_grid_opening = true;
    } else {
        const char *gridnames = pj_param(P->ctx, P->params, "sgrids").s;
        gMutexVGridShift.lock();
        const bool isKnownGrid = gKnownGridsVGridShift.find(gridnames) !=
                                 gKnownGridsVGridShift.end();
        gMutexVGridShift.unlock();

        if (isKnownGrid) {
            Q->defer_grid_opening = true;
        } else {
            /* Build gridlist. P->vgridlist_geoid can be empty if +grids only
             * ask for optional grids. */
            Q->grids = pj_vgrid_init(P, "grids");

            /* Was gridlist compiled properly? */
            if (proj_errno(P)) {
                proj_log_error(P, _("could not find required grid(s)."));
                return pj_vgridshift_destructor(
                    P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
            }

            gMutexVGridShift.lock();
            gKnownGridsVGridShift.insert(gridnames);
            gMutexVGridShift.unlock();
        }
    }

    P->fwd4d = pj_vgridshift_forward_4d;
    P->inv4d = pj_vgridshift_reverse_4d;
    P->fwd3d = pj_vgridshift_forward_3d;
    P->inv3d = pj_vgridshift_reverse_3d;
    P->fwd = nullptr;
    P->inv = nullptr;

    P->left = PJ_IO_UNITS_RADIANS;
    P->right = PJ_IO_UNITS_RADIANS;

    return P;
}

void pj_clear_vgridshift_knowngrids_cache() {
    std::lock_guard<std::mutex> lock(gMutexVGridShift);
    gKnownGridsVGridShift.clear();
}

Coverage Report

Created: 2025-06-13 06:18

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