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 * API functions to support the grdedit application.
19 *
20 * Brief synopsis: grdedit reads an existing grid file and takes command
21 * line arguments to redefine some of the grdheader parameters:
22 *
23 * x_min/x_max OR x_inc (the other is recomputed)
24 * y_min/y_max OR y_inc (the other is recomputed)
25 * z_scale_factor/z_add_offset
26 * x_units/y_units/z_units
27 * title/command/remark
28 *
29 * Author: Paul Wessel
30 * Date: 1-JAN-2010
31 * Version: 6 API
32 */
33
34 #include "gmt_dev.h"
35
36 #define THIS_MODULE_CLASSIC_NAME "grdedit"
37 #define THIS_MODULE_MODERN_NAME "grdedit"
38 #define THIS_MODULE_LIB "core"
39 #define THIS_MODULE_PURPOSE "Modify header or content of a grid"
40 #define THIS_MODULE_KEYS "<G{,ND(,GG}"
41 #define THIS_MODULE_NEEDS ""
42 #define THIS_MODULE_OPTIONS "-:JRVbdefhiw" GMT_OPT("H")
43
44 struct GRDEDIT_CTRL {
45 struct GRDEDIT_In {
46 bool active;
47 char *file;
48 } In;
49 struct GRDEDIT_A { /* -A */
50 bool active;
51 } A;
52 struct GRDEDIT_C { /* -C */
53 bool active;
54 } C;
55 struct GRDEDIT_D { /* -D[+x<xname>][+yyname>][+z<zname>][+s<scale>][+ooffset>][+n<invalid>][+t<title>][+r<remark>] */
56 bool active;
57 char *information;
58 } D;
59 struct GRDEDIT_E { /* -E[a|h|l|r|t|v] */
60 bool active;
61 char mode; /* l rotate 90 degrees left (CCW), t = transpose, r = rotate 90 degrees right (CW) */
62 /* a rotate around (180), h flip grid horizontally (FLIPLR), v flip grid vertically (FLIPUD) */
63 } E;
64 struct GRDEDIT_G {
65 bool active;
66 char *file;
67 } G;
68 struct GRDEDIT_L { /* -L[+n|p] */
69 bool active;
70 int mode;
71 } L;
72 struct GRDEDIT_N { /* N<xyzfile> */
73 bool active;
74 char *file;
75 } N;
76 struct GRDEDIT_S { /* -S */
77 bool active;
78 } S;
79 struct GRDEDIT_T { /* -T */
80 bool active;
81 } T;
82 };
83
New_Ctrl(struct GMT_CTRL * GMT)84 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
85 struct GRDEDIT_CTRL *C;
86
87 C = gmt_M_memory (GMT, NULL, 1, struct GRDEDIT_CTRL);
88
89 /* Initialize values whose defaults are not 0/false/NULL */
90
91 return (C);
92 }
93
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDEDIT_CTRL * C)94 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDEDIT_CTRL *C) { /* Deallocate control structure */
95 if (!C) return;
96 gmt_M_str_free (C->In.file);
97 gmt_M_str_free (C->D.information);
98 gmt_M_str_free (C->G.file);
99 gmt_M_str_free (C->N.file);
100 gmt_M_free (GMT, C);
101 }
102
usage(struct GMTAPI_CTRL * API,int level)103 static int usage (struct GMTAPI_CTRL *API, int level) {
104 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
105 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
106 GMT_Usage (API, 0, "usage: %s %s [-A] [-C] [%s] [-E[a|e|h|l|r|t|v]] [-G%s] [%s] [-L[+n|p]] "
107 "[-N<table>] [%s] [-S] [-T] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n", name, GMT_INGRID, GMT_GRDEDIT2D,
108 GMT_OUTGRID, GMT_J_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT,
109 GMT_i_OPT, GMT_w_OPT, GMT_colon_OPT, GMT_PAR_OPT);
110
111 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
112
113 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
114 gmt_ingrid_syntax (API, 0, "Name of grid to be modified");
115 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
116 GMT_Usage (API, 1, "\n-A Adjust dx/dy to be compatible with the file domain (or new -R).");
117 GMT_Usage (API, 1, "\n-C Remove the command history from the header.");
118 gmt_grd_info_syntax (API->GMT, 'D');
119 GMT_Usage (API, 1, "\n-E[a|e|h|l|r|t|v]");
120 GMT_Usage (API, -2, "Transform the entire grid (this may exchange x and y). Append operation:");
121 GMT_Usage (API, 3, "a: Flip grid horizontally and vertically (h + v).");
122 GMT_Usage (API, 3, "e: Exchange x(lon) and y(lat).");
123 GMT_Usage (API, 3, "h: Flip grid left-to-right (as grdmath FLIPLR).");
124 GMT_Usage (API, 3, "l: Rotate grid 90 degrees left (counter-clockwise).");
125 GMT_Usage (API, 3, "r: Rotate grid 90 degrees right (clockwise).");
126 GMT_Usage (API, 3, "t: Transpose grid [Default].");
127 GMT_Usage (API, 3, "v: Flip grid top-to-bottom (as grdmath FLIPUD).");
128 gmt_outgrid_syntax (API, 'G', "Specify new output grid file [Default updates given grid file]");
129 GMT_Option (API, "J");
130 GMT_Usage (API, 1, "\n-L[+n|p]");
131 GMT_Usage (API, -2, "Shift the grid\'s longitude range (geographic grids only):");
132 GMT_Usage (API, 3, "+n Adjust <west>/<east> so <east> <= 0.");
133 GMT_Usage (API, 3, "+p Adjust <west>/<east> so <west> >= 0.");
134 GMT_Usage (API, -2, "Default adjusts <west>/<east> so <west> >= -180 or <east> <= +180.");
135 GMT_Usage (API, 1, "\n-N<table>");
136 GMT_Usage (API, -2, "Give <table> with new xyz values to replace the corresponding, existing grid nodes.");
137 GMT_Option (API, "R");
138 GMT_Usage (API, 1, "\n-S For global grids of 360 degree longitude range: "
139 "Will rotate entire grid in longitude to coincide with new borders in -R.");
140 GMT_Usage (API, 1, "\n-T Toggle header from grid-line to pixel-registered grid or vice versa. "
141 "Note: This shrinks the grid region by 0.5*{dx,dy} going from pixel to grid-line registration "
142 "and expands it by 0.5*{dx,dy} going from grid-line to pixel registration. "
143 "No grid values will be changed (it is a non-destructive change to the grid).");
144 GMT_Option (API, "V,bi3,di,e,f,h,i,w,:,.");
145
146 return (GMT_MODULE_USAGE);
147 }
148
parse(struct GMT_CTRL * GMT,struct GRDEDIT_CTRL * Ctrl,struct GMT_OPTION * options)149 static int parse (struct GMT_CTRL *GMT, struct GRDEDIT_CTRL *Ctrl, struct GMT_OPTION *options) {
150
151 /* This parses the options provided to grdedit and sets parameters in Ctrl.
152 * Note Ctrl has already been initialized and non-zero default values set.
153 * Any GMT common options will override values set previously by other commands.
154 * It also replaces any file names specified as input or output with the data ID
155 * returned when registering these sources/destinations with the API.
156 */
157
158 unsigned int n_errors = 0, n_files = 0;
159 struct GMT_OPTION *opt = NULL;
160 struct GMTAPI_CTRL *API = GMT->parent;
161
162 for (opt = options; opt; opt = opt->next) { /* Process all the options given */
163 switch (opt->option) {
164 /* Common parameters */
165
166 case '<': /* Input file (only one is accepted) */
167 if (n_files++ > 0) {n_errors++; continue; }
168 Ctrl->In.active = true;
169 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
170 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
171 break;
172
173 /* Processes program-specific parameters */
174
175 case 'A': /* Adjust increments */
176 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
177 Ctrl->A.active = true;
178 break;
179 case 'C': /* Clear history */
180 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
181 Ctrl->C.active = true;
182 break;
183 case 'D': /* Give grid information */
184 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
185 Ctrl->D.active = true;
186 Ctrl->D.information = strdup (opt->arg);
187 break;
188 case 'E': /* Transpose or rotate grid */
189 n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
190 Ctrl->E.active = true;
191 if (opt->arg[0] == '\0') /* Default transpose */
192 Ctrl->E.mode = 't';
193 else if (strchr ("aehlrtv", opt->arg[0]))
194 Ctrl->E.mode = opt->arg[0];
195 else {
196 n_errors++;
197 GMT_Report (API, GMT_MSG_ERROR, "Option -E: Unrecognized modifier %c\n", opt->arg[0]);
198 }
199 break;
200 case 'G': /* Separate output grid file */
201 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
202 Ctrl->G.active = true;
203 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
204 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
205 break;
206 case 'L': /* Rotate w/e */
207 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
208 Ctrl->L.active = true;
209 if (strstr (opt->arg, "+n")) Ctrl->L.mode = -1;
210 else if (strstr (opt->arg, "+p")) Ctrl->L.mode = +1;
211 else if (opt->arg[0])
212 n_errors++;
213 break;
214 case 'N': /* Replace nodes */
215 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
216 Ctrl->N.active = true;
217 if (opt->arg[0]) Ctrl->N.file = strdup (opt->arg);
218 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->N.file))) n_errors++;
219 break;
220 case 'S': /* Rotate global grid */
221 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
222 Ctrl->S.active = true;
223 break;
224 case 'T': /* Toggle registration */
225 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
226 Ctrl->T.active = true;
227 break;
228
229 default: /* Report bad options */
230 n_errors += gmt_default_error (GMT, opt->option);
231 break;
232 }
233 }
234
235 n_errors += gmt_M_check_condition (GMT, Ctrl->G.active && !Ctrl->G.file, "Option -G: Must specify an output grid file\n");
236 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->A.active,
237 "Option -S: Incompatible with -A\n");
238 n_errors += gmt_M_check_condition (GMT, Ctrl->E.active &&
239 (Ctrl->A.active || Ctrl->D.active || Ctrl->N.active || Ctrl->S.active || Ctrl->T.active),
240 "Option -E: Incompatible with -A, -D, -N, -S, and -T\n");
241 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->T.active,
242 "Option -S: Incompatible with -T\n");
243 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && Ctrl->N.active,
244 "Option -S: Incompatible with -N\n");
245 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && !GMT->common.R.active[RSET],
246 "Option -S: Must also specify -R\n");
247 n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && !gmt_M_360_range (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]),
248 "Option -S: -R longitudes must span exactly 360 degrees\n");
249 n_errors += gmt_M_check_condition (GMT, n_files != 1, "Must specify a single grid file\n");
250 if (Ctrl->N.active) {
251 n_errors += gmt_check_binary_io (GMT, 3);
252 }
253
254 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
255 }
256
257 #define bailout(code) {gmt_M_free_options (mode); return (code);}
258 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
259
GMT_grdedit(void * V_API,int mode,void * args)260 EXTERN_MSC int GMT_grdedit (void *V_API, int mode, void *args) {
261 /* High-level function that implements the grdedit task */
262 bool grid_was_read = false, do_J = false;
263
264 unsigned int row, col;
265 int error;
266
267 uint64_t ij, n_data, n_use;
268
269 double shift_amount = 0.0;
270
271 char *registration[2] = {"gridline", "pixel"}, *out_file = NULL, *projstring = NULL;
272
273 struct GRDEDIT_CTRL *Ctrl = NULL;
274 struct GMT_GRID *G = NULL;
275 struct GMT_GRID_HEADER_HIDDEN *HH = NULL;
276 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
277 struct GMT_OPTION *options = NULL;
278 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
279
280 /*----------------------- Standard module initialization and parsing ----------------------*/
281
282 if (API == NULL) return (GMT_NOT_A_SESSION);
283 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
284 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
285
286 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
287
288 /* Parse the command-line arguments */
289
290 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 */
291 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
292 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
293 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
294
295 /*---------------------------- This is the grdedit main code ----------------------------*/
296
297 if ((do_J = GMT->common.J.active)) { /* Only gave -J to set the proj4 flag, so grid is already projected (i.e., Cartesian) */
298 projstring = gmt_export2proj4 (GMT); /* Convert the GMT -J<...> into a proj4 string */
299 if (strstr (projstring, "+proj=longlat") || strstr (projstring, "+proj=latlong")) { /* These means we have a geographic grid */
300 gmt_set_geographic (GMT, GMT_IN); /* Force Geographic */
301 gmt_set_geographic (GMT, GMT_OUT);
302 }
303 else {
304 gmt_set_cartesian (GMT, GMT_IN); /* Force Cartesian since processing -J changed that */
305 gmt_set_cartesian (GMT, GMT_OUT);
306 gmt_parse_R_option (GMT, GMT->common.R.string); /* Since parsing under -J would have imposed checks. */
307 }
308 GMT->current.proj.projection = GMT->current.proj.projection_GMT = GMT_NO_PROJ;
309 GMT->common.J.active = false; /* Leave no trace of this, except in the history... */
310 }
311
312 if ((G = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, Ctrl->G.active ? GMT_CONTAINER_AND_DATA : GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get header only */
313 free(projstring);
314 Return (API->error);
315 }
316 grid_was_read = Ctrl->G.active;
317 HH = gmt_get_H_hidden (G->header);
318 if (gmt_M_is_geographic (GMT, GMT_IN)) gmt_set_geographic (GMT, GMT_OUT); /* Out same as in */
319
320 if ((G->header->type == GMT_GRID_IS_SF || G->header->type == GMT_GRID_IS_SD) && Ctrl->T.active) {
321 GMT_Report (API, GMT_MSG_ERROR, "Toggling registrations not possible for Surfer grid formats\n");
322 GMT_Report (API, GMT_MSG_ERROR, "(Use grdconvert to convert to GMT default format and work on that file)\n");
323 free(projstring);
324 Return (GMT_RUNTIME_ERROR);
325 }
326
327 if (Ctrl->S.active && !gmt_grd_is_global (GMT, G->header)) {
328 GMT_Report (API, GMT_MSG_ERROR, "Shift only allowed for global grids\n");
329 free(projstring);
330 Return (GMT_RUNTIME_ERROR);
331 }
332
333 out_file = (Ctrl->G.active) ? Ctrl->G.file : Ctrl->In.file; /* Where to write the modified grid */
334
335 GMT_Report (API, GMT_MSG_INFORMATION, "Editing parameters for grid %s:\n", out_file);
336
337 /* Decode grd information given, if any */
338
339 if (do_J) /* Save proj4 string in the header */
340 G->header->ProjRefPROJ4 = projstring;
341
342 if (Ctrl->D.active) {
343 double scale_factor, add_offset;
344 gmt_grdfloat nan_value;
345 GMT_Report (API, GMT_MSG_INFORMATION, "Decode and change attributes in file %s\n", out_file);
346 scale_factor = G->header->z_scale_factor;
347 add_offset = G->header->z_add_offset;
348 nan_value = G->header->nan_value;
349 if (gmt_decode_grd_h_info (GMT, Ctrl->D.information, G->header))
350 Return (GMT_PARSE_ERROR);
351
352 if (nan_value != G->header->nan_value) {
353 /* Must read data */
354 if (!grid_was_read && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, G) == NULL)
355 Return (API->error);
356 grid_was_read = true;
357 /* Recalculate z_min/z_max */
358 gmt_grd_zminmax (GMT, G->header, G->data);
359 }
360 if (scale_factor != G->header->z_scale_factor || add_offset != G->header->z_add_offset) {
361 G->header->z_min = (G->header->z_min - add_offset) / scale_factor * G->header->z_scale_factor + G->header->z_add_offset;
362 G->header->z_max = (G->header->z_max - add_offset) / scale_factor * G->header->z_scale_factor + G->header->z_add_offset;
363 }
364 }
365
366 if (Ctrl->C.active) { /* Wipe history */
367 gmt_M_memset (G->header->command, GMT_GRID_COMMAND_LEN320, char);
368 }
369
370 if (Ctrl->S.active) {
371 shift_amount = GMT->common.R.wesn[XLO] - G->header->wesn[XLO];
372 GMT_Report (API, GMT_MSG_INFORMATION, "Shifting longitudes in file %s by %g degrees\n", out_file, shift_amount);
373 if (!grid_was_read && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, G) == NULL) { /* Get data */
374 Return (API->error);
375 }
376 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
377 Return (API->error);
378 }
379 gmt_grd_shift (GMT, G, shift_amount);
380 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, out_file, G) != GMT_NOERROR) {
381 Return (API->error);
382 }
383 }
384 else if (Ctrl->N.active) {
385 int in_ID;
386 struct GMT_RECORD *In = NULL;
387 double *in = NULL;
388 GMT_Report (API, GMT_MSG_INFORMATION, "Replacing nodes using xyz values from file %s\n", Ctrl->N.file);
389
390 if (!grid_was_read && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, G) == NULL) { /* Get data */
391 Return (API->error);
392 }
393 /* Must register Ctrl->N.file first since we are going to read rec-by-rec from all available sources */
394 if ((in_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_IN, NULL, Ctrl->N.file)) == GMT_NOTSET) {
395 GMT_Report (API, GMT_MSG_ERROR, "Unable to register file %s\n", Ctrl->N.file);
396 Return (GMT_RUNTIME_ERROR);
397 }
398 /* Initialize the i/o since we are doing record-by-record reading/writing */
399 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_EXISTING, 0, options) != GMT_NOERROR) {
400 Return (API->error); /* Establishes data input */
401 }
402 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data input and sets access mode */
403 Return (API->error);
404 }
405 GMT_Set_Columns (API, GMT_IN, 3, GMT_COL_FIX_NO_TEXT);
406
407 n_data = n_use = 0;
408 do { /* Keep returning records until we reach EOF */
409 if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
410 if (gmt_M_rec_is_error (GMT)) { /* Bail if there are any read errors */
411 Return (GMT_RUNTIME_ERROR);
412 }
413 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
414 break;
415 continue; /* Go back and read the next record */
416 }
417 if (In->data == NULL) {
418 gmt_quit_bad_record (API, In);
419 Return (API->error);
420 }
421
422 in = In->data; /* Only need to process numerical part here */
423
424 /* Data record to process */
425
426 n_data++;
427
428 if (gmt_M_y_is_outside (GMT, in[GMT_Y], G->header->wesn[YLO], G->header->wesn[YHI])) continue; /* Outside y-range */
429 if (gmt_x_is_outside (GMT, &in[GMT_X], G->header->wesn[XLO], G->header->wesn[XHI])) continue; /* Outside x-range */
430 if (gmt_row_col_out_of_bounds (GMT, in, G->header, &row, &col)) continue; /* Outside grid node range */
431
432 ij = gmt_M_ijp (G->header, row, col);
433 G->data[ij] = (gmt_grdfloat)in[GMT_Z];
434 n_use++;
435 if (gmt_M_grd_duplicate_column (GMT, G->header, GMT_IN)) { /* Make sure longitudes got replicated */
436 /* Possibly need to replicate e/w value */
437 if (col == 0) {ij = gmt_M_ijp (G->header, row, G->header->n_columns-1); G->data[ij] = (gmt_grdfloat)in[GMT_Z]; n_use++; }
438 else if (col == (G->header->n_columns-1)) {ij = gmt_M_ijp (G->header, row, 0); G->data[ij] = (gmt_grdfloat)in[GMT_Z]; n_use++; }
439 }
440 } while (true);
441
442 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
443 Return (API->error);
444 }
445
446 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, out_file, G) != GMT_NOERROR) {
447 Return (API->error);
448 }
449 GMT_Report (API, GMT_MSG_INFORMATION, "Read %" PRIu64 " new data points, updated %" PRIu64 ".\n", n_data, n_use);
450 }
451 else if (Ctrl->E.active) { /* Transpose, flip, or rotate the matrix and possibly exchange x and y info */
452 struct GMT_GRID_HEADER *h_tr = NULL;
453 uint64_t ij, ij_tr = 0;
454 gmt_grdfloat *a_tr = NULL, *save_grid_pointer = NULL;
455
456 if (!grid_was_read && GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, G) == NULL) { /* Get data */
457 Return (API->error);
458 }
459
460 switch (Ctrl->E.mode) {
461 case 'a': /* Rotate grid around 180 degrees */
462 GMT_Report (API, GMT_MSG_INFORMATION, "Flip grid horizontally and vertically\n");
463 break;
464 case 'e': /* Exchange lon and lat */
465 GMT_Report (API, GMT_MSG_INFORMATION, "Exchange x|longitude and y|latitude\n");
466 break;
467 case 'h': /* Flip grid horizontally */
468 GMT_Report (API, GMT_MSG_INFORMATION, "Flip grid horizontally (FLIPLR)\n");
469 break;
470 case 'l': /* Rotate grid 90 CW */
471 GMT_Report (API, GMT_MSG_INFORMATION, "Rotate grid 90 degrees left (Counter-Clockwise)\n");
472 break;
473 case 't': /* Transpose grid */
474 GMT_Report (API, GMT_MSG_INFORMATION, "Transpose grid\n");
475 break;
476 case 'r': /* Rotate grid 90 CCW */
477 GMT_Report (API, GMT_MSG_INFORMATION, "Rotate grid 90 degrees right (Clockwise)\n");
478 break;
479 case 'v': /* Flip grid vertically */
480 GMT_Report (API, GMT_MSG_INFORMATION, "Flip grid vertically (FLIPUD)\n");
481 break;
482 }
483
484 h_tr = gmt_get_header (GMT);
485 gmt_copy_gridheader (GMT, h_tr, G->header); /* First make a copy of header */
486 if (strchr ("eltr", Ctrl->E.mode)) { /* These operators interchange x and y */
487 h_tr->wesn[XLO] = G->header->wesn[YLO];
488 h_tr->wesn[XHI] = G->header->wesn[YHI];
489 h_tr->inc[GMT_X] = G->header->inc[GMT_Y];
490 strncpy (h_tr->x_units, G->header->y_units, GMT_GRID_UNIT_LEN80);
491 h_tr->wesn[YLO] = G->header->wesn[XLO];
492 h_tr->wesn[YHI] = G->header->wesn[XHI];
493 h_tr->inc[GMT_Y] = G->header->inc[GMT_X];
494 strncpy (h_tr->y_units, G->header->x_units, GMT_GRID_UNIT_LEN80);
495 gmt_set_grddim (GMT, h_tr); /* Recompute n_columns, n_rows, mx, size, etc */
496 gmt_M_doublep_swap (G->x, G->y);
497 }
498
499 /* Now transpose the matrix */
500
501 a_tr = gmt_M_memory (GMT, NULL, G->header->size, gmt_grdfloat);
502 gmt_M_grd_loop (GMT, G, row, col, ij) {
503 switch (Ctrl->E.mode) {
504 case 'a': /* Rotate grid around 180 degrees */
505 ij_tr = gmt_M_ijp (h_tr, G->header->n_rows-1-row, G->header->n_columns-1-col);
506 break;
507 case 'e': /* Exchange coordinates */
508 ij_tr = gmt_M_ijp (h_tr, h_tr->n_rows-1-col, h_tr->n_columns-1-row);
509 break;
510 case 'h': /* Flip horizontally (FLIPLR) */
511 ij_tr = gmt_M_ijp (h_tr, row, G->header->n_columns-1-col);
512 break;
513 case 'l': /* Rotate 90 CCW */
514 ij_tr = gmt_M_ijp (h_tr, G->header->n_columns-1-col, row);
515 break;
516 case 'r': /* Rotate 90 CW */
517 ij_tr = gmt_M_ijp (h_tr, col, G->header->n_rows-1-row);
518 break;
519 case 't': /* Transpose */
520 ij_tr = gmt_M_ijp (h_tr, col, row);
521 break;
522 case 'v': /* Flip vertically (FLIPUD) */
523 ij_tr = gmt_M_ijp (h_tr, G->header->n_rows-1-row, col);
524 break;
525 }
526 a_tr[ij_tr] = G->data[ij];
527 }
528 save_grid_pointer = G->data; /* Save original grid pointer and hook on the modified grid instead */
529 G->data = a_tr;
530 gmt_copy_gridheader (GMT, G->header, h_tr); /* Update to the new header */
531 gmt_M_free (GMT, h_tr->hidden);
532 gmt_M_free (GMT, h_tr);
533 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, out_file, G) != GMT_NOERROR) {
534 Return (API->error);
535 }
536 G->data = save_grid_pointer;
537 gmt_M_free (GMT, a_tr);
538 }
539 else if (Ctrl->L.active) { /* Wrap the longitude boundaries */
540 double wesn[4];
541
542 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
543 Return (API->error);
544 }
545 gmt_M_memcpy (wesn, G->header->wesn, 4, double); /* Copy */
546 if (Ctrl->L.mode == -1) {
547 while (G->header->wesn[XHI] > 0.0) G->header->wesn[XLO] -= 360.0, G->header->wesn[XHI] -= 360.0;
548 }
549 else if (Ctrl->L.mode == +1) {
550 while (G->header->wesn[XLO] < 0.0) G->header->wesn[XLO] += 360.0, G->header->wesn[XHI] += 360.0;
551 }
552 else { /* Aim for a reasonable w/e range */
553 if (G->header->wesn[XHI] < -180.0) G->header->wesn[XLO] += 360.0, G->header->wesn[XHI] += 360.0;
554 if (G->header->wesn[XLO] < -180.0 && (G->header->wesn[XHI] + 360.0) <= 180.0) G->header->wesn[XLO] += 360.0, G->header->wesn[XHI] += 360.0;
555 if (G->header->wesn[XLO] > 360.0) G->header->wesn[XLO] -= 360.0, G->header->wesn[XHI] -= 360.0;
556 if (G->header->wesn[XHI] > 180.0 && (G->header->wesn[XLO] - 360.0) >= -180.0) G->header->wesn[XLO] -= 360.0, G->header->wesn[XHI] -= 360.0;
557 }
558 if (doubleAlmostEqual (G->header->wesn[XLO], wesn[XLO])) {
559 GMT_Report (API, GMT_MSG_INFORMATION, "Change region in file %s from %g/%g/%g/%g to %g/%g/%g/%g\n", out_file,
560 G->header->wesn[XLO], G->header->wesn[XHI], G->header->wesn[YLO], G->header->wesn[YHI],
561 wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
562 }
563 else {
564 GMT_Report (API, GMT_MSG_INFORMATION, "Region in file remains unchanged: %g/%g/%g/%g\n", out_file,
565 G->header->wesn[XLO], G->header->wesn[XHI], G->header->wesn[YLO], G->header->wesn[YHI]);
566 }
567 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, Ctrl->G.active ? GMT_CONTAINER_AND_DATA : GMT_CONTAINER_ONLY, NULL, out_file, G) != GMT_NOERROR) {
568 Return (API->error);
569 }
570 }
571 else { /* Change the domain boundaries */
572 if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) { /* Disables further data input */
573 Return (API->error);
574 }
575 if (Ctrl->T.active) { /* Grid-line <---> Pixel toggling of the header */
576 gmt_change_grdreg (GMT, G->header, 1 - G->header->registration);
577 GMT_Report (API, GMT_MSG_INFORMATION, "Toggled registration mode in file %s from %s to %s\n",
578 out_file, registration[1-G->header->registration], registration[G->header->registration]);
579 GMT_Report (API, GMT_MSG_INFORMATION, "Reset region in file %s to %g/%g/%g/%g\n",
580 out_file, G->header->wesn[XLO], G->header->wesn[XHI], G->header->wesn[YLO], G->header->wesn[YHI]);
581 }
582 if (GMT->common.R.active[RSET]) {
583 GMT_Report (API, GMT_MSG_INFORMATION, "Reset region in file %s to %g/%g/%g/%g\n",
584 out_file, GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI], GMT->common.R.wesn[YLO], GMT->common.R.wesn[YHI]);
585 gmt_M_memcpy (G->header->wesn, GMT->common.R.wesn, 4, double);
586 Ctrl->A.active = true; /* Must ensure -R -I compatibility */
587 }
588 if (Ctrl->A.active) {
589 G->header->inc[GMT_X] = gmt_M_get_inc (GMT, G->header->wesn[XLO], G->header->wesn[XHI], G->header->n_columns, G->header->registration);
590 G->header->inc[GMT_Y] = gmt_M_get_inc (GMT, G->header->wesn[YLO], G->header->wesn[YHI], G->header->n_rows, G->header->registration);
591 GMT_Report (API, GMT_MSG_INFORMATION, "Reset grid-spacing in file %s to %g/%g\n",
592 out_file, G->header->inc[GMT_X], G->header->inc[GMT_Y]);
593 }
594 if (gmt_M_is_geographic (GMT, GMT_IN) && gmt_M_is_cartesian (GMT, GMT_OUT)) { /* Force a switch from geographic to Cartesian */
595 HH->grdtype = GMT_GRID_CARTESIAN;
596 strcpy (G->header->x_units, "x_units");
597 strcpy (G->header->y_units, "y_units");
598 }
599 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, Ctrl->G.active ? GMT_CONTAINER_AND_DATA : GMT_CONTAINER_ONLY, NULL, out_file, G) != GMT_NOERROR) {
600 Return (API->error);
601 }
602 }
603
604 GMT_Report (API, GMT_MSG_INFORMATION, Ctrl->G.active ? "Modified grid written to file %s.\n" : "File %s updated.\n", out_file);
605
606 Return (GMT_NOERROR);
607 }
608