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: grdfill.c reads a grid that has one or more holes in
19 * it flagged by NaNs and interpolates across the holes given a choice
20 * of algorithm.
21 *
22 * Author: Paul Wessel
23 * Date: 7-JAN-2018
24 * Version: 6 API
25 */
26
27 #include "gmt_dev.h"
28
29 #define THIS_MODULE_CLASSIC_NAME "grdfill"
30 #define THIS_MODULE_MODERN_NAME "grdfill"
31 #define THIS_MODULE_LIB "core"
32 #define THIS_MODULE_PURPOSE "Interpolate across holes in a grid"
33 #define THIS_MODULE_KEYS "<G{,>?}"
34 #define THIS_MODULE_NEEDS ""
35 #define THIS_MODULE_OPTIONS "-RVf"
36
37 enum GRDFILL_mode {
38 ALG_CONSTANT,
39 ALG_NN,
40 ALG_SPLINE
41 };
42
43 struct GRDFILL_CTRL {
44 struct GRDFILL_In {
45 bool active;
46 char *file;
47 } In;
48 struct GRDFILL_A { /* -A<algo>[<options>] */
49 bool active;
50 unsigned int mode;
51 double value;
52 } A;
53 struct GRDFILL_G { /* -G<outgrid> */
54 bool active;
55 char *file;
56 } G;
57 struct GRDFILL_L { /* -L[p] */
58 bool active;
59 unsigned int mode; /* 0 = region, 1 = polygons */
60 } L;
61 struct GRDFILL_N { /* -N<value> */
62 bool active;
63 float value;
64 } N;
65 };
66
New_Ctrl(struct GMT_CTRL * GMT)67 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
68 struct GRDFILL_CTRL *C;
69
70 C = gmt_M_memory (GMT, NULL, 1, struct GRDFILL_CTRL);
71
72 /* Initialize values whose defaults are not 0/false/NULL */
73
74 return (C);
75 }
76
Free_Ctrl(struct GMT_CTRL * GMT,struct GRDFILL_CTRL * C)77 static void Free_Ctrl (struct GMT_CTRL *GMT, struct GRDFILL_CTRL *C) { /* Deallocate control structure */
78 if (!C) return;
79 gmt_M_str_free (C->In.file);
80 gmt_M_str_free (C->G.file);
81 gmt_M_free (GMT, C);
82 }
83
usage(struct GMTAPI_CTRL * API,int level)84 static int usage (struct GMTAPI_CTRL *API, int level) {
85 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
86 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
87 GMT_Usage (API, 0, "usage: %s %s [-Ac|n|s[<arg>]] [-G%s] [-L[p]] [-N<value>] [%s] [%s] [%s] [%s]\n",
88 name, GMT_INGRID, GMT_OUTGRID, GMT_Rgeo_OPT, GMT_V_OPT, GMT_f_OPT, GMT_PAR_OPT);
89
90 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
91
92 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
93 gmt_ingrid_syntax (API, 0, "Name of grid with NaN holes");
94 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
95 GMT_Usage (API, 1, "\n-Ac|n|s[<arg>]");
96 GMT_Usage (API, -2, "Specify algorithm and any required parameters for in-fill:");
97 GMT_Usage (API, 3, "c: Fill in NaN holes with the given constant <value>.");
98 GMT_Usage (API, 3, "n: Fill in NaN holes with nearest neighbor values; "
99 "append <max_radius> nodes for the outward search "
100 "[Default radius is sqrt(nx^2+ny^2), with (nx,ny) the dimensions of the grid].");
101 GMT_Usage (API, 3, "s: Fill in NaN holes with a spline (optionally append tension).");
102 GMT_Usage (API, -2, "Note: -A is required unless -L is used.");
103 gmt_outgrid_syntax (API, 'G', "Give filename for where to write the filled-in grid");
104 GMT_Usage (API, 1, "\n-L[p]");
105 GMT_Usage (API, -2, "Just list the sub-regions w/e/s/n of each hole. "
106 "No grid fill takes place and -G is ignored. "
107 "Optionally, append p to write polygons corresponding to these regions.");
108 GMT_Usage (API, 1, "\n-N<value>");
109 GMT_Usage (API, -2, "Set alternate node <value> to indicate a hole [Default looks for NaN-nodes].");
110 GMT_Option (API, "R,V,f,.");
111
112 return (GMT_MODULE_USAGE);
113 }
114
parse(struct GMT_CTRL * GMT,struct GRDFILL_CTRL * Ctrl,struct GMT_OPTION * options)115 static int parse (struct GMT_CTRL *GMT, struct GRDFILL_CTRL *Ctrl, struct GMT_OPTION *options) {
116 /* This parses the options provided to grdfill and sets parameters in CTRL.
117 * Any GMT common options will override values set previously by other commands.
118 * It also replaces any file names specified as input or output with the data ID
119 * returned when registering these sources/destinations with the API.
120 */
121
122 unsigned int n_errors = 0, n_files = 0;
123 struct GMT_OPTION *opt = NULL;
124 struct GMTAPI_CTRL *API = GMT->parent;
125
126 for (opt = options; opt; opt = opt->next) {
127 switch (opt->option) {
128
129 case '<': /* Input and Output files */
130 if (n_files++ > 0) {n_errors++; continue; }
131 Ctrl->In.active = true;
132 if (opt->arg[0]) Ctrl->In.file = strdup (opt->arg);
133 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file))) n_errors++;
134 break;
135 case '>': /* Output file */
136 Ctrl->G.active = true;
137 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
138 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
139 break;
140
141 /* Processes program-specific parameters */
142
143 case 'A':
144 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
145 Ctrl->A.active = true;
146 switch (opt->arg[0]) {
147 case 'c': /* Constant fill */
148 Ctrl->A.mode = ALG_CONSTANT;
149 Ctrl->A.value = atof (&opt->arg[1]);
150 break;
151 case 'n': /* Nearest neighbor fill */
152 Ctrl->A.mode = ALG_NN;
153 Ctrl->A.value = (opt->arg[1]) ? atof (&opt->arg[1]) : -1;
154 break;
155 case 's': /* Spline fill */
156 Ctrl->A.mode = ALG_SPLINE;
157 if (opt->arg[1]) Ctrl->A.value = atof (&opt->arg[1]);
158 break;
159 default:
160 GMT_Report (API, GMT_MSG_ERROR, "Option -A: Unrecognized algorithm (%c)\n", opt->arg[0]);
161 n_errors++;
162 }
163 break;
164
165 case 'G':
166 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
167 Ctrl->G.active = true;
168 if (Ctrl->G.file) {
169 GMT_Report (API, GMT_MSG_ERROR, "Specify only one output file\n");
170 n_errors++;
171 }
172 else {
173 Ctrl->G.file = strdup (opt->arg);
174 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
175 }
176 break;
177
178 case 'L':
179 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active);
180 Ctrl->L.active = true;
181 if (opt->arg[0] == 'p')
182 Ctrl->L.mode = 1;
183 break;
184
185 case 'N':
186 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
187 Ctrl->N.active = true;
188 if (opt->arg[0] && !strchr ("Nn", opt->arg[0]))
189 Ctrl->N.value = (float) atof (opt->arg);
190 else {
191 GMT_Report (API, GMT_MSG_ERROR, "Option -N: No value (or NaN) given\n");
192 n_errors++;
193 }
194 break;
195
196 default: /* Report bad options */
197 n_errors += gmt_default_error (GMT, opt->option);
198 break;
199 }
200 }
201
202 n_errors += gmt_M_check_condition (GMT, !Ctrl->In.file, "Must specify input grid file\n");
203 n_errors += gmt_M_check_condition (GMT, !(Ctrl->A.active || Ctrl->L.active), "Must specify an algorithm with -A unless -L is used\n");
204 n_errors += gmt_M_check_condition (GMT, !(Ctrl->L.active || Ctrl->G.file), "Must specify output grid file\n");
205 n_errors += gmt_M_check_condition (GMT, !(Ctrl->A.active || Ctrl->L.active), "Must specify either -A or -L\n");
206
207 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
208 }
209
grdfill_do_constant_fill(struct GMT_GRID * G,unsigned int limit[],gmt_grdfloat value)210 GMT_LOCAL int grdfill_do_constant_fill (struct GMT_GRID *G, unsigned int limit[], gmt_grdfloat value) {
211 /* Algorithm 1: Replace NaNs with a constant value */
212 uint64_t node;
213
214 for (unsigned int row = limit[YLO]; row <= limit[YHI]; row++) {
215 for (unsigned int col = limit[XLO]; col <= limit[XHI]; col++) {
216 node = gmt_M_ijp (G->header, row, col);
217 if (gmt_M_is_fnan (G->data[node]))
218 G->data[node] = value;
219 }
220 }
221 return GMT_NOERROR;
222 }
223
224 #if 1
grdfill_do_splinefill(struct GMTAPI_CTRL * API,struct GMT_GRID * G,double wesn[],unsigned int limit[],unsigned int n_in_hole,double value)225 GMT_LOCAL int grdfill_do_splinefill (struct GMTAPI_CTRL *API, struct GMT_GRID *G, double wesn[], unsigned int limit[], unsigned int n_in_hole, double value) {
226 /* Algorithm 2: Replace NaNs with a spline */
227 char input[GMT_VF_LEN] = {""}, output[GMT_VF_LEN] = {""}, args[GMT_LEN256] = {""}, method[GMT_LEN32] = {""};
228 unsigned int row, col, row_hole, col_hole, mode, d_limit[4], n_constraints;
229 uint64_t node, node_hole, k = 0, dim[GMT_DIM_SIZE] = {0, 0, 0, 0};
230 double *x = NULL, *y = NULL;
231 gmt_grdfloat *z = NULL;
232 struct GMT_VECTOR *V = NULL;
233 struct GMT_GRID *G_hole = NULL;
234 struct GMT_CTRL *GMT = API->GMT;
235
236 /* Allocate a vector container for input to greenspline */
237 dim[0] = 3; /* Want three input columns but let length be 0 - this signals that no vector allocations should take place */
238 if ((V = GMT_Create_Data (API, GMT_IS_VECTOR, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
239 return (API->error);
240 }
241 /* Create a virtual file to hold the resampled grid */
242 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, output) == GMT_NOTSET) {
243 return (API->error);
244 }
245 /* Add up to 2 rows/cols around hole, but watch for grid edges */
246 gmt_M_memcpy (d_limit, limit, 4, unsigned int); /* d_limit will be used to set the grid domain */
247 /* Undo the half grid inc padding done in the main */
248 wesn[XLO] += 0.5 * G->header->inc[GMT_X];
249 wesn[XHI] -= 0.5 * G->header->inc[GMT_X];
250 wesn[YLO] += 0.5 * G->header->inc[GMT_Y];
251 wesn[YHI] -= 0.5 * G->header->inc[GMT_Y];
252 for (k = 1; k <= 2; k++) {
253 if (d_limit[XLO]) d_limit[XLO]--, wesn[XLO] -= G->header->inc[GMT_X]; /* Move one column westward */
254 if (d_limit[XHI] < (G->header->n_columns-1)) d_limit[XHI]++, wesn[XHI] += G->header->inc[GMT_X]; /* Move one column eastward */
255 if (d_limit[YLO]) d_limit[YLO]--, wesn[YHI] += G->header->inc[GMT_Y]; /* Move one row northward */
256 if (d_limit[YHI] < (G->header->n_rows-1)) d_limit[YHI]++, wesn[YLO] -= G->header->inc[GMT_Y]; /* Move one row southward */
257 }
258 n_constraints = (d_limit[YHI] - d_limit[YLO] + 1) * (d_limit[XHI] - d_limit[XLO] + 1) - n_in_hole;
259 x = gmt_M_memory (GMT, NULL, n_constraints, double);
260 y = gmt_M_memory (GMT, NULL, n_constraints, double);
261 z = gmt_M_memory (GMT, NULL, n_constraints, gmt_grdfloat);
262 for (row = d_limit[YLO], k = 0; row <= d_limit[YHI]; row++) {
263 for (col = d_limit[XLO]; col <= d_limit[XHI]; col++) {
264 node = gmt_M_ijp (G->header, row, col);
265 if (gmt_M_is_fnan (G->data[node])) continue;
266 x[k] = gmt_M_grd_col_to_x (GMT, col, G->header);
267 y[k] = gmt_M_grd_row_to_y (GMT, row, G->header);
268 z[k] = G->data[node];
269 k++;
270 }
271 }
272 V->n_rows = n_constraints; /* Must specify how many input points we have */
273 GMT_Put_Vector (API, V, GMT_X, GMT_DOUBLE, x);
274 GMT_Put_Vector (API, V, GMT_Y, GMT_DOUBLE, y);
275 GMT_Put_Vector (API, V, GMT_Z, GMT_FLOAT, z);
276 /* Associate our input data vectors with a virtual input file */
277 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET|GMT_VIA_VECTOR, GMT_IS_POINT, GMT_IN|GMT_IS_REFERENCE, V, input) == GMT_NOTSET)
278 return (API->error);
279 /* Prepare the greenspline command-line arguments */
280 mode = (gmt_M_is_geographic (GMT, GMT_IN)) ? 2 : 1;
281 if (value > 0.0)
282 sprintf (method, "t%g", value);
283 else
284 sprintf (method, "c");
285 sprintf (args, "%s -G%s -S%s -R%.16g/%.16g/%.16g/%.16g -I%.16g/%.16g -D%d", input, output, method, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], G->header->inc[GMT_X], G->header->inc[GMT_Y], mode);
286 if (G->header->registration == GMT_GRID_PIXEL_REG) strcat (args, " -r");
287 strcat (args, " --GMT_HISTORY=readonly");
288 /* Run the greenspline module */
289 GMT_Report (API, GMT_MSG_INFORMATION, "Calling greenspline with args %s\n", args);
290 if (GMT_Call_Module (API, "greenspline", GMT_MODULE_CMD, args)) {
291 return (API->error); /* Some sort of failure */
292 }
293 if ((G_hole = GMT_Read_VirtualFile (API, output)) == NULL) { /* Load in the resampled grid */
294 return (API->error); /* Some sort of failure */
295 }
296 /* Use 0-nx,0-ny for the hole index and the large grid G nodes for the holes */
297 for (row_hole = 0, row = d_limit[YLO]; row_hole < G_hole->header->n_rows; row_hole++, row++) {
298 for (col_hole = 0, col = d_limit[XLO]; col_hole < G_hole->header->n_columns; col_hole++, col++) {
299 node = gmt_M_ijp (G->header, row, col);
300 if (!gmt_M_is_fnan (G->data[node])) continue; /* Had data value */
301 node_hole = gmt_M_ijp (G_hole->header, row_hole, col_hole);
302 G->data[node] = G_hole->data[node_hole]; /* Replace the NaN with splined value */
303 }
304 }
305
306 /* Close the two virtual files */
307 if (GMT_Close_VirtualFile (API, output)) {
308 GMT_Report (API, GMT_MSG_ERROR, "Failed to close virtual output file %s\n", output);
309 return (API->error);
310 }
311 if (GMT_Close_VirtualFile (API, input)) {
312 GMT_Report (API, GMT_MSG_ERROR, "Failed to close virtual input file %s\n", input);
313 return (API->error);
314 }
315
316 /* Free our custom vectors */
317 gmt_M_free (GMT, x); gmt_M_free (GMT, y); gmt_M_free (GMT, z);
318 if (GMT_Destroy_Data (API, &G_hole) != GMT_NOERROR) {
319 GMT_Report (API, GMT_MSG_ERROR, "Failed to destroy temporary hold grid\n");
320 return (API->error);
321 }
322 return (GMT_NOERROR);
323 }
324 #endif
325
grdfill_trace_the_hole(struct GMT_GRID * G,uint64_t node,unsigned int row,unsigned int col,int64_t * step,char * ID,unsigned int limit[])326 GMT_LOCAL unsigned int grdfill_trace_the_hole (struct GMT_GRID *G, uint64_t node, unsigned int row, unsigned int col, int64_t *step, char *ID, unsigned int limit[]) {
327 /* Determine all the direct neighbor nodes in the W/E/S/N directions that are also NaN, recursively.
328 * This is a limited form of Moore neighborhood called Von Neumann neighborhood since we only do the
329 * four cardinal directions and ignore the diagonals. Ignoring the diagonal means two holes that
330 * touch at a diagonal will be encoded as two separate holes whereas a full Moore neighborhood
331 * calculation would find a single hole. For our purposes (filling in the hole), smaller holes
332 * are less computationally intensive, hence we limit ourselves to good old Von Neumann. */
333 static int drow[4] = {1, 0, -1, 0}, dcol[4] = {0, 1, 0, -1}; /* Change in row,col per cardinal direction */
334 unsigned int side, next_row, next_rcol, n_nodes = 0;
335 int64_t ij;
336
337 for (side = 0; side < 4; side++) { /* For each of the 4 cardinal directions */
338 ij = node + step[side]; /* This is the grid node (and ID node) of this cardinal point */
339 if (ID[ij] == 0 && gmt_M_is_fnan (G->data[ij])) { /* Hole extends in this direction, follow it to its next neighbor */
340 next_row = (unsigned int)(row + drow[side]); /* Set col and row of next point */
341 next_rcol = (unsigned int)(col + dcol[side]);
342 ID[ij] = 1; /* Mark this node as part of the current hole */
343 if (next_rcol < limit[XLO]) limit[XLO] = next_rcol; /* Update min/max row and col values */
344 else if (next_rcol > limit[XHI]) limit[XHI] = next_rcol;
345 if (next_row < limit[YLO]) limit[YLO] = next_row;
346 else if (next_row > limit[YHI]) limit[YHI] = next_row;
347 /* Recursively trace this nodes next neighbors as well */
348 n_nodes = n_nodes + 1 + grdfill_trace_the_hole (G, ij, next_row, next_rcol, step, ID, limit);
349 }
350 }
351 return (n_nodes);
352 }
353
grdfill_find_nearest(int64_t i,int64_t j,int64_t * r2,int64_t * is,int64_t * js,int64_t * xs,int64_t * ys)354 GMT_LOCAL int64_t grdfill_find_nearest (int64_t i, int64_t j, int64_t *r2, int64_t *is, int64_t *js, int64_t *xs, int64_t *ys) {
355 /* function to find the nearest point based on previous search, smallest distance outside a radius */
356 int64_t ct = 0, nx, ny, nx1, ii, k = 0, rr, nx1_2, ny_2;
357
358 rr = INTMAX_MAX; /* Ensure we reset this the first time */
359
360 /* starting with nx = ny */
361 nx = (int64_t)(sqrt((double)(*r2) / 2.0));
362 /* loop over possible nx, find smallest rr */
363 for (nx1 = nx; nx1 <= (int64_t)sqrt((double)(*r2)) + 1; nx1++) {
364 if ((nx1_2 = (nx1 * nx1)) < (*r2))
365 ny = (int64_t)(sqrt((double)((*r2) - nx1_2)));
366 else
367 ny = 0;
368 while ((nx1_2 + ny * ny ) <= (*r2) && ny <= nx1)
369 ny++;
370 ny_2 = ny * ny;
371 if (ny <= nx1) {
372 if (rr > (nx1_2 + ny_2)) {
373 k = 0;
374 rr = nx1_2 + ny_2;
375 xs[0] = nx1;
376 ys[0] = ny;
377 }
378 else if (rr == (nx1_2 + ny_2)) {
379 k++;
380 xs[k] = nx1;
381 ys[k] = ny;
382 }
383 }
384 }
385
386 /* Return the grid index, changing the order may lead to different solution, but mainly because nearest-neighbor is non-unique */
387 for (ii = 0; ii <= k; ii++) {
388 nx = xs[ii];
389 ny = ys[ii];
390
391 if (ny == 0) {
392 js[ct] = 0; is[ct++] = nx;
393 js[ct] = 0; is[ct++] = -nx;
394 js[ct] = nx; is[ct++] = 0;
395 js[ct] = -nx; is[ct++] = 0;
396 }
397 else if (nx != ny){
398 js[ct] = ny; is[ct++] = nx;
399 js[ct] = ny; is[ct++] = -nx;
400 js[ct] = -ny; is[ct++] = nx;
401 js[ct] = -ny; is[ct++] = -nx;
402 js[ct] = nx; is[ct++] = ny;
403 js[ct] = -nx; is[ct++] = ny;
404 js[ct] = nx; is[ct++] = -ny;
405 js[ct] = -nx; is[ct++] = -ny;
406 }
407 else {
408 js[ct] = nx; is[ct++] = nx;
409 js[ct] = nx; is[ct++] = -nx;
410 js[ct] = -nx; is[ct++] = nx;
411 js[ct] = -nx; is[ct++] = -nx;
412 }
413 }
414
415 for (ii = 0; ii < ct; ii++) { /* Set absolute row/col values */
416 is[ii] = is[ii] + i;
417 js[ii] = js[ii] + j;
418 }
419
420 *r2 = rr; /* Update radius */
421
422 /* return the number of indexes of the nearest neighbor */
423 return (ct);
424 }
425
grdfill_nearest_interp(struct GMT_CTRL * GMT,struct GMT_GRID * In,struct GMT_GRID * Out,int64_t radius)426 GMT_LOCAL void grdfill_nearest_interp (struct GMT_CTRL *GMT, struct GMT_GRID *In, struct GMT_GRID *Out, int64_t radius) {
427 uint64_t ij, node;
428 int64_t nx = In->header->n_columns, ny = In->header->n_rows;
429 int64_t i, j, flag, ct, k, recx = 1, recy = 1, cs = 0, rr;
430 int64_t *is = NULL, *js = NULL, *xs = NULL, *ys = NULL;
431 gmt_grdfloat *m = In->data, *m_interp = Out->data; /* Short-hand for input and output grids */
432 double rad2;
433
434 /* Allocate memory for temporary indexes */
435 is = gmt_M_memory (GMT, NULL, 2048, int64_t);
436 js = gmt_M_memory (GMT, NULL, 2048, int64_t);
437 xs = gmt_M_memory (GMT, NULL, 512, int64_t);
438 ys = gmt_M_memory (GMT, NULL, 512, int64_t);
439
440 if (radius == -1) /* Set default radius */
441 radius = (int64_t)floor (sqrt ((double)(nx*nx + ny*ny)));
442 rad2 = (double)(radius * radius);
443
444 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Interpolating to nearest neighbor...\n");
445 gmt_M_row_loop (GMT, In, i) { /* Loop over each row in grid */
446 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Working on row %" PRIi64 "\n", i);
447 rr = 0;
448 gmt_M_col_loop (GMT, In, i, j, ij) { /* Loop over all columns */
449 if (!gmt_M_is_fnan (m[ij])) /* Already duplicated in the calling program */
450 rr = 0;
451 else { /* search nearest neighbor, use previous nearest distance to exclude certain search area. */
452 flag = 0;
453 /* set the starting search radius based on last nearest distance */
454 if (rr >= 4) {
455 if (recy > 0 && recx > 0)
456 rr = (recx-1)*(recx-1)+(recy-1)*(recy-1)-1;
457 else if (recy == 0 && recx > 0)
458 rr = (recx-1)*(recx-1)-1;
459 else if (recy > 0 && recx == 0)
460 rr = (recy-1)*(recy-1)-1;
461 }
462 else
463 rr = 0;
464
465 while (flag == 0 && rr <= rad2) {
466 ct = grdfill_find_nearest (i, j, &rr, is, js, xs, ys);
467 cs++;
468 for (k = 0; k < ct; k++) {
469 if (is[k] >= 0 && is[k] < ny && js[k] >=0 && js[k] < nx) {
470 node = gmt_M_ijp (In->header, is[k], js[k]);
471 if (!gmt_M_is_fnan (m[node])) {
472 m_interp[ij] = m[node];
473 flag = 1;
474 recx = int64_abs (is[k]-i);
475 recy = int64_abs (js[k]-j);
476 break;
477 }
478 }
479 }
480 }
481 }
482 if (i == 0 && rr == -1)
483 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "(%d %d %d %d)\n", j, rr, recx, recy);
484 }
485 }
486
487 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%" PRIi64 " number of searches used\n", cs);
488
489 gmt_M_free (GMT, is);
490 gmt_M_free (GMT, js);
491 gmt_M_free (GMT, xs);
492 gmt_M_free (GMT, ys);
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
GMT_grdfill(void * V_API,int mode,void * args)498 EXTERN_MSC int GMT_grdfill (void *V_API, int mode, void *args) {
499 char *ID = NULL, *RG_orig_hist = NULL;
500 int error = 0, RG_id = 0;
501 unsigned int hole_number = 0, row, col, limit[4], n_nodes;
502 uint64_t node, offset;
503 int64_t off[4];
504 double wesn[4];
505 struct GMT_GRID *Grid = NULL;
506 struct GMT_RECORD *Out = NULL;
507 struct GRDFILL_CTRL *Ctrl = NULL;
508 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
509 struct GMT_OPTION *options = NULL;
510 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
511
512 /*----------------------- Standard module initialization and parsing ----------------------*/
513
514 if (API == NULL) return (GMT_NOT_A_SESSION);
515 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
516 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
517
518 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
519
520 /* Parse the command-line arguments */
521
522 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 */
523 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
524 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
525 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
526
527 /*---------------------------- This is the grdfill main code ----------------------------*/
528
529 if ((Grid = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, Ctrl->In.file, NULL)) == NULL) { /* Get header only */
530 Return (API->error);
531 }
532
533 if (GMT->common.R.active[RSET]) { /* Specified a subset */
534 bool global = false;
535 global = gmt_grd_is_global (GMT, Grid->header);
536 if (!global && (GMT->common.R.wesn[XLO] < Grid->header->wesn[XLO] || GMT->common.R.wesn[XHI] > Grid->header->wesn[XHI])) error++;
537 if (GMT->common.R.wesn[YLO] < Grid->header->wesn[YLO] || GMT->common.R.wesn[YHI] > Grid->header->wesn[YHI]) error++;
538 if (error) {
539 GMT_Report (API, GMT_MSG_ERROR, "Subset exceeds data domain!\n");
540 Return (GMT_RUNTIME_ERROR);
541 }
542 if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, GMT->common.R.wesn, Ctrl->In.file, Grid) == NULL) {
543 Return (API->error); /* Get subset */
544 }
545 }
546 else if (GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, Ctrl->In.file, Grid) == NULL) {
547 Return (API->error); /* Get all */
548 }
549
550 if (Ctrl->N.active) { /* User wants a specific value to indicate a hole instead of NaN; replace those values with NaNs since that is what is expected below */
551 gmt_M_grd_loop (GMT, Grid, row, col, node) { /* Loop over all grid nodes */
552 if (floatAlmostEqualZero (Grid->data[node], Ctrl->N.value))
553 Grid->data[node] = GMT->session.f_NaN;
554 }
555 }
556 /* Here any hole is identified as a patch of NaNs */
557
558 if (Ctrl->A.mode == ALG_NN) { /* Do Eric Xu's NN algorithm and bail */
559 int64_t radius = lrint (Ctrl->A.value);
560 struct GMT_GRID *New = NULL;
561 if ((New = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_DATA, Grid)) == NULL) {
562 GMT_Report (API, GMT_MSG_ERROR, "Unable to duplicate input grid!\n");
563 Return (API->error);
564 }
565 grdfill_nearest_interp (GMT, Grid, New, radius); /* Perform the NN replacements */
566
567 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_GRID_ALL, NULL, Ctrl->G.file, New)) {
568 GMT_Report (API, GMT_MSG_ERROR, "Failed to write output grid!\n");
569 Return (API->error);
570 }
571 Return (GMT_NOERROR);
572 }
573
574 if (Ctrl->L.active) {
575 if (GMT_Init_IO (API, GMT_IS_DATASET, (Ctrl->L.mode) ? GMT_IS_POLYGON : GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) { /* Registers default output destination, unless already set */
576 Return (API->error);
577 }
578 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_OFF) != GMT_NOERROR) { /* Enables data output and sets access mode */
579 Return (API->error);
580 }
581 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) { /* Sets output geometry */
582 Return (API->error);
583 }
584 if ((error = GMT_Set_Columns (API, GMT_OUT, 2 + 2*(1-Ctrl->L.mode), GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
585 Return (API->error);
586 }
587 if (Ctrl->L.mode) gmt_set_segmentheader (GMT, GMT_OUT, true);
588 }
589
590 if (Ctrl->A.mode == ALG_SPLINE) { /* Must preserve the original -RG history which greenspline messes with */
591 RG_id = gmt_get_option_id (0, "R") + 1;
592 if (GMT->init.history[RG_id]) RG_orig_hist = strdup (GMT->init.history[RG_id]);
593 }
594
595 /* To avoid having to check every row,col for being inside the grid we set
596 * the boundary row/cols in the ID grid to 1. */
597
598 ID = gmt_M_memory_aligned (GMT, NULL, Grid->header->size, char);
599 /* Set the top and bottom boundary rows to UINT_MAX */
600 offset = (uint64_t)(Grid->header->pad[YHI] + Grid->header->n_rows) * Grid->header->mx;
601 for (node = 0; node < (uint64_t)Grid->header->pad[YHI]*Grid->header->mx; node++) ID[node] = ID[node+offset] = 1;
602 /* Set the left and right boundary columnss to UINT_MAX */
603 offset = Grid->header->pad[XLO] + Grid->header->n_columns;
604 for (row = 0; row < Grid->header->my; row++) {
605 for (col = 0; col < Grid->header->pad[XLO]; col++)
606 ID[row*Grid->header->mx+col] = 1;
607 for (col = 0; col < Grid->header->pad[XHI]; col++)
608 ID[row*Grid->header->mx+offset+col] = 1;
609 }
610 /* Initiate the node offsets in the cardinal directions */
611 off[0] = Grid->header->mx; off[1] = 1; off[2] = -off[0]; off[3] = -off[1];
612
613 Out = gmt_new_record (GMT, wesn, NULL); /* Since we only need to worry about numerics in this module */
614 gmt_M_grd_loop (GMT, Grid, row, col, node) { /* Loop over all grid nodes */
615 if (ID[node]) continue; /* Already identified as part of a hole */
616 if (gmt_M_is_fnan (Grid->data[node])) { /* Node is part of a new hole */
617 limit[XLO] = limit[XHI] = col; /* Initiate min/max col to this single node */
618 limit[YLO] = limit[YHI] = row; /* Initiate min/max row to this single node */
619 ID[node] = 1; /* Flag this node as part of a hole */
620 ++hole_number; /* Increase the current hole number */
621 /* Trace all the contiguous neighbors, updating the bounding box as we go along */
622 n_nodes = 1 + grdfill_trace_the_hole (Grid, node, row, col, off, ID, limit);
623 wesn[XLO] = gmt_M_grd_col_to_x (GMT, limit[XLO], Grid->header) - 0.5 * Grid->header->inc[GMT_X];
624 wesn[XHI] = gmt_M_grd_col_to_x (GMT, limit[XHI], Grid->header) + 0.5 * Grid->header->inc[GMT_X];
625 wesn[YLO] = gmt_M_grd_row_to_y (GMT, limit[YHI], Grid->header) - 0.5 * Grid->header->inc[GMT_Y];
626 wesn[YHI] = gmt_M_grd_row_to_y (GMT, limit[YLO], Grid->header) + 0.5 * Grid->header->inc[GMT_Y];
627 GMT_Report (API, GMT_MSG_INFORMATION, "Hole BB %u: -R: %g/%g/%g/%g [%u nodes]\n", hole_number, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], n_nodes);
628 if (Ctrl->L.active) {
629 if (Ctrl->L.mode) { /* Write a closed polygon */
630 double tmp[4];
631 gmt_M_memcpy (tmp, wesn, 4, double);
632 GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL);
633 wesn[GMT_X] = tmp[XLO]; wesn[GMT_Y] = tmp[YLO];
634 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
635 wesn[GMT_X] = tmp[XHI];
636 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
637 wesn[GMT_Y] = tmp[YHI];
638 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
639 wesn[GMT_X] = tmp[XLO];
640 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
641 wesn[GMT_Y] = tmp[YLO];
642 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
643 }
644 else /* Write west east south north limits */
645 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
646 }
647 else {
648 switch (Ctrl->A.mode) {
649 case ALG_CONSTANT: /* Fill in using a constant value */
650 error = grdfill_do_constant_fill (Grid, limit, (gmt_grdfloat)Ctrl->A.value);
651 break;
652 case ALG_SPLINE: /* Fill in using a spline */
653 error = grdfill_do_splinefill (API, Grid, wesn, limit, n_nodes, Ctrl->A.value);
654 break;
655 }
656 if (error) {
657 GMT_Report (API, GMT_MSG_INFORMATION, "Failed to fill hole %u\n", hole_number);
658 continue;
659 }
660 }
661 }
662 }
663 if (hole_number) GMT_Report (API, GMT_MSG_INFORMATION, "Found %u holes\n", hole_number);
664 gmt_M_free_aligned (GMT, ID);
665 gmt_M_free (GMT, Out);
666
667 if (Ctrl->L.active) {
668 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
669 free(RG_orig_hist);
670 Return (API->error);
671 }
672 }
673 else if (hole_number) { /* Must write the revised grid if there were any holes*/
674 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) {
675 free(RG_orig_hist);
676 Return (API->error);
677 }
678 if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
679 free(RG_orig_hist);
680 Return (API->error);
681 }
682 }
683 else {
684 GMT_Report (API, GMT_MSG_WARNING, "No holes detected in grid - grid was not updated\n");
685 }
686
687 if (Ctrl->A.mode == ALG_SPLINE) {
688 if (GMT->init.history[RG_id]) gmt_M_str_free (GMT->init.history[RG_id]);
689 if (RG_orig_hist) GMT->init.history[RG_id] = RG_orig_hist;
690 }
691
692 Return (GMT_NOERROR);
693 }
694