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: grd2cpt reads a 2d binary gridded grdfile and creates a continuous-color-
19 * palette CPT, with a non-linear histogram-equalized mapping between
20 * hue and data value. (The linear mapping can be made with grd2cpt.)
21 *
22 * Creates a cumulative distribution function f(z) describing the data
23 * in the grdfile. f(z) is sampled at z values supplied by the user
24 * [with -S option] or guessed from the sample mean and standard deviation.
25 * f(z) is then found by looping over the grd array for each z and counting
26 * data values <= z. Once f(z) is found then a master CPT is resampled
27 * based on a normalized f(z).
28 *
29 * Author: Walter H. F. Smith
30 * Date: 1-JAN-2010
31 * Version: 6 API
32 *
33 */
34
35 #include "gmt_dev.h"
36
37 #define THIS_MODULE_CLASSIC_NAME "grd2cpt"
38 #define THIS_MODULE_MODERN_NAME "grd2cpt"
39 #define THIS_MODULE_LIB "core"
40 #define THIS_MODULE_PURPOSE "Make linear or histogram-equalized color palette table from grid"
41 #define THIS_MODULE_KEYS "<G{+,>C},ED)=f"
42 #define THIS_MODULE_NEEDS ""
43 #define THIS_MODULE_OPTIONS "->RVhbo"
44
45 #define GRD2CPT_N_LEVELS 11 /* The default number of levels if nothing is specified */
46
47 struct GRD2CPT_CTRL {
48 struct GRD2CPT_In {
49 bool active;
50 } In;
51 struct GRD2CPT_Out { /* -> */
52 bool active;
53 char *file;
54 } Out;
55 struct GRD2CPT_A { /* -A<transp>[+a] */
56 bool active;
57 unsigned int mode;
58 double value;
59 } A;
60 struct GRD2CPT_C { /* -C<cpt> or -C<color1>,<color2>[,<color3>,...] */
61 bool active;
62 char *file;
63 } C;
64 struct GRD2CPT_D { /* -D[i|o] */
65 bool active;
66 unsigned int mode;
67 } D;
68 struct GRD2CPT_E { /* -E<nlevels>[+c][+f<file>] */
69 bool active;
70 bool cdf;
71 char *file;
72 unsigned int levels;
73 } E;
74 struct GRD2CPT_F { /* -F[r|R|h|c][+c[<label>]] */
75 bool active;
76 bool cat;
77 unsigned int model;
78 char *label;
79 } F;
80 struct GRD2CPT_G { /* -Glow/high for input CPT truncation */
81 bool active;
82 double z_low, z_high;
83 } G;
84 struct GRD2CPT_H { /* -H */
85 bool active;
86 } H;
87 struct GRD2CPT_I { /* -I[z][c] */
88 bool active;
89 unsigned int mode;
90 } I;
91 struct GRD2CPT_L { /* -L<min_limit>/<max_limit> */
92 bool active;
93 bool minimum_given, maximum_given;
94 double min, max;
95 } L;
96 struct GRD2CPT_M { /* -M */
97 bool active;
98 } M;
99 struct GRD2CPT_N { /* -N */
100 bool active;
101 } N;
102 struct GRD2CPT_Q { /* -Q[i|o] */
103 bool active;
104 unsigned int mode;
105 } Q;
106 struct GRD2CPT_T { /* -T<start>/<stop>/<inc> or -T<n_levels> */
107 bool active;
108 bool interpolate;
109 unsigned int mode; /* 0 or 1 (-Tn) */
110 unsigned int n_levels;
111 double low, high, inc;
112 char *file;
113 } T;
114 struct GRD2CPT_S { /* -S<kind> */
115 bool active;
116 int kind; /* -1 symmetric +-zmin, +1 +-zmax, -2 = +-Minx(|zmin|,|zmax|), +2 = +-Max(|zmin|,|zmax|), 0 = min to max [Default] */
117 } S;
118 struct GRD2CPT_W { /* -W[w] */
119 bool active;
120 bool wrap;
121 } W;
122 struct GRD2CPT_Z { /* -Z */
123 bool active;
124 } Z;
125 };
126
New_Ctrl(struct GMT_CTRL * GMT)127 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
128 struct GRD2CPT_CTRL *C;
129
130 C = gmt_M_memory (GMT, NULL, 1, struct GRD2CPT_CTRL);
131
132 /* Initialize values whose defaults are not 0/false/NULL */
133 C->G.z_low = C->G.z_high = GMT->session.d_NaN; /* No truncation */
134 return (C);
135 }
136
Free_Ctrl(struct GMT_CTRL * GMT,struct GRD2CPT_CTRL * C)137 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRD2CPT_CTRL *C) { /* Deallocate control structure */
138 if (!C) return;
139 gmt_M_str_free (C->Out.file);
140 gmt_M_str_free (C->C.file);
141 gmt_M_str_free (C->E.file);
142 gmt_M_str_free (C->F.label);
143 gmt_M_str_free (C->T.file);
144 gmt_M_free (GMT, C);
145 }
146
usage(struct GMTAPI_CTRL * API,int level)147 static int usage (struct GMTAPI_CTRL *API, int level) {
148 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
149 const char *H_OPT = (API->GMT->current.setting.run_mode == GMT_MODERN) ? " [-H]" : "";
150 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
151 GMT_Usage (API, 0, "usage: %s %s [-A<transparency>[+a]] [-C<cpt>] [-D[i|o]] [-E[<nlevels>][+c][+f<file>]] "
152 "[-F[R|r|h|c][+c[<label>]]] [-G<zlo>/<zhi>]%s [-I[c][z]] [-L<min_limit>/<max_limit>] [-M] [-N] [-Q[i|o]] "
153 "[%s] [-Sh|l|m|u] [-T<start>/<stop>/<inc>|<n>] [%s] [-W[w]] [-Z] [%s] [%s] [%s] [%s]\n",
154 name, GMT_INGRID, H_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bo_OPT, GMT_ho_OPT, GMT_o_OPT, GMT_PAR_OPT);
155
156 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
157
158 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
159 gmt_ingrid_syntax (API, 0, "Name of one or more grid files");
160 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
161 GMT_Usage (API, 1, "\n-A<transparency>[+a]");
162 GMT_Usage (API, -2, "Set constant transparency for all colors; append +a to also include back-, for-, and nan-colors [0].");
163 if (gmt_list_cpt (API->GMT, 'C')) return (GMT_CPT_READ_ERROR); /* Display list of available color tables */
164 GMT_Usage (API, 1, "\n-D[i|o]");
165 GMT_Usage (API, -2, "Set back- and foreground color to match the bottom/top limits "
166 "in the output CPT [Default (-D or -Do) uses color output table]. Append i "
167 "to match the bottom/top values in the input CPT instead.");
168 GMT_Usage (API, 1, "\n-E[<nlevels>][+c][+f<file>]");
169 GMT_Usage (API, -2, "Set CPT to go from grid zmin to zmax (i.e., a linear scale). "
170 "Alternatively, append <nlevels> to sample equidistant color levels from zmin to zmax. "
171 "Instead, append +c to use the cumulative density function (cdf) to set color bounds. "
172 "Append +f<file> to save the CDF table to a file.");
173 GMT_Usage (API, 1, "\n-F[R|r|h|c][+c[<label>]]");
174 GMT_Usage (API, -2, "Select the color model for output [Default uses the input model]:");
175 GMT_Usage (API, 3, "R: Output r/g/b or grayscale or colorname.");
176 GMT_Usage (API, 3, "r: Output r/g/b only.");
177 GMT_Usage (API, 3, "h: Output h-s-v.");
178 GMT_Usage (API, 3, "c: Output c/m/y/k.");
179 GMT_Usage (API, -2, "One modifier controls generation of categorical labels:");
180 GMT_Usage (API, 3, "+c Output a discrete CPT in categorical CPT format. "
181 "The <label>, if appended, sets the labels for each category. It may be a "
182 "comma-separated list of category names, or <start>[-] where we automatically build "
183 "labels from <start> (a letter or an integer). Append - to build range labels <start>-<start+1>.");
184 GMT_Usage (API, 1, "\n-G<zlo>/<zhi>");
185 GMT_Usage (API, -2, "Truncate incoming CPT to be limited to the z-range <zlo>/<zhi>. "
186 "To accept one of the incoming limits, set that limit to NaN.");
187 GMT_Usage (API, 1, "\n-H Modern mode only: Also write CPT to standard output [Default just saves as current CPT].");
188 GMT_Usage (API, 1, "\n-I[c][z]");
189 GMT_Usage (API, -2, "Invert sense of CPT in one or two ways:");
190 GMT_Usage (API, 3, "c: Invert sense of color table as well as back- and foreground color [Default].");
191 GMT_Usage (API, 3, "z: Invert sign of z-values in the color table (takes affect before -G, T are consulted).");
192 GMT_Usage (API, 1, "\n-L<min_limit>/<max_limit>");
193 GMT_Usage (API, -2, "Limit the range of the data. Node values outside this range are set to NaN. "
194 "To only give min or max limit, set the other to - [Default uses actual min,max of data].");
195 GMT_Usage (API, 1, "\n-M Use GMT defaults to set back-, foreground, and NaN colors [Default uses color table].");
196 GMT_Usage (API, 1, "\n-N Do not write back-, foreground, and NaN colors [Default will].");
197 GMT_Usage (API, 1, "\n-Q[i|o]");
198 GMT_Usage (API, -2, "Assign a logarithmic colortable [Default is linear]. Append mode:");
199 GMT_Usage (API, 3, "i: z-values are actually log10(z). Assign colors and write z [Default].");
200 GMT_Usage (API, 3, "o: z-values are z, but take log10(z), assign colors and write z.");
201 GMT_Option (API, "R");
202 GMT_Usage (API, 1, "\n-Sh|l|m|u");
203 GMT_Usage (API, -2, "Force color tables to be symmetric about 0. Append one modifier:");
204 GMT_Usage (API, 3, "l: (low) for values symmetric about zero from -|zmin| to +|zmin|.");
205 GMT_Usage (API, 3, "u: (upper) for values symmetric about zero from -|zmax| to +|zmax|.");
206 GMT_Usage (API, 3, "m: (min) for values symmetric about zero -+min(|zmin|,|zmax|).");
207 GMT_Usage (API, 3, "h: (high) for values symmetric about zero -+max(|zmin|,|zmax|).");
208 GMT_Usage (API, 1, "\n-T<start>/<stop>/<inc>|<n>");
209 GMT_Usage (API, -2, "CDF sample points should range from <start> to <stop> by <inc>. "
210 "Alternatively, use -T<n> to select <n> points from a cumulative normal distribution [%d]. "
211 "Here, <start> maps to data min and <stop> maps to data max (but see -L) "
212 "[Default uses equidistant steps for a Gaussian CDF].", GRD2CPT_N_LEVELS);
213 GMT_Option (API, "V");
214 GMT_Usage (API, 1, "\n-W[w]");
215 GMT_Usage (API, -2, "Do not interpolate color palette. Alternatively, append w for a wrapped CPT.");
216 GMT_Usage (API, 1, "\n-Z Force a continuous color palette [Default is discontinuous, i.e., constant color intervals].");
217 GMT_Option (API, "bo,h,o,.");
218
219 return (GMT_MODULE_USAGE);
220 }
221
parse(struct GMT_CTRL * GMT,struct GRD2CPT_CTRL * Ctrl,struct GMT_OPTION * options)222 static int parse (struct GMT_CTRL *GMT, struct GRD2CPT_CTRL *Ctrl, struct GMT_OPTION *options) {
223 /* This parses the options provided to grdcut and sets parameters in CTRL.
224 * Any GMT common options will override values set previously by other commands.
225 * It also replaces any file names specified as input or output with the data ID
226 * returned when registering these sources/destinations with the API.
227 */
228
229 int n;
230 unsigned int n_errors = 0, pos = 0, n_files[2] = {0, 0};
231 char txt_a[GMT_LEN512] = {""}, txt_b[GMT_LEN32] = {""}, *c = NULL;
232 char *T_arg = NULL, *S_arg = NULL;
233 struct GMT_OPTION *opt = NULL;
234 struct GMTAPI_CTRL *API = GMT->parent;
235
236 for (opt = options; opt; opt = opt->next) {
237 switch (opt->option) {
238
239 case '<': /* Input files */
240 Ctrl->In.active = true;
241 n_files[GMT_IN]++;
242 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
243 break;
244 case '>': /* Got named output file */
245 if (n_files[GMT_OUT]++ == 0) Ctrl->Out.file = strdup (opt->arg);
246 break;
247
248 /* Processes program-specific parameters */
249
250 case 'A': /* Sets transparency [-A<transp>[+a]] */
251 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
252 Ctrl->A.active = true;
253 if (opt->arg[0] == '+') { /* Old syntax */
254 Ctrl->A.mode = 1;
255 Ctrl->A.value = 0.01 * atof (&opt->arg[1]);
256 }
257 else if ((c = strstr (opt->arg, "+a"))) {
258 Ctrl->A.mode = 1;
259 c[0] = '\0';
260 Ctrl->A.value = 0.01 * atof (opt->arg);
261 c[0] = '+';
262 }
263 else
264 Ctrl->A.value = 0.01 * atof (opt->arg);
265 break;
266 case 'C': /* Get CPT */
267 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
268 Ctrl->C.active = true;
269 if (opt->arg[0]) Ctrl->C.file = strdup (opt->arg);
270 break;
271 case 'D': /* Set fore/back-ground to match end-colors */
272 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
273 Ctrl->D.active = true;
274 Ctrl->D.mode = 1;
275 if (opt->arg[0] == 'i') Ctrl->D.mode = 2;
276 break;
277 case 'E': /* Use n levels */
278 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
279 Ctrl->E.active = true;
280 if (opt->arg[0]) { /* Got an argument */
281 if (gmt_validate_modifiers (GMT, opt->arg, 'E', "cf", GMT_MSG_ERROR)) n_errors++;
282 if ((c = gmt_first_modifier (GMT, opt->arg, "cf")) ) {
283 while (gmt_getmodopt (GMT, 'E', c, "cf", &pos, txt_a, &n_errors) && n_errors == 0) {
284 switch (txt_a[0]) {
285 case 'c': Ctrl->E.cdf = true; break; /* Determine Cumulative Density Function */
286 case 'f':
287 if (txt_a[1])
288 Ctrl->E.file = strdup (&txt_a[1]);
289 else {
290 GMT_Report (API, GMT_MSG_ERROR, "Option -E: No filename given via +f\n");
291 n_errors++;
292 }
293 break; /* Incremental distance */
294 default: break; /* These are caught in gmt_getmodopt so break is just for Coverity */
295 }
296 }
297 c[0] = '\0'; /* Chop off the modifiers */
298 }
299 if (sscanf (opt->arg, "%d", &Ctrl->E.levels) != 1) {
300 GMT_Report (API, GMT_MSG_ERROR, "Option -E: Cannot decode value given\n");
301 n_errors++;
302 }
303 if (c) c[0] = '+'; /* Restore modifiers */
304 }
305 break;
306 case 'F': /* Set color model for output */
307 n_errors += gmt_M_repeated_module_option (API, Ctrl->F.active);
308 if (gmt_validate_modifiers (GMT, opt->arg, 'F', "c", GMT_MSG_ERROR)) n_errors++;
309 if (gmt_get_modifier (opt->arg, 'c', txt_a)) {
310 Ctrl->F.cat = true;
311 if (txt_a[0]) Ctrl->F.label = strdup (txt_a);
312 }
313 Ctrl->F.active = true;
314 switch (txt_a[0]) {
315 case 'r': Ctrl->F.model = GMT_RGB + GMT_NO_COLORNAMES; break;
316 case 'h': Ctrl->F.model = GMT_HSV; break;
317 case 'c': Ctrl->F.model = GMT_CMYK; break;
318 default: Ctrl->F.model = GMT_RGB; break;
319 }
320 break;
321 case 'G': /* truncate incoming CPT */
322 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
323 Ctrl->G.active = true;
324 n_errors += gmt_get_limits (GMT, 'G', opt->arg, 0, &Ctrl->G.z_low, &Ctrl->G.z_high);
325 break;
326 case 'H': /* Modern mode only: write CPT to stdout */
327 n_errors += gmt_M_repeated_module_option (API, Ctrl->H.active);
328 Ctrl->H.active = true;
329 break;
330 case 'I': /* Invert table */
331 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
332 Ctrl->I.active = true;
333 if ((Ctrl->I.mode = gmt_parse_inv_cpt (GMT, opt->arg)) == UINT_MAX)
334 n_errors++;
335 break;
336 case 'L': /* Limit data range */
337 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
338 Ctrl->L.active = true;
339 if ((n = sscanf (opt->arg, "%[^/]/%s", txt_a, txt_b)) != 2) {
340 GMT_Report (API, GMT_MSG_ERROR, "Option -L: Cannot decode two limits\n");
341 n_errors++;
342 }
343 else { /* Assign limits unless give as "-" which means to skip that limit */
344 if (strcmp (txt_a, "-")) Ctrl->L.min = atof (txt_a), Ctrl->L.minimum_given = true;
345 if (strcmp (txt_b, "-")) Ctrl->L.max = atof (txt_b), Ctrl->L.maximum_given = true;
346 }
347 break;
348 case 'M': /* Override fore/back/NaN using GMT defaults */
349 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
350 Ctrl->M.active = true;
351 break;
352 case 'N': /* Do not write F/B/N colors */
353 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
354 Ctrl->N.active = true;
355 break;
356 case 'Q': /* Logarithmic data */
357 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
358 Ctrl->Q.active = true;
359 if (opt->arg[0] == 'o') /* Input data is z, but take log10(z) before interpolation colors */
360 Ctrl->Q.mode = 2;
361 else /* Input is log10(z) */
362 Ctrl->Q.mode = 1;
363 break;
364 case 'T': /* Sets sample range */
365 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
366 Ctrl->T.active = true;
367 T_arg = opt->arg;
368 break;
369 case 'S': /* Force symmetry */
370 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
371 Ctrl->S.active = true;
372 S_arg = opt->arg;
373 break;
374 case 'W': /* Do not interpolate colors */
375 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
376 if (opt->arg[0] == 'w')
377 Ctrl->W.wrap = true;
378 else
379 Ctrl->W.active = true;
380 break;
381 case 'Z': /* Continuous colors */
382 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
383 Ctrl->Z.active = true;
384 break;
385
386 default: /* Report bad options */
387 n_errors += gmt_default_error (GMT, opt->option);
388 break;
389 }
390 }
391
392 if (T_arg) { /* Must handle old or new syntax */
393 if (strchr ("-+_=", T_arg[0]) && T_arg[1] == '\0') { /* Old -Targs for symmetry given */
394 if (gmt_M_compat_check (GMT, 5)) { /* OK to parse that? */
395 GMT_Report (API, GMT_MSG_COMPAT, "Option -T-|+|_|= is deprecated; use -Sl|u|m|h instead.\n");
396 Ctrl->S.active = true;
397 if (!S_arg) Ctrl->T.active = false;
398 switch (T_arg[0]) {
399 case '-': Ctrl->S.kind = -1; break; /* Symmetric with |zmin| range */
400 case '+': Ctrl->S.kind = +1; break; /* Symmetric with |zmax| range */
401 case '_': Ctrl->S.kind = -2; break; /* Symmetric with min(|zmin|,|zmax|) range */
402 case '=': Ctrl->S.kind = +2; break; /* Symmetric with max(|zmin|,|zmax|) range */
403 default: break;
404 }
405 }
406 else {
407 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Cannot decode values %s\n", T_arg);
408 n_errors++;
409 }
410 }
411 else { /* Got correct modern args */
412 if (strchr (T_arg, '/')) { /* Gave low/high/inc */
413 if (sscanf (T_arg, "%lf/%lf/%lf", &Ctrl->T.low, &Ctrl->T.high, &Ctrl->T.inc) != 3) {
414 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Cannot decode values\n");
415 n_errors++;
416 }
417 Ctrl->T.mode = 0;
418 }
419 else if (T_arg[0]) { /* Gave -T<nlevels> */
420 Ctrl->T.n_levels = atoi (T_arg);
421 Ctrl->T.mode = 1;
422 }
423 }
424 }
425 if (S_arg) { /* Must handle old or new syntax */
426 if (strchr ("hlmu", S_arg[0]) && S_arg[1] == '\0') { /* New -S for symmetry */
427 switch (S_arg[0]) {
428 case 'l': Ctrl->S.kind = -1; break; /* Symmetric with |zmin| range */
429 case 'u': Ctrl->S.kind = +1; break; /* Symmetric with |zmax| range */
430 case 'm': Ctrl->S.kind = -2; break; /* Symmetric with min(|zmin|,|zmax|) range */
431 case 'h': Ctrl->S.kind = +2; break; /* Symmetric with max(|zmin|,|zmax|) range */
432 default: break;
433 }
434 }
435 else if (gmt_M_compat_check (GMT, 5)) { /* Old-style -S range */
436 GMT_Report (API, GMT_MSG_COMPAT, "Option -S<start>/<stop>/<inc> or -S<n> is deprecated; use -T instead.\n");
437 if (strchr (S_arg, '/')) { /* Gave low/high/inc */
438 if (sscanf (S_arg, "%lf/%lf/%lf", &Ctrl->T.low, &Ctrl->T.high, &Ctrl->T.inc) != 3) {
439 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Cannot decode values %s\n", S_arg);
440 n_errors++;
441 }
442 Ctrl->T.mode = 0;
443 }
444 else if (S_arg[0]) { /* Gave -S<nlevels> */
445 Ctrl->T.n_levels = atoi (S_arg);
446 Ctrl->T.mode = 1;
447 }
448 if (!T_arg) Ctrl->S.active = false;
449 Ctrl->T.active = true;
450 }
451 else {
452 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Cannot decode values %s\n", S_arg);
453 n_errors++;
454 }
455 }
456
457 if (Ctrl->H.active && GMT->current.setting.run_mode == GMT_CLASSIC) {
458 GMT_Report (API, GMT_MSG_WARNING, "Option -H: Only available in modern mode - ignored in classic mode\n");
459 Ctrl->H.active = false;
460 }
461 n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->C.file == NULL,
462 "Options -C: No CPT argument given\n");
463 n_errors += gmt_M_check_condition (GMT, n_files[GMT_IN] < 1, "No grid name(s) specified.\n");
464 n_errors += gmt_M_check_condition (GMT, Ctrl->W.active && Ctrl->Z.active,
465 "Options -W and -Z cannot be used simultaneously\n");
466 n_errors += gmt_M_check_condition (GMT, Ctrl->F.cat && Ctrl->Z.active,
467 "Options -F+c and -Z cannot be used simultaneously\n");
468 n_errors += gmt_M_check_condition (GMT, Ctrl->L.active && Ctrl->L.minimum_given && Ctrl->L.maximum_given && Ctrl->L.min >= Ctrl->L.max,
469 "Option -L: min_limit must be less than max_limit.\n");
470 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && Ctrl->T.mode == 0 && (Ctrl->T.high <= Ctrl->T.low || Ctrl->T.inc <= 0.0),
471 "Option -S: Bad arguments\n");
472 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && Ctrl->T.mode == 1 && Ctrl->T.n_levels == 0,
473 "Option -S: Bad arguments\n");
474 n_errors += gmt_M_check_condition (GMT, Ctrl->T.active && (Ctrl->S.active || Ctrl->E.active),
475 "Option -T: Cannot be combined with -E nor -S option.\n");
476 n_errors += gmt_M_check_condition (GMT, Ctrl->E.file && !Ctrl->E.cdf, "Option -E modifier +f is only valid if +c is used\n");
477 n_errors += gmt_M_check_condition (GMT, n_files[GMT_OUT] > 1, "Only one output destination can be specified\n");
478 n_errors += gmt_M_check_condition (GMT, Ctrl->A.active && (Ctrl->A.value < 0.0 || Ctrl->A.value > 1.0),
479 "Option -A: Transparency must be n 0-100 range [0 or opaque]\n");
480
481 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
482 }
483
grd2cpt_free_the_grids(struct GMTAPI_CTRL * API,struct GMT_GRID ** G,char ** grdfile,uint64_t n)484 GMT_LOCAL int grd2cpt_free_the_grids (struct GMTAPI_CTRL *API, struct GMT_GRID **G, char **grdfile, uint64_t n) {
485 /* Free what we are pointing to */
486 uint64_t k;
487 for (k = 0; k < n; k++) {
488 gmt_M_str_free (grdfile[k]);
489 if (GMT_Destroy_Data (API, &G[k]) != GMT_NOERROR)
490 return (API->error);
491 }
492 return (GMT_NOERROR);
493 }
494
495 #define bailout(code) {gmt_M_free_options (mode); return (code);}
496 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
497
498 EXTERN_MSC int gmtlib_compare_observation (const void *a, const void *b);
499
GMT_grd2cpt(void * V_API,int mode,void * args)500 EXTERN_MSC int GMT_grd2cpt (void *V_API, int mode, void *args) {
501 uint64_t ij, k, ngrd = 0, nxyg, nxy = 0, nfound, ngood;
502 unsigned int row, col, j, cpt_flags = 0;
503 int signed_levels, error = 0;
504 size_t n_alloc = GMT_TINY_CHUNK;
505 bool write = false, interpolate = true;
506
507 char format[GMT_BUFSIZ] = {""}, **grdfile = NULL;
508
509 double *z = NULL, wesn[4], mean, sd, wsum = 0.0, scale;
510
511 struct CDF_CPT {
512 double z; /* Data value */
513 double f; /* Cumulative distribution function f(z) */
514 } *cdf_cpt = NULL;
515
516 struct GMT_OPTION *opt = NULL;
517 struct GMT_PALETTE *Pin = NULL, *Pout = NULL;
518 struct GMT_GRID **G, *W = NULL;
519 struct GMT_OBSERVATION *pair = NULL;
520 struct GRD2CPT_CTRL *Ctrl = NULL;
521 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
522 struct GMT_OPTION *options = NULL;
523 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
524
525 /*----------------------- Standard module initialization and parsing ----------------------*/
526
527 if (API == NULL) return (GMT_NOT_A_SESSION);
528 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
529 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
530
531 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
532
533 /* Parse the command-line arguments */
534
535 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 */
536 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
537 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
538 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
539
540 /*---------------------------- This is the grd2cpt main code ----------------------------*/
541
542 if (!Ctrl->C.active) { /* No table specified; set GMT->current.setting.cpt table */
543 Ctrl->C.active = true;
544 Ctrl->C.file = strdup (GMT->current.setting.cpt);
545 }
546
547 if (!Ctrl->E.active) Ctrl->E.levels = (Ctrl->T.n_levels > 0) ? Ctrl->T.n_levels : GRD2CPT_N_LEVELS; /* Default number of levels */
548 if (Ctrl->E.cdf && Ctrl->E.levels == 0) Ctrl->E.levels = GRD2CPT_N_LEVELS;
549 if (Ctrl->M.active) cpt_flags |= GMT_CPT_NO_BNF; /* bit 0 controls if BFN is determined by parameters */
550 if (Ctrl->D.mode == 2) cpt_flags |= GMT_CPT_EXTEND_BNF; /* bit 1 controls if BF will be set to equal bottom/top rgb value */
551 if (Ctrl->Z.active) cpt_flags |= GMT_CPT_CONTINUOUS; /* Controls if final CPT should be continuous in case input is a list of colors */
552
553 if ((Pin = GMT_Read_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, cpt_flags, NULL, Ctrl->C.file, NULL)) == NULL) {
554 Return (API->error);
555 }
556 if ((API->error = gmt_validate_cpt_parameters (GMT, Pin, Ctrl->C.file, &interpolate, &(Ctrl->Z.active))))
557 Return (API->error)
558
559 if (Ctrl->I.mode & GMT_CPT_Z_REVERSE) /* Must reverse the z-values before anything else */
560 gmt_scale_cpt (GMT, Pin, -1.0);
561
562 if (Ctrl->G.active) { /* Attempt truncation */
563 struct GMT_PALETTE *Ptrunc = gmt_truncate_cpt (GMT, Pin, Ctrl->G.z_low, Ctrl->G.z_high); /* Possibly truncate the CPT */
564 if (Ptrunc == NULL)
565 Return (EXIT_FAILURE);
566 Pin = Ptrunc;
567 }
568 if (Ctrl->W.wrap) Pin->is_wrapping = true; /* A cyclic CPT has been requested */
569
570 write = (GMT->current.setting.run_mode == GMT_CLASSIC || Ctrl->H.active); /* Only output to stdout in classic mode and with -H in modern mode */
571
572 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grid(s)\n");
573
574 gmt_M_memset (wesn, 4, double);
575 if (GMT->common.R.active[RSET]) gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Subset */
576
577 G = gmt_M_memory (GMT, NULL, n_alloc, struct GMT_GRID *); /* Potentially an array of grids */
578 grdfile = gmt_M_memory (GMT, NULL, n_alloc, char *); /* Potentially an array of gridfile names */
579
580 for (opt = options, k = 0; opt; opt = opt->next) {
581 if (opt->option != '<') continue; /* We are only processing input files here */
582
583 if ((G[k] = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, wesn, opt->arg, NULL)) == NULL) {
584 error = grd2cpt_free_the_grids (API, G, grdfile, k);
585 gmt_M_free (GMT, G);
586 gmt_M_free (GMT, grdfile);
587 Return ((error) ? error : API->error);
588 }
589 grdfile[k] = strdup (opt->arg);
590 if (k && !(G[k]->header->n_columns == G[k-1]->header->n_columns && G[k]->header->n_rows == G[k-1]->header->n_rows)) {
591 GMT_Report (API, GMT_MSG_ERROR, "Grids do not have the same domain!\n");
592 error = grd2cpt_free_the_grids (API, G, grdfile, k);
593 gmt_M_free (GMT, G);
594 gmt_M_free (GMT, grdfile);
595 Return ((error) ? error : API->error);
596 }
597
598 k++;
599 if (k == n_alloc) {
600 size_t old_n_alloc = n_alloc;
601 n_alloc += GMT_TINY_CHUNK;
602 G = gmt_M_memory (GMT, G, n_alloc, struct GMT_GRID *);
603 gmt_M_memset (&(G[old_n_alloc]), n_alloc - old_n_alloc, struct GMT_GRID *); /* Set to NULL */
604 grdfile = gmt_M_memory (GMT, grdfile, n_alloc, char *);
605 gmt_M_memset (&(grdfile[old_n_alloc]), n_alloc - old_n_alloc, char *); /* Set to NULL */
606 }
607 }
608
609 ngrd = k;
610 if (ngrd < n_alloc) {
611 G = gmt_M_memory (GMT, G, ngrd, struct GMT_GRID *);
612 grdfile = gmt_M_memory (GMT, grdfile, ngrd, char *);
613 }
614
615 nxyg = G[0]->header->nm * ngrd;
616
617 if (Ctrl->E.cdf) {
618 pair = gmt_M_memory (GMT, NULL, nxyg, struct GMT_OBSERVATION);
619 W = gmt_duplicate_grid (GMT, G[0], GMT_DUPLICATE_ALLOC);
620 /* Determine the area weights */
621 gmt_get_cellarea (GMT, W);
622 }
623
624 /* Loop over the files and find NaNs. If set limits, may create more NaNs */
625 /* We use the G[0]->header to keep limits representing all the grids */
626
627 nfound = 0;
628 mean = sd = 0.0;
629 if (Ctrl->L.active) { /* Loop over the grdfiles, and set anything outside the limiting values to NaN. */
630 G[0]->header->z_min = (Ctrl->L.minimum_given) ? Ctrl->L.min : G[0]->header->z_max;
631 G[0]->header->z_max = (Ctrl->L.maximum_given) ? Ctrl->L.max : G[0]->header->z_min;
632 for (k = 0; k < ngrd; k++) { /* For each grid */
633 gmt_M_grd_loop (GMT, G[k], row, col, ij) {
634 if (gmt_M_is_fnan (G[k]->data[ij]))
635 nfound++;
636 else {
637 if (Ctrl->L.minimum_given && G[k]->data[ij] < Ctrl->L.min) {
638 nfound++;
639 G[k]->data[ij] = GMT->session.f_NaN;
640 }
641 else if (Ctrl->L.maximum_given && G[k]->data[ij] > Ctrl->L.max) {
642 nfound++;
643 G[k]->data[ij] = GMT->session.f_NaN;
644 }
645 else if (Ctrl->E.cdf) {
646 pair[nxy].value = G[k]->data[ij];
647 pair[nxy++].weight = W->data[ij];
648 wsum += W->data[ij];
649 }
650 else {
651 mean += G[k]->data[ij];
652 sd += G[k]->data[ij] * G[k]->data[ij];
653 if (!Ctrl->L.minimum_given && G[k]->data[ij] < Ctrl->L.min) Ctrl->L.min = G[k]->data[ij];
654 if (!Ctrl->L.maximum_given && G[k]->data[ij] > Ctrl->L.max) Ctrl->L.max = G[k]->data[ij];
655 }
656 }
657 }
658 }
659 if (!Ctrl->L.minimum_given) G[0]->header->z_min = Ctrl->L.min;
660 if (!Ctrl->L.maximum_given) G[0]->header->z_max = Ctrl->L.max;
661 }
662 else {
663 Ctrl->L.min = G[0]->header->z_max; /* This is just to double check G[k]->header->z_min, G[k]->header->z_max */
664 Ctrl->L.max = G[0]->header->z_min;
665 for (k = 0; k < ngrd; k++) { /* For each grid */
666 gmt_M_grd_loop (GMT, G[k], row, col, ij) {
667 if (gmt_M_is_fnan (G[k]->data[ij]))
668 nfound++;
669 else if (Ctrl->E.cdf) {
670 pair[nxy].value = G[k]->data[ij];
671 pair[nxy++].weight = W->data[ij];
672 wsum += W->data[ij];
673 }
674 else {
675 if (G[k]->data[ij] < Ctrl->L.min) Ctrl->L.min = G[k]->data[ij];
676 if (G[k]->data[ij] > Ctrl->L.max) Ctrl->L.max = G[k]->data[ij];
677 mean += G[k]->data[ij];
678 sd += G[k]->data[ij] * G[k]->data[ij];
679 }
680 }
681 }
682 G[0]->header->z_min = Ctrl->L.min;
683 G[0]->header->z_max = Ctrl->L.max;
684 }
685
686 if (Ctrl->E.cdf) {
687 gmt_free_grid (GMT, &W, true); /* Done with the area weights grid */
688 /* Sort observations on z */
689 qsort (pair, nxy, sizeof (struct GMT_OBSERVATION), gmtlib_compare_observation);
690 /* Compute normalized cumulative weights */
691 scale = 1.0 / wsum; /* Do avoid division later */
692 wsum = 0.0; /* Do this in double precision since GMT_OBSERVATION is just float */
693 for (k = 0; k < nxy; k++) { /* Build CDF from tiny to 1 */
694 pair[k].weight *= scale;
695 wsum += pair[k].weight;
696 pair[k].weight = (gmt_grdfloat)wsum;
697 }
698 if (Ctrl->E.file) { /* Save the CDF to file */
699 uint64_t dim[4] = {1, 1, nxy, 2};
700 struct GMT_DATASET *CDF = NULL;
701 struct GMT_DATASEGMENT *S = NULL;
702 if ((CDF = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
703 GMT_Report (API, GMT_MSG_ERROR, "Unable to allocate memory for CDF output.\n");
704 Return (GMT_MEMORY_ERROR);
705 }
706 S = CDF->table[0]->segment[0];
707 for (k = 0; k < nxy; k++) {
708 S->data[GMT_X][k] = pair[k].value;
709 S->data[GMT_Y][k] = pair[k].weight;
710 }
711 if (GMT_Write_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_WRITE_NORMAL, NULL, Ctrl->E.file, CDF) != GMT_NOERROR) {
712 GMT_Report (API, GMT_MSG_ERROR, "Unable to write CDF output.\n");
713 Return (API->error);
714 }
715 }
716 }
717
718 if (Ctrl->E.active && Ctrl->E.levels == 0) { /* Use existing CPT structure, just linearly change z */
719 if ((Pout = GMT_Duplicate_Data (API, GMT_IS_PALETTE, GMT_DUPLICATE_ALLOC, Pin)) == NULL) return (API->error);
720 gmt_stretch_cpt (GMT, Pout, Ctrl->L.min, Ctrl->L.max);
721 if (Ctrl->I.mode & GMT_CPT_C_REVERSE)
722 gmt_invert_cpt (GMT, Pout); /* Also flip the colors */
723 cpt_flags = 0;
724 if (Ctrl->N.active) cpt_flags |= GMT_CPT_NO_BNF; /* bit 0 controls if BFN will be written out */
725 if (Ctrl->D.mode == 1) cpt_flags |= GMT_CPT_EXTEND_BNF; /* bit 1 controls if BF will be set to equal bottom/top rgb value */
726 if (Ctrl->F.active) Pout->model = Ctrl->F.model;
727 if (Ctrl->F.cat) Pout->categorical = GMT_CPT_CATEGORICAL_VAL;
728 if (write && GMT_Write_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, cpt_flags, NULL, Ctrl->Out.file, Pout) != GMT_NOERROR) {
729 Return (API->error);
730 }
731 if (!write)
732 gmt_save_current_cpt (GMT, Pout, cpt_flags); /* Save for use by session, if modern */
733 grd2cpt_free_the_grids (API, G, grdfile, ngrd);
734 gmt_M_free (GMT, G);
735 gmt_M_free (GMT, grdfile);
736
737 Return (GMT_NOERROR);
738 }
739
740 ngood = nxyg - nfound; /* This is the number of non-NaN points for the cdf function */
741 mean /= ngood;
742 sd /= ngood;
743 sd = sqrt (sd - mean * mean);
744 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) {
745 sprintf (format, "Mean and S.D. of data are %s %s\n",
746 GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
747 GMT_Report (API, GMT_MSG_INFORMATION, format, mean, sd);
748 }
749
750 /* Decide how to make steps in z. */
751 if (Ctrl->T.active && Ctrl->T.mode == 0) { /* Use predefined levels and interval */
752 unsigned int i, j;
753
754 Ctrl->E.levels = (G[0]->header->z_min < Ctrl->T.low) ? 1 : 0;
755 Ctrl->E.levels += urint (floor((Ctrl->T.high - Ctrl->T.low)/Ctrl->T.inc)) + 1;
756 if (G[0]->header->z_max > Ctrl->T.high) Ctrl->E.levels++;
757 cdf_cpt = gmt_M_memory (GMT, NULL, Ctrl->E.levels, struct CDF_CPT);
758 if (G[0]->header->z_min < Ctrl->T.low) {
759 cdf_cpt[0].z = G[0]->header->z_min;
760 cdf_cpt[1].z = Ctrl->T.low;
761 i = 2;
762 }
763 else {
764 cdf_cpt[0].z = Ctrl->T.low;
765 i = 1;
766 }
767 j = (G[0]->header->z_max > Ctrl->T.high) ? Ctrl->E.levels - 1 : Ctrl->E.levels;
768 while (i < j) {
769 cdf_cpt[i].z = cdf_cpt[i-1].z + Ctrl->T.inc;
770 i++;
771 }
772 if (j == Ctrl->E.levels-1) cdf_cpt[j].z = G[0]->header->z_max;
773 }
774 else if (Ctrl->S.active || (Ctrl->E.active && !Ctrl->E.cdf)) { /* Make an equaldistant color map from G[k]->header->z_min to G[k]->header->z_max */
775 double start, range;
776
777 switch (Ctrl->S.kind) {
778 case -1:
779 start = -fabs (G[0]->header->z_min);
780 break;
781 case 1:
782 start = -fabs (G[0]->header->z_max);
783 break;
784 case -2:
785 start = -MIN (fabs (G[0]->header->z_min), fabs (G[0]->header->z_max));
786 break;
787 case 2:
788 start = -MAX (fabs (G[0]->header->z_min), fabs (G[0]->header->z_max));
789 break;
790 default:
791 start = G[0]->header->z_min;
792 break;
793 }
794 range = (Ctrl->S.kind) ? 2.0 * fabs (start) : G[0]->header->z_max - G[0]->header->z_min;
795 range *= (1.0 + GMT_CONV8_LIMIT); /* To ensure the max grid values do not exceed the CPT limit due to round-off issues */
796 start -= fabs (start) * GMT_CONV8_LIMIT; /* To ensure the start of cpt is less than min value due to roundoff */
797 Ctrl->T.inc = range / (double)(Ctrl->E.levels - 1);
798 cdf_cpt = gmt_M_memory (GMT, NULL, Ctrl->E.levels, struct CDF_CPT);
799 for (j = 0; j < Ctrl->E.levels; j++) cdf_cpt[j].z = start + j * Ctrl->T.inc;
800 }
801 else if (Ctrl->E.cdf) { /* Use the cumulative weighted distribution to issue the desired equal-area level boundaries */
802 double dw = pair[nxy-1].weight / (Ctrl->E.levels - 1), p = 0.0, dp = 1.0 / (Ctrl->E.levels - 1);
803
804 cdf_cpt = gmt_M_memory (GMT, NULL, Ctrl->E.levels, struct CDF_CPT);
805 cdf_cpt[0].z = pair[0].value;
806 wsum = dw; p = dp;
807 GMT_Report (API, GMT_MSG_INFORMATION, "Evaluated %d equidistant points on the cumulative density function:\n", Ctrl->E.levels);
808 GMT_Report (API, GMT_MSG_INFORMATION, "z = %16g cdf(z) = %6.4f\n", cdf_cpt[0].z, 0.0);
809 for (j = 1, k = 0; j < (Ctrl->E.levels - 1); j++) {
810 while (k < nxy && pair[k].weight < wsum) k++; /* k is the point with weight >= wsum; so a linear interpolation */
811 cdf_cpt[j].z = pair[k-1].value + (wsum - pair[k-1].weight) * (pair[k].value - pair[k-1].value) / (pair[k].weight - pair[k-1].weight);
812 GMT_Report (API, GMT_MSG_INFORMATION, "z = %16g cdf(z) = %6.4f\n", cdf_cpt[j].z, p);
813 wsum += dw; /* Next area boundary */
814 p += dp; /* Next CDF value */
815 }
816 cdf_cpt[Ctrl->E.levels-1].z = pair[nxy-1].value;
817 GMT_Report (API, GMT_MSG_INFORMATION, "z = %16g cdf(z) = %6.4f\n", cdf_cpt[Ctrl->E.levels-1].z, 1.0);
818 gmt_M_free (GMT, pair);
819 /* Make sure we do not have slices with no z-range */
820 p = 0.0;
821 dp = GMT_CONV8_LIMIT * (cdf_cpt[Ctrl->E.levels-1].z - cdf_cpt[0].z);
822 for (j = 1; j < (Ctrl->E.levels - 1); j++) {
823 if (doubleAlmostEqualZero (cdf_cpt[j-1].z, cdf_cpt[j].z+p)) {
824 GMT_Report (API, GMT_MSG_WARNING, "CDF is vertical, adding %g to have monotonic increasing z-values:\n", dp);
825 p += dp;
826 cdf_cpt[j].z += p;
827 }
828 }
829 }
830 else { /* This is completely ad-hoc. It chooses z based on equidistant steps [of 0.1 unless -Sn set] for a Gaussian CDF: */
831 double z_inc = 1.0 / (Ctrl->E.levels - 1); /* Increment between selected points [0.1] */
832 double zcrit_tail = gmt_zcrit (GMT, 1.0 - z_inc); /* Get the +/- z-value containing bulk of distribution, with z_inc in each tail */
833 cdf_cpt = gmt_M_memory (GMT, NULL, Ctrl->E.levels, struct CDF_CPT);
834 if ((mean - zcrit_tail*sd) <= G[0]->header->z_min || (mean + zcrit_tail*sd) >= G[0]->header->z_max) {
835 /* Adjust mean/std so that our critical locations are still inside the min/max of the data */
836 mean = 0.5 * (G[0]->header->z_min + G[0]->header->z_max);
837 sd = (G[0]->header->z_max - mean) / 1.5; /* This factor of 1.5 probably needs to change since z_inc is no longer fixed at 0.1 */
838 if (sd <= 0.0) {
839 GMT_Report (API, GMT_MSG_ERROR, "Min and Max data values are equal.\n");
840 gmt_M_free (GMT, cdf_cpt);
841 Return (GMT_RUNTIME_ERROR);
842 }
843 } /* End of stupid bug fix */
844
845 /* So we go in steps of z_inc in the Gaussian CDF except we start and stop at actual min/max */
846 cdf_cpt[0].z = G[0]->header->z_min;
847 for (j = 1; j < (Ctrl->E.levels - 1); j++) cdf_cpt[j].z = mean + gmt_zcrit (GMT, j *z_inc) * sd;
848 cdf_cpt[Ctrl->E.levels-1].z = G[0]->header->z_max;
849 }
850
851 /* Get here when we are ready to go. cdf_cpt[].z contains the sample points. */
852
853 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) sprintf (format, "z = %s and CDF(z) = %s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
854 for (j = 0; j < Ctrl->E.levels; j++) {
855 if (cdf_cpt[j].z <= G[0]->header->z_min)
856 cdf_cpt[j].f = 0.0;
857 else if (cdf_cpt[j].z >= G[0]->header->z_max)
858 cdf_cpt[j].f = 1.0;
859 else {
860 nfound = 0;
861 for (k = 0; k < ngrd; k++) { /* For each grid */
862 gmt_M_grd_loop (GMT, G[k], row, col, ij) {
863 if (!gmt_M_is_fnan (G[k]->data[ij]) && G[k]->data[ij] <= cdf_cpt[j].z) nfound++;
864 }
865 }
866 cdf_cpt[j].f = (double)(nfound-1)/(double)(ngood-1);
867 }
868 GMT_Report (API, GMT_MSG_INFORMATION, format, cdf_cpt[j].z, cdf_cpt[j].f);
869 }
870
871 /* Now the cdf function has been found. We now resample the chosen CPT */
872
873 z = gmt_M_memory (GMT, NULL, Ctrl->E.levels, double);
874 for (j = 0; j < Ctrl->E.levels; j++) z[j] = cdf_cpt[j].z;
875 if (Ctrl->Q.mode == 2) for (j = 0; j < Ctrl->E.levels; j++) z[j] = d_log10 (GMT, z[j]); /* Make log10(z) values for interpolation step */
876
877 signed_levels = Ctrl->E.levels;
878 Pout = gmt_sample_cpt (GMT, Pin, z, -signed_levels, Ctrl->Z.active, Ctrl->I.active, Ctrl->Q.mode, Ctrl->W.active); /* -ve to keep original colors */
879
880 /* Determine mode flags for output */
881 cpt_flags = 0;
882 if (Ctrl->N.active) cpt_flags |= GMT_CPT_NO_BNF; /* bit 0 controls if BFN will be written out */
883 if (Ctrl->D.mode == 1) cpt_flags |= GMT_CPT_EXTEND_BNF; /* bit 1 controls if BF will be set to equal bottom/top rgb value */
884 if (Ctrl->F.active) Pout->model = Ctrl->F.model;
885 if (Ctrl->F.cat) { /* Flag as a categorical CPT */
886 Pout->categorical = GMT_CPT_CATEGORICAL_VAL;
887 if (Ctrl->F.label) { /* Want categorical labels */
888 unsigned int ns = 0;
889 char **label = gmt_cat_cpt_strings (GMT, Ctrl->F.label, Pout->n_colors, &ns);
890 for (unsigned int k = 0; k < MIN (Pout->n_colors, ns); k++) {
891 if (Pout->data[k].label) gmt_M_str_free (Pout->data[k].label);
892 Pout->data[k].label = label[k]; /* Now the job of the CPT to free these strings */
893 }
894 gmt_M_free (GMT, label);
895 }
896 }
897
898 if (Ctrl->A.active) gmt_cpt_transparency (GMT, Pout, Ctrl->A.value, Ctrl->A.mode); /* Set transparency */
899
900 if (write && GMT_Write_Data (API, GMT_IS_PALETTE, GMT_IS_FILE, GMT_IS_NONE, cpt_flags, NULL, Ctrl->Out.file, Pout) != GMT_NOERROR)
901 error = API->error;
902
903 if (!write)
904 gmt_save_current_cpt (GMT, Pout, cpt_flags); /* Save for use by session, if modern */
905
906 gmt_M_free (GMT, cdf_cpt);
907 gmt_M_free (GMT, z);
908 if (error == GMT_NOERROR)
909 error = grd2cpt_free_the_grids (API, G, grdfile, ngrd);
910 else
911 grd2cpt_free_the_grids (API, G, grdfile, ngrd);
912 gmt_M_free (GMT, G);
913 gmt_M_free (GMT, grdfile);
914
915 Return ((error) ? error : GMT_NOERROR);
916 }
917