1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
#define PJ_LIB__
#include "proj_internal.h"
#include "projects.h"
PROJ_HEAD(hgridshift, "Horizontal grid shift");
static XYZ forward_3d(LPZ lpz, PJ *P) {
PJ_TRIPLET point;
point.lpz = lpz;
if (P->gridlist != NULL) {
/* Only try the gridshift if at least one grid is loaded,
* otherwise just pass the coordinate through unchanged. */
point.lp = proj_hgrid_apply(P, point.lp, PJ_FWD);
}
return point.xyz;
}
static LPZ reverse_3d(XYZ xyz, PJ *P) {
PJ_TRIPLET point;
point.xyz = xyz;
if (P->gridlist != NULL) {
/* Only try the gridshift if at least one grid is loaded,
* otherwise just pass the coordinate through unchanged. */
point.lp = proj_hgrid_apply(P, point.lp, PJ_INV);
}
return point.lpz;
}
static PJ_COORD forward_4d (PJ_COORD obs, PJ *P) {
obs.xyz = forward_3d (obs.lpz, P);
return obs;
}
static PJ_COORD reverse_4d (PJ_COORD obs, PJ *P) {
obs.lpz = reverse_3d (obs.xyz, P);
return obs;
}
PJ *TRANSFORMATION(hgridshift,0) {
P->fwd4d = forward_4d;
P->inv4d = reverse_4d;
P->fwd3d = forward_3d;
P->inv3d = reverse_3d;
P->fwd = 0;
P->inv = 0;
P->left = PJ_IO_UNITS_RADIANS;
P->right = PJ_IO_UNITS_RADIANS;
if (0==pj_param(P->ctx, P->params, "tgrids").i) {
proj_log_error(P, "hgridshift: +grids parameter missing.");
return pj_default_destructor (P, PJD_ERR_NO_ARGS);
}
proj_hgrid_init(P, "grids");
/* Was gridlist compiled properly? */
if ( proj_errno(P) ) {
proj_log_error(P, "hgridshift: could not find required grid(s).");
return pj_default_destructor(P, PJD_ERR_FAILED_TO_LOAD_GRID);
}
return P;
}
#ifndef PJ_SELFTEST
/* selftest stub */
int pj_hgridshift_selftest (void) {return 0;}
#else
int pj_hgridshift_selftest (void) {
PJ *P;
PJ_COORD expect, a, b;
double dist;
/* fail on purpose: +grids parameter is mandatory*/
P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift");
if (0!=P) {
proj_destroy (P);
return 99;
}
/* fail on purpose: open non-existing grid */
P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=@nonexistinggrid.gsb,anothernonexistinggrid.gsb");
if (0!=P) {
proj_destroy (P);
return 999;
}
/* Failure most likely means the grid is missing */
P = proj_create(PJ_DEFAULT_CTX, "+proj=hgridshift +grids=nzgd2kgrid0005.gsb +ellps=GRS80");
if (0==P)
return 10;
a = proj_coord (0,0,0,0);
a.lpz.lam = PJ_TORAD(173);
a.lpz.phi = PJ_TORAD(-45);
b = a;
dist = proj_roundtrip (P, PJ_FWD, 1, &b);
if (dist > 0.00000001) {
printf("dist: %f\n",dist);
return 1;
}
expect.lpz.lam = PJ_TORAD(172.999892181021551);
expect.lpz.phi = PJ_TORAD(-45.001620431954613);
b = proj_trans(P, PJ_FWD, a);
if (proj_xy_dist(expect.xy, b.xy) > 1e-4)
return 2;
proj_destroy(P);
return 0;
}
#endif
|
