1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /*
18 * Brief synopsis: project.c reads (x,y,[z]) data and writes some combination of (x,y,z,p,q,u,v),
19 * where p,q is the distance along,across the track of the projection of (x,y),
20 * and u,v are the un-transformed (x,y) coordinates of the projected position.
21 * Can also create (x,y) along track.
22 *
23 * Author: Walter H. F. Smith
24 * Date: 1 JAN 2010
25 * Version: 6 API
26 */
27
28 #include "gmt_dev.h"
29
30 #define THIS_MODULE_CLASSIC_NAME "project"
31 #define THIS_MODULE_MODERN_NAME "project"
32 #define THIS_MODULE_LIB "core"
33 #define THIS_MODULE_PURPOSE "Project data onto lines or great circles, or generate tracks"
34 #define THIS_MODULE_KEYS "<D{,>D},G-("
35 #define THIS_MODULE_NEEDS ""
36 #define THIS_MODULE_OPTIONS "-:>Vbdefghioqs" GMT_OPT("HMm")
37
38 #define PROJECT_N_FARGS 7
39
40 struct PROJECT_CTRL { /* All control options for this program (except common args) */
41 /* active is true if the option has been activated */
42 struct PROJECT_A { /* -A<azimuth> */
43 bool active;
44 double azimuth;
45 } A;
46 struct PROJECT_C { /* -C<ox>/<oy> */
47 bool active;
48 double x, y;
49 } C;
50 struct PROJECT_E { /* -E<bx>/<by> */
51 bool active;
52 double x, y;
53 } E;
54 struct PROJECT_F { /* -F<flags> */
55 bool active;
56 char col[PROJECT_N_FARGS]; /* Character codes for desired output in the right order */
57 } F;
58 struct PROJECT_G { /* -G<inc>[<unit>[/<colat>][+c|h][+n] */
59 bool active;
60 bool header;
61 bool number;
62 bool through_C;
63 unsigned int mode;
64 double inc;
65 double colat;
66 char unit;
67 int d_mode; /* Could be negative */
68 } G;
69 struct PROJECT_L { /* -L[w][<l_min>/<l_max>] */
70 bool active;
71 bool constrain;
72 double min, max;
73 } L;
74 struct PROJECT_N { /* -N */
75 bool active;
76 } N;
77 struct PROJECT_Q { /* -Q */
78 bool active;
79 } Q;
80 struct PROJECT_S { /* -S */
81 bool active;
82 } S;
83 struct PROJECT_T { /* -T<px>/<py> */
84 bool active;
85 double x, y;
86 } T;
87 struct PROJECT_W { /* -W<w_min>/<w_max> */
88 bool active;
89 double min, max;
90 } W;
91 struct PROJECT_Z { /* -Z<major/minor/azimuth>[+e][+n] */
92 bool active;
93 bool exact;
94 bool number;
95 char unit;
96 int mode; /* Could be negative */
97 double major, minor, azimuth;
98 } Z;
99 };
100
101 struct PROJECT_DATA {
102 double a[6]; /* xypqrs */
103 double *z; /* all z's */
104 char *t; /* Text endings of each record */
105 };
106
107 struct PROJECT_INFO {
108 uint64_t n_used;
109 uint64_t n_outputs;
110 uint64_t n_z;
111 int output_choice[PROJECT_N_FARGS];
112 bool find_new_point;
113 bool want_z_output;
114 bool first;
115 double pole[3];
116 double plon, plat; /* Pole location */
117 };
118
project_compare_distances(const void * point_1,const void * point_2)119 GMT_LOCAL int project_compare_distances (const void *point_1, const void *point_2) {
120 double d_1, d_2;
121
122 d_1 = ((struct PROJECT_DATA *)point_1)->a[2];
123 d_2 = ((struct PROJECT_DATA *)point_2)->a[2];
124
125 if (d_1 < d_2) return (-1);
126 if (d_1 > d_2) return (1);
127 return (0);
128 }
129
project_oblique_setup(struct GMT_CTRL * GMT,double plat,double plon,double * p,double * clat,double * clon,double * c,bool c_given,bool generate)130 GMT_LOCAL double project_oblique_setup (struct GMT_CTRL *GMT, double plat, double plon, double *p, double *clat, double *clon, double *c, bool c_given, bool generate) {
131 /* routine sets up a unit 3-vector p, the pole of an
132 oblique projection, given plat, plon, the position
133 of this pole in the usual coordinate frame.
134 c_given = true means that clat, clon are to be used
135 as the usual coordinates of a point through which the
136 user wants the central meridian of the oblique
137 projection to go. If such a point is not given, then
138 the central meridian will go through p and the usual
139 N pole. In either case, a unit 3-vector c is created
140 which is the directed normal to the plane of the central
141 meridian, pointing in the positive normal (east) sense.
142 Latitudes and longitudes are in degrees. */
143
144 double s[3]; /* s points to the south pole */
145 double x[3]; /* tmp vector */
146 double cp, sin_lat_to_pole;
147
148 s[0] = s[1] = 0.0; s[2] = -1.0;
149
150 gmt_geo_to_cart (GMT, plat, plon, p, true);
151
152 if (c_given) gmt_geo_to_cart (GMT, *clat, *clon, s, true); /* s points to user's clat, clon */
153 gmt_cross3v (GMT, p, s, x);
154 gmt_normalize3v (GMT, x);
155 gmt_cross3v (GMT, x, p, c);
156 gmt_normalize3v (GMT, c);
157 if (!generate) gmt_M_memcpy (c, x, 3, double);
158 if (!c_given) gmt_cart_to_geo (GMT, clat, clon, c, true); /* return the adjusted center */
159 cp = gmt_dot3v (GMT, p, c);
160 sin_lat_to_pole = d_sqrt (1.0 - cp * cp);
161 return (sin_lat_to_pole);
162 }
163
project_oblique_transform(struct GMT_CTRL * GMT,double xlat,double xlon,double * x_t_lat,double * x_t_lon,double * p,double * c)164 GMT_LOCAL void project_oblique_transform (struct GMT_CTRL *GMT, double xlat, double xlon, double *x_t_lat, double *x_t_lon, double *p, double *c) {
165 /* routine takes the point x at conventional (xlat, xlon) and
166 computes the transformed coordinates (x_t_lat, x_t_lon) in
167 an oblique reference frame specified by the unit 3-vectors
168 p (the pole) and c (the directed normal to the oblique
169 central meridian). p and c have been computed earlier by
170 the routine project_oblique_setup().
171 Latitudes and longitudes are in degrees. */
172
173 double x[3], p_cross_x[3], temp1, temp2;
174
175 gmt_geo_to_cart (GMT, xlat, xlon, x, true);
176
177 temp1 = gmt_dot3v (GMT, x,p);
178 *x_t_lat = d_asind(temp1);
179
180 gmt_cross3v (GMT, p,x,p_cross_x);
181 gmt_normalize3v (GMT, p_cross_x);
182
183 temp1 = gmt_dot3v (GMT, p_cross_x, c);
184 temp2 = gmt_dot3v (GMT, x, c);
185 *x_t_lon = copysign(d_acosd(temp1), temp2);
186 }
187
project_make_euler_matrix(double * p,double * e,double theta)188 GMT_LOCAL void project_make_euler_matrix (double *p, double *e, double theta) {
189 /* Routine to fill an euler matrix e with the elements
190 needed to rotate a 3-vector about the pole p through
191 an angle theta (in degrees). p is a unit 3-vector.
192 Latitudes and longitudes are in degrees. */
193
194 double cos_theta, sin_theta, one_minus_cos_theta;
195 double pxsin, pysin, pzsin, temp;
196
197 sincosd (theta, &sin_theta, &cos_theta);
198 one_minus_cos_theta = 1.0 - cos_theta;
199
200 pxsin = p[0] * sin_theta;
201 pysin = p[1] * sin_theta;
202 pzsin = p[2] * sin_theta;
203
204 temp = p[0] * one_minus_cos_theta;
205 e[0] = temp * p[0] + cos_theta;
206 e[1] = temp * p[1] - pzsin;
207 e[2] = temp * p[2] + pysin;
208
209 temp = p[1] * one_minus_cos_theta;
210 e[3] = temp * p[0] + pzsin;
211 e[4] = temp * p[1] + cos_theta;
212 e[5] = temp * p[2] - pxsin;
213
214 temp = p[2] * one_minus_cos_theta;
215 e[6] = temp * p[0] - pysin;
216 e[7] = temp * p[1] + pxsin;
217 e[8] = temp * p[2] + cos_theta;
218 }
219
project_matrix_3v(double * a,double * x,double * b)220 GMT_LOCAL void project_matrix_3v (double *a, double *x, double *b) {
221 /* routine to find b, where Ax = b, A is a 3 by 3 square matrix,
222 and x and b are 3-vectors. A is stored row wise, that is:
223
224 A = { a11, a12, a13, a21, a22, a23, a31, a32, a33 } */
225
226 b[0] = x[0]*a[0] + x[1]*a[1] + x[2]*a[2];
227 b[1] = x[0]*a[3] + x[1]*a[4] + x[2]*a[5];
228 b[2] = x[0]*a[6] + x[1]*a[7] + x[2]*a[8];
229 }
230
project_matrix_2v(double * a,double * x,double * b)231 GMT_LOCAL void project_matrix_2v (double *a, double *x, double *b) {
232 /* routine to find b, where Ax = b, A is a 2 by 2 square matrix,
233 and x and b are 2-vectors. A is stored row wise, that is:
234
235 A = { a11, a12, a21, a22 } */
236
237 b[0] = x[0]*a[0] + x[1]*a[1];
238 b[1] = x[0]*a[2] + x[1]*a[3];
239 }
240
project_sphere_setup(struct GMT_CTRL * GMT,double alat,double alon,double * a,double blat,double blon,double * b,double azim,double * p,double * c,bool two_pts)241 GMT_LOCAL void project_sphere_setup (struct GMT_CTRL *GMT, double alat, double alon, double *a, double blat, double blon, double *b, double azim, double *p, double *c, bool two_pts) {
242 /* routine to initialize a pole vector, p, and a central meridian
243 normal vector, c, for use in projecting points onto a great circle.
244
245 The great circle is specified in either one of two ways:
246 if (two_pts), then the user has given two points, a and b,
247 which specify the great circle (directed from a to b);
248 if !(two_pts), then the user has given one point, a, and an azimuth,
249 azim, clockwise from north, which defines the projection.
250
251 The strategy is to use the project_oblique_transform operations above,
252 in such a way that the great circle of the projection is the
253 equator of an oblique transform, and the central meridian goes
254 through a. Then the transformed longitude gives the distance
255 along the projection circle, and the transformed latitude gives
256 the distance normal to the projection circle.
257
258 If (two_pts), then p = normalized(a X b). If not, we temporarily
259 create p_temp = normalized(a X n), where n is the north pole.
260 p_temp is then rotated about a through the angle azim to give p.
261 After p is found, then c = normalized(p X a).
262
263 Latitudes and longitudes are in degrees.
264 */
265
266 double e[9]; /* Euler rotation matrix, if needed */
267
268 /* First find p vector */
269
270 if (two_pts) {
271 gmt_geo_to_cart (GMT, alat, alon, a, true);
272 gmt_geo_to_cart (GMT, blat, blon, b, true);
273 gmt_cross3v (GMT, a, b, p);
274 gmt_normalize3v (GMT, p);
275 }
276 else {
277 gmt_geo_to_cart (GMT, alat, alon, a, true);
278 b[0] = b[1] = 0.0; /* set b to north pole */
279 b[2] = 1.0;
280 gmt_cross3v (GMT, a, b, c); /* use c for p_temp */
281 gmt_normalize3v (GMT, c);
282 project_make_euler_matrix(a, e, -azim);
283 project_matrix_3v(e, c, p); /* c (p_temp) rotates to p */
284 }
285
286 /* Now set c vector */
287
288 gmt_cross3v (GMT, p, a, c);
289 gmt_normalize3v (GMT, c);
290 }
291
project_flat_setup(double alat,double alon,double blat,double blon,double plat,double plon,double * azim,double * e,bool two_pts,bool pole_set,bool azim_set)292 GMT_LOCAL void project_flat_setup (double alat, double alon, double blat, double blon, double plat, double plon, double *azim, double *e, bool two_pts, bool pole_set, bool azim_set) {
293 /* Sets up stuff for rotation of Cartesian 2-vectors, analogous
294 to the spherical three vector stuff above.
295 Output is the Cartesian azimuth in degrees.
296 Latitudes and longitudes are in degrees. */
297
298 if (two_pts)
299 *azim = d_atan2d (blat - alat, blon - alon);
300 else if (pole_set)
301 *azim = d_atan2d (plat - alat, plon - alon);
302 else if (azim_set)
303 *azim = 90 - *azim;
304
305 e[0] = e[3] = cosd (*azim);
306 e[1] = sind (*azim);
307 e[2] = -e[1];
308 }
309
New_Ctrl(struct GMT_CTRL * GMT)310 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
311 struct PROJECT_CTRL *C;
312
313 C = gmt_M_memory (GMT, NULL, 1U, struct PROJECT_CTRL);
314
315 /* Initialize values whose defaults are not 0/false/NULL */
316
317 C->G.colat = 90.0; /* Great circle path */
318 return (C);
319 }
320
Free_Ctrl(struct GMT_CTRL * GMT,struct PROJECT_CTRL * C)321 static void Free_Ctrl (struct GMT_CTRL *GMT, struct PROJECT_CTRL *C) { /* Deallocate control structure */
322 if (!C) return;
323 gmt_M_free (GMT, C);
324 }
325
usage(struct GMTAPI_CTRL * API,int level)326 static int usage (struct GMTAPI_CTRL *API, int level) {
327 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
328 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
329 GMT_Usage (API, 0, "usage: %s [<table>] -C<ox>/<oy> [-A<azimuth>] [-E<bx>/<by>] [-F<flags>] [-G<dist>[unit][/<colat>][+c][+h][+n]] "
330 "[-L[w|<l_min>/<l_max>]] [-N] [-Q] [-S] [-T<px>/<py>] [%s] [-W<w_min>/<w_max>] [-Z<major>[unit][/<minor>/<azimuth>][+e]] "
331 "[%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n",
332 name, GMT_V_OPT, GMT_b_OPT, GMT_d_OPT, GMT_e_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT,
333 GMT_q_OPT, GMT_s_OPT, GMT_colon_OPT, GMT_PAR_OPT);
334
335 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
336
337 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
338 GMT_Usage (API, 1, "\nNote: project does not want input if -G option is given. "
339 "The projection may be defined in (only) one of three ways:");
340 GMT_Usage (API, 3, "1. By a center -C and an azimuth -A,");
341 GMT_Usage (API, 3, "2. By a center -C and end point of the path -E,");
342 GMT_Usage (API, 3, "3. By a center -C and a roTation pole position -T.");
343 GMT_Usage (API, -2, "In a spherical projection [default], all cases place the central meridian "
344 "of the transformed coordinates (p,q) through -C (p = 0 at -C). The equator "
345 "of the (p,q) system (line q = 0) passes through -C and makes an angle "
346 "<azimuth> with North (case 1), or passes through -E (case 2), or is "
347 "determined by the pole -T (case 3). In (3), point -C need not be on equator. "
348 "In a Cartesian [-N option] projection, p = q = 0 at -O in all cases; "
349 "(1) and (2) orient the p axis, while (3) orients the q axis.");
350 GMT_Option (API, "<");
351 GMT_Usage (API, 1, "\n-C<ox>/<oy>");
352 GMT_Usage (API, -2, "Set the location of the center to be <ox>/<oy>.");
353 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
354 GMT_Usage (API, 1, "\n-A<azimuth>");
355 GMT_Usage (API, -2, "Set the option (1) Azimuth, (<azimuth> in degrees CW from North).");
356 GMT_Usage (API, 1, "\n-E<bx>/<by>");
357 GMT_Usage (API, -2, "Set the option (2) location of end point E to be <bx>/<by>.");
358 GMT_Usage (API, 1, "\n-F<flags>");
359 GMT_Usage (API, -2, "Indicate what output you want as one or more of xyzpqrs in any order:");
360 GMT_Usage (API, 3, "%s x,y,[z] refer to input data locations and optional values.", GMT_LINE_BULLET);
361 GMT_Usage (API, 3, "%s p,q are the coordinates of x,y in the projection's coordinate system.", GMT_LINE_BULLET);
362 GMT_Usage (API, 3, "%s r,s is the projected position of x,y (taking q = 0) in the (x,y) coordinate system.", GMT_LINE_BULLET);
363 GMT_Usage (API, -2, "Note 1: p,q may be scaled from degrees into kilometers by the -Q option. See -L, -Q, -W. "
364 "Note 2: z refers to all input data columns beyond the required x,y. "
365 "Note 3: If -G is set, -F is not available and output defaults to rsp "
366 "[Default is all fields, i.e., -Fxyzpqrs].");
367 GMT_Usage (API, 1, "\n-G<dist>[<unit>][/<colat>][+c|h]");
368 GMT_Usage (API, -2, "Generate (r,s,p) points along profile every <dist> units. (No input data used.) "
369 "If E given, will generate from C to E; else must give -L<l_min>/<l_max> for length. "
370 "Optionally, append /<colat> for a small circle path through C and E (requires -C -E). Some modifiers:");
371 GMT_Usage (API, 3, "+c If given -T and -C, compute and use <colat> of C [90].");
372 GMT_Usage (API, 3, "+h Place information about the pole in a segment header [no header].");
373 GMT_Usage (API, 3, "+n Increment specifies the number of points instead (requires -C -E or -Z).");
374 GMT_Usage (API, -2, "For geographic profile, append e (meter), f (foot), k (km), M (mile), n (nautical mile), u (survey foot), d (arc degree), "
375 "m (arc minute), or s (arc second) to <inc> [k]. ");
376 GMT_Usage (API, 1, "\n-L[w|<l_min>/<l_max>]");
377 GMT_Usage (API, -2, "Check the Length along the projected track and use only certain points:n");
378 GMT_Usage (API, 3, "%s Append w to only use those points Within the span from C to E (Must have set -E).", GMT_LINE_BULLET);
379 GMT_Usage (API, 3, "%s Append <l_min>/<l_max> will only use points whose p is [l_min <= p <= l_max].", GMT_LINE_BULLET);
380 GMT_Usage (API, -2, "Default uses all points. Note p = 0 at C and increases toward E in <azimuth> direction.");
381 GMT_Usage (API, 1, "\n-N Flat Earth mode; a Cartesian projection is made. Default is spherical.");
382 GMT_Usage (API, 1, "\n-Q Convert to Map units, so x,y,r,s are degrees, "
383 "while p,q,dist,l_min,l_max,w_min,w_max are km. "
384 "If not set, then p,q,dist,l_min,l_max,w_min,w_max are assumed to be in same units as x,y,r,s.");
385 GMT_Usage (API, 1, "\n-S Output should be Sorted into increasing p value.");
386 GMT_Usage (API, 1, "\n-T<px>/<py>");
387 GMT_Usage (API, -2, "Set the option (3) location of the roTation pole to the projection to be <px>/<py>.");
388 GMT_Option (API, "V");
389 GMT_Usage (API, 1, "\n-W<w_min>/<w_max>");
390 GMT_Usage (API, -2, "Check the width across the projected track and use only certain points. "
391 "This will use only those points whose q is [w_min <= q <= w_max]. "
392 "Note: q is positive to your LEFT as you walk from C toward E in the <azimuth> direction.");
393 GMT_Usage (API, 1, "\n-Z<major>[unit][/<minor>/<azimuth>][+e|n]");
394 GMT_Usage (API, -2, "With -G and -C, generate an ellipse of given major and minor axes and the azimuth of the major axis.");
395 GMT_Usage (API, 3, "+e Adjust increment to fix perimeter exactly [use increment as given in -G].");
396 GMT_Usage (API, -2, "For a degenerate ellipse, i.e., circle, you may just give <major> only as the <diameter>. "
397 "Append e (meter), f (foot), k (km), M (mile), n (nautical mile), u (survey foot), d (arc degree), "
398 "m (arc minute), or s (arc second) to <major> [k]; we assume -G is in the same units (unless +n is used). "
399 "For Cartesian ellipses, add -N and provide <direction> (CCW from horizontal) instead of <azimuth>.");
400 GMT_Option (API, "bi2,bo,d,e,f,g,h,i,o,q,s,:,.");
401
402 return (GMT_MODULE_USAGE);
403 }
404
parse(struct GMT_CTRL * GMT,struct PROJECT_CTRL * Ctrl,struct GMT_OPTION * options)405 static int parse (struct GMT_CTRL *GMT, struct PROJECT_CTRL *Ctrl, struct GMT_OPTION *options) {
406 /* This parses the options provided to project and sets parameters in CTRL.
407 * Any GMT common options will override values set previously by other commands.
408 * It also replaces any file names specified as input or output with the data ID
409 * returned when registering these sources/destinations with the API.
410 */
411
412 bool n_active = false;
413 unsigned int n_errors = 0, j, k, pos;
414 size_t len;
415 char txt_a[GMT_LEN64] = {""}, txt_b[GMT_LEN64] = {""}, p[GMT_LEN256] = {""}, *ce = NULL, *ch = NULL, *c = NULL, dummy;
416 struct GMT_OPTION *opt = NULL;
417 struct GMTAPI_CTRL *API = GMT->parent;
418
419 for (opt = options; opt; opt = opt->next) {
420 if (opt->option == 'N') /* Must find -N first, if given */
421 Ctrl->N.active = true;
422 else if (opt->option == 'Q') { /* Geographic coordinates */
423 gmt_set_geographic (GMT, GMT_IN);
424 gmt_set_geographic (GMT, GMT_OUT);
425 }
426 }
427
428 for (opt = options; opt; opt = opt->next) {
429 switch (opt->option) {
430
431 case '<': /* Input files are OK */
432 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
433 break;
434
435 /* Processes program-specific parameters */
436
437 case 'A':
438 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
439 Ctrl->A.active = true;
440 Ctrl->A.azimuth = atof (opt->arg);
441 break;
442 case 'C':
443 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
444 Ctrl->C.active = true;
445 if (sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b) != 2) {
446 GMT_Report (API, GMT_MSG_ERROR, "Option -C: Expected -C<lon0>/<lat0>\n");
447 n_errors++;
448 }
449 else {
450 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->C.x), txt_a);
451 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->C.y), txt_b);
452 if (n_errors) GMT_Report (API, GMT_MSG_ERROR, "Option -C: Undecipherable argument %s\n", opt->arg);
453 }
454 break;
455 case 'D':
456 if (gmt_M_compat_check (GMT, 4)) {
457 GMT_Report (API, GMT_MSG_COMPAT, "Option -D is deprecated; use --FORMAT_GEO_OUT instead\n");
458 if (opt->arg[0] == 'g') GMT->current.io.geo.range = GMT_IS_0_TO_P360_RANGE;
459 if (opt->arg[0] == 'd') GMT->current.io.geo.range = GMT_IS_M180_TO_P180_RANGE;
460 }
461 else
462 n_errors += gmt_default_error (GMT, opt->option);
463 break;
464 case 'E':
465 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
466 Ctrl->E.active = true;
467 if (sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b) != 2) {
468 GMT_Report (API, GMT_MSG_ERROR, "Option -E: Expected -E<lon1>/<lat1>\n");
469 n_errors++;
470 }
471 else {
472 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->E.x), txt_a);
473 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->E.y), txt_b);
474 if (n_errors) GMT_Report (API, GMT_MSG_ERROR, "Option -E: Undecipherable argument %s\n", opt->arg);
475 }
476 break;
477 case 'F':
478 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
479 Ctrl->F.active = true;
480 for (j = 0, k = 0; opt->arg[j]; j++, k++) {
481 if (k < PROJECT_N_FARGS) {
482 Ctrl->F.col[k] = opt->arg[j];
483 if (!strchr ("xyzpqrs", opt->arg[j])) {
484 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Choose from -Fxyzpqrs\n");
485 n_errors++;
486 }
487 }
488 else {
489 n_errors++;
490 GMT_Report (API, GMT_MSG_ERROR, "Option -F: Too many output columns selected\n");
491 }
492 }
493 break;
494 case 'G':
495 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
496 Ctrl->G.active = true;
497 len = strlen (opt->arg) - 1;
498 if (len > 0 && opt->arg[len] == '+') { /* Obsolete way to say +h */
499 Ctrl->G.header = true; /* Wish to place a segment header on output */
500 opt->arg[len] = 0; /* Temporarily remove the trailing + sign */
501 ch = &opt->arg[len];
502 }
503 if ((c = gmt_first_modifier (GMT, opt->arg, "chn"))) { /* Process any modifiers */
504 pos = 0; /* Reset to start of new word */
505 while (gmt_getmodopt (GMT, 'G', c, "chn", &pos, p, &n_errors) && n_errors == 0) {
506 switch (p[0]) {
507 case 'c': Ctrl->G.through_C = true; break; /* Compute required colatitude for small circle to go through C */
508 case 'h': Ctrl->G.header = true; break; /* Output segment header */
509 case 'n': Ctrl->G.number = true; break; /* Gave number of points rather than increment */
510 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
511 }
512 }
513 c[0] = '\0'; /* Hide modifiers */
514 }
515 if (sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b) == 2) { /* Got dist/colat */
516 Ctrl->G.mode = 1;
517 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->G.colat), txt_b);
518 }
519 else
520 strcpy (txt_a, opt->arg);
521 if (Ctrl->G.number) {
522 Ctrl->Z.number = true; /* Since -Z +n needs backwards compatibility support */
523 Ctrl->G.inc = atof (txt_a);
524 }
525 else {
526 if (gmt_M_is_geographic (GMT, GMT_IN) && strchr (GMT_LEN_UNITS, txt_a[strlen(txt_a)-1]) == NULL) /* No unit given, append k for default km */
527 strcat (txt_a, "k");
528 Ctrl->G.d_mode = gmt_get_distance (GMT, txt_a, &(Ctrl->G.inc), &(Ctrl->G.unit));
529 }
530 if (ch) ch[0] = '+'; /* Restore the obsolete plus-sign */
531 if (c) c[0] = '+'; /* Restore modifier */
532 break;
533 case 'L':
534 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
535 Ctrl->L.active = true;
536 if (opt->arg[0] == 'W' || opt->arg[0] == 'w')
537 Ctrl->L.constrain = true;
538 else {
539 n_errors += gmt_M_check_condition (GMT, sscanf(opt->arg, "%lf/%lf", &Ctrl->L.min, &Ctrl->L.max) != 2, "Option -L: Expected -L[w | <min>/<max>]\n");
540 }
541 break;
542 case 'N': /* Handled above but still in argv */
543 n_errors += gmt_M_repeated_module_option (API, n_active);
544 n_active = true;
545 break;
546 case 'Q':
547 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
548 Ctrl->Q.active = true;
549 break;
550 case 'S':
551 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
552 Ctrl->S.active = true;
553 break;
554 case 'T':
555 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
556 Ctrl->T.active = true;
557 if (sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b) != 2) {
558 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Expected -T<lonp>/<latp>\n");
559 n_errors++;
560 }
561 else {
562 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, txt_a, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->T.x), txt_a);
563 n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, txt_b, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->T.y), txt_b);
564 if (n_errors) GMT_Report (API, GMT_MSG_ERROR, "Option -T: Undecipherable argument %s\n", opt->arg);
565 }
566 break;
567 case 'W':
568 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
569 Ctrl->W.active = true;
570 n_errors += gmt_M_check_condition (GMT, sscanf (opt->arg, "%lf/%lf", &Ctrl->W.min, &Ctrl->W.max) != 2,
571 "Option -W: Expected -W<min>/<max>\n");
572 break;
573 case 'Z': /* Parameters of ellipse */
574 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
575 Ctrl->Z.active = true;
576 if ((ce = strstr (opt->arg, "+e"))) {
577 Ctrl->Z.exact = true;
578 ce[0] = '\0'; /* Chop off +e */
579 }
580 else if ((ce = strstr (opt->arg, "+n"))) {
581 Ctrl->Z.number = true;
582 ce[0] = '\0'; /* Chop off +n */
583 }
584 if (strchr (opt->arg, '/')) { /* Ellipse setting */
585 k = sscanf (opt->arg, "%[^/]/%[^/]/%lf", txt_a, txt_b, &Ctrl->Z.azimuth);
586 if (k == 3) { /* Ellipse, check that major >= minor */
587 if (Ctrl->N.active) {
588 Ctrl->Z.major = atof (txt_a);
589 Ctrl->Z.minor = atof (txt_b);
590 }
591 else { /* Watch out for units and convert to km */
592 if (gmt_M_is_geographic (GMT, GMT_IN) && strchr (GMT_LEN_UNITS, txt_a[strlen(txt_a)-1]) == NULL) /* No unit given, append k for default km */
593 strcat (txt_a, "k");
594 Ctrl->Z.mode = gmt_get_distance (GMT, txt_a, &(Ctrl->Z.major), &(Ctrl->Z.unit));
595 (void) gmt_get_distance (GMT, txt_b, &(Ctrl->Z.minor), &dummy);
596 }
597 if (Ctrl->Z.major < Ctrl->Z.minor) {
598 GMT_Report (API, GMT_MSG_ERROR, "Option -Z: Major axis must be equal to or larger than the minor axis\n");
599 n_errors++;
600 }
601 }
602 else { /* Bad number of arguments */
603 if (Ctrl->N.active)
604 GMT_Report (API, GMT_MSG_ERROR, "Option -Z: Expected -Z<major/minor/direction>[+e|n] or -Z<diameter>[+e|n]\n");
605 else
606 GMT_Report (API, GMT_MSG_ERROR, "Option -Z: Expected -Z<major[unit]/minor/azimuth>[+e|n] or -Z<diameter>[unit][+e|n]\n");
607 n_errors++;
608 }
609 }
610 else { /* Circle as degenerated ellipse */
611 if (Ctrl->N.active)
612 Ctrl->Z.minor = Ctrl->Z.major = atof (opt->arg);
613 else {
614 strcpy (txt_a, opt->arg);
615 if (gmt_M_is_geographic (GMT, GMT_IN) && strchr (GMT_LEN_UNITS, txt_a[strlen(txt_a)-1]) == NULL) /* No unit given, append k for default km */
616 strcat (txt_a, "k");
617 Ctrl->Z.mode = gmt_get_distance (GMT, opt->arg, &(Ctrl->Z.minor), &(Ctrl->Z.unit));
618 Ctrl->Z.major = Ctrl->Z.minor;
619 }
620 }
621
622 if (ce) ce[0] = '+'; /* Restore the plus-sign */
623 if (Ctrl->N.active) { /* Convert to semi-axis or radius and azimuth since that is now the Cartesian code operates */
624 Ctrl->Z.minor /= 2.0;
625 Ctrl->Z.major /= 2.0;
626 if (k == 3) Ctrl->Z.azimuth = 90.0 - Ctrl->Z.azimuth;
627 }
628 break;
629 default: /* Report bad options */
630 n_errors += gmt_default_error (GMT, opt->option);
631 break;
632 }
633 }
634
635 if (!Ctrl->N.active && ((Ctrl->C.active && (Ctrl->C.x < -360 || Ctrl->C.x > 360) && (Ctrl->C.y < -90 || Ctrl->C.y > 90)) || (Ctrl->E.active && (Ctrl->E.x < -360 || Ctrl->E.x > 360) && (Ctrl->E.y < -90 || Ctrl->E.y > 90)))) {
636 GMT_Report (API, GMT_MSG_ERROR, "Your -C or -E options suggest Cartesian coordinates. Please see -N\n");
637 n_errors++;
638 }
639
640 n_errors += gmt_M_check_condition (GMT, Ctrl->G.number && !(Ctrl->C.active && Ctrl->E.active),
641 "Option -G: Can only use +n if both -C and -E are specified\n");
642 n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && !Ctrl->L.constrain && Ctrl->L.min >= Ctrl->L.max,
643 "Option -L: w_min must be < w_max\n");
644 n_errors += gmt_M_check_condition (GMT, Ctrl->W.active && Ctrl->W.min >= Ctrl->W.max,
645 "Option -W: w_min must be < w_max\n");
646 n_errors += gmt_M_check_condition (GMT, (Ctrl->A.active + Ctrl->E.active + Ctrl->T.active) > 1,
647 "Specify only one of -A, -E, and -T\n");
648 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->T.active,
649 "Option -T: Cannot be used with -N; specify using -E or -A instead\n");
650 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active && (Ctrl->C.x == Ctrl->E.x) && (Ctrl->C.y == Ctrl->E.y),
651 "Option -E: Second point must differ from origin!\n");
652 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && Ctrl->L.min == Ctrl->L.max && !(Ctrl->E.active || Ctrl->Z.active),
653 "Option -G: Must also specify -Lmin/max or use -E instead\n");
654 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && Ctrl->F.active,
655 "Option -G: -F not allowed [Defaults to rsp]\n");
656 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && Ctrl->G.inc <= 0.0,
657 "Option -G: Must specify a positive increment\n");
658 n_errors += gmt_M_check_condition (GMT, Ctrl->L.constrain && !Ctrl->E.active,
659 "Option -L: Must specify -Lmin/max or use -E instead\n");
660 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && (gmt_M_is_geographic (GMT, GMT_IN) || gmt_M_is_geographic (GMT, GMT_OUT)),
661 "Option -N: Cannot be used with -fg\n");
662 n_errors += gmt_M_check_condition (GMT, Ctrl->N.active && Ctrl->G.mode,
663 "Option -N: Cannot be used with -G<dist>/<colat>\n");
664 n_errors += gmt_M_check_condition (GMT, Ctrl->G.through_C && !Ctrl->C.active,
665 "Option -G: Modifier +c requires both -C and -T\n");
666 n_errors += gmt_M_check_condition (GMT, Ctrl->Z.active && !Ctrl->G.active,
667 "Option -Z: Requires option -G as well\n");
668 n_errors += gmt_M_check_condition (GMT, Ctrl->G.through_C && !Ctrl->T.active,
669 "Option -G: Modifier +c requires both -C and -T\n");
670 n_errors += gmt_check_binary_io (GMT, 2);
671
672 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
673 }
674
project_write_one_segment(struct GMT_CTRL * GMT,struct PROJECT_CTRL * Ctrl,double theta,struct PROJECT_DATA * p_data,struct PROJECT_INFO * P)675 GMT_LOCAL int project_write_one_segment (struct GMT_CTRL *GMT, struct PROJECT_CTRL *Ctrl, double theta, struct PROJECT_DATA *p_data, struct PROJECT_INFO *P) {
676 uint64_t col, n_items, rec, k;
677 double sin_theta, cos_theta, e[9], x[3], xt[3], *out = NULL;
678 struct GMT_RECORD *Out = NULL;
679
680 if (Ctrl->S.active) qsort (p_data, P->n_used, sizeof (struct PROJECT_DATA), project_compare_distances);
681
682 /* Get here when all data are loaded with p,q and p is in increasing order if desired. */
683
684 if (P->find_new_point) { /* We need to find r,s */
685
686 if (Ctrl->N.active) {
687 sincosd (theta, &sin_theta, &cos_theta);
688 for (rec = 0; rec < P->n_used; rec++) {
689 p_data[rec].a[4] = Ctrl->C.x + p_data[rec].a[2] * cos_theta;
690 p_data[rec].a[5] = Ctrl->C.y + p_data[rec].a[2] * sin_theta;
691 }
692 }
693 else {
694 gmt_geo_to_cart (GMT, Ctrl->C.y, Ctrl->C.x, x, true);
695 for (rec = 0; rec < P->n_used; rec++) {
696 project_make_euler_matrix (P->pole, e, p_data[rec].a[2]);
697 project_matrix_3v (e,x,xt);
698 gmt_cart_to_geo (GMT, &(p_data[rec].a[5]), &(p_data[rec].a[4]), xt, true);
699 }
700 }
701 }
702
703 /* At this stage, all values are still in degrees. */
704
705 if (Ctrl->Q.active) {
706 for (rec = 0; rec < P->n_used; rec++) {
707 p_data[rec].a[2] *= GMT->current.proj.DIST_KM_PR_DEG;
708 p_data[rec].a[3] *= GMT->current.proj.DIST_KM_PR_DEG;
709 }
710 }
711
712 n_items = gmt_get_cols (GMT, GMT_OUT);
713 out = gmt_M_memory (GMT, NULL, n_items, double);
714 Out = gmt_new_record (GMT, out, NULL); /* text will be set below */
715
716 /* Now output this segment */
717
718 for (rec = 0; rec < P->n_used; rec++) {
719 for (col = k = 0; col < P->n_outputs; col++) {
720 if (P->output_choice[col] == -1) { /* Copy over all z columns */
721 gmt_M_memcpy (&out[k], p_data[rec].z, P->n_z, double);
722 k += P->n_z;
723 }
724 else
725 out[k++] = p_data[rec].a[P->output_choice[col]];
726 }
727 Out->text = p_data[rec].t; /* The trailing text (may be NULL) */
728 GMT_Put_Record (GMT->parent, GMT_WRITE_DATA, Out); /* Write this to output */
729 }
730 gmt_M_free (GMT, out);
731 gmt_M_free (GMT, Out);
732 return (0);
733 }
734
735 #define bailout(code) {gmt_M_free_options (mode); return (code);}
736 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
737
GMT_project(void * V_API,int mode,void * args)738 EXTERN_MSC int GMT_project (void *V_API, int mode, void *args) {
739 uint64_t rec, n_total_read, col, n_total_used = 0;
740 bool skip, z_first = true, z_set_auto = false;
741 int error = 0;
742
743 size_t n_alloc = GMT_CHUNK;
744
745 double xx, yy, cos_theta, sin_theta, sin_lat_to_pole = 1.0;
746 double theta = 0.0, d_along, sign = 1.0, s, c, *in = NULL;
747 double a[3], b[3], x[3], xt[3], center[3], e[9];
748
749 struct PROJECT_DATA *p_data = NULL;
750 struct PROJECT_INFO P;
751 struct PROJECT_CTRL *Ctrl = NULL;
752 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
753 struct GMT_OPTION *options = NULL;
754 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
755
756 /*----------------------- Standard module initialization and parsing ----------------------*/
757
758 if (API == NULL) return (GMT_NOT_A_SESSION);
759 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
760 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
761
762 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
763
764 /* Parse the command-line arguments */
765
766 if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, NULL, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
767 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
768 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
769 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
770
771 /*---------------------------- This is the project main code ----------------------------*/
772
773 gmt_M_memset (&P, 1, struct PROJECT_INFO);
774 gmt_M_memset (a, 3, double);
775 gmt_M_memset (b, 3, double);
776 gmt_M_memset (x, 3, double);
777 gmt_M_memset (xt, 3, double);
778 gmt_M_memset (center, 3, double);
779 gmt_M_memset (e, 9, double);
780 P.first = true;
781 if (Ctrl->N.active) { /* Must undo an optional -fg that was set before */
782 gmt_set_cartesian (GMT, GMT_IN);
783 gmt_set_cartesian (GMT, GMT_OUT);
784 }
785 else { /* Make sure we set -fg */
786 gmt_set_geographic (GMT, GMT_IN);
787 gmt_set_geographic (GMT, GMT_OUT);
788 if (Ctrl->Z.mode && Ctrl->Z.unit != 'k') { /* Degenerate ellipse with diameter given in units */
789 if (gmt_init_distaz (GMT, Ctrl->Z.unit, Ctrl->Z.mode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE)
790 Return (GMT_NOT_A_VALID_TYPE);
791 if (GMT->current.map.dist[GMT_MAP_DIST].arc) { /* Got angular measures, convert to km */
792 Ctrl->Z.minor *= (GMT->current.proj.KM_PR_DEG / GMT->current.map.dist[GMT_MAP_DIST].scale); /* Now in km */
793 Ctrl->Z.major *= (GMT->current.proj.KM_PR_DEG / GMT->current.map.dist[GMT_MAP_DIST].scale); /* Now in km */
794 }
795 else {
796 Ctrl->Z.minor /= (GMT->current.map.dist[GMT_MAP_DIST].scale * METERS_IN_A_KM); /* Now in km */
797 Ctrl->Z.major /= (GMT->current.map.dist[GMT_MAP_DIST].scale * METERS_IN_A_KM); /* Now in km */
798 }
799 }
800 }
801
802 /* Convert user's -F choices to internal parameters */
803 for (col = P.n_outputs = 0; col < PROJECT_N_FARGS && Ctrl->F.col[col]; col++) {
804 switch (Ctrl->F.col[col]) {
805 case 'z': /* Special flag, can mean any number of z columns and trailing text */
806 P.output_choice[col] = -1;
807 P.want_z_output = true;
808 break;
809 case 'x':
810 P.output_choice[col] = 0;
811 break;
812 case 'y':
813 P.output_choice[col] = 1;
814 break;
815 case 'p':
816 P.output_choice[col] = 2;
817 break;
818 case 'q':
819 P.output_choice[col] = 3;
820 break;
821 case 'r':
822 P.output_choice[col] = 4;
823 P.find_new_point = true;
824 break;
825 case 's':
826 P.output_choice[col] = 5;
827 P.find_new_point = true;
828 break;
829 default: /* Already checked in parser that this cannot happen */
830 break;
831 }
832 P.n_outputs++; /* Since we count the extra z's later and add them */
833 }
834
835 if (P.n_outputs == 0 && !Ctrl->G.active) { /* Generate default -F setting (xyzpqrs) */
836 P.n_outputs = PROJECT_N_FARGS;
837 for (col = 0; col < 2; col++) P.output_choice[col] = (int)col; /* Do xy */
838 P.output_choice[2] = -1; /* Do z as col 2 */
839 P.want_z_output = z_set_auto = true;
840 for (col = 3; col < P.n_outputs; col++) P.output_choice[col] = (int)col - 1; /* Do pqrs */
841 P.find_new_point = true;
842 }
843 if (Ctrl->G.active) { /* Hardwire 3 output columns and set their types */
844 if (gmt_M_is_geographic (GMT, GMT_IN) && Ctrl->G.d_mode) {
845 if (gmt_init_distaz (GMT, Ctrl->G.unit, Ctrl->G.d_mode, GMT_MAP_DIST) == GMT_NOT_A_VALID_TYPE)
846 Return (GMT_NOT_A_VALID_TYPE);
847 }
848 P.n_outputs = 3;
849 gmt_set_column_type (GMT, GMT_OUT, GMT_X, (Ctrl->N.active) ? GMT_IS_FLOAT : GMT_IS_LON);
850 gmt_set_column_type (GMT, GMT_OUT, GMT_Y, (Ctrl->N.active) ? GMT_IS_FLOAT : GMT_IS_LAT);
851 gmt_set_column_type (GMT, GMT_OUT, GMT_Z, GMT_IS_FLOAT);
852 if (Ctrl->Q.active && Ctrl->G.unit != 'k' && !Ctrl->G.number)
853 Ctrl->G.inc *= 0.001 / GMT->current.map.dist[GMT_MAP_DIST].scale; /* Now in km */
854 }
855 else if (!Ctrl->N.active) { /* Decode and set the various output column types in the geographic case */
856 for (col = 0; col < P.n_outputs; col++) {
857 switch (P.output_choice[col]) {
858 case 0: case 4: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_LON); break;
859 case 1: case 5: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_LAT); break;
860 default: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_FLOAT); break;
861 }
862 }
863 }
864 p_data = gmt_M_memory (GMT, NULL, n_alloc, struct PROJECT_DATA);
865
866 if (Ctrl->G.active && Ctrl->E.active && (Ctrl->L.min == Ctrl->L.max)) Ctrl->L.constrain = true; /* Default generate from A to B */
867
868 /* Set up rotation matrix e for flat earth, or pole and center for spherical; get Ctrl->L.min, Ctrl->L.max if stay_within */
869
870 if (Ctrl->N.active) { /* Flat Earth mode */
871 theta = Ctrl->A.azimuth;
872 project_flat_setup (Ctrl->C.y, Ctrl->C.x, Ctrl->E.y, Ctrl->E.x, Ctrl->T.y, Ctrl->T.x, &theta, e, Ctrl->E.active, Ctrl->T.active, Ctrl->A.active);
873 /* Azimuth (theta) is now Cartesian in degrees */
874 if (Ctrl->L.constrain) {
875 Ctrl->L.min = 0.0;
876 xx = Ctrl->E.x - Ctrl->C.x;
877 yy = Ctrl->E.y - Ctrl->C.y;
878 Ctrl->L.max = hypot (xx, yy);
879 if (Ctrl->Q.active) Ctrl->L.max *= GMT->current.proj.DIST_KM_PR_DEG;
880 }
881 }
882 else { /* Spherical Earth mode */
883 if (Ctrl->T.active) { /* Gave the pole */
884 sin_lat_to_pole = project_oblique_setup (GMT, Ctrl->T.y, Ctrl->T.x, P.pole, &Ctrl->C.y, &Ctrl->C.x, center, Ctrl->C.active, Ctrl->G.active);
885 gmt_cart_to_geo (GMT, &P.plat, &P.plon, x, true); /* Save lon, lat of the pole */
886 if (Ctrl->G.through_C) { /* Must compute the colatitude from P to C and use that */
887 gmt_geo_to_cart (GMT, Ctrl->C.y, Ctrl->C.x, x, true); /* Cartesian vector to C */
888 gmt_geo_to_cart (GMT, Ctrl->T.y, Ctrl->T.x, xt, true); /* Cartesian vector to T */
889 Ctrl->G.colat = d_acosd (gmt_dot3v (GMT, x, xt)); /* Compute the angle between C & T which equals the small circle's colatitude */
890 GMT_Report (API, GMT_MSG_INFORMATION, "Colatitude of point C w.r.t. pole P is %g.\n", Ctrl->G.colat);
891 }
892 }
893 else { /* Using -C, -E or -A */
894 double s_hi, s_lo, s_mid, radius, m[3], ap[3], bp[3];
895 int done, n_iter = 0;
896
897 project_sphere_setup (GMT, Ctrl->C.y, Ctrl->C.x, a, Ctrl->E.y, Ctrl->E.x, b, Ctrl->A.azimuth, P.pole, center, Ctrl->E.active);
898 gmt_cart_to_geo (GMT, &P.plat, &P.plon, x, true); /* Save lon, lat of the pole */
899 radius = 0.5 * d_acosd (gmt_dot3v (GMT, a, b));
900 if (radius > fabs (Ctrl->G.colat)) {
901 GMT_Report (API, GMT_MSG_ERROR, "Center [-C] and end point [-E] are too far apart (%g) to define a small-circle with colatitude %g. Revert to great-circle.\n", radius, Ctrl->G.colat);
902 Ctrl->G.mode = 0;
903 }
904 else if (doubleAlmostEqual (Ctrl->G.colat, 90.0)) { /* Great circle pole needed */
905 if (Ctrl->L.constrain) Ctrl->L.max = d_acosd (gmt_dot3v (GMT, a, b));
906 }
907 else { /* Find small circle pole so C and E are |lat| degrees from it. */
908 for (col = 0; col < 3; col++) m[col] = a[col] + b[col]; /* Mid point along A-B */
909 gmt_normalize3v (GMT, m);
910 sign = copysign (1.0, Ctrl->G.colat);
911 s_hi = 90.0; s_lo = 0.0;
912 done = false;
913 do { /* Trial for finding pole S */
914 n_iter++;
915 s_mid = 0.5 * (s_lo + s_hi);
916 sincosd (sign * s_mid, &s, &c);
917 for (col = 0; col < 3; col++) x[col] = P.pole[col] * s + m[col] * c;
918 gmt_normalize3v (GMT, x);
919 radius = d_acosd (gmt_dot3v (GMT, a, x));
920 if (fabs (radius - fabs (Ctrl->G.colat)) < GMT_CONV8_LIMIT)
921 done = true;
922 else if (radius > fabs (Ctrl->G.colat))
923 s_hi = s_mid;
924 else
925 s_lo = s_mid;
926 if (n_iter > 500) done = true; /* Safety valve */
927 } while (!done);
928 gmt_cart_to_geo (GMT, &P.plat, &P.plon, x, true); /* Save lon, lat of the new pole */
929 GMT_Report (API, GMT_MSG_INFORMATION, "Pole for small circle located at %g %g\n", radius, P.plon, P.plat);
930 gmt_M_memcpy (P.pole, x, 3, double); /* Replace great circle pole with small circle pole */
931 sin_lat_to_pole = s;
932 gmt_cross3v (GMT, P.pole, a, x);
933 gmt_normalize3v (GMT, x);
934 gmt_cross3v (GMT, x, P.pole, ap);
935 gmt_normalize3v (GMT, ap);
936 gmt_cross3v (GMT, P.pole, b, x);
937 gmt_normalize3v (GMT, x);
938 gmt_cross3v (GMT, x, P.pole, bp);
939 gmt_normalize3v (GMT, bp);
940 if (Ctrl->L.constrain) Ctrl->L.max = d_acosd (gmt_dot3v (GMT, ap, bp)) * sin_lat_to_pole;
941 }
942 }
943 if (Ctrl->L.constrain) {
944 Ctrl->L.min = 0.0;
945 if (Ctrl->Q.active) Ctrl->L.max *= GMT->current.proj.DIST_KM_PR_DEG;
946 }
947 }
948 if (Ctrl->G.number) Ctrl->G.inc = (Ctrl->L.max - Ctrl->L.min) / (Ctrl->G.inc - 1.0);
949
950 /* Now things are initialized. We will work in degrees for awhile, so we convert things */
951
952 if (Ctrl->Q.active) {
953 Ctrl->G.inc /= GMT->current.proj.DIST_KM_PR_DEG;
954 Ctrl->L.min /= GMT->current.proj.DIST_KM_PR_DEG;
955 Ctrl->L.max /= GMT->current.proj.DIST_KM_PR_DEG;
956 Ctrl->W.min /= GMT->current.proj.DIST_KM_PR_DEG;
957 Ctrl->W.max /= GMT->current.proj.DIST_KM_PR_DEG;
958 }
959
960 /* Now we are ready to work */
961
962 P.n_used = n_total_read = 0;
963
964 if (Ctrl->G.active) { /* No input data expected, just generate x,y,d track from arguments given */
965 double out[3];
966 struct GMT_RECORD *Out = gmt_new_record (GMT, out, NULL);
967 char *z_header = NULL;
968 P.output_choice[0] = 4;
969 P.output_choice[1] = 5;
970 P.output_choice[2] = 2;
971 if (Ctrl->Z.active) { /* Do the full ellipse */
972 double h = pow (Ctrl->Z.major - Ctrl->Z.minor, 2.0) / pow (Ctrl->Z.major + Ctrl->Z.minor, 2.0);
973 Ctrl->L.min = 0.0;
974 Ctrl->L.max = M_PI * (Ctrl->Z.major + Ctrl->Z.minor) * (1.0 + (3.0 * h)/(10.0 + sqrt (4.0 - 3.0 * h))); /* Ramanujan approximation of ellipse circumference */
975 if (gmt_M_is_geographic (GMT, GMT_IN)) Ctrl->L.max /= GMT->current.proj.DIST_KM_PR_DEG;
976 if (Ctrl->Z.number) /* Want a specific number of points */
977 Ctrl->G.inc = Ctrl->L.max / rint (Ctrl->G.inc);
978 else if (Ctrl->Z.exact) { /* Adjust inc to fit the ellipse perimeter exactly */
979 double f = rint (Ctrl->L.max / Ctrl->G.inc);
980 Ctrl->G.inc = Ctrl->L.max / f;
981 }
982 }
983
984 GMT_Report (API, GMT_MSG_INFORMATION, "Generate table data\n");
985 GMT_Report (API, GMT_MSG_INFORMATION, "Go from min dist = %g to max dist = %g\n", Ctrl->L.min, Ctrl->L.max);
986 d_along = Ctrl->L.min;
987 while ((Ctrl->L.max - d_along) > (GMT_CONV8_LIMIT*Ctrl->G.inc)) {
988 p_data[P.n_used].a[2] = d_along;
989 p_data[P.n_used].t = NULL; /* Initialize since that is not done by realloc */
990 p_data[P.n_used].z = NULL; /* Initialize since that is not done by realloc */
991 P.n_used++;
992 d_along = Ctrl->L.min + P.n_used * Ctrl->G.inc;
993 if (P.n_used == (n_alloc-1)) {
994 size_t old_n_alloc = n_alloc;
995 n_alloc <<= 1;
996 p_data = gmt_M_memory (GMT, p_data, n_alloc, struct PROJECT_DATA);
997 gmt_M_memset (&(p_data[old_n_alloc]), n_alloc - old_n_alloc, struct PROJECT_DATA); /* Set to NULL/0 */
998 }
999 }
1000 p_data[P.n_used].a[2] = Ctrl->L.max;
1001 p_data[P.n_used].t = NULL; /* Initialize since that is not done by realloc */
1002 p_data[P.n_used].z = NULL; /* Initialize since that is not done by realloc */
1003 P.n_used ++;
1004
1005 /* We need to find r,s */
1006
1007 if (Ctrl->Z.active) {
1008 uint64_t ne = P.n_used - 1;
1009 if (Ctrl->N.active) { /* Cartesian ellipse */
1010 double r, ca, sa, d_azim = TWO_PI / ne, e2 = 1.0 - pow (Ctrl->Z.minor/Ctrl->Z.major, 2.0);
1011 sincosd (Ctrl->Z.azimuth - 90.0, &sin_theta, &cos_theta);
1012 e[0] = e[3] = cos_theta;
1013 e[1] = sin_theta;
1014 e[2] = -e[1];
1015 for (rec = 0; rec < P.n_used; rec++) {
1016 sincos (d_azim*rec, &sa, &ca);
1017 r = Ctrl->Z.minor / sqrt (1.0 - e2 * pow (ca, 2.0));
1018 x[0] = r * ca; x[1] = r * sa;
1019 project_matrix_2v (e, x, xt);
1020 p_data[rec].a[4] = Ctrl->C.x + xt[GMT_X];
1021 p_data[rec].a[5] = Ctrl->C.y + xt[GMT_Y];
1022 }
1023 z_header = strdup ("Testing Cartesian ellipse");
1024 }
1025 else { /* Geographic ellipse */
1026 struct GMT_DATASEGMENT *S = gmt_get_geo_ellipse (GMT, Ctrl->C.x, Ctrl->C.y, Ctrl->Z.major, Ctrl->Z.minor, Ctrl->Z.azimuth, ne);
1027 for (rec = 0; rec < P.n_used; rec++) {
1028 p_data[rec].a[4] = S->data[GMT_X][rec];
1029 p_data[rec].a[5] = S->data[GMT_Y][rec];
1030 }
1031 z_header = strdup (S->header);
1032 gmt_free_segment (GMT, &S);
1033 }
1034 }
1035 else if (Ctrl->N.active) {
1036 sincosd (theta, &sin_theta, &cos_theta);
1037 for (rec = 0; rec < P.n_used; rec++) {
1038 p_data[rec].a[4] = Ctrl->C.x + p_data[rec].a[2] * cos_theta;
1039 p_data[rec].a[5] = Ctrl->C.y + p_data[rec].a[2] * sin_theta;
1040 }
1041 }
1042 else {
1043 /* Must set generating vector to point along zero-meridian so it is the desired number of degrees [Ctrl->G.colat]
1044 * from the pole. */
1045 double C[3], N[3], counteract = 1.0;
1046 gmt_geo_to_cart (GMT, Ctrl->C.y, Ctrl->C.x, C, true); /* User origin C */
1047 gmt_cross3v (GMT, P.pole, C, N); /* This is vector normal to meridian plan */
1048 gmt_normalize3v (GMT, N); /* Make it a unit vector */
1049 project_make_euler_matrix (N, e, Ctrl->G.colat); /* Rotation matrix about N */
1050 project_matrix_3v (e, P.pole, x); /* This is the generating vector for our circle */
1051 if (Ctrl->L.constrain) counteract = 1.0 / sin_lat_to_pole; /* Increase angle to counteract effect of small circle settings */
1052 for (rec = 0; rec < P.n_used; rec++) {
1053 project_make_euler_matrix (P.pole, e, sign * p_data[rec].a[2] * counteract);
1054 project_matrix_3v (e,x,xt);
1055 gmt_cart_to_geo (GMT, &(p_data[rec].a[5]), &(p_data[rec].a[4]), xt, true);
1056 }
1057 }
1058
1059 /* At this stage, all values are still in degrees. */
1060
1061 if (Ctrl->Q.active) {
1062 for (rec = 0; rec < P.n_used; rec++) {
1063 p_data[rec].a[2] *= GMT->current.proj.DIST_KM_PR_DEG;
1064 p_data[rec].a[3] *= GMT->current.proj.DIST_KM_PR_DEG;
1065 }
1066 }
1067
1068 /* Now output generated track */
1069
1070 if ((error = GMT_Set_Columns (API, GMT_OUT, (unsigned int)P.n_outputs, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
1071 gmt_M_free (GMT, p_data);
1072 gmt_M_free (GMT, Out);
1073 if (Ctrl->Z.active) gmt_M_str_free (z_header);
1074 Return (error);
1075 }
1076 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers data output failed */
1077 for (rec = 0; rec < P.n_used; rec++) {
1078 gmt_M_str_free (p_data[rec].t); gmt_M_free (GMT, p_data[rec].z);
1079 }
1080 gmt_M_free (GMT, p_data);
1081 gmt_M_free (GMT, Out);
1082 if (Ctrl->Z.active) gmt_M_str_free (z_header);
1083 Return (API->error);
1084 }
1085 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Failed to enable data output and set access mode */
1086 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_LINE) != GMT_NOERROR) { /* Sets output geometry */
1087 gmt_M_free (GMT, p_data);
1088 gmt_M_free (GMT, Out);
1089 if (Ctrl->Z.active) gmt_M_str_free (z_header);
1090 Return (API->error);
1091 }
1092 for (rec = 0; rec < P.n_used; rec++) {
1093 gmt_M_str_free (p_data[rec].t); gmt_M_free (GMT, p_data[rec].z);
1094 }
1095 gmt_M_free (GMT, p_data);
1096 gmt_M_free (GMT, Out);
1097 if (Ctrl->Z.active) gmt_M_str_free (z_header);
1098 Return (API->error);
1099 }
1100
1101 if (!GMT->common.b.active[GMT_OUT] && Ctrl->G.header) { /* Want segment header on output */
1102 int kind = (doubleAlmostEqualZero (Ctrl->G.colat, 90.0)) ? 0 : 1;
1103 char *type[2] = {"Great", "Small"};
1104 gmt_set_segmentheader (GMT, GMT_OUT, true); /* Turn on segment headers on output */
1105 if (Ctrl->Z.active)
1106 sprintf (GMT->current.io.segment_header, "%s", z_header);
1107 else
1108 sprintf (GMT->current.io.segment_header, "%s-circle Pole at %g %g", type[kind], P.plon, P.plat);
1109 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); /* Write segment header */
1110 }
1111 if (Ctrl->Z.active) {
1112 gmt_M_str_free (z_header);
1113 if (Ctrl->Z.number) P.n_used--; /* To avoid repeated point if +n is used */
1114 }
1115 for (rec = 0; rec < P.n_used; rec++) {
1116 for (col = 0; col < P.n_outputs; col++) out[col] = p_data[rec].a[P.output_choice[col]];
1117 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
1118 }
1119 gmt_M_free (GMT, Out);
1120 }
1121 else { /* Must read input file */
1122 struct GMT_RECORD *In = NULL;
1123
1124 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");
1125 /* Specify input and output expected columns */
1126 if ((error = GMT_Set_Columns (API, GMT_IN, 0, GMT_COL_FIX)) != GMT_NOERROR) {
1127 Return (error);
1128 }
1129
1130 /* Initialize the i/o since we are doing record-by-record reading/writing */
1131 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Establishes data input */
1132 Return (API->error);
1133 }
1134 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
1135 Return (API->error);
1136 }
1137
1138 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers data output */
1139 Return (API->error);
1140 }
1141 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
1142 Return (API->error);
1143 }
1144 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
1145 Return (API->error);
1146 }
1147
1148 do { /* Keep returning records until we reach EOF */
1149 if ((In = GMT_Get_Record (API, GMT_READ_MIXED, NULL)) == NULL) { /* Read next record, get NULL if special case */
1150 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
1151 Return (GMT_RUNTIME_ERROR);
1152 }
1153 else if (gmt_M_rec_is_table_header (GMT)) /* Echo table headers */
1154 GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, NULL);
1155 else if (gmt_M_rec_is_segment_header (GMT)) { /* Echo segment headers */
1156 if (P.n_used) { /* Write out previous segment */
1157 if ((error = project_write_one_segment (GMT, Ctrl, theta, p_data, &P)) != 0) Return (error);
1158 n_total_used += P.n_used;
1159 P.n_used = 0;
1160 }
1161 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
1162 }
1163 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
1164 break;
1165 continue; /* Go back and read the next record */
1166 }
1167 if (In->data == NULL) {
1168 gmt_quit_bad_record (API, In);
1169 Return (API->error);
1170 }
1171
1172 /* Data record to process */
1173 in = In->data; /* Only need to process numerical part here */
1174
1175 if (z_first) {
1176 uint64_t n_cols = gmt_get_cols (GMT, GMT_IN), n_tot_cols;
1177 if (n_cols == 2 && P.want_z_output && In->text == NULL) {
1178 if (z_set_auto) { /* Implicitly set -Fxyzpqrs earlier but input file has no z values so roll back to -Fxypqrs */
1179 P.want_z_output = z_set_auto = false;
1180 for (col = 3; col < P.n_outputs; col++) P.output_choice[col-1] = P.output_choice[col]; /* Shuffle pqrs to the left */
1181 P.n_outputs--;
1182 if (!Ctrl->N.active) {
1183 for (col = 0; col < P.n_outputs; col++) {
1184 switch (P.output_choice[col]) {
1185 case 0: case 4: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_LON); break;
1186 case 1: case 5: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_LAT); break;
1187 default: gmt_set_column_type (GMT, GMT_OUT, (unsigned int)col, GMT_IS_FLOAT); break;
1188 }
1189 }
1190 }
1191 }
1192 else {
1193 GMT_Report (API, GMT_MSG_ERROR, "No data columns or trailing text after leading coordinates, cannot use z flag in -F\n");
1194 Return (GMT_RUNTIME_ERROR);
1195 }
1196 }
1197 else if (n_cols < 2) {
1198 GMT_Report (API, GMT_MSG_ERROR, "Input file must at least have x,y or lon,lat columns\n");
1199 gmt_M_free (GMT, p_data);
1200 Return (GMT_RUNTIME_ERROR);
1201 }
1202 else if (!Ctrl->N.active) { /* Must update col types since # of z values will need to be considered */
1203 unsigned int k, kk;
1204 P.n_z = n_cols - 2;
1205 for (col = kk = 0; col < P.n_outputs; col++, kk++) {
1206 switch (P.output_choice[col]) {
1207 case -1: /* Need to set float type for all the nz columns based on what the input z column types are */
1208 for (k = 0; k < P.n_z; k++, kk++) gmt_set_column_type (GMT, GMT_OUT, kk, GMT->current.io.col_type[GMT_IN][k+GMT_Z]);
1209 kk--; /* Since we also do this at end of loop */
1210 break;
1211 case 0: case 4: gmt_set_column_type (GMT, GMT_OUT, kk, GMT_IS_LON); break;
1212 case 1: case 5: gmt_set_column_type (GMT, GMT_OUT, kk, GMT_IS_LAT); break;
1213 default: gmt_set_column_type (GMT, GMT_OUT, kk, GMT_IS_FLOAT); break;
1214 }
1215 }
1216 }
1217 else /* Just need to know how many */
1218 P.n_z = n_cols - 2;
1219 z_first = false;
1220 n_tot_cols = (P.want_z_output) ? P.n_outputs - 1 + P.n_z : P.n_outputs;
1221 if ((error = GMT_Set_Columns (API, GMT_OUT, (unsigned int)n_tot_cols, gmt_M_colmode (In->text))) != GMT_NOERROR) {
1222 Return (error);
1223 }
1224 }
1225
1226 xx = in[GMT_X];
1227 yy = in[GMT_Y];
1228
1229 n_total_read ++;
1230
1231 if (Ctrl->N.active) {
1232 x[0] = xx - Ctrl->C.x;
1233 x[1] = yy - Ctrl->C.y;
1234 project_matrix_2v (e, x, xt);
1235 }
1236 else
1237 project_oblique_transform (GMT, yy, xx, &xt[1], &xt[0], P.pole, center);
1238
1239 skip = ((Ctrl->L.active && (xt[0] < Ctrl->L.min || xt[0] > Ctrl->L.max)) || (Ctrl->W.active && (xt[1] < Ctrl->W.min || xt[1] > Ctrl->W.max)));
1240
1241 if (skip) continue;
1242
1243 p_data[P.n_used].a[0] = xx;
1244 p_data[P.n_used].a[1] = yy;
1245 p_data[P.n_used].a[2] = xt[0];
1246 p_data[P.n_used].a[3] = xt[1];
1247 p_data[P.n_used].t = NULL; /* Initialize since that is not done by realloc */
1248 p_data[P.n_used].z = NULL; /* Initialize since that is not done by realloc */
1249 if (P.want_z_output && In->text) /* Must store all trailing text */
1250 p_data[P.n_used].t = strdup (In->text);
1251 if (P.n_z) { /* Copy over z column(s) */
1252 p_data[P.n_used].z = gmt_M_memory (GMT, NULL, P.n_z, double);
1253 gmt_M_memcpy (p_data[P.n_used].z, &in[GMT_Z], P.n_z, double);
1254 }
1255 P.n_used++;
1256 if (P.n_used == n_alloc) {
1257 size_t old_n_alloc = n_alloc;
1258 n_alloc <<= 1;
1259 p_data = gmt_M_memory (GMT, p_data, n_alloc, struct PROJECT_DATA);
1260 gmt_M_memset (&(p_data[old_n_alloc]), n_alloc - old_n_alloc, struct PROJECT_DATA); /* Set to NULL/0 */
1261 }
1262 } while (true);
1263
1264 if (P.n_used < n_alloc) p_data = gmt_M_memory (GMT, p_data, P.n_used, struct PROJECT_DATA);
1265
1266 if (P.n_used) { /* Finish last segment output */
1267 if ((error = project_write_one_segment (GMT, Ctrl, theta, p_data, &P)) != 0) Return (error);
1268 n_total_used += P.n_used;
1269 }
1270
1271 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
1272 Return (API->error);
1273 }
1274 }
1275
1276 for (rec = 0; rec < P.n_used; rec++) {
1277 gmt_M_str_free (p_data[rec].t);
1278 gmt_M_free (GMT, p_data[rec].z);
1279 }
1280
1281 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
1282 Return (API->error);
1283 }
1284
1285 GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " read, %" PRIu64 " used\n", n_total_read, n_total_used);
1286
1287 gmt_M_free (GMT, p_data);
1288
1289 Return (GMT_NOERROR);
1290 }
1291