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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
#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. */
pj_apply_gridshift_3( P->ctx, P->gridlist,
P->gridlist_count, 1, 1, 0,
&point.xyz.x, &point.xyz.y, &point.xyz.z );
}
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. */
pj_apply_gridshift_3( P->ctx, P->gridlist,
P->gridlist_count, 0, 1, 0,
&point.xyz.x, &point.xyz.y, &point.xyz.z );
}
return point.lpz;
}
static PJ_OBS forward_obs(PJ_OBS obs, PJ *P) {
PJ_OBS point;
point.coo.xyz = forward_3d (obs.coo.lpz, P);
return point;
}
static PJ_OBS reverse_obs(PJ_OBS obs, PJ *P) {
PJ_OBS point;
point.coo.lpz = reverse_3d (obs.coo.xyz, P);
return point;
}
#if 0
static XY forward_xy(LP lp, PJ *P) {
PJ_TRIPLET point;
point.lp = lp;
point.lpz.z = 0;
point.xyz = forward_3d (point.lpz, P);
return point.xy;
}
static LP reverse_lp(XY xy, PJ *P) {
PJ_TRIPLET point;
point.xy = xy;
point.xyz.z = 0;
point.lpz = reverse_3d (point.xyz, P);
return point.lp;
}
#endif
PJ *PROJECTION(hgridshift) {
P->fwdobs = forward_obs;
P->invobs = reverse_obs;
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);
}
/* Build gridlist. P->gridlist can be empty if +grids only ask for optional grids. */
P->gridlist = pj_gridlist_from_nadgrids( P->ctx, pj_param(P->ctx, P->params, "sgrids").s,
&(P->gridlist_count) );
/* Was gridlist compiled properly? */
if ( pj_ctx_get_errno(pj_get_ctx(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_OBS 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_obs_null;
a.coo.lpz.lam = PJ_TORAD(173);
a.coo.lpz.phi = PJ_TORAD(-45);
dist = proj_roundtrip (P, PJ_FWD, 1, a);
if (dist > 0.00000001)
return 1;
expect.coo.lpz.lam = PJ_TORAD(172.999892181021551);
expect.coo.lpz.phi = PJ_TORAD(-45.001620431954613);
b = proj_trans_obs(P, PJ_FWD, a);
if (proj_xy_dist(expect.coo.xy, b.coo.xy) > 1e-4)
return 2;
proj_destroy(P);
return 0;
}
#endif
|
