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: grdpaste.c reads two grid files and writes a new file with
19 * the first two pasted together along their common edge.
20 *
21 * Author: Walter Smith
22 * Date: 1-JAN-2010
23 * Version: 6 API
24 */
25
26 #include "gmt_dev.h"
27
28 #define THIS_MODULE_CLASSIC_NAME "grdpaste"
29 #define THIS_MODULE_MODERN_NAME "grdpaste"
30 #define THIS_MODULE_LIB "core"
31 #define THIS_MODULE_PURPOSE "Join two grids along their common edge"
32 #define THIS_MODULE_KEYS "<G{2,GG}"
33 #define THIS_MODULE_NEEDS ""
34 #define THIS_MODULE_OPTIONS "-Vf"
35
36 struct GRDPASTE_CTRL {
37 struct GRDPASTE_In {
38 bool active;
39 char *file[2];
40 } In;
41 struct GRDPASTE_G { /* -G<output_grdfile> */
42 bool active;
43 char *file;
44 } G;
45 };
46
New_Ctrl(struct GMT_CTRL * GMT)47 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
48 struct GRDPASTE_CTRL *C = NULL;
49
50 C = gmt_M_memory (GMT, NULL, 1, struct GRDPASTE_CTRL);
51
52 /* Initialize values whose defaults are not 0/false/NULL */
53
54 return (C);
55 }
56
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDPASTE_CTRL * C)57 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDPASTE_CTRL *C) { /* Deallocate control structure */
58 if (!C) return;
59 gmt_M_str_free (C->G.file);
60 gmt_M_str_free (C->In.file[GMT_IN]);
61 gmt_M_str_free (C->In.file[GMT_OUT]);
62 gmt_M_free (GMT, C);
63 }
64
usage(struct GMTAPI_CTRL * API,int level)65 static int usage (struct GMTAPI_CTRL *API, int level) {
66 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
67 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
68 GMT_Usage (API, 0, "usage: %s <grid1> <grid2> -G%s [%s] [%s] [%s]\n", name, GMT_OUTGRID, GMT_V_OPT, GMT_f_OPT, GMT_PAR_OPT);
69
70 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
71
72 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
73 GMT_Usage (API, 1, "\n<grid1> and <grid2> are to be combined into <outgrid>. "
74 "They must have same increments and registration and one edge in common. "
75 "If in doubt, run grdinfo first and check your files. "
76 "Use grdpaste and/or grdsample to adjust files as necessary. "
77 "If grids are geographic and adds to full 360-degree range then <grid1> "
78 "determines west. Use grdedit -S to rotate grid to another -Rw/e/s/n.");
79 gmt_outgrid_syntax (API, 'G', "Set name of the output grid file");
80 if (gmt_M_showusage (API)) GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
81 GMT_Option (API, "V,f,.");
82
83 return (GMT_MODULE_USAGE);
84 }
85
parse(struct GMT_CTRL * GMT,struct GRDPASTE_CTRL * Ctrl,struct GMT_OPTION * options)86 static int parse (struct GMT_CTRL *GMT, struct GRDPASTE_CTRL *Ctrl, struct GMT_OPTION *options) {
87 /* This parses the options provided to grdpaste and sets parameters in CTRL.
88 * Any GMT common options will override values set previously by other commands.
89 * It also replaces any file names specified as input or output with the data ID
90 * returned when registering these sources/destinations with the API.
91 */
92
93 unsigned int n_errors = 0, n_in = 0;
94 struct GMT_OPTION *opt = NULL;
95 struct GMTAPI_CTRL *API = GMT->parent;
96
97 for (opt = options; opt; opt = opt->next) {
98 switch (opt->option) {
99
100 case '<': /* Input files */
101 if (n_in == 2) {
102 n_errors++;
103 GMT_Report (API, GMT_MSG_ERROR, "Only two files may be pasted\n");
104 }
105 else {
106 Ctrl->In.file[n_in] = strdup (opt->arg);
107 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file[n_in]))) n_errors++;
108 n_in++;
109 }
110 break;
111
112 /* Processes program-specific parameters */
113
114 case 'G':
115 Ctrl->G.active = true;
116 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
117 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
118 break;
119
120 default: /* Report bad options */
121 n_errors += gmt_default_error (GMT, opt->option);
122 break;
123 }
124 }
125
126 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file[0] || !Ctrl->In.file[1], "Must specify two input files\n");
127 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file, "Option -G: Must specify output file\n");
128
129 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
130 }
131
132 #define bailout(code) {gmt_M_free_options (mode); return (code);}
133 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
134
135 /* True if grid is a COARDS/CF netCDF file */
grdpaste_is_nc_grid(struct GMT_GRID * grid)136 GMT_LOCAL inline bool grdpaste_is_nc_grid (struct GMT_GRID *grid) {
137 return
138 grid->header->type == GMT_GRID_IS_NB ||
139 grid->header->type == GMT_GRID_IS_NS ||
140 grid->header->type == GMT_GRID_IS_NI ||
141 grid->header->type == GMT_GRID_IS_NF ||
142 grid->header->type == GMT_GRID_IS_ND;
143 }
144
GMT_grdpaste(void * V_API,int mode,void * args)145 EXTERN_MSC int GMT_grdpaste (void *V_API, int mode, void *args) {
146 int error = 0, way = 0;
147 unsigned int one_or_zero;
148 bool common_y = false;
149
150 char format[GMT_BUFSIZ];
151
152 double x_noise, y_noise;
153
154 struct GMT_GRID *A = NULL, *B = NULL, *C = NULL;
155 struct GMT_GRID_HEADER_HIDDEN *AH = NULL, *BH = NULL;
156 struct GRDPASTE_CTRL *Ctrl = NULL;
157 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
158 struct GMT_OPTION *options = NULL;
159 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
160
161 /*----------------------- Standard module initialization and parsing ----------------------*/
162
163 if (API == NULL) return (GMT_NOT_A_SESSION);
164 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
165 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
166
167 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
168
169 /* Parse the command-line arguments */
170
171 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 */
172 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
173 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
174 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
175
176 /*---------------------------- This is the grdpaste main code ----------------------------*/
177
178 GMT_Report (API, GMT_MSG_INFORMATION, "Processing input grids\n");
179 gmt_set_pad (GMT, 0); /* No padding */
180
181 /* Try to find a common side to join on */
182
183 if ((A = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file[0], NULL)) == NULL) { /* Get header only */
184 Return (API->error);
185 }
186 if ((B = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file[1], NULL)) == NULL) { /* Get header only */
187 Return (API->error);
188 }
189
190 if (A->header->registration != B->header->registration)
191 error++;
192 if ((A->header->z_scale_factor != B->header->z_scale_factor) || (A->header->z_add_offset != B->header->z_add_offset)) {
193 GMT_Report (API, GMT_MSG_ERROR, "Scale/offset not compatible!\n");
194 Return (GMT_RUNTIME_ERROR);
195 }
196
197 if (! (fabs (A->header->inc[GMT_X] - B->header->inc[GMT_X]) < 1.0e-6 && fabs (A->header->inc[GMT_Y] - B->header->inc[GMT_Y]) < 1.0e-6)) {
198 GMT_Report (API, GMT_MSG_ERROR, "Grid intervals do not match!\n");
199 Return (GMT_RUNTIME_ERROR);
200 }
201
202 if ((C = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_NONE, A)) == NULL) Return (API->error); /* Just to get a header */
203
204 one_or_zero = A->header->registration == GMT_GRID_NODE_REG;
205 x_noise = GMT_CONV4_LIMIT * C->header->inc[GMT_X] * 10;
206 y_noise = GMT_CONV4_LIMIT * C->header->inc[GMT_Y] * 10;
207
208 AH = gmt_get_H_hidden (A->header);
209 BH = gmt_get_H_hidden (B->header);
210
211 common_y = (fabs (A->header->wesn[YLO] - B->header->wesn[YLO]) < y_noise && fabs (A->header->wesn[YHI] - B->header->wesn[YHI]) < y_noise);
212
213 if (gmt_M_x_is_lon (GMT, GMT_IN)) { /* Must be careful in determining a match since grids may differ by +/-360 in x */
214 double del;
215 if (common_y) { /* A and B are side-by-side, may differ by +-360 +- 1 pixel width */
216 del = A->header->wesn[XLO] - B->header->wesn[XHI]; /* Test if B left of A */
217 if (doubleAlmostEqual (fabs (del), 360.0)) { /* Let A remain left of B */
218 /* Since new range is a full 360 we always let the first grid decide west boundary.
219 * If this is not desirable the user should switch A and B or sue grdedit -S */
220 way = 4;
221 }
222 else if (del < (360.0 + C->header->inc[GMT_X] + x_noise) && del > (360.0 - C->header->inc[GMT_X] - x_noise)) {
223 B->header->wesn[XLO] += 360.0; B->header->wesn[XHI] += 360.0;
224 }
225 else if (del < (-360.0 + C->header->inc[GMT_X] + x_noise) && del > (-360.0 - C->header->inc[GMT_X] - x_noise)) {
226 A->header->wesn[XLO] += 360.0; A->header->wesn[XHI] += 360.0;
227 }
228 else { /* Neither, check if A is left of B */
229 del = B->header->wesn[XLO] - A->header->wesn[XHI];
230 if (del < (360.0 + C->header->inc[GMT_X] + x_noise) && del > (360.0 - C->header->inc[GMT_X] - x_noise)) {
231 A->header->wesn[XLO] += 360.0; A->header->wesn[XHI] += 360.0;
232 }
233 else if (del < (-360.0 + C->header->inc[GMT_X] + x_noise) && del > (-360.0 - C->header->inc[GMT_X] - x_noise)) {
234 B->header->wesn[XLO] += 360.0; B->header->wesn[XHI] += 360.0;
235 }
236 }
237 }
238 else { /* A and B are on top of each other, may differ by +-360 */
239 del = A->header->wesn[XLO] - B->header->wesn[XLO]; /* Test if B left of A */
240 if (del < (360.0 + x_noise) && del > (360.0 - x_noise)) {
241 B->header->wesn[XLO] += 360.0; B->header->wesn[XHI] += 360.0;
242 }
243 else if (del < (-360.0 + x_noise) && del > (-360.0 - x_noise)) {
244 A->header->wesn[XLO] += 360.0; A->header->wesn[XHI] += 360.0;
245 }
246 }
247 }
248
249 gmt_M_memcpy (C->header->wesn, A->header->wesn, 4, double); /* Output region is set as the same as A... */
250 if (common_y) {
251
252 C->header->n_rows = A->header->n_rows;
253
254 if (way == 4 || fabs (A->header->wesn[XHI] - B->header->wesn[XLO]) < x_noise) { /* A is on the left of B */
255 way = 4;
256 C->header->n_columns = A->header->n_columns + B->header->n_columns - one_or_zero;
257 C->header->wesn[XHI] = B->header->wesn[XHI]; /* ...but not for east */
258 }
259 else if (fabs (A->header->wesn[XLO] - B->header->wesn[XHI]) < x_noise) { /* A is on the right of B */
260 way = 3;
261 C->header->n_columns = A->header->n_columns + B->header->n_columns - one_or_zero;
262 C->header->wesn[XLO] = B->header->wesn[XLO]; /* ...but not for west */
263 }
264 else if ((fabs (A->header->wesn[XLO] - B->header->wesn[XHI]) < (C->header->inc[GMT_X] + x_noise)) ) {
265 /* A is on right of B but their pixel|grid reg limits under|overlap by one cell */
266 if (one_or_zero) /* Grid registration - underlap */
267 way = 32;
268 else /* Pixel registration - overlap */
269 way = 33;
270 C->header->n_columns = A->header->n_columns + B->header->n_columns - !one_or_zero;
271 C->header->wesn[XLO] = B->header->wesn[XLO]; /* ...but not for west */
272 }
273 else if ((fabs (A->header->wesn[XHI] - B->header->wesn[XLO]) < (C->header->inc[GMT_X] + x_noise)) ) {
274 /* A is on left of B but their pixel|grid reg limits under|overlap by one cell */
275 if (one_or_zero) /* Grid registration - underlap */
276 way = 43;
277 else /* Pixel registration - overlap */
278 way = 44;
279 C->header->n_columns = A->header->n_columns + B->header->n_columns - !one_or_zero;
280 C->header->wesn[XHI] = B->header->wesn[XHI]; /* ...but not for east */
281 }
282 else {
283 GMT_Report (API, GMT_MSG_ERROR, "Grids do not share a common edge!\n");
284 Return (GMT_RUNTIME_ERROR);
285 }
286 }
287 else if (fabs (A->header->wesn[XLO] - B->header->wesn[XLO]) < x_noise && fabs (A->header->wesn[XHI] - B->header->wesn[XHI]) < x_noise) {
288
289 C->header->n_columns = A->header->n_columns;
290
291 if (fabs (A->header->wesn[YHI] - B->header->wesn[YLO]) < y_noise) { /* B is exactly on top of A */
292 way = 1;
293 C->header->n_rows = A->header->n_rows + B->header->n_rows - one_or_zero;
294 C->header->wesn[YHI] = B->header->wesn[YHI]; /* ...but not for north */
295 }
296 else if (fabs (A->header->wesn[YLO] - B->header->wesn[YHI]) < y_noise) { /* A is exactly on top of B */
297 way = 2;
298 C->header->n_rows = A->header->n_rows + B->header->n_rows - one_or_zero;
299 C->header->wesn[YLO] = B->header->wesn[YLO]; /* ...but not for south */
300 }
301 else if ((fabs (A->header->wesn[YHI] - B->header->wesn[YLO]) < (C->header->inc[GMT_Y] + y_noise)) ) {
302 /* B is on top of A but their pixel|grid reg limits under|overlap by one cell */
303 if (one_or_zero) /* Grid registration - underlap */
304 way = 10;
305 else /* Pixel registration - overlap */
306 way = 11;
307 C->header->n_rows = A->header->n_rows + B->header->n_rows - !one_or_zero;
308 C->header->wesn[YHI] = B->header->wesn[YHI]; /* ...but not for north */
309 }
310 else if ((fabs (A->header->wesn[YLO] - B->header->wesn[YHI]) < (C->header->inc[GMT_Y] + y_noise)) ) {
311 /* A is on top of B but their pixel|grid reg limits under|overlap by one cell */
312 if (one_or_zero) /* Grid registration - underlap */
313 way = 21;
314 else /* Pixel registration - overlap */
315 way = 22;
316 C->header->n_rows = A->header->n_rows + B->header->n_rows - !one_or_zero;
317 C->header->wesn[YLO] = B->header->wesn[YLO]; /* ...but not for south */
318 }
319 else {
320 GMT_Report (API, GMT_MSG_ERROR, "Grids do not share a common edge!\n");
321 Return (GMT_RUNTIME_ERROR);
322 }
323 }
324 else {
325 GMT_Report (API, GMT_MSG_ERROR, "Grids do not share a common edge!\n");
326 Return (GMT_RUNTIME_ERROR);
327 }
328 if (gmt_M_x_is_lon (GMT, GMT_IN) && C->header->wesn[XHI] > 360.0) { /* Take out 360 */
329 C->header->wesn[XLO] -= 360.0;
330 C->header->wesn[XHI] -= 360.0;
331 }
332
333 /* Now we can do it */
334
335 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
336 sprintf (format, "%%s\t%s\t%s\t%s\t%s\t%s\t%s\t%%d\t%%d\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
337 GMT_Report (API, GMT_MSG_INFORMATION, "File\tW\tE\tS\tN\tdx\tdy\tnx\tny\n");
338 GMT_Report (API, GMT_MSG_INFORMATION, format, Ctrl->In.file[0], A->header->wesn[XLO], A->header->wesn[XHI], A->header->wesn[YLO], A->header->wesn[YHI], A->header->inc[GMT_X], A->header->inc[GMT_Y], A->header->n_columns, A->header->n_rows);
339 GMT_Report (API, GMT_MSG_INFORMATION, format, Ctrl->In.file[1], B->header->wesn[XLO], B->header->wesn[XHI], B->header->wesn[YLO], B->header->wesn[YHI], B->header->inc[GMT_X], B->header->inc[GMT_Y], B->header->n_columns, B->header->n_rows);
340 GMT_Report (API, GMT_MSG_INFORMATION, format, Ctrl->G.file, C->header->wesn[XLO], C->header->wesn[XHI], C->header->wesn[YLO], C->header->wesn[YHI], C->header->inc[GMT_X], C->header->inc[GMT_Y], C->header->n_columns, C->header->n_rows);
341 }
342
343 gmt_set_grddim (GMT, C->header);
344 if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL,
345 NULL, 0, 0, C) == NULL) Return (API->error); /* Note: 0 for pad since no BC work needed */
346 A->data = B->data = C->data; /* A and B share the same final matrix declared for C */
347 A->header->size = B->header->size = C->header->size; /* Set A & B's size to the same as C */
348 AH->no_BC = BH->no_BC = true; /* We must disable the BC machinery */
349
350 switch (way) { /* How A and B are positioned relative to each other */
351 case 1: /* B is on top of A */
352 case 10: /* B is on top of A but their grid reg limits underlap by one cell */
353 case 11: /* B is on top of A but their pixel reg limits overlap by one cell */
354 if (grdpaste_is_nc_grid(A)) {
355 AH->data_offset = B->header->n_columns * (B->header->n_rows - one_or_zero);
356 if (way == 11)
357 AH->data_offset -= B->header->n_columns;
358 else if (way == 10)
359 AH->data_offset += B->header->n_columns;
360 }
361 else {
362 GMT->current.io.pad[YHI] = B->header->n_rows - one_or_zero;
363 if (way == 11)
364 GMT->current.io.pad[YHI]--;
365 else if (way == 10)
366 GMT->current.io.pad[YHI]++;
367 gmt_M_grd_setpad (GMT, A->header, GMT->current.io.pad);
368 }
369 A->header->my = C->header->my; /* Needed if grid is read by gmt_gdal_read_grd */
370 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[0], A) == NULL) { /* Get data from A */
371 Return (API->error);
372 }
373 if (grdpaste_is_nc_grid(B)) {
374 gmt_set_pad (GMT, 0U); /* Reset padding */
375 }
376 else {
377 GMT->current.io.pad[YHI] = 0;
378 GMT->current.io.pad[YLO] = A->header->n_rows - one_or_zero;
379 if (way == 11)
380 GMT->current.io.pad[YLO]--;
381 else if (way == 10)
382 GMT->current.io.pad[YLO]++;
383 gmt_M_grd_setpad (GMT, B->header, GMT->current.io.pad);
384 }
385 B->header->my = C->header->my; /* Needed if grid is read by gmt_gdal_read_grd */
386 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[1], B) == NULL) { /* Get data from B */
387 Return (API->error);
388 }
389 break;
390 case 2: /* A is on top of B */
391 case 21: /* A is on top of B but their grid reg limits underlap by one cell */
392 case 22: /* A is on top of B but their pixel reg limits overlap by one cell */
393 if (!grdpaste_is_nc_grid(A)) {
394 GMT->current.io.pad[YLO] = B->header->n_rows - one_or_zero;
395 if (way == 22)
396 GMT->current.io.pad[YLO]--;
397 else if (way == 21)
398 GMT->current.io.pad[YLO]++;
399 gmt_M_grd_setpad (GMT, A->header, GMT->current.io.pad);
400 }
401 A->header->my = C->header->my; /* Needed if grid is read by gmt_gdal_read_grd */
402 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[0], A) == NULL) { /* Get data from A */
403 Return (API->error);
404 }
405 if (grdpaste_is_nc_grid(B)) {
406 gmt_set_pad (GMT, 0U); /* Reset padding */
407 BH->data_offset = A->header->n_columns * (A->header->n_rows - one_or_zero);
408 if (way == 22)
409 BH->data_offset -= A->header->n_columns;
410 else if (way == 21)
411 BH->data_offset += A->header->n_columns;
412 }
413 else {
414 GMT->current.io.pad[YLO] = 0;
415 GMT->current.io.pad[YHI] = A->header->n_rows - one_or_zero;
416 if (way == 22)
417 GMT->current.io.pad[YHI]--;
418 else if (way == 21)
419 GMT->current.io.pad[YHI]++;
420 gmt_M_grd_setpad (GMT, B->header, GMT->current.io.pad);
421 }
422 B->header->my = C->header->my; /* Needed if grid is read by gmt_gdal_read_grd */
423 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[1], B) == NULL) { /* Get data from B */
424 Return (API->error);
425 }
426 break;
427 case 3: /* A is on the right of B */
428 case 32: /* A is on right of B but their grid reg limits underlap by one cell */
429 case 33: /* A is on right of B but their pixel reg limits overlap by one cell */
430 if (grdpaste_is_nc_grid(A)) {
431 AH->stride = C->header->n_columns;
432 AH->data_offset = B->header->n_columns - one_or_zero;
433 if (way == 33)
434 AH->data_offset--;
435 else if (way == 32)
436 AH->data_offset++;
437 }
438 else {
439 GMT->current.io.pad[XLO] = B->header->n_columns - one_or_zero;
440 if (way == 33)
441 GMT->current.io.pad[XLO]--;
442 else if (way == 32)
443 GMT->current.io.pad[XLO]++;
444 gmt_M_grd_setpad (GMT, A->header, GMT->current.io.pad);
445 }
446 A->header->mx = C->header->mx; /* Needed if grid is read by gmt_gdal_read_grd */
447 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[0], A) == NULL) { /* Get data from A */
448 Return (API->error);
449 }
450 if (grdpaste_is_nc_grid(B)) {
451 gmt_set_pad (GMT, 0U); /* Reset padding */
452 BH->stride = C->header->n_columns;
453 }
454 else {
455 GMT->current.io.pad[XLO] = 0; GMT->current.io.pad[XHI] = A->header->n_columns - one_or_zero;
456 if (way == 33)
457 GMT->current.io.pad[XHI]--;
458 else if (way == 32)
459 GMT->current.io.pad[XHI]++;
460 gmt_M_grd_setpad (GMT, B->header, GMT->current.io.pad);
461 }
462 B->header->mx = C->header->mx; /* Needed if grid is read by gmt_gdal_read_grd */
463 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[1], B) == NULL) { /* Get data from B */
464 Return (API->error);
465 }
466 break;
467 case 4: /* A is on the left of B */
468 case 43: /* A is on left of B but their grid reg limits underlap by one cell */
469 case 44: /* A is on left of B but their pixel reg limits overlap by one cell */
470 if (grdpaste_is_nc_grid(A)) {
471 AH->stride = C->header->n_columns;
472 }
473 else {
474 GMT->current.io.pad[XHI] = B->header->n_columns - one_or_zero;
475 if (way == 44)
476 GMT->current.io.pad[XHI]--;
477 else if (way == 43)
478 GMT->current.io.pad[XHI]++;
479 gmt_M_grd_setpad (GMT, A->header, GMT->current.io.pad);
480 }
481 A->header->mx = C->header->mx; /* Needed if grid is read by gmt_gdal_read_grd */
482 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[0], A) == NULL) { /* Get data from A */
483 Return (API->error);
484 }
485 if (grdpaste_is_nc_grid(B)) {
486 gmt_set_pad (GMT, 0U); /* Reset padding */
487 BH->stride = C->header->n_columns;
488 BH->data_offset = A->header->n_columns - one_or_zero;
489 if (way == 44)
490 BH->data_offset--;
491 else if (way == 43)
492 BH->data_offset++;
493 }
494 else {
495 GMT->current.io.pad[XHI] = 0;
496 GMT->current.io.pad[XLO] = A->header->n_columns - one_or_zero;
497 if (way == 44)
498 GMT->current.io.pad[XLO]--;
499 else if (way == 43)
500 GMT->current.io.pad[XLO]++;
501 gmt_M_grd_setpad (GMT, B->header, GMT->current.io.pad);
502 }
503 B->header->mx = C->header->mx; /* Needed if grid is read by gmt_gdal_read_grd */
504 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file[1], B) == NULL) { /* Get data from B */
505 Return (API->error);
506 }
507 break;
508 }
509
510 gmt_set_pad (GMT, 0U); /* Reset padding */
511 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, C)) Return (API->error);
512 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, C) != GMT_NOERROR) {
513 Return (API->error);
514 }
515 A->data = B->data = NULL; /* Since these were never actually allocated */
516
517 gmt_set_pad (GMT, API->pad); /* Restore to GMT Defaults */
518 Return (GMT_NOERROR);
519 }
520