/src/proj/src/transformations/xyzgridshift.cpp

Source (jump to first uncovered line)
/******************************************************************************
 * Project:  PROJ
 * Purpose:  Geocentric translation using a grid
 * Author:   Even Rouault, <even.rouault at spatialys.com>
 *
 ******************************************************************************
 * Copyright (c) 2019, Even Rouault, <even.rouault at spatialys.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included
 * in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

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

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

#include <algorithm>

PROJ_HEAD(xyzgridshift, "Geocentric grid shift");

using namespace NS_PROJ;

namespace { // anonymous namespace
struct xyzgridshiftData {
    PJ *cart = nullptr;
    bool grid_ref_is_input = true;
    ListOfGenericGrids grids{};
    bool defer_grid_opening = false;
    int error_code_in_defer_grid_opening = 0;
    double multiplier = 1.0;
};
} // anonymous namespace

// ---------------------------------------------------------------------------

static bool get_grid_values(PJ *P, xyzgridshiftData *Q, const PJ_LP &lp,
                            double &dx, double &dy, double &dz) {
    if (Q->defer_grid_opening) {
        Q->defer_grid_opening = false;
        Q->grids = pj_generic_grid_init(P, "grids");
        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 false;
    }

    GenericShiftGridSet *gridset = nullptr;
    auto grid = pj_find_generic_grid(Q->grids, lp, gridset);
    if (!grid) {
        return false;
    }
    if (grid->isNullGrid()) {
        dx = 0;
        dy = 0;
        dz = 0;
        return true;
    }
    const auto samplesPerPixel = grid->samplesPerPixel();
    if (samplesPerPixel < 3) {
        proj_log_error(P, "xyzgridshift: grid has not enough samples");
        return false;
    }
    int sampleX = 0;
    int sampleY = 1;
    int sampleZ = 2;
    for (int i = 0; i < samplesPerPixel; i++) {
        const auto desc = grid->description(i);
        if (desc == "x_translation") {
            sampleX = i;
        } else if (desc == "y_translation") {
            sampleY = i;
        } else if (desc == "z_translation") {
            sampleZ = i;
        }
    }
    const auto unit = grid->unit(sampleX);
    if (!unit.empty() && unit != "metre") {
        proj_log_error(P, "xyzgridshift: Only unit=metre currently handled");
        return false;
    }

    bool must_retry = false;
    if (!pj_bilinear_interpolation_three_samples(P->ctx, grid, lp, sampleX,
                                                 sampleY, sampleZ, dx, dy, dz,
                                                 must_retry)) {
        if (must_retry)
            return get_grid_values(P, Q, lp, dx, dy, dz);
        return false;
    }

    dx *= Q->multiplier;
    dy *= Q->multiplier;
    dz *= Q->multiplier;
    return true;
}

// ---------------------------------------------------------------------------

#define SQUARE(x) ((x) * (x))

// ---------------------------------------------------------------------------

static PJ_COORD iterative_adjustment(PJ *P, xyzgridshiftData *Q,
                                     const PJ_COORD &pointInit, double factor) {
    PJ_COORD point = pointInit;
    for (int i = 0; i < 10; i++) {
        PJ_COORD geodetic;
        geodetic.lpz = pj_inv3d(point.xyz, Q->cart);

        double dx, dy, dz;
        if (!get_grid_values(P, Q, geodetic.lp, dx, dy, dz)) {
            return proj_coord_error();
        }

        dx *= factor;
        dy *= factor;
        dz *= factor;

        const double err = SQUARE((point.xyz.x - pointInit.xyz.x) - dx) +
                           SQUARE((point.xyz.y - pointInit.xyz.y) - dy) +
                           SQUARE((point.xyz.z - pointInit.xyz.z) - dz);

        point.xyz.x = pointInit.xyz.x + dx;
        point.xyz.y = pointInit.xyz.y + dy;
        point.xyz.z = pointInit.xyz.z + dz;
        if (err < 1e-10) {
            break;
        }
    }
    return point;
}

// ---------------------------------------------------------------------------

static PJ_COORD direct_adjustment(PJ *P, xyzgridshiftData *Q, PJ_COORD point,
                                  double factor) {
    PJ_COORD geodetic;
    geodetic.lpz = pj_inv3d(point.xyz, Q->cart);

    double dx, dy, dz;
    if (!get_grid_values(P, Q, geodetic.lp, dx, dy, dz)) {
        return proj_coord_error();
    }
    point.xyz.x += factor * dx;
    point.xyz.y += factor * dy;
    point.xyz.z += factor * dz;
    return point;
}

// ---------------------------------------------------------------------------

static PJ_XYZ pj_xyzgridshift_forward_3d(PJ_LPZ lpz, PJ *P) {
    auto Q = static_cast<xyzgridshiftData *>(P->opaque);
    PJ_COORD point = {{0, 0, 0, 0}};
    point.lpz = lpz;

    if (Q->grid_ref_is_input) {
        point = direct_adjustment(P, Q, point, 1.0);
    } else {
        point = iterative_adjustment(P, Q, point, 1.0);
    }

    return point.xyz;
}

static PJ_LPZ pj_xyzgridshift_reverse_3d(PJ_XYZ xyz, PJ *P) {
    auto Q = static_cast<xyzgridshiftData *>(P->opaque);
    PJ_COORD point = {{0, 0, 0, 0}};
    point.xyz = xyz;

    if (Q->grid_ref_is_input) {
        point = iterative_adjustment(P, Q, point, -1.0);
    } else {
        point = direct_adjustment(P, Q, point, -1.0);
    }

    return point.lpz;
}

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

    auto Q = static_cast<struct xyzgridshiftData *>(P->opaque);
    if (Q) {
        if (Q->cart)
            Q->cart->destructor(Q->cart, errlev);
        delete Q;
    }
    P->opaque = nullptr;

    return pj_default_destructor(P, errlev);
}

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

PJ *PJ_TRANSFORMATION(xyzgridshift, 0) {
    auto Q = new xyzgridshiftData;
    P->opaque = (void *)Q;
    P->destructor = pj_xyzgridshift_destructor;
    P->reassign_context = pj_xyzgridshift_reassign_context;

    P->fwd4d = nullptr;
    P->inv4d = nullptr;
    P->fwd3d = pj_xyzgridshift_forward_3d;
    P->inv3d = pj_xyzgridshift_reverse_3d;
    P->fwd = nullptr;
    P->inv = nullptr;

    P->left = PJ_IO_UNITS_CARTESIAN;
    P->right = PJ_IO_UNITS_CARTESIAN;

    // Pass a dummy ellipsoid definition that will be overridden just afterwards
    Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
    if (Q->cart == nullptr)
        return pj_xyzgridshift_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);

    /* inherit ellipsoid definition from P to Q->cart */
    pj_inherit_ellipsoid_def(P, Q->cart);

    const char *grid_ref = pj_param(P->ctx, P->params, "sgrid_ref").s;
    if (grid_ref) {
        if (strcmp(grid_ref, "input_crs") == 0) {
            // default
        } else if (strcmp(grid_ref, "output_crs") == 0) {
            // Convention use for example for NTF->RGF93 grid that contains
            // delta x,y,z from NTF to RGF93, but the grid itself is referenced
            // in RGF93
            Q->grid_ref_is_input = false;
        } else {
            proj_log_error(P, _("unusupported value for grid_ref"));
            return pj_xyzgridshift_destructor(
                P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
        }
    }

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

    /* multiplier for delta x,y,z */
    if (pj_param(P->ctx, P->params, "tmultiplier").i) {
        Q->multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
    }

    if (P->ctx->defer_grid_opening) {
        Q->defer_grid_opening = true;
    } else {
        Q->grids = pj_generic_grid_init(P, "grids");
        /* Was gridlist compiled properly? */
        if (proj_errno(P)) {
            proj_log_error(P, _("could not find required grid(s)."));
            return pj_xyzgridshift_destructor(
                P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
        }
    }

    return P;
}

Coverage Report

Created: 2025-06-22 06:59

Line	Count	Source (jump to first uncovered line)
1		/******************************************************************************
2		* Project: PROJ
3		* Purpose: Geocentric translation using a grid
4		* Author: Even Rouault, <even.rouault at spatialys.com>
5		*
6		******************************************************************************
7		* Copyright (c) 2019, 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		#include <errno.h>
29		#include <stddef.h>
30		#include <string.h>
31		#include <time.h>
32
33		#include "grids.hpp"
34		#include "proj_internal.h"
35
36		#include <algorithm>
37
38		PROJ_HEAD(xyzgridshift, "Geocentric grid shift");
39
40		using namespace NS_PROJ;
41
42		namespace { // anonymous namespace
43		struct xyzgridshiftData {
44		PJ *cart = nullptr;
45		bool grid_ref_is_input = true;
46		ListOfGenericGrids grids{};
47		bool defer_grid_opening = false;
48		int error_code_in_defer_grid_opening = 0;
49		double multiplier = 1.0;
50		};
51		} // anonymous namespace
52
53		// ---------------------------------------------------------------------------
54
55		static bool get_grid_values(PJ P, xyzgridshiftData Q, const PJ_LP &lp,
56	0	double &dx, double &dy, double &dz) {
57	0	if (Q->defer_grid_opening) {
58	0	Q->defer_grid_opening = false;
59	0	Q->grids = pj_generic_grid_init(P, "grids");
60	0	Q->error_code_in_defer_grid_opening = proj_errno(P);
61	0	}
62	0	if (Q->error_code_in_defer_grid_opening) {
63	0	proj_errno_set(P, Q->error_code_in_defer_grid_opening);
64	0	return false;
65	0	}
66
67	0	GenericShiftGridSet *gridset = nullptr;
68	0	auto grid = pj_find_generic_grid(Q->grids, lp, gridset);
69	0	if (!grid) {
70	0	return false;
71	0	}
72	0	if (grid->isNullGrid()) {
73	0	dx = 0;
74	0	dy = 0;
75	0	dz = 0;
76	0	return true;
77	0	}
78	0	const auto samplesPerPixel = grid->samplesPerPixel();
79	0	if (samplesPerPixel < 3) {
80	0	proj_log_error(P, "xyzgridshift: grid has not enough samples");
81	0	return false;
82	0	}
83	0	int sampleX = 0;
84	0	int sampleY = 1;
85	0	int sampleZ = 2;
86	0	for (int i = 0; i < samplesPerPixel; i++) {
87	0	const auto desc = grid->description(i);
88	0	if (desc == "x_translation") {
89	0	sampleX = i;
90	0	} else if (desc == "y_translation") {
91	0	sampleY = i;
92	0	} else if (desc == "z_translation") {
93	0	sampleZ = i;
94	0	}
95	0	}
96	0	const auto unit = grid->unit(sampleX);
97	0	if (!unit.empty() && unit != "metre") {
98	0	proj_log_error(P, "xyzgridshift: Only unit=metre currently handled");
99	0	return false;
100	0	}
101
102	0	bool must_retry = false;
103	0	if (!pj_bilinear_interpolation_three_samples(P->ctx, grid, lp, sampleX,
104	0	sampleY, sampleZ, dx, dy, dz,
105	0	must_retry)) {
106	0	if (must_retry)
107	0	return get_grid_values(P, Q, lp, dx, dy, dz);
108	0	return false;
109	0	}
110
111	0	dx *= Q->multiplier;
112	0	dy *= Q->multiplier;
113	0	dz *= Q->multiplier;
114	0	return true;
115	0	}
116
117		// ---------------------------------------------------------------------------
118
119	0	#define SQUARE(x) ((x) * (x))
120
121		// ---------------------------------------------------------------------------
122
123		static PJ_COORD iterative_adjustment(PJ P, xyzgridshiftData Q,
124	0	const PJ_COORD &pointInit, double factor) {
125	0	PJ_COORD point = pointInit;
126	0	for (int i = 0; i < 10; i++) {
127	0	PJ_COORD geodetic;
128	0	geodetic.lpz = pj_inv3d(point.xyz, Q->cart);
129
130	0	double dx, dy, dz;
131	0	if (!get_grid_values(P, Q, geodetic.lp, dx, dy, dz)) {
132	0	return proj_coord_error();
133	0	}
134
135	0	dx *= factor;
136	0	dy *= factor;
137	0	dz *= factor;
138
139	0	const double err = SQUARE((point.xyz.x - pointInit.xyz.x) - dx) +
140	0	SQUARE((point.xyz.y - pointInit.xyz.y) - dy) +
141	0	SQUARE((point.xyz.z - pointInit.xyz.z) - dz);
142
143	0	point.xyz.x = pointInit.xyz.x + dx;
144	0	point.xyz.y = pointInit.xyz.y + dy;
145	0	point.xyz.z = pointInit.xyz.z + dz;
146	0	if (err < 1e-10) {
147	0	break;
148	0	}
149	0	}
150	0	return point;
151	0	}
152
153		// ---------------------------------------------------------------------------
154
155		static PJ_COORD direct_adjustment(PJ P, xyzgridshiftData Q, PJ_COORD point,
156	0	double factor) {
157	0	PJ_COORD geodetic;
158	0	geodetic.lpz = pj_inv3d(point.xyz, Q->cart);
159
160	0	double dx, dy, dz;
161	0	if (!get_grid_values(P, Q, geodetic.lp, dx, dy, dz)) {
162	0	return proj_coord_error();
163	0	}
164	0	point.xyz.x += factor * dx;
165	0	point.xyz.y += factor * dy;
166	0	point.xyz.z += factor * dz;
167	0	return point;
168	0	}
169
170		// ---------------------------------------------------------------------------
171
172	0	static PJ_XYZ pj_xyzgridshift_forward_3d(PJ_LPZ lpz, PJ *P) {
173	0	auto Q = static_cast<xyzgridshiftData *>(P->opaque);
174	0	PJ_COORD point = {{0, 0, 0, 0}};
175	0	point.lpz = lpz;
176
177	0	if (Q->grid_ref_is_input) {
178	0	point = direct_adjustment(P, Q, point, 1.0);
179	0	} else {
180	0	point = iterative_adjustment(P, Q, point, 1.0);
181	0	}
182
183	0	return point.xyz;
184	0	}
185
186	0	static PJ_LPZ pj_xyzgridshift_reverse_3d(PJ_XYZ xyz, PJ *P) {
187	0	auto Q = static_cast<xyzgridshiftData *>(P->opaque);
188	0	PJ_COORD point = {{0, 0, 0, 0}};
189	0	point.xyz = xyz;
190
191	0	if (Q->grid_ref_is_input) {
192	0	point = iterative_adjustment(P, Q, point, -1.0);
193	0	} else {
194	0	point = direct_adjustment(P, Q, point, -1.0);
195	0	}
196
197	0	return point.lpz;
198	0	}
199
200	397	static PJ pj_xyzgridshift_destructor(PJ P, int errlev) {
201	397	if (nullptr == P)
202	0	return nullptr;
203
204	397	auto Q = static_cast<struct xyzgridshiftData *>(P->opaque);
205	397	if (Q) {
206	397	if (Q->cart)
207	397	Q->cart->destructor(Q->cart, errlev);
208	397	delete Q;
209	397	}
210	397	P->opaque = nullptr;
211
212	397	return pj_default_destructor(P, errlev);
213	397	}
214
215	0	static void pj_xyzgridshift_reassign_context(PJ P, PJ_CONTEXT ctx) {
216	0	auto Q = (struct xyzgridshiftData *)P->opaque;
217	0	for (auto &grid : Q->grids) {
218	0	grid->reassign_context(ctx);
219	0	}
220	0	}
221
222	397	PJ *PJ_TRANSFORMATION(xyzgridshift, 0) {
223	397	auto Q = new xyzgridshiftData;
224	397	P->opaque = (void *)Q;
225	397	P->destructor = pj_xyzgridshift_destructor;
226	397	P->reassign_context = pj_xyzgridshift_reassign_context;
227
228	397	P->fwd4d = nullptr;
229	397	P->inv4d = nullptr;
230	397	P->fwd3d = pj_xyzgridshift_forward_3d;
231	397	P->inv3d = pj_xyzgridshift_reverse_3d;
232	397	P->fwd = nullptr;
233	397	P->inv = nullptr;
234
235	397	P->left = PJ_IO_UNITS_CARTESIAN;
236	397	P->right = PJ_IO_UNITS_CARTESIAN;
237
238		// Pass a dummy ellipsoid definition that will be overridden just afterwards
239	397	Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
240	397	if (Q->cart == nullptr)
241	0	return pj_xyzgridshift_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
242
243		/* inherit ellipsoid definition from P to Q->cart */
244	397	pj_inherit_ellipsoid_def(P, Q->cart);
245
246	397	const char *grid_ref = pj_param(P->ctx, P->params, "sgrid_ref").s;
247	397	if (grid_ref) {
248	0	if (strcmp(grid_ref, "input_crs") == 0) {
249		// default
250	0	} else if (strcmp(grid_ref, "output_crs") == 0) {
251		// Convention use for example for NTF->RGF93 grid that contains
252		// delta x,y,z from NTF to RGF93, but the grid itself is referenced
253		// in RGF93
254	0	Q->grid_ref_is_input = false;
255	0	} else {
256	0	proj_log_error(P, _("unusupported value for grid_ref"));
257	0	return pj_xyzgridshift_destructor(
258	0	P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
259	0	}
260	0	}
261
262	397	if (0 == pj_param(P->ctx, P->params, "tgrids").i) {
263	1	proj_log_error(P, _("+grids parameter missing."));
264	1	return pj_xyzgridshift_destructor(P, PROJ_ERR_INVALID_OP_MISSING_ARG);
265	1	}
266
267		/* multiplier for delta x,y,z */
268	396	if (pj_param(P->ctx, P->params, "tmultiplier").i) {
269	0	Q->multiplier = pj_param(P->ctx, P->params, "dmultiplier").f;
270	0	}
271
272	396	if (P->ctx->defer_grid_opening) {
273	0	Q->defer_grid_opening = true;
274	396	} else {
275	396	Q->grids = pj_generic_grid_init(P, "grids");
276		/* Was gridlist compiled properly? */
277	396	if (proj_errno(P)) {
278	3	proj_log_error(P, _("could not find required grid(s)."));
279	3	return pj_xyzgridshift_destructor(
280	3	P, PROJ_ERR_INVALID_OP_FILE_NOT_FOUND_OR_INVALID);
281	3	}
282	396	}
283
284	393	return P;
285	396	}