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
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
|
/*******************************************************************************
Tiny test of an evolving new API, demonstrating simple examples of
2D and 3D transformations.
The main transformation setup object is PJ, well known from the two
former proj APIs (projects.h and proj_api.h)
The main data object PJ_OBS is new, but encapsulates the older
LP, XY etc. objects in a framework for storing a 3D observation taken at
a 4D point in space-time, and including some additional metadata (e.g.
a serial number or an epsg code).
PJ_OBS uses unions to enforce explicit statement of what kind of
coordinates are expected at a given spot in the code.
The primary elements of the API are:
pj_create (char *desc):
Create a new PJ object from a text description.
pj_trans (PJ *P, int direction, PJ_OBS obs):
Forward or inverse transformation of obs.
pj_error (PJ *P):
Check error status of P.
pj_free (PJ *P):
Clean up when done with the transformation object.
These 4 functions should cover most of the common use cases. Additional
functionality for handling thread contexts are provided by the functions
pj_debug_set, pj_error_set, pj_log_set, pj_context_renew,
pj_context_inherit, pj_context_free, pj_fileapi_set.
The proj thread contexts have not seen widespread use, so one of the
intentions with this new API is to make them less visible on the API
surface: Contexts do not have a life by themselves, they are only
visible through their associated PJs, and the number of functions
supporting them is limited.
Compilation example:
gcc -pedantic -W -Wall -Wextra -I../src -o pj_proj_test pj_proj_test.c -L../../../build/proj.4/lib -lproj_4_9
Thomas Knudsen, 2016-10-30/11-03
*******************************************************************************/
#include <proj.h>
int pj_pipeline_selftest (void);
int main (void) {
PJ *P;
PJ_OBS a, b;
char *args[3] = {"proj=utm", "zone=32", "ellps=GRS80"};
int err;
double dist;
XY cph_utm32;
/* Log everything libproj offers to log for you */
pj_log_level (0, PJ_LOG_TRACE);
/* An utm projection on the GRS80 ellipsoid */
P = pj_create ("+proj=utm +zone=32 +ellps=GRS80");
if (0==P)
return puts ("Oops"), 0;
/* Clean up */
pj_free (P);
/* Same projection, now using argc/argv style initialization */
P = pj_create_argv (3, args);
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude) to (12, 55) */
a = pj_obs_null;
/* a.coo.lp: The coordinate part of a, interpreted as a classic LP pair */
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(55);
/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);
cph_utm32 = b.coo.xy;
/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
/* Null projection */
a = pj_trans (P, PJ_IDENT, a);
printf ("IDENT: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
/* Forward again, to get two linear items for comparison */
a = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.enh.e, b.coo.enh.n);
dist = pj_xy_dist (a.coo.xy, b.coo.xy);
printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);
/* should be identical - checking whether null-init is done */
dist = pj_xyz_dist (a.coo.xyz, b.coo.xyz);
printf ("Roundtrip deviation, (nm): %15.9f\n", dist*1e9);
/* Invalid projection */
a = pj_trans (P, 42, a);
if (a.coo.lpz.lam!=DBL_MAX)
printf ("%15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
err = pj_error (P);
printf ("pj_error: %d\n", err);
/* Clean up */
pj_free (P);
/* Now do some 3D transformations */
P = pj_create ("+proj=cart +ellps=GRS80");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
a = b = pj_obs_null;
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
a.coo.lpz.z = 100;
/* Forward projection: 3D-Cartesian-to-Ellipsoidal */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);
/* Check roundtrip precision for 10000 iterations each way */
dist = pj_roundtrip (P, PJ_FWD, 10000, a);
printf ("Roundtrip deviation, fwd (nm): %15.9f\n", dist*1e9*1e5);
dist = pj_roundtrip (P, PJ_INV, 10000, b);
printf ("Roundtrip deviation, inv (nm): %15.9f\n", dist*1e9);
/* Inverse projection: Ellipsoidal-to-3D-Cartesian */
b = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);
/* Move p to another context */
pj_context_renew (P);
b = pj_trans (P, PJ_FWD, b);
printf ("CTX1: %15.9f %15.9f %15.9f\n", b.coo.xyz.x, b.coo.xyz.y, b.coo.xyz.z);
/* Move it back to the default context */
pj_context_free (P);
b = pj_trans (P, PJ_INV, b);
printf ("CTX0: %15.9f %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi), b.coo.lpz.z);
pj_free (P);
/***************************************************************************
P I P E L I N E T E S T S
****************************************************************************
***************************************************************************/
/* forward-reverse geo->utm->geo */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +step +proj=utm +inv");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (longitude, latitude, height) to (12, 55, 100) */
a = b = pj_obs_null;
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
a.coo.lpz.z = 100;
printf ("PRE: %15.9f %15.9f\n", a.coo.lpz.lam, a.coo.lpz.phi);
/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", TODEG(b.coo.lpz.lam), TODEG(b.coo.lpz.phi));
/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lpz.lam), TODEG(a.coo.lpz.phi));
pj_free (P);
/* And now the back-to-back situation utm->geo->utm */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm +inv +step +proj=utm");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
a = b = pj_obs_null;
a.coo.xy = cph_utm32;
printf ("PRE: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);
/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);
/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", a.coo.xy.x, a.coo.xy.y);
pj_free (P);
/* Finally testing a corner case: A rather pointless one-step pipeline geo->utm */
P = pj_create ("+proj=pipeline +ellps=GRS80 +zone=32 +step +proj=utm");
if (0==P)
return puts ("Oops"), 0;
/* zero initialize everything, then set (easting, northing) to utm(12, 55) */
a = b = pj_obs_null;
a.coo.lpz.lam = TORAD(12);
a.coo.lpz.phi = TORAD(55);
printf ("PRE: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
printf ("EXP: %15.9f %15.9f\n", cph_utm32.x, cph_utm32.y);
/* Forward projection */
b = pj_trans (P, PJ_FWD, a);
printf ("FWD: %15.9f %15.9f\n", b.coo.xy.x, b.coo.xy.y);
/* Inverse projection */
a = pj_trans (P, PJ_INV, b);
printf ("INV: %15.9f %15.9f\n", TODEG(a.coo.lp.lam), TODEG(a.coo.lp.phi));
/* Geodesic distance between two points with angular 2D coordinates */
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(60);
b.coo.lp.lam = TORAD(12);
b.coo.lp.phi = TORAD(61);
dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
printf ("1 deg at 60N: %15.9f\n", dist);
a.coo.lp.lam = TORAD(12);
a.coo.lp.phi = TORAD(0.);
b.coo.lp.lam = TORAD(12);
b.coo.lp.phi = TORAD(1.);
dist = pj_lp_dist (P, a.coo.lp, b.coo.lp);
printf ("1 deg at Equator: %15.9f\n", dist);
pj_free (P);
pj_pipeline_selftest ();
return 0;
}
|
