/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;
};
} // 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);
        if (proj_errno(P)) {
            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);
        if (proj_errno(P)) {
            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: 2024-02-25 06:14

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