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 * surface.c: a gridding program using splines in tension.
19 * Reads xyz Cartesian triples and fits a surface to the data.
20 * The surface satisfies (1 - T) D4 z - T D2 z = 0,
21 * where D4 is the 2-D biharmonic operator, D2 is the
22 * 2-D Laplacian, and T is a "tension factor" between 0 and 1.
23 * End member T = 0 is the classical minimum curvature
24 * surface. T = 1 gives a harmonic surface. Use T = 0.25
25 * or so for potential data; something more for topography.
26 *
27 * Program includes overrelaxation for fast convergence and
28 * automatic optimal grid factorization.
29 *
30 * See reference Smith & Wessel (Geophysics, 3, 293-305, 1990) for details.
31 *
32 * Authors: Walter H. F. Smith and Paul Wessel
33 * Date: 1-JAN-2010
34 * Version: 6 API
35 *
36 * This 6.0 version is a complete re-write that differs from the original code:
37 * 1. It uses scan-line grid structures, so we no longer need to transpose grids
38 * 2. It keeps node spacing at 1, thus we no longer need complicated strides between
39 * active nodes. That spacing is now always 1 and we expand the grid as we
40 * go to larger grids (i.e., adding more nodes).
41 * 3. It relies more on functions and macros from GMT to handle rows/cols/node calculations.
42 */
43
44 #include "gmt_dev.h"
45
46 #define THIS_MODULE_CLASSIC_NAME "surface"
47 #define THIS_MODULE_MODERN_NAME "surface"
48 #define THIS_MODULE_LIB "core"
49 #define THIS_MODULE_PURPOSE "Grid table data using adjustable tension continuous curvature splines"
50 #define THIS_MODULE_KEYS "<D{,DD(=,LG(,GG}"
51 #define THIS_MODULE_NEEDS "R"
52 #define THIS_MODULE_OPTIONS "-:RVabdefhiqrw" GMT_ADD_x_OPT GMT_OPT("FH")
53
54 struct SURFACE_CTRL {
55 struct SURFACE_A { /* -A<aspect_ratio> */
56 bool active;
57 unsigned int mode; /* 1 if given as fraction */
58 double value;
59 } A;
60 struct SURFACE_C { /* -C<converge_limit> */
61 bool active;
62 unsigned int mode; /* 1 if given as fraction */
63 double value;
64 } C;
65 struct SURFACE_D { /* -D<line.xyz>[+d][+z[<zval>]] */
66 bool active;
67 bool debug;
68 bool fix_z;
69 double z;
70 char *file; /* Name of file with breaklines */
71 } D;
72 struct SURFACE_G { /* -G<file> */
73 bool active;
74 char *file;
75 } G;
76 struct SURFACE_I { /* -I (for checking only) */
77 bool active;
78 } I;
79 struct SURFACE_J { /* -G<file> */
80 bool active;
81 char *projstring;
82 } J;
83 struct SURFACE_L { /* -Ll|u<limit> */
84 bool active[2];
85 char *file[2];
86 double limit[2];
87 unsigned int mode[2];
88 } L;
89 struct SURFACE_M { /* -M<radius> */
90 bool active;
91 char *arg;
92 } M;
93 struct SURFACE_N { /* -N<max_iterations> */
94 bool active;
95 unsigned int value;
96 } N;
97 struct SURFACE_Q { /* -Q */
98 bool active;
99 } Q;
100 struct SURFACE_S { /* -S<radius>[m|s] */
101 bool active;
102 double radius;
103 char unit;
104 } S;
105 struct SURFACE_T { /* -T<tension>[i][b] */
106 bool active;
107 double b_tension, i_tension;
108 } T;
109 struct SURFACE_W { /* -W[<logfile>] */
110 bool active;
111 char *file;
112 } W;
113 struct SURFACE_Z { /* -Z<over_relaxation> */
114 bool active;
115 double value;
116 } Z;
117 };
118
119 /* Various constants used in surface */
120
121 #define SURFACE_OUTSIDE LONG_MAX /* Index number indicating data is outside usable area */
122 #define SURFACE_CONV_LIMIT 0.0001 /* Default is 100 ppm of data range as convergence criterion */
123 #define SURFACE_MAX_ITERATIONS 500 /* Default iterations at final grid size */
124 #define SURFACE_OVERRELAXATION 1.4 /* Default over-relaxation value */
125 #define SURFACE_CLOSENESS_FACTOR 0.05 /* A node is considered known if the nearest data is within 0.05 of a gridspacing of the node */
126 #define SURFACE_IS_UNCONSTRAINED 0 /* Node has no data constraint within its bin box */
127 #define SURFACE_DATA_IS_IN_QUAD1 1 /* Nearnest data constraint is in quadrant 1 relative to current node */
128 #define SURFACE_DATA_IS_IN_QUAD2 2 /* Nearnest data constraint is in quadrant 2 relative to current node */
129 #define SURFACE_DATA_IS_IN_QUAD3 3 /* Nearnest data constraint is in quadrant 3 relative to current node */
130 #define SURFACE_DATA_IS_IN_QUAD4 4 /* Nearnest data constraint is in quadrant 4 relative to current node */
131 #define SURFACE_IS_CONSTRAINED 5 /* Node has already been set (either data constraint < 5% of grid size or during filling) */
132 #define SURFACE_UNCONSTRAINED 0 /* Use coefficients set for unconstrained node */
133 #define SURFACE_CONSTRAINED 1 /* Use coefficients set for constrained node */
134 #define SURFACE_BREAKLINE 1 /* Flag for breakline constraints that should overrule data constraints */
135
136 /* Misc. macros used to get row, cols, index, node, x, y, plane trend etc. */
137
138 /* Go from row, col to grid node location, accounting for the 2 boundary rows and columns: */
139 #define row_col_to_node(row,col,mx) ((uint64_t)(((int64_t)(row)+(int64_t)2)*((int64_t)(mx))+(int64_t)(col)+(int64_t)2))
140 /* Go from row, col to index array position, where no boundary rows and columns involved: */
141 #define row_col_to_index(row,col,n_columns) ((uint64_t)((int64_t)(row)*((int64_t)(n_columns))+(int64_t)(col)))
142 /* Go from data x to fractional column x: */
143 #define x_to_fcol(x,x0,idx) (((x) - (x0)) * (idx))
144 /* Go from x to grid integer column knowing it is a gridline-registered grid: */
145 #define x_to_col(x,x0,idx) ((int64_t)floor(x_to_fcol(x,x0,idx)+0.5))
146 /* Go from data y to fractional row y_up measured from south (y_up = 0) towards north (y_up = n_rows-1): */
147 #define y_to_frow(y,y0,idy) (((y) - (y0)) * (idy))
148 /* Go from y to row (row = 0 is north) knowing it is a gridline-registered grid: */
149 #define y_to_row(y,y0,idy,n_rows) ((n_rows) - 1 - x_to_col(y,y0,idy))
150 /* Go from col to x knowing it is a gridline-registered grid: */
151 #define col_to_x(col,x0,x1,dx,n_columns) (((int)(col) == (int)((n_columns)-1)) ? (x1) : (x0) + (col) * (dx))
152 /* Go from row to y knowing it is a gridline-registered grid: */
153 #define row_to_y(row,y0,y1,dy,n_rows) (((int)(row) == (int)((n_rows)-1)) ? (y0) : (y1) - (row) * (dy))
154 /* Evaluate the change in LS plane from (0,0) to (xx,y_up) (so in intercept involved): */
155 #define evaluate_trend(C,xx,y_up) (C->plane_sx * (xx) + C->plane_sy * (y_up))
156 /* Evaluate the LS plane at location (xx,y_up) (this includes the intercept): */
157 #define evaluate_plane(C,xx,y_up) (C->plane_icept + evaluate_trend (C, xx, y_up))
158 /* Extract col from index: */
159 #define index_to_col(index,n_columns) ((index) % (n_columns))
160 /* Extract row from index: */
161 #define index_to_row(index,n_columns) ((index) / (n_columns))
162
163 enum surface_nodes { /* Node locations relative to current node, using compass directions */
164 N2 = 0, NW, N1, NE, W2, W1, E1, E2, SW, S1, SE, S2 }; /* I.e., 0-11 */
165
166 /* The 4 indices per quadrant refer to points A-D in Figure A-1 in the reference for quadrant 1 */
167
168 static unsigned int p[5][4] = { /* Indices into C->offset for each of the 4 quadrants, i.e., C->offset[p[quadrant][k]]], k = 0-3 */
169 { 0, 0, 0, 0}, /* This row is never used, so allows us to use quadrant 1-4 as first array index directly */
170 { NW, W1, S1, SE}, /* Indices for 1st quadrant */
171 { SW, S1, E1, NE}, /* Indices for 2nd quadrant */
172 { SE, E1, N1, NW}, /* Indices for 3rd quadrant */
173 { NE, N1, W1, SW} /* Indices for 4th quadrant */
174 };
175
176 enum surface_bound { LO = 0, HI = 1 };
177
178 enum surface_limit { NONE = 0, DATA = 1, VALUE = 2, SURFACE = 3 };
179
180 enum surface_conv { BY_VALUE = 0, BY_PERCENT = 1 };
181
182 enum surface_iter { GRID_NODES = 0, GRID_DATA = 1 };
183
184 struct SURFACE_DATA { /* Data point and index to node it currently constrains */
185 gmt_grdfloat x, y, z;
186 unsigned int kind;
187 uint64_t index;
188 };
189
190 struct SURFACE_BRIGGS { /* Coefficients in Taylor series for Laplacian(z) a la I. C. Briggs (1974) */
191 gmt_grdfloat b[6];
192 };
193
194 struct SURFACE_SEARCH { /* Things needed inside compare function will be passed to qsort_r */
195 int current_nx; /* Number of nodes in y-dir for a given grid factor */
196 int current_ny; /* Number of nodes in y-dir for a given grid factor */
197 double inc[2]; /* Size of each grid cell for a given grid factor */
198 double wesn[4]; /* Grid domain */
199 };
200
201 struct SURFACE_INFO { /* Control structure for surface setup and execution */
202 size_t n_alloc; /* Number of data point positions allocated */
203 uint64_t npoints; /* Number of data points */
204 uint64_t node_sw_corner; /* Node index of southwest interior grid corner for current stride */
205 uint64_t node_se_corner; /* Node index of southeast interior grid corner for current stride */
206 uint64_t node_nw_corner; /* Node index of northwest interior grid corner for current stride */
207 uint64_t node_ne_corner; /* Node index of northeast interior grid corner for current stride */
208 uint64_t n_empty; /* No of unconstrained nodes at initialization */
209 uint64_t nxny; /* Total number of grid nodes without boundaries */
210 uint64_t mxmy; /* Total number of grid nodes with padding */
211 uint64_t total_iterations; /* Total iterations so far. */
212 FILE *fp_log; /* File pointer to log file, if -W is selected */
213 struct SURFACE_DATA *data; /* All the data constraints */
214 struct SURFACE_BRIGGS *Briggs; /* Array with Briggs 6-coefficients per nearest active data constraint */
215 struct GMT_GRID *Grid; /* The final grid */
216 struct GMT_GRID *Bound[2]; /* Optional grids for lower and upper limits on the solution */
217 struct SURFACE_SEARCH info; /* Information needed by the compare function passed to qsort_r */
218 unsigned int n_factors; /* Number of factors in common for the dimensions (n_rows-1, n_columns-1) */
219 unsigned int factors[32]; /* Array of these ommon factors */
220 unsigned int set_limit[2]; /* For low and high: NONE = unconstrained, DATA = by min data value, VALUE = by user value, SURFACE by a grid */
221 unsigned int max_iterations; /* Max iterations per call to iterate */
222 unsigned int converge_mode; /* BY_PERCENT if -C set fractional convergence limit [BY_VALUE] */
223 unsigned int p[5][4]; /* Arrays with four nodes as function of quadrant in constrained fit */
224 int current_stride; /* Current node spacings relative to final spacing */
225 int previous_stride; /* Previous node spacings relative to final spacing */
226 int n_columns; /* Number of nodes in x-dir. (Final grid) */
227 int n_rows; /* Number of nodes in y-dir. (Final grid) */
228 int mx; /* Width of final grid including padding */
229 int my; /* Height of final grid including padding */
230 int current_nx; /* Number of nodes in x-dir for current stride */
231 int current_ny; /* Number of nodes in y-dir for current stride */
232 int current_mx; /* Number of current nodes in x-dir plus 4 extra columns */
233 int previous_nx; /* Number of nodes in x-dir for previous stride */
234 int previous_ny; /* Number of nodes in y-dir for previous stride */
235 int previous_mx; /* Number of current nodes in x-dir plus 4 extra columns */
236 int current_mxmy; /* Total number of grid nodes with padding */
237 int offset[12]; /* Node-indices shifts of 12 nearby points relative center node */
238 unsigned char *status; /* Array with node status or quadrants */
239 char mode_type[2]; /* D = include data points when iterating, I = just interpolate from larger grid */
240 char format[GMT_BUFSIZ]; /* Format statement used in some messages */
241 char *limit_file[2]; /* Pointers to grids with low and high limits, if selected */
242 bool periodic; /* true if geographic grid and west-east == 360 */
243 bool constrained; /* true if set_limit[LO] or set_limit[HI] is true */
244 bool logging; /* true if -W was specified */
245 double limit[2]; /* Low and high constrains on range of solution */
246 double inc[2]; /* Size of each grid cell for current grid factor */
247 double r_inc[2]; /* Reciprocal grid spacings */
248 double converge_limit; /* Convergence limit */
249 double radius; /* Search radius for initializing grid */
250 double tension; /* Tension parameter on the surface */
251 double boundary_tension; /* Tension parameter at the boundary */
252 double interior_tension; /* Tension parameter in the interior */
253 double z_mean; /* Mean value of the data constraints z */
254 double z_rms; /* Root mean square range of z after removing planar trend */
255 double r_z_rms; /* Reciprocal of z_rms (to avoid dividing) */
256 double plane_icept; /* Intercept of best fitting plane to data */
257 double plane_sx; /* Slope of best fitting plane to data in x-direction */
258 double plane_sy; /* Slope of best fitting plane to data in y-direction */
259 double *fraction; /* Hold fractional increments of row and column used in fill_in_forecast */
260 double coeff[2][12]; /* Coefficients for 12 nearby nodes, for constrained [0] and unconstrained [1] nodes */
261 double relax_old, relax_new; /* Coefficients for relaxation factor to speed up convergence */
262 double wesn_orig[4]; /* Original -R domain as we might have shifted it due to -r */
263 double alpha; /* Aspect ratio dy/dx (1 for square pixels) */
264 double a0_const_1, a0_const_2; /* Various constants for off gridnode point equations */
265 double alpha2, e_m2, one_plus_e2;
266 double eps_p2, eps_m2, two_plus_ep2;
267 double two_plus_em2;
268 };
269
surface_set_coefficients(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)270 GMT_LOCAL void surface_set_coefficients (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
271 /* These are the coefficients in the finite-difference expressions given
272 * by equations (A-4) [SURFACE_UNCONSTRAINED=0] and (A-7) [SURFACE_CONSTRAINED=1] in the reference.
273 * Note that the SURFACE_UNCONSTRAINED coefficients are normalized by a0 (20 for no tension/aspects)
274 * whereas the SURFACE_CONSTRAINED is used for a partial sum hence the normalization is done when the
275 * sum over the Briggs coefficients have been included in iterate. */
276 double alpha4, loose, a0;
277
278 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Set finite-difference coefficients [stride = %d]\n", C->current_stride);
279
280 loose = 1.0 - C->interior_tension;
281 C->alpha2 = C->alpha * C->alpha;
282 alpha4 = C->alpha2 * C->alpha2;
283 C->eps_p2 = C->alpha2;
284 C->eps_m2 = 1.0 / C->alpha2;
285 C->one_plus_e2 = 1.0 + C->alpha2;
286 C->two_plus_ep2 = 2.0 + 2.0 * C->eps_p2;
287 C->two_plus_em2 = 2.0 + 2.0 * C->eps_m2;
288
289 C->e_m2 = 1.0 / C->alpha2;
290
291 a0 = 1.0 / ( (6 * alpha4 * loose + 10 * C->alpha2 * loose + 8 * loose - 2 * C->one_plus_e2) + 4 * C->interior_tension * C->one_plus_e2);
292 C->a0_const_1 = 2.0 * loose * (1.0 + alpha4);
293 C->a0_const_2 = 2.0 - C->interior_tension + 2 * loose * C->alpha2;
294
295 C->coeff[SURFACE_CONSTRAINED][W2] = C->coeff[SURFACE_CONSTRAINED][E2] = -loose;
296 C->coeff[SURFACE_CONSTRAINED][N2] = C->coeff[SURFACE_CONSTRAINED][S2] = -loose * alpha4;
297 C->coeff[SURFACE_UNCONSTRAINED][W2] = C->coeff[SURFACE_UNCONSTRAINED][E2] = -loose * a0;
298 C->coeff[SURFACE_UNCONSTRAINED][N2] = C->coeff[SURFACE_UNCONSTRAINED][S2] = -loose * alpha4 * a0;
299 C->coeff[SURFACE_CONSTRAINED][W1] = C->coeff[SURFACE_CONSTRAINED][E1] = 2 * loose * C->one_plus_e2;
300 C->coeff[SURFACE_UNCONSTRAINED][W1] = C->coeff[SURFACE_UNCONSTRAINED][E1] = (2 * C->coeff[SURFACE_CONSTRAINED][W1] + C->interior_tension) * a0;
301 C->coeff[SURFACE_CONSTRAINED][N1] = C->coeff[SURFACE_CONSTRAINED][S1] = C->coeff[SURFACE_CONSTRAINED][W1] * C->alpha2;
302 C->coeff[SURFACE_UNCONSTRAINED][N1] = C->coeff[SURFACE_UNCONSTRAINED][S1] = C->coeff[SURFACE_UNCONSTRAINED][W1] * C->alpha2;
303 C->coeff[SURFACE_CONSTRAINED][NW] = C->coeff[SURFACE_CONSTRAINED][NE] = C->coeff[SURFACE_CONSTRAINED][SW] =
304 C->coeff[SURFACE_CONSTRAINED][SE] = -2 * loose * C->alpha2;
305 C->coeff[SURFACE_UNCONSTRAINED][NW] = C->coeff[SURFACE_UNCONSTRAINED][NE] = C->coeff[SURFACE_UNCONSTRAINED][SW] =
306 C->coeff[SURFACE_UNCONSTRAINED][SE] = C->coeff[SURFACE_CONSTRAINED][NW] * a0;
307
308 C->alpha2 *= 2; /* We will need these coefficients times two in the boundary conditions; do the doubling here */
309 C->e_m2 *= 2;
310 }
311
surface_set_offset(struct SURFACE_INFO * C)312 GMT_LOCAL void surface_set_offset (struct SURFACE_INFO *C) {
313 /* The offset array holds the offset in 1-D index relative
314 * to the current node. For movement along a row this is
315 * always -2, -1, 0, +1, +2 but along a column we move in
316 * multiples of current_mx, the extended grid row width,
317 * which is current_mx = current_nx + 4.
318 */
319 C->offset[N2] = -2 * C->current_mx; /* N2: 2 rows above */
320 C->offset[NW] = -C->current_mx - 1; /* NW: 1 row above and one column left */
321 C->offset[N1] = -C->current_mx; /* N1: 1 row above */
322 C->offset[NE] = -C->current_mx + 1; /* NE: 1 row above and one column right */
323 C->offset[W2] = -2; /* W2: 2 columns left */
324 C->offset[W1] = -1; /* W1 : 1 column left */
325 C->offset[E1] = +1; /* E1 : 1 column right */
326 C->offset[E2] = +2; /* E2 : 2 columns right */
327 C->offset[SW] = C->current_mx - 1; /* SW : 1 row below and one column left */
328 C->offset[S1] = C->current_mx; /* S1 : 1 row below */
329 C->offset[SE] = C->current_mx + 1; /* SE : 1 row below and one column right */
330 C->offset[S2] = 2 * C->current_mx; /* S2 : 2 rows below */
331 }
332
fill_in_forecast(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)333 GMT_LOCAL void fill_in_forecast (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
334
335 /* Fills in bilinear estimates into new node locations after grid is expanded.
336 These new nodes are marked as unconstrained while the coarser data are considered
337 constraints in the next iteration. We do this in two steps:
338 a) We sweep through the grid from last to first node and copy each node to the
339 new location due to increased grid dimensions.
340 b) Once nodes are in place we sweep through and apply the bilinear interpolation.
341 */
342
343 uint64_t index_00, index_10, index_11, index_01, index_new, current_node, previous_node;
344 int previous_row, previous_col, i, j, col, row, expand, first;
345 unsigned char *status = C->status;
346 double c, sx, sy, sxy, r_prev_size, c_plus_sy_dy, sx_plus_sxy_dy;
347 gmt_grdfloat *u = C->Grid->data;
348
349 /* First we expand the active grid to allow for more nodes. We do this by
350 * looping backwards from last node to first so that the locations we copy
351 * the old node values to will always have higher node number than their source.
352 * The previous grid solution has dimensions previous_nx x previous_ny while
353 * the new grid has dimensions current_nx x current_ny. We thus loop over
354 * the old grid and place these nodes into the new grid. */
355
356 expand = C->previous_stride / C->current_stride; /* Multiplicity of new nodes in both x and y dimensions */
357 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Expand grid by factor of %d when going from stride = %d to %d\n", expand, C->previous_stride, C->current_stride);
358
359 for (previous_row = C->previous_ny - 1; previous_row >= 0; previous_row--) { /* Loop backward over the previous grid rows */
360 row = previous_row * expand; /* Corresponding row in the new extended grid */
361 for (previous_col = C->previous_nx - 1; previous_col >= 0; previous_col--) { /* Loop backward over previous grid cols */
362 col = previous_col * expand; /* Corresponding col in the new extended grid */
363 current_node = row_col_to_node (row, col, C->current_mx); /* Current node index */
364 previous_node = row_col_to_node (previous_row, previous_col, C->previous_mx); /* Previous node index */
365 C->Grid->data[current_node] = C->Grid->data[previous_node]; /* Copy the value over */
366 }
367 }
368
369 /* The active grid has now increased in size and the previous values have been copied to their new nodes.
370 * The grid nodes in-between these new "constrained" nodes are partly filled with old values (since
371 * we just copied, not moved, the nodes) or zeros (since we expanded the grid into new unused memory).
372 * This does not matter since we will now fill in those in-between nodes with a bilinear interpolation
373 * based on the coarser (previous) nodes. At the end all nodes in the active grid are valid, except
374 * in the boundary rows/cols. These are reset by set_BC before the iteration starts.
375 */
376
377 /* Precalculate the fractional increments of rows and cols in-between the old constrained rows and cols.
378 * These are all fractions between 0 and 1. E.g., if we quadruple the grid dimensions in x and y then
379 * expand == 4 and we need 4 fractions = {0, 0.25, 0.5, 0.75}. */
380
381 r_prev_size = 1.0 / (double)C->previous_stride;
382 for (i = 0; i < expand; i++) C->fraction[i] = i * r_prev_size;
383
384 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Fill in expanded grid by bilinear interpolation [stride = %d]\n", C->current_stride);
385
386 /* Loop over 4-point "bin squares" from the first northwest bin to the last southeast bin. The bin vertices are the expanded previous nodes */
387
388 for (previous_row = 1; previous_row < C->previous_ny; previous_row++) { /* Starts at row 1 since it is the baseline for the bin extending up to row 0 (north) */
389 row = previous_row * expand; /* Corresponding row in the new extended grid */
390
391 for (previous_col = 0; previous_col < (C->previous_nx-1); previous_col++) { /* Stop 1 column short of east since east is the right boundary of last bin */
392 col = previous_col * expand; /* Corresponding col in the new extended grid */
393
394 /* Get the indices of the bilinear square defined by nodes {00, 10, 11, 01}, with 00 referring to the current (lower left) node */
395 index_00 = row_col_to_node (row, col, C->current_mx); /* Lower left corner of square bin and our origin */
396 index_01 = index_00 - expand * C->current_mx; /* Upper left corner of square bin */
397 index_10 = index_00 + expand; /* Lower right corner of square bin */
398 index_11 = index_01 + expand; /* Upper right corner of square bin */
399
400 /* Get bilinear coefficients for interpolation z = c + sx * delta_x + sy * delta_y + sxy * delta_x * delta_y,
401 * which we will use as z = (c + sy * delta_y) + delta_x * (sx + sxy * delta_y).
402 * Below, delta_x and delta_y are obtained via C->fraction[i|j] that we pre-calculated above. */
403 c = u[index_00]; sx = u[index_10] - c;
404 sy = u[index_01] - c; sxy = u[index_11] - u[index_10] - sy;
405
406 /* Fill in all the denser nodes except the lower-left starting point */
407
408 for (j = 0, first = 1; j < expand; j++) { /* Set first = 1 so we skip the first column when j = 0 */
409 c_plus_sy_dy = c + sy * C->fraction[j]; /* Compute terms that remain constant for this j */
410 sx_plus_sxy_dy = sx + sxy * C->fraction[j];
411 index_new = index_00 - j * C->current_mx + first; /* Start node on this intermediate row */
412 for (i = first; i < expand; i++, index_new++) { /* Sweep across this row and interpolate */
413 u[index_new] = (gmt_grdfloat)(c_plus_sy_dy + C->fraction[i] * sx_plus_sxy_dy);
414 status[index_new] = SURFACE_IS_UNCONSTRAINED; /* These are considered temporary estimates */
415 }
416 first = 0; /* Reset to 0 for the remainder of the j loop */
417 }
418 status[index_00] = SURFACE_IS_CONSTRAINED; /* The previous node values will be kept fixed in the next iterate call */
419 }
420 }
421
422 /* The loops above exclude the north and east boundaries. First do linear interpolation along the east edge */
423
424 index_00 = C->node_ne_corner; /* Upper NE node */
425 for (previous_row = 1; previous_row < C->previous_ny; previous_row++) { /* So first edge is from row = 1 up to row = 0 on eastern edge */
426 index_01 = index_00; /* Previous lower becomes current upper node */
427 index_00 += expand * C->current_mx; /* Lower node after striding down */
428 sy = u[index_01] - u[index_00]; /* Vertical gradient in u toward ymax (for increasing j) */
429 index_new = index_00 - C->current_mx; /* Since we start at j = 1 we skip up one row here */
430 for (j = 1; j < expand; j++, index_new -= C->current_mx) { /* Start at 1 since we skip the constrained index_00 node */
431 u[index_new] = u[index_00] + (gmt_grdfloat)(C->fraction[j] * sy);
432 status[index_new] = SURFACE_IS_UNCONSTRAINED; /* These are considered temporary estimates */
433 }
434 status[index_00] = SURFACE_IS_CONSTRAINED; /* The previous node values will be kept fixed in the next iterate call */
435 }
436 /* Next do linear interpolation along the north edge */
437 index_10 = C->node_nw_corner; /* Left NW node */
438 for (previous_col = 0; previous_col < (C->previous_nx-1); previous_col++) { /* To ensure last edge ends at col = C->previous_nx-1 */
439 index_00 = index_10; /* Previous right node becmes current left node */
440 index_10 = index_00 + expand; /* Right node after striding to the right */
441 sx = u[index_10] - u[index_00]; /* Horizontal gradient in u toward xmax (for increasing i) */
442 index_new = index_00 + 1; /* Start at 1 since we skip the constrained index_00 node */
443 for (i = 1; i < expand; i++, index_new++) {
444 u[index_new] = u[index_00] + (gmt_grdfloat)(C->fraction[i] * sx);
445 status[index_new] = SURFACE_IS_UNCONSTRAINED; /* These are considered temporary estimates */
446 }
447 status[index_00] = SURFACE_IS_CONSTRAINED; /* The previous node values will be kept fixed in the next iterate call */
448 }
449 /* Finally set the northeast corner to be considered fixed in the next iterate call and our work here is done */
450 status[C->node_ne_corner] = SURFACE_IS_CONSTRAINED;
451 }
452
surface_compare_points(const void * point_1v,const void * point_2v,void * arg)453 GMT_LOCAL int surface_compare_points (const void *point_1v, const void *point_2v, void *arg) {
454 /* Routine for qsort_r to sort data structure for fast access to data by node location.
455 Sorts on index first, then on radius to node corresponding to index, so that index
456 goes from low to high, and so does radius. Note: These are simple Cartesian distance
457 * calculations. The metadata needed to do the calculations are passed via *arg.
458 */
459 uint64_t col, row, index_1, index_2;
460 double x0, y0, dist_1, dist_2;
461 const struct SURFACE_DATA *point_1 = point_1v, *point_2 = point_2v;
462 struct SURFACE_SEARCH *info;
463 index_1 = point_1->index;
464 index_2 = point_2->index;
465 if (index_1 < index_2) return (-1);
466 if (index_1 > index_2) return (+1);
467 if (index_1 == SURFACE_OUTSIDE) return (0);
468 /* Points are in same grid cell. First check for breakline points to sort those ahead of data points */
469 if (point_1->kind == SURFACE_BREAKLINE && point_2->kind == 0) return (-1);
470 if (point_2->kind == SURFACE_BREAKLINE && point_1->kind == 0) return (+1);
471 /* Now find the one who is nearest to grid point */
472 /* Note: index calculations do not include boundary pad */
473 info = arg; /* Get the needed metadata for distance calculations */
474 row = index_to_row (index_1, info->current_nx);
475 col = index_to_col (index_1, info->current_nx);
476 x0 = col_to_x (col, info->wesn[XLO], info->wesn[XHI], info->inc[GMT_X], info->current_nx);
477 y0 = row_to_y (row, info->wesn[YLO], info->wesn[YHI], info->inc[GMT_Y], info->current_ny);
478 dist_1 = (point_1->x - x0) * (point_1->x - x0) + (point_1->y - y0) * (point_1->y - y0);
479 /* Try to speed things up by first checking if point_2 x-distance from x0 alone exceeds point_1's radial distance */
480 dist_2 = (point_2->x - x0) * (point_2->x - x0); /* Just dx^2 */
481 if (dist_1 < dist_2) return (-1); /* Don't need to consider the y-distance */
482 /* Did not exceed, so now we must finalize the dist_2 calculation by including the y-separation */
483 dist_2 += (point_2->y - y0) * (point_2->y - y0);
484 if (dist_1 < dist_2) return (-1);
485 if (dist_1 > dist_2) return (+1);
486 return (0);
487 }
488
surface_smart_divide(struct SURFACE_INFO * C)489 GMT_LOCAL void surface_smart_divide (struct SURFACE_INFO *C) {
490 /* Divide grid by its next largest prime factor and shift that setting by one */
491 C->current_stride /= C->factors[C->n_factors - 1];
492 C->n_factors--;
493 }
494
surface_set_index(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)495 GMT_LOCAL void surface_set_index (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
496 /* Recomputes data[k].index for the new value of the stride,
497 sorts the data again on index and radii, and throws away
498 data which are now outside the usable limits.
499 Note: These indices exclude the padding. */
500 int col, row;
501 uint64_t k, k_skipped = 0;
502 struct GMT_GRID_HEADER *h = C->Grid->header;
503
504 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Recompute data index for next iteration [stride = %d]\n", C->current_stride);
505
506 for (k = 0; k < C->npoints; k++) {
507 col = (int)x_to_col (C->data[k].x, h->wesn[XLO], C->r_inc[GMT_X]);
508 row = (int)y_to_row (C->data[k].y, h->wesn[YLO], C->r_inc[GMT_Y], C->current_ny);
509 if (col < 0 || col >= C->current_nx || row < 0 || row >= C->current_ny) {
510 C->data[k].index = SURFACE_OUTSIDE;
511 k_skipped++;
512 }
513 else
514 C->data[k].index = row_col_to_index (row, col, C->current_nx);
515 }
516
517 qsort_r (C->data, C->npoints, sizeof (struct SURFACE_DATA), surface_compare_points, &(C->info));
518
519 C->npoints -= k_skipped;
520 }
521
surface_solve_Briggs_coefficients(struct SURFACE_INFO * C,gmt_grdfloat * b,double xx,double yy,gmt_grdfloat z)522 GMT_LOCAL void surface_solve_Briggs_coefficients (struct SURFACE_INFO *C, gmt_grdfloat *b, double xx, double yy, gmt_grdfloat z) {
523 /* Given the normalized offset (xx,yy) from current node (value z) we determine the
524 * Briggs coefficients b_k, k = 1,5 [Equation (A-6) in the reference]
525 * Here, xx, yy are the fractional distances, accounting for any anisotropy.
526 * Note b[5] initially contains the sum of the 5 Briggs coefficients but
527 * we actually need to divide by it so we do that change here as well.
528 * Finally, b[4] will be multiplied with the off-node constraint so we do that here.
529 */
530 double xx2, yy2, xx_plus_yy, xx_plus_yy_plus_one, inv_xx_plus_yy_plus_one, inv_delta, b_4;
531
532 xx_plus_yy = xx + yy;
533 xx_plus_yy_plus_one = 1.0 + xx_plus_yy;
534 inv_xx_plus_yy_plus_one = 1.0 / xx_plus_yy_plus_one;
535 xx2 = xx * xx; yy2 = yy * yy;
536 inv_delta = inv_xx_plus_yy_plus_one / xx_plus_yy;
537 b[0] = (gmt_grdfloat)((xx2 + 2.0 * xx * yy + xx - yy2 - yy) * inv_delta);
538 b[1] = (gmt_grdfloat)(2.0 * (yy - xx + 1.0) * inv_xx_plus_yy_plus_one);
539 b[2] = (gmt_grdfloat)(2.0 * (xx - yy + 1.0) * inv_xx_plus_yy_plus_one);
540 b[3] = (gmt_grdfloat)((-xx2 + 2.0 * xx * yy - xx + yy2 + yy) * inv_delta);
541 b_4 = 4.0 * inv_delta;
542 /* We also need to normalize by the sum of the b[k] values, so sum them here */
543 b[5] = b[0] + b[1] + b[2] + b[3] + (gmt_grdfloat)b_4;
544 /* We need to sum k = 0<5 of u[k]*b[k], where u[k] are the nodes of the points A-D,
545 * but the k = 4 point (E) is our data constraint. We multiply that in here, once,
546 * add add b[4] to the rest of the sum inside the iteration loop. */
547 b[4] = (gmt_grdfloat)(b_4 * z);
548
549 /* b[5] is part of a denominator so we do the division here instead of inside iterate loop */
550 b[5] = (gmt_grdfloat)(1.0 / (C->a0_const_1 + C->a0_const_2 * b[5]));
551 }
552
surface_find_nearest_constraint(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)553 GMT_LOCAL void surface_find_nearest_constraint (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
554 /* Determines the nearest data point per bin and sets the
555 * Briggs parameters or, if really close, fixes the node value */
556 uint64_t k, last_index, node, briggs_index;
557 int row, col;
558 double xx, yy, x0, y0, dx, dy;
559 gmt_grdfloat z_at_node, *u = C->Grid->data;
560 unsigned char *status = C->status;
561 struct GMT_GRID_HEADER *h = C->Grid->header;
562
563 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Determine nearest point and set Briggs coefficients [stride = %d]\n", C->current_stride);
564
565 gmt_M_grd_loop (GMT, C->Grid, row, col, node) { /* Reset status of all interior grid nodes */
566 status[node] = SURFACE_IS_UNCONSTRAINED;
567 }
568
569 last_index = UINTMAX_MAX; briggs_index = 0U;
570
571 for (k = 0; k < C->npoints; k++) { /* Find constraining value */
572 if (C->data[k].index != last_index) { /* Moving to the next node to address its nearest data constraint */
573 /* Note: Index calculations do not consider the boundary padding */
574 row = (int)index_to_row (C->data[k].index, C->current_nx);
575 col = (int)index_to_col (C->data[k].index, C->current_nx);
576 last_index = C->data[k].index; /* Now this is the last unique index we worked on */
577 node = row_col_to_node (row, col, C->current_mx);
578 /* Get coordinates of this node */
579 x0 = col_to_x (col, h->wesn[XLO], h->wesn[XHI], C->inc[GMT_X], C->current_nx);
580 y0 = row_to_y (row, h->wesn[YLO], h->wesn[YHI], C->inc[GMT_Y], C->current_ny);
581 /* Get offsets dx,dy of data point location relative to this node (dy is positive up) in fraction of grid increments */
582 dx = x_to_fcol (C->data[k].x, x0, C->r_inc[GMT_X]);
583 dy = y_to_frow (C->data[k].y, y0, C->r_inc[GMT_Y]);
584 /* So dx, dy are here fractions of C->inc[GMT_X] and C-inc[GMT_Y] */
585 /* "Really close" will mean within 5% of the current grid spacing from the center node */
586
587 if (fabs (dx) < SURFACE_CLOSENESS_FACTOR && fabs (dy) < SURFACE_CLOSENESS_FACTOR) { /* Considered close enough to assign fixed value to node */
588 status[node] = SURFACE_IS_CONSTRAINED;
589 /* Since data constraint is forcibly moved from (dx, dy) to (0,0) we must adjust for
590 * the small change in the planar trend between the two locations, and then
591 * possibly clip the value if constraining surfaces were given. Note that
592 * dx, dy is in -1/1 range normalized by (current_x|y_inc) so to recover the
593 * corresponding dx,dy in units of current grid fractions we must scale both
594 * dx and dy by current_stride; this is equivalent to scaling the trend.
595 * This trend then is normalized by dividing by the z rms.*/
596
597 z_at_node = C->data[k].z + (gmt_grdfloat) (C->r_z_rms * C->current_stride * evaluate_trend (C, dx, dy));
598 if (C->constrained) {
599 if (C->set_limit[LO] && !gmt_M_is_fnan (C->Bound[LO]->data[node]) && z_at_node < C->Bound[LO]->data[node])
600 z_at_node = C->Bound[LO]->data[node];
601 else if (C->set_limit[HI] && !gmt_M_is_fnan (C->Bound[HI]->data[node]) && z_at_node > C->Bound[HI]->data[node])
602 z_at_node = C->Bound[HI]->data[node];
603 }
604 u[node] = z_at_node;
605 }
606 else { /* We have a nearby data point in one of the quadrants */
607 /* Note: We must swap dx,dy for 2nd and 4th quadrants and always use absolute values since we are
608 rotating other cases (quadrants 2-4) to look like quadrant 1 */
609 if (dy >= 0.0) { /* Upper two quadrants */
610 if (dx >= 0.0) { /* Both positive, use as is */
611 status[node] = SURFACE_DATA_IS_IN_QUAD1;
612 xx = dx; yy = dy;
613 }
614 else { /* dx negative, so remove sign, and swap */
615 status[node] = SURFACE_DATA_IS_IN_QUAD2;
616 yy = -dx; xx = dy;
617 }
618 }
619 else { /* Lower two quadrants where we need to remove sign from dy */
620 if (dx >= 0.0) { /* Also swap x and y */
621 status[node] = SURFACE_DATA_IS_IN_QUAD4;
622 yy = dx; xx = -dy;
623 }
624 else { /* Just remove both signs */
625 status[node] = SURFACE_DATA_IS_IN_QUAD3;
626 xx = -dx; yy = -dy;
627 }
628 }
629 /* Evaluate the Briggs coefficients */
630 surface_solve_Briggs_coefficients (C, C->Briggs[briggs_index].b, xx, yy, C->data[k].z);
631 briggs_index++;
632 }
633 }
634 }
635 }
636
surface_set_grid_parameters(struct SURFACE_INFO * C)637 GMT_LOCAL void surface_set_grid_parameters (struct SURFACE_INFO *C) {
638 /* Set the previous settings to the current settings */
639 C->previous_nx = C->current_nx;
640 C->previous_mx = C->current_mx;
641 C->previous_ny = C->current_ny;
642 /* Update the current parameters given the new C->current_stride setting */
643 C->info.current_nx = C->current_nx = (C->n_columns - 1) / C->current_stride + 1;
644 C->info.current_ny = C->current_ny = (C->n_rows - 1) / C->current_stride + 1;
645 C->current_mx = C->current_nx + 4;
646 C->current_mxmy = C->current_mx * (C->current_ny + 4); /* Only place where "my" is used */
647 C->info.inc[GMT_X] = C->inc[GMT_X] = C->current_stride * C->Grid->header->inc[GMT_X];
648 C->info.inc[GMT_Y] = C->inc[GMT_Y] = C->current_stride * C->Grid->header->inc[GMT_Y];
649 C->r_inc[GMT_X] = 1.0 / C->inc[GMT_X];
650 C->r_inc[GMT_Y] = 1.0 / C->inc[GMT_Y];
651 /* Update the grid node indices of the 4 corners */
652 C->node_nw_corner = 2 * C->current_mx + 2;
653 C->node_sw_corner = C->node_nw_corner + (C->current_ny - 1) * C->current_mx;
654 C->node_se_corner = C->node_sw_corner + C->current_nx - 1;
655 C->node_ne_corner = C->node_nw_corner + C->current_nx - 1;
656 }
657
surface_initialize_grid(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)658 GMT_LOCAL void surface_initialize_grid (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
659 /* For the initial gridsize, compute weighted averages of data inside the search radius
660 * and assign the values to u[col,row], where col,row are multiples of gridsize.
661 * Weights are Gaussian, i.e., this is a MA Gaussian filter operation.
662 */
663 uint64_t index_1, index_2, k, k_index, node;
664 int del_col, del_row, col, row, col_min, col_max, row_min, row_max, ki, kj;
665 double r, rfact, sum_w, sum_zw, weight, x0, y0;
666 gmt_grdfloat *u = C->Grid->data;
667 struct GMT_GRID_HEADER *h = C->Grid->header;
668
669 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Initialize grid using moving average scheme [stride = %d]\n", C->current_stride);
670
671 del_col = irint (ceil (C->radius / C->inc[GMT_X]));
672 del_row = irint (ceil (C->radius / C->inc[GMT_Y]));
673 rfact = -4.5 / (C->radius*C->radius);
674 for (row = 0; row < C->current_ny; row++) {
675 y0 = row_to_y (row, h->wesn[YLO], h->wesn[YHI], C->inc[GMT_Y], C->current_ny);
676 for (col = 0; col < C->current_nx; col++) {
677 /* For this node on the grid, find all data points within the radius */
678 x0 = col_to_x (col, h->wesn[XLO], h->wesn[XHI], C->inc[GMT_X], C->current_nx);
679 col_min = col - del_col;
680 if (col_min < 0) col_min = 0;
681 col_max = col + del_col;
682 if (col_max >= C->current_nx) col_max = C->current_nx - 1;
683 row_min = row - del_row;
684 if (row_min < 0) row_min = 0;
685 row_max = row + del_row;
686 if (row_max >= C->current_ny) row_max = C->current_ny - 1;
687 index_1 = row_col_to_index (row_min, col_min, C->current_nx);
688 index_2 = row_col_to_index (row_max, col_max+1, C->current_nx);
689 sum_w = sum_zw = 0.0;
690 k = 0;
691 while (k < C->npoints && C->data[k].index < index_1) k++;
692 /* This double loop visits all nodes within the rectangle of dimensions (2*del_col by 2*del_row) centered on x0,y0 */
693 for (kj = row_min; k < C->npoints && kj <= row_max && C->data[k].index < index_2; kj++) {
694 for (ki = col_min; k < C->npoints && ki <= col_max && C->data[k].index < index_2; ki++) {
695 k_index = row_col_to_index (kj, ki, C->current_nx);
696 while (k < C->npoints && C->data[k].index < k_index) k++;
697 while (k < C->npoints && C->data[k].index == k_index) {
698 r = (C->data[k].x-x0)*(C->data[k].x-x0) + (C->data[k].y-y0)*(C->data[k].y-y0);
699 if (r > C->radius) continue; /* Outside the circle */
700 weight = exp (rfact * r);
701 sum_w += weight;
702 sum_zw += weight * C->data[k].z;
703 k++;
704 }
705 }
706 }
707 node = row_col_to_node (row, col, C->current_mx);
708 if (sum_w == 0.0) {
709 sprintf (C->format, "No data inside search radius at: %s %s [node set to data mean]\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
710 GMT_Report (GMT->parent, GMT_MSG_WARNING, C->format, x0, y0);
711 u[node] = (gmt_grdfloat)C->z_mean;
712 }
713 else
714 u[node] = (gmt_grdfloat)(sum_zw/sum_w);
715 }
716 }
717 }
718
surface_read_data(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,struct GMT_OPTION * options)719 GMT_LOCAL int surface_read_data (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, struct GMT_OPTION *options) {
720 /* Procdss input data into data structure */
721 int col, row, error;
722 uint64_t k = 0, kmax = 0, kmin = 0, n_dup = 0;
723 double *in, half_dx, zmin = DBL_MAX, zmax = -DBL_MAX, wesn_lim[4];
724 struct GMT_GRID_HEADER *h = C->Grid->header;
725 struct GMT_RECORD *In = NULL;
726
727 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Processing input table data\n");
728 C->data = gmt_M_memory (GMT, NULL, C->n_alloc, struct SURFACE_DATA);
729
730 /* Read in xyz data and computes index no and store it in a structure */
731
732 if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, 3, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR)
733 return (error);
734 if (GMT_Init_IO (GMT->parent, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) /* Establishes data input */
735 return (GMT->parent->error);
736
737 C->z_mean = 0.0;
738 /* Initially allow points to be within 1 grid spacing of the grid */
739 wesn_lim[XLO] = h->wesn[XLO] - C->inc[GMT_X]; wesn_lim[XHI] = h->wesn[XHI] + C->inc[GMT_X];
740 wesn_lim[YLO] = h->wesn[YLO] - C->inc[GMT_Y]; wesn_lim[YHI] = h->wesn[YHI] + C->inc[GMT_Y];
741 half_dx = 0.5 * C->inc[GMT_X];
742
743 if (GMT_Begin_IO (GMT->parent, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) /* Enables data input and sets access mode */
744 return (GMT->parent->error);
745
746 do { /* Keep returning records until we reach EOF */
747 if ((In = GMT_Get_Record (GMT->parent, GMT_READ_DATA, NULL)) == NULL) { /* Read next record, get NULL if special case */
748 if (gmt_M_rec_is_error (GMT)) /* Bail if there are any read errors */
749 return (GMT_RUNTIME_ERROR);
750 else if (gmt_M_rec_is_eof (GMT)) /* Reached end of file */
751 break;
752 continue; /* Go back and read the next record */
753 }
754
755 if (In->data == NULL) {
756 gmt_quit_bad_record (GMT->parent, In);
757 return (GMT->parent->error);
758 }
759
760 /* Data record to process */
761 in = In->data; /* Only need to process numerical part here */
762
763 if (gmt_M_is_dnan (in[GMT_Z])) continue; /* Cannot use NaN values */
764 if (gmt_M_y_is_outside (GMT, in[GMT_Y], wesn_lim[YLO], wesn_lim[YHI])) continue; /* Outside y-range (or latitude) */
765 if (gmt_x_is_outside (GMT, &in[GMT_X], wesn_lim[XLO], wesn_lim[XHI])) continue; /* Outside x-range (or longitude) */
766 row = (int)y_to_row (in[GMT_Y], h->wesn[YLO], C->r_inc[GMT_Y], C->current_ny);
767 if (row < 0 || row >= C->current_ny) continue;
768 if (C->periodic && ((h->wesn[XHI]-in[GMT_X]) < half_dx)) { /* Push all values to the western nodes */
769 in[GMT_X] -= 360.0; /* Make this point constrain the western node value and then duplicate to east later */
770 col = 0;
771 }
772 else /* Regular point not at the periodic boundary */
773 col = (int)x_to_col (in[GMT_X], h->wesn[XLO], C->r_inc[GMT_X]);
774 if (col < 0 || col >= C->current_nx) continue;
775
776 C->data[k].index = row_col_to_index (row, col, C->current_nx);
777 C->data[k].x = (gmt_grdfloat)in[GMT_X];
778 C->data[k].y = (gmt_grdfloat)in[GMT_Y];
779 C->data[k].z = (gmt_grdfloat)in[GMT_Z];
780 /* Determine the mean, min and max z-values */
781 if (zmin > in[GMT_Z]) zmin = in[GMT_Z], kmin = k;
782 if (zmax < in[GMT_Z]) zmax = in[GMT_Z], kmax = k;
783 C->z_mean += in[GMT_Z];
784 if (++k == C->n_alloc) {
785 C->n_alloc <<= 1;
786 C->data = gmt_M_memory (GMT, C->data, C->n_alloc, struct SURFACE_DATA);
787 }
788 if (C->periodic && col == 0) { /* Now we must replicate information from the western to the eastern boundary */
789 col = C->current_nx - 1;
790 C->data[k].index = row_col_to_index (row, col, C->current_nx);
791 C->data[k].x = (gmt_grdfloat)(in[GMT_X] + 360.0);
792 C->data[k].y = (gmt_grdfloat)in[GMT_Y];
793 C->data[k].z = (gmt_grdfloat)in[GMT_Z];
794 C->z_mean += in[GMT_Z];
795 if (++k == C->n_alloc) {
796 C->n_alloc <<= 1;
797 C->data = gmt_M_memory (GMT, C->data, C->n_alloc, struct SURFACE_DATA);
798 }
799 n_dup++;
800 }
801 } while (true);
802
803
804 if (GMT_End_IO (GMT->parent, GMT_IN, 0) != GMT_NOERROR) /* Disables further data input */
805 return (GMT->parent->error);
806
807 C->npoints = k; /* Number of data points that passed being "inside" the grid region */
808
809 if (C->npoints == 0) {
810 GMT_Report (GMT->parent, GMT_MSG_ERROR, "No datapoints inside region, aborting\n");
811 return (GMT_RUNTIME_ERROR);
812 }
813
814 C->z_mean /= C->npoints; /* Estimate mean data value */
815 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
816 char msg[GMT_LEN256] = {""};
817 sprintf (C->format, "%s %s %s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
818 sprintf (msg, C->format, (double)C->data[kmin].x, (double)C->data[kmin].y, (double)C->data[kmin].z);
819 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Minimum value of your dataset x,y,z at: %s", msg);
820 sprintf (msg, C->format, (double)C->data[kmax].x, (double)C->data[kmax].y, (double)C->data[kmax].z);
821 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Maximum value of your dataset x,y,z at: %s", msg);
822 if (C->periodic && n_dup) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Number of input values shared between repeating west and east column nodes: %" PRIu64 "\n", n_dup);
823 }
824 C->data = gmt_M_memory (GMT, C->data, C->npoints, struct SURFACE_DATA);
825
826 if (C->set_limit[LO] == DATA) /* Wanted to set lower limit based on minimum observed z value */
827 C->limit[LO] = C->data[kmin].z;
828 else if (C->set_limit[LO] == VALUE && C->limit[LO] > C->data[kmin].z)
829 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your lower value is > than min data value.\n");
830 if (C->set_limit[HI] == DATA) /* Wanted to set upper limit based on maximum observed z value */
831 C->limit[HI] = C->data[kmax].z;
832 else if (C->set_limit[HI] == VALUE && C->limit[HI] < C->data[kmax].z)
833 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Your upper value is < than max data value.\n");
834 return (0);
835 }
836
surface_load_constraints(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,int transform)837 GMT_LOCAL int surface_load_constraints (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, int transform) {
838 /* Deal with the constants or grids supplied via -L. Note: Because we remove a
839 * best-fitting plane from the data, even a simple constant constraint will become
840 * a plane and thus must be represented on a grid. */
841 unsigned int end, col, row;
842 uint64_t node;
843 char *limit[2] = {"Lower", "Upper"};
844 double y_up;
845 struct GMTAPI_CTRL *API = GMT->parent;
846
847 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Load any data constraint limit grids\n");
848
849 /* Load lower/upper limits, verify range, deplane, and rescale */
850
851 for (end = LO; end <= HI; end++) {
852 if (C->set_limit[end] == NONE) continue; /* Nothing desired */
853 if (C->set_limit[end] < SURFACE) { /* Got a constant level for this end */
854 if ((C->Bound[end] = GMT_Duplicate_Data (API, GMT_IS_GRID, GMT_DUPLICATE_ALLOC, C->Grid)) == NULL) return (API->error);
855 for (node = 0; node < C->mxmy; node++) C->Bound[end]->data[node] = (gmt_grdfloat)C->limit[end];
856 }
857 else { /* Got a grid with a surface */
858 if ((C->Bound[end] = GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, C->limit_file[end], NULL)) == NULL) return (API->error); /* Get header only */
859 if (C->Bound[end]->header->n_columns != C->Grid->header->n_columns || C->Bound[end]->header->n_rows != C->Grid->header->n_rows) {
860 GMT_Report (API, GMT_MSG_ERROR, "%s limit file not of proper dimensions!\n", limit[end]);
861 return (GMT_RUNTIME_ERROR);
862 }
863 if (GMT_Read_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, C->limit_file[end], C->Bound[end]) == NULL) return (API->error);
864 }
865 if (transform) { /* Remove best-fitting plane and normalize the bounding values */
866 for (row = 0; row < C->Grid->header->n_rows; row++) {
867 y_up = (double)(C->Grid->header->n_rows - row - 1); /* Require y_up = 0 at south and positive toward north */
868 node = row_col_to_node (row, 0, C->current_mx);
869 for (col = 0; col < C->Grid->header->n_columns; col++, node++) {
870 if (gmt_M_is_fnan (C->Bound[end]->data[node])) continue;
871 C->Bound[end]->data[node] -= (gmt_grdfloat)evaluate_plane (C, col, y_up); /* Remove plane */
872 C->Bound[end]->data[node] *= (gmt_grdfloat)C->r_z_rms; /* Normalize residuals */
873 }
874 }
875 }
876 C->constrained = true; /* At least one of the limits will be constrained */
877 }
878
879 return (0);
880 }
881
surface_write_grid(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,char * grdfile)882 GMT_LOCAL int surface_write_grid (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, char *grdfile) {
883 /* Write output grid to file */
884 uint64_t node;
885 int row, col, err, end;
886 char *limit[2] = {"lower", "upper"};
887 gmt_grdfloat *u = C->Grid->data;
888
889 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Prepare final output grid [stride = %d]\n", C->current_stride);
890
891 strcpy (C->Grid->header->title, "Data gridded with continuous surface splines in tension");
892
893 if (GMT->common.R.registration == GMT_GRID_PIXEL_REG) { /* Pixel registration request. Reset region to the original extent */
894 gmt_M_memcpy (C->Grid->header->wesn, C->wesn_orig, 4, double);
895 C->Grid->header->registration = GMT->common.R.registration;
896 /* Must reduce both n_columns,n_rows by 1 and make the eastmost column and northernmost row part of the grid pad */
897 C->Grid->header->n_columns--; C->n_columns--;
898 C->Grid->header->n_rows--; C->n_rows--;
899 C->Grid->header->pad[XHI]++; /* Presumably increase pad from 2 to 3 */
900 C->Grid->header->pad[YHI]++; /* Presumably increase pad from 2 to 3 */
901 gmt_set_grddim (GMT, C->Grid->header); /* Reset all integer dimensions and xy_off */
902 }
903 if (C->constrained) { /* Must check that we don't exceed any imposed limits. */
904 /* Reload the constraints, but this time do not transform the data */
905 if ((err = surface_load_constraints (GMT, C, false)) != 0) return (err);
906
907 gmt_M_grd_loop (GMT, C->Grid, row, col, node) { /* Make sure we clip to the specified bounds */
908 if (C->set_limit[LO] && !gmt_M_is_fnan (C->Bound[LO]->data[node]) && u[node] < C->Bound[LO]->data[node]) u[node] = C->Bound[LO]->data[node];
909 if (C->set_limit[HI] && !gmt_M_is_fnan (C->Bound[HI]->data[node]) && u[node] > C->Bound[HI]->data[node]) u[node] = C->Bound[HI]->data[node];
910 }
911 /* Free any bounding surfaces */
912 for (end = LO; end <= HI; end++) {
913 if ((C->set_limit[end] > NONE && C->set_limit[end] < SURFACE) && GMT_Destroy_Data (GMT->parent, &C->Bound[end]) != GMT_NOERROR) {
914 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failed to free %s boundary\n", limit[end]);
915 }
916 }
917 }
918 if (C->periodic) { /* Ensure exact periodicity at E-W boundaries */
919 for (row = 0; row < C->current_ny; row++) {
920 node = row_col_to_node (row, 0, C->current_mx);
921 u[node] = u[node+C->current_nx-1] = (gmt_grdfloat)(0.5 * (u[node] + u[node+C->current_nx-1])); /* Set these to the same as their average */
922 }
923 }
924 /* Time to write our final grid */
925 if (GMT_Write_Data (GMT->parent, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, grdfile, C->Grid) != GMT_NOERROR)
926 return (GMT->parent->error);
927
928 return (0);
929 }
930
surface_set_BCs(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,gmt_grdfloat * u)931 GMT_LOCAL void surface_set_BCs (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, gmt_grdfloat *u) {
932 /* Fill in auxiliary boundary rows and columns; see equations (A-8,9,10) in the reference */
933 uint64_t n, n_s, n_n, n_w, n_e; /* Node variables */
934 int col, row, *d_n = C->offset; /* Relative changes in node index from present node n */
935 double x_0_const = 4.0 * (1.0 - C->boundary_tension) / (2.0 - C->boundary_tension);
936 double x_1_const = (3 * C->boundary_tension - 2.0) / (2.0 - C->boundary_tension);
937 double y_denom = 2 * C->alpha * (1.0 - C->boundary_tension) + C->boundary_tension;
938 double y_0_const = 4 * C->alpha * (1.0 - C->boundary_tension) / y_denom;
939 double y_1_const = (C->boundary_tension - 2 * C->alpha * (1.0 - C->boundary_tension) ) / y_denom;
940
941 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Apply all boundary conditions [stride = %d]\n", C->current_stride);
942
943 /* First set (1-T)d2[]/dn2 + Td[]/dn = 0 along edges */
944
945 for (col = 0, n_s = C->node_sw_corner, n_n = C->node_nw_corner; col < C->current_nx; col++, n_s++, n_n++) { /* set BC1 along south and north side */
946 u[n_s+d_n[S1]] = (gmt_grdfloat)(y_0_const * u[n_s] + y_1_const * u[n_s+d_n[N1]]); /* South: u_{0,-1} = 2 * u_{0,0} - u_{0,+1} */
947 u[n_n+d_n[N1]] = (gmt_grdfloat)(y_0_const * u[n_n] + y_1_const * u[n_n+d_n[S1]]); /* North: u_{0,+1} = 2 * u_{0,0} - u_{0,-1} */
948 }
949 if (C->periodic) { /* Set periodic boundary conditions in longitude at west and east boundaries */
950 for (row = 0, n_w = C->node_nw_corner, n_e = C->node_ne_corner; row < C->current_ny; row++, n_w += C->current_mx, n_e += C->current_mx) {
951 u[n_w+d_n[W1]] = u[n_e+d_n[W1]]; /* West side */
952 u[n_e+d_n[E1]] = u[n_w+d_n[E1]]; /* East side */
953 u[n_e] = u[n_w] = 0.5f * (u[n_e] + u[n_w]); /* Set to average of east and west */
954 }
955 }
956 else { /* Regular natural BC */
957 for (row = 0, n_w = C->node_nw_corner, n_e = C->node_ne_corner; row < C->current_ny; row++, n_w += C->current_mx, n_e += C->current_mx) {
958 /* West: u_{-10} = 2 * u_{00} - u_{10} */
959 u[n_w+d_n[W1]] = (gmt_grdfloat)(x_1_const * u[n_w+d_n[E1]] + x_0_const * u[n_w]);
960 /* East: u_{10} = 2 * u_{00} - u_{-10} */
961 u[n_e+d_n[E1]] = (gmt_grdfloat)(x_1_const * u[n_e+d_n[W1]] + x_0_const * u[n_e]);
962 }
963 }
964
965 /* Now set d2[]/dxdy = 0 at each of the 4 corners */
966
967 n = C->node_sw_corner; /* Just use shorthand in each expression */
968 u[n+d_n[SW]] = u[n+d_n[SE]] + u[n+d_n[NW]] - u[n+d_n[NE]];
969 n = C->node_nw_corner;
970 u[n+d_n[NW]] = u[n+d_n[NE]] + u[n+d_n[SW]] - u[n+d_n[SE]];
971 n = C->node_se_corner;
972 u[n+d_n[SE]] = u[n+d_n[SW]] + u[n+d_n[NE]] - u[n+d_n[NW]];
973 n = C->node_ne_corner;
974 u[n+d_n[NE]] = u[n+d_n[NW]] + u[n+d_n[SE]] - u[n+d_n[SW]];
975
976 /* Now set dC/dn = 0 at each edge */
977
978 for (col = 0, n_s = C->node_sw_corner, n_n = C->node_nw_corner; col < C->current_nx; col++, n_s++, n_n++) { /* set BC2 along south and north side */
979 /* South side */
980 u[n_s+d_n[S2]] = (gmt_grdfloat)(u[n_s+d_n[N2]] + C->eps_m2 * (u[n_s+d_n[NW]] + u[n_s+d_n[NE]]
981 - u[n_s+d_n[SW]] - u[n_s+d_n[SE]]) + C->two_plus_em2 * (u[n_s+d_n[S1]] - u[n_s+d_n[N1]]));
982 /* North side */
983 u[n_n+d_n[N2]] = (gmt_grdfloat)(u[n_n+d_n[S2]] + C->eps_m2 * (u[n_n+d_n[SW]] + u[n_n+d_n[SE]]
984 - u[n_n+d_n[NW]] - u[n_n+d_n[NE]]) + C->two_plus_em2 * (u[n_n+d_n[N1]] - u[n_n+d_n[S1]]));
985 }
986
987 for (row = 0, n_w = C->node_nw_corner, n_e = C->node_ne_corner; row < C->current_ny; row++, n_w += C->current_mx, n_e += C->current_mx) { /* set BC2 along west and east side */
988 if (C->periodic) { /* Set periodic boundary conditions in longitude */
989 u[n_w+d_n[W2]] = u[n_e+d_n[W2]]; /* West side */
990 u[n_e+d_n[E2]] = u[n_w+d_n[E2]]; /* East side */
991 }
992 else { /* Natural BCs */
993 /* West side */
994 u[n_w+d_n[W2]] = (gmt_grdfloat)(u[n_w+d_n[E2]] + C->eps_p2 * (u[n_w+d_n[NE]] + u[n_w+d_n[SE]]
995 - u[n_w+d_n[NW]] - u[n_w+d_n[SW]]) + C->two_plus_ep2 * (u[n_w+d_n[W1]] - u[n_w+d_n[E1]]));
996 /* East side */
997 u[n_e+d_n[E2]] = (gmt_grdfloat)(u[n_e+d_n[W2]] + C->eps_p2 * (u[n_e+d_n[NW]] + u[n_e+d_n[SW]]
998 - u[n_e+d_n[NE]] - u[n_e+d_n[SE]]) + C->two_plus_ep2 * (u[n_e+d_n[E1]] - u[n_e+d_n[W1]]));
999 }
1000 }
1001 }
1002
surface_iterate(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,int mode)1003 GMT_LOCAL uint64_t surface_iterate (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, int mode) {
1004 /* Main finite difference solver */
1005 uint64_t node, briggs_index, iteration_count = 0;
1006 unsigned int set, quadrant, current_max_iterations = C->max_iterations * C->current_stride;
1007 int col, row, k, *d_node = C->offset; /* Relative changes in node index from present node */
1008 unsigned char *status = C->status; /* Quadrant or status information for each node */
1009 char *mode_name[2] = {"node", "data"};
1010 bool finished;
1011 double current_limit = C->converge_limit / C->current_stride;
1012 double u_change, max_u_change, max_z_change, sum_bk_uk, u_00;
1013 gmt_grdfloat *b = NULL;
1014 gmt_grdfloat *u_new = C->Grid->data, *u_old = C->Grid->data;
1015
1016 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Starting iterations, mode = %s Max iterations = %d [stride = %d]\n", mode_name[mode], current_max_iterations, C->current_stride);
1017
1018 sprintf (C->format, "%%4ld\t%%c\t%%8" PRIu64 "\t%s\t%s\t%%10" PRIu64 "\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
1019 if (C->logging) fprintf (C->fp_log, "%c Grid size = %d Mode = %c Convergence limit = %g -Z%d\n",
1020 GMT->current.setting.io_seg_marker[GMT_OUT], C->current_stride, C->mode_type[mode], current_limit, C->current_stride);
1021
1022 /* We need to do an even number of iterations so that the final result for this iteration resides in C->Grid->data */
1023 do {
1024
1025 surface_set_BCs (GMT, C, u_old); /* Set the boundary rows and columns */
1026
1027 briggs_index = 0; /* Reset the Briggs constraint table index */
1028 max_u_change = -1.0; /* Ensure max_u_change is < 0 for starters */
1029
1030 /* Now loop over all interior data nodes */
1031 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Iteration %d\n", iteration_count);
1032
1033 for (row = 0; row < C->current_ny; row++) { /* Loop over rows */
1034 node = C->node_nw_corner + row * C->current_mx; /* Node at left side of this row */
1035 for (col = 0; col < C->current_nx; col++, node++) { /* Loop over all columns */
1036 if (status[node] == SURFACE_IS_CONSTRAINED) { /* Data constraint fell exactly on the node, keep it as is */
1037 continue;
1038 }
1039
1040 /* Here we must estimate a solution via equations (A-4) [SURFACE_UNCONSTRAINED] or (A-7) [SURFACE_CONSTRAINED] */
1041 u_00 = 0.0; /* Start with zero, build updated solution for central node */
1042 set = (status[node] == SURFACE_IS_UNCONSTRAINED) ? SURFACE_UNCONSTRAINED : SURFACE_CONSTRAINED; /* Index to C->coeff set to use */
1043 for (k = 0; k < 12; k++) { /* This is either equation (A-4) or the corresponding part of (A-7), depending on the value of set */
1044 u_00 += (u_old[node+d_node[k]] * C->coeff[set][k]);
1045 }
1046 if (set == SURFACE_CONSTRAINED) { /* Solution is (A-7) and modifications depend on which quadrant the point lies in */
1047 b = C->Briggs[briggs_index].b; /* Shorthand to this node's Briggs b-array */
1048 quadrant = status[node]; /* Which quadrant did the point fall in? */
1049 for (k = 0, sum_bk_uk = 0.0; k < 4; k++) { /* Sum over b[k]*u[k] for nodes A-D in Fig A-1 */
1050 sum_bk_uk += b[k] * u_old[node+d_node[C->p[quadrant][k]]];
1051 }
1052 u_00 = (u_00 + C->a0_const_2 * (sum_bk_uk + b[4])) * b[5]; /* Add point E in Fig A-1 to sum_bk_uk and normalize */
1053 briggs_index++; /* Got to next sequential Briggs array index */
1054 }
1055 /* We now apply the over-relaxation: */
1056 u_00 = u_old[node] * C->relax_old + u_00 * C->relax_new;
1057
1058 if (C->constrained) { /* Must check that we don't exceed any imposed limits. */
1059 if (C->set_limit[LO] && !gmt_M_is_fnan (C->Bound[LO]->data[node]) && u_00 < C->Bound[LO]->data[node])
1060 u_00 = C->Bound[LO]->data[node];
1061 else if (C->set_limit[HI] && !gmt_M_is_fnan (C->Bound[HI]->data[node]) && u_00 > C->Bound[HI]->data[node])
1062 u_00 = C->Bound[HI]->data[node];
1063 }
1064 u_change = fabs (u_00 - u_old[node]); /* Change in node value between iterations */
1065 u_new[node] = (gmt_grdfloat)u_00; /* Our updated estimate at this node */
1066 if (u_change > max_u_change) max_u_change = u_change; /* Keep track of max u_change across all nodes */
1067 } /* End of loop over columns */
1068 } /* End of loop over rows [and possibly threads via OpenMP] */
1069
1070 iteration_count++; C->total_iterations++; /* Update iteration counts for this stride and for total */
1071 max_z_change = max_u_change * C->z_rms; /* Scale max_u_change back into original z units -> max_z_change */
1072 GMT_Report (GMT->parent, GMT_MSG_DEBUG, C->format,
1073 C->current_stride, C->mode_type[mode], iteration_count, max_z_change, current_limit, C->total_iterations);
1074 if (C->logging) fprintf (C->fp_log, "%d\t%c\t%" PRIu64 "\t%.8g\t%.8g\t%" PRIu64 "\n", C->current_stride, C->mode_type[mode], iteration_count, max_z_change, current_limit, C->total_iterations);
1075 finished = (max_z_change <= current_limit || iteration_count >= current_max_iterations);
1076
1077 } while (!finished);
1078
1079 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, C->format,
1080 C->current_stride, C->mode_type[mode], iteration_count, max_z_change, current_limit, C->total_iterations);
1081
1082 return (iteration_count);
1083 }
1084
surface_check_errors(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)1085 GMT_LOCAL void surface_check_errors (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
1086 /* Compute misfits at original data locations, This is only done at the
1087 * final grid resolution, hence current_stride == 1. */
1088
1089 uint64_t node, k;
1090 int row, col, *d_node = C->offset;
1091 unsigned char *status = C->status;
1092
1093 double x0, y0, dx, dy, mean_error = 0.0, mean_squared_error = 0.0, z_est, z_err, curvature, c;
1094 double du_dx, du_dy, d2u_dx2, d2u_dxdy, d2u_dy2, d3u_dx3, d3u_dx2dy, d3u_dxdy2, d3u_dy3;
1095
1096 gmt_grdfloat *u = C->Grid->data;
1097 struct GMT_GRID_HEADER *h = C->Grid->header;
1098 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
1099
1100 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Compute rms misfit and curvature.\n");
1101
1102 surface_set_BCs (GMT, C, u); /* First update the boundary values */
1103
1104 /* Estimate solution at all data constraints using 3rd order Taylor expansion from nearest node */
1105
1106 for (k = 0; k < C->npoints; k++) {
1107 row = (int)index_to_row (C->data[k].index, C->n_columns);
1108 col = (int)index_to_col (C->data[k].index, C->n_columns);
1109 node = row_col_to_node (row, col, C->mx);
1110 if (status[node] == SURFACE_IS_CONSTRAINED) continue; /* Since misfit by definition is zero so no point adding it */
1111 /* Get coordinates of this node */
1112 x0 = col_to_x (col, h->wesn[XLO], h->wesn[XHI], h->inc[GMT_X], h->n_columns);
1113 y0 = row_to_y (row, h->wesn[YLO], h->wesn[YHI], h->inc[GMT_Y], h->n_rows);
1114 /* Get dx,dy of data point away from this node */
1115 dx = x_to_fcol (C->data[k].x, x0, HH->r_inc[GMT_X]);
1116 dy = y_to_frow (C->data[k].y, y0, HH->r_inc[GMT_Y]);
1117
1118 du_dx = 0.5 * (u[node+d_node[E1]] - u[node+d_node[W1]]);
1119 du_dy = 0.5 * (u[node+d_node[N1]] - u[node+d_node[S1]]);
1120 d2u_dx2 = u[node+d_node[E1]] + u[node+d_node[W1]] - 2 * u[node];
1121 d2u_dy2 = u[node+d_node[N1]] + u[node+d_node[S1]] - 2 * u[node];
1122 d2u_dxdy = 0.25 * (u[node+d_node[NE]] - u[node+d_node[NW]]
1123 - u[node+d_node[SE]] + u[node+d_node[SW]]);
1124 d3u_dx3 = 0.5 * (u[node+d_node[E2]] - 2 * u[node+d_node[E1]]
1125 + 2 * u[node+d_node[W1]] - u[node+d_node[W2]]);
1126 d3u_dy3 = 0.5 * (u[node+d_node[N2]] - 2 * u[node+d_node[N1]]
1127 + 2 * u[node+d_node[S1]] - u[node+d_node[S2]]);
1128 d3u_dx2dy = 0.5 * ((u[node+d_node[NE]] + u[node+d_node[NW]] - 2 * u[node+d_node[N1]])
1129 - (u[node+d_node[SE]] + u[node+d_node[SW]] - 2 * u[node+d_node[S1]]));
1130 d3u_dxdy2 = 0.5 * ((u[node+d_node[NE]] + u[node+d_node[SE]] - 2 * u[node+d_node[E1]])
1131 - (u[node+d_node[NW]] + u[node+d_node[SW]] - 2 * u[node+d_node[W1]]));
1132
1133 /* Compute the 3rd order Taylor approximation from current node */
1134
1135 z_est = u[node] + dx * (du_dx + dx * ((0.5 * d2u_dx2) + dx * (d3u_dx3 / 6.0)))
1136 + dy * (du_dy + dy * ((0.5 * d2u_dy2) + dy * (d3u_dy3 / 6.0)))
1137 + dx * dy * (d2u_dxdy) + (0.5 * dx * d3u_dx2dy) + (0.5 * dy * d3u_dxdy2);
1138
1139 z_err = z_est - C->data[k].z; /* Misfit between surface estimate and observation */
1140 mean_error += z_err;
1141 mean_squared_error += (z_err * z_err);
1142 }
1143 mean_error /= C->npoints;
1144 mean_squared_error = sqrt (mean_squared_error / C->npoints);
1145
1146 /* Compute the total curvature of the grid */
1147
1148 curvature = 0.0;
1149 gmt_M_grd_loop (GMT, C->Grid, row, col, node) {
1150 c = u[node+d_node[E1]] + u[node+d_node[W1]] + u[node+d_node[N1]] + u[node+d_node[S1]] - 4.0 * u[node+d_node[E1]];
1151 curvature += (c * c);
1152 }
1153
1154 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Fit info: N data points N nodes\tmean error\trms error\tcurvature\n");
1155 sprintf (C->format,"\t%%8ld\t%%8ld\t%s\t%s\t%s\n", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
1156 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, C->format, C->npoints, C->nxny, mean_error, mean_squared_error, curvature);
1157 }
1158
surface_remove_planar_trend(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)1159 GMT_LOCAL void surface_remove_planar_trend (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
1160 /* Fit LS plane and remove trend from our (x,y,z) input data; we add trend to grid before output.
1161 * Note: Here, x and y are first converted to fractional grid spacings from 0 to {n_columns,n_rows}-1.
1162 * Hence the same scheme is used by replace_planar trend (i.e., just use row,col as coordinates).
1163 * Note: The plane is fit to the original data z-values before normalizing by rms. */
1164 uint64_t k;
1165 double a, b, c, d, xx, y_up, zz, sx, sy, sz, sxx, sxy, sxz, syy, syz;
1166 struct GMT_GRID_HEADER *h = C->Grid->header;
1167 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (h);
1168
1169 sx = sy = sz = sxx = sxy = sxz = syy = syz = 0.0;
1170
1171 for (k = 0; k < C->npoints; k++) { /* Sum up normal equation terms */
1172 xx = x_to_fcol (C->data[k].x, h->wesn[XLO], HH->r_inc[GMT_X]); /* Distance from west to this point */
1173 y_up = y_to_frow (C->data[k].y, h->wesn[YLO], HH->r_inc[GMT_Y]); /* Distance from south to this point */
1174 zz = C->data[k].z;
1175 sx += xx; sy += y_up; sz += zz; sxx += (xx * xx);
1176 sxy += (xx * y_up); sxz += (xx * zz); syy += (y_up * y_up); syz += (y_up * zz);
1177 }
1178
1179 d = C->npoints*sxx*syy + 2*sx*sy*sxy - C->npoints*sxy*sxy - sx*sx*syy - sy*sy*sxx;
1180
1181 if (d == 0.0) { /* When denominator is zero we have a horizontal plane */
1182 C->plane_icept = C->plane_sx = C->plane_sy = 0.0;
1183 return;
1184 }
1185
1186 a = sz*sxx*syy + sx*sxy*syz + sy*sxy*sxz - sz*sxy*sxy - sx*sxz*syy - sy*syz*sxx;
1187 b = C->npoints*sxz*syy + sz*sy*sxy + sy*sx*syz - C->npoints*sxy*syz - sz*sx*syy - sy*sy*sxz;
1188 c = C->npoints*sxx*syz + sx*sy*sxz + sz*sx*sxy - C->npoints*sxy*sxz - sx*sx*syz - sz*sy*sxx;
1189
1190 C->plane_icept = a / d;
1191 C->plane_sx = b / d;
1192 C->plane_sy = c / d;
1193 if (C->periodic) C->plane_sx = 0.0; /* Cannot have x-trend for periodic geographic data */
1194
1195 for (k = 0; k < C->npoints; k++) { /* Now remove this plane from the data constraints */
1196 xx = x_to_fcol (C->data[k].x, h->wesn[XLO], HH->r_inc[GMT_X]);
1197 y_up = y_to_frow (C->data[k].y, h->wesn[YLO], HH->r_inc[GMT_Y]);
1198 C->data[k].z -= (gmt_grdfloat)evaluate_plane(C, xx, y_up);
1199 }
1200
1201 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Plane fit z = %g + (%g * col) + (%g * row)\n", C->plane_icept, C->plane_sx, C->plane_sy);
1202 }
1203
surface_restore_planar_trend(struct SURFACE_INFO * C)1204 GMT_LOCAL void surface_restore_planar_trend (struct SURFACE_INFO *C) {
1205 /* Scale grid back up by the data rms and restore the least-square plane.
1206 * Note: In determining the plane and in evaluating it, remember that the
1207 * x and y coordinates needed are not data coordinates but fractional col
1208 * and row distances from an origin at the lower left (southwest corner).
1209 * This means the y-values are positive up and increase in the opposite
1210 * direction than how rows increase. Hence the use of y_up below. */
1211 unsigned int row, col;
1212 uint64_t node;
1213 gmt_grdfloat *u = C->Grid->data;
1214 double y_up; /* Measure y up from south in fractional rows */
1215
1216 for (row = 0; row < C->Grid->header->n_rows; row++) {
1217 y_up = (double)(C->Grid->header->n_rows - row - 1); /* # of rows from south (where y_up = 0) to this node */
1218 node = row_col_to_node (row, 0, C->current_mx); /* Node index at left end of interior row */
1219 for (col = 0; col < C->Grid->header->n_columns; col++, node++) /* March across this row */
1220 u[node] = (gmt_grdfloat)((u[node] * C->z_rms) + (evaluate_plane (C, col, y_up)));
1221 }
1222 }
1223
surface_throw_away_unusables(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)1224 GMT_LOCAL void surface_throw_away_unusables (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
1225 /* We eliminate data which will become unusable on the final iteration, when current_stride = 1.
1226 It assumes current_stride = 1 and that surface_set_grid_parameters has been called.
1227 We sort, mark redundant data as SURFACE_OUTSIDE, and sort again, chopping off the excess.
1228 */
1229
1230 uint64_t last_index = UINTMAX_MAX, n_outside = 0, k;
1231
1232 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Eliminate data points that are not nearest a node.\n");
1233
1234 /* Sort the data */
1235
1236 qsort_r (C->data, C->npoints, sizeof (struct SURFACE_DATA), surface_compare_points, &(C->info));
1237
1238 /* If more than one datum is indexed to the same node, only the first should be kept.
1239 Mark the additional ones as SURFACE_OUTSIDE.
1240 */
1241 for (k = 0; k < C->npoints; k++) {
1242 if (C->data[k].index == last_index) { /* Same node but further away than our guy */
1243 C->data[k].index = SURFACE_OUTSIDE;
1244 n_outside++;
1245 }
1246 else { /* New index, just update last_index */
1247 last_index = C->data[k].index;
1248 }
1249 }
1250
1251 if (n_outside) { /* Sort again; this time the SURFACE_OUTSIDE points will be sorted to end of the array */
1252 qsort_r (C->data, C->npoints, sizeof (struct SURFACE_DATA), surface_compare_points, &(C->info));
1253 C->npoints -= n_outside; /* Effectively chopping off the eliminated points */
1254 C->data = gmt_M_memory (GMT, C->data, C->npoints, struct SURFACE_DATA); /* Adjust memory accordingly */
1255 GMT_Report (GMT->parent, GMT_MSG_WARNING, "%" PRIu64 " unusable points were supplied; these will be ignored.\n", n_outside);
1256 GMT_Report (GMT->parent, GMT_MSG_WARNING, "You should have pre-processed the data with block-mean, -median, or -mode.\n");
1257 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Check that previous processing steps write results with enough decimals.\n");
1258 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Possibly some data were half-way between nodes and subject to IEEE 754 rounding.\n");
1259 }
1260 }
1261
surface_rescale_z_values(struct GMT_CTRL * GMT,struct SURFACE_INFO * C)1262 GMT_LOCAL int surface_rescale_z_values (struct GMT_CTRL *GMT, struct SURFACE_INFO *C) {
1263 /* Find and normalize data by their rms value */
1264 uint64_t k;
1265 double ssz = 0.0;
1266
1267 for (k = 0; k < C->npoints; k++) ssz += (C->data[k].z * C->data[k].z);
1268 C->z_rms = sqrt (ssz / C->npoints);
1269 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Normalize detrended data constraints by z rms = %g\n", C->z_rms);
1270
1271 if (C->z_rms < GMT_CONV8_LIMIT) {
1272 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Input data lie exactly on a plane.\n");
1273 C->r_z_rms = C->z_rms = 1.0;
1274 return (1); /* Flag to tell the main to just write out the plane */
1275 }
1276 else
1277 C->r_z_rms = 1.0 / C->z_rms;
1278
1279 for (k = 0; k < C->npoints; k++) C->data[k].z *= (gmt_grdfloat)C->r_z_rms;
1280
1281 if (C->converge_limit == 0.0 || C->converge_mode == BY_PERCENT) { /* Set default values for convergence criteria */
1282 unsigned int ppm;
1283 double limit = (C->converge_mode == BY_PERCENT) ? C->converge_limit : SURFACE_CONV_LIMIT;
1284 ppm = urint (limit / 1.0e-6);
1285 C->converge_limit = limit * C->z_rms; /* i.e., 100 ppm of L2 scale */
1286 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Select default convergence limit of %g (%u ppm of L2 scale)\n", C->converge_limit, ppm);
1287 }
1288 return (0);
1289 }
1290
surface_suggest_sizes(struct GMT_CTRL * GMT,struct GMT_GRID * G,unsigned int factors[],unsigned int n_columns,unsigned int n_rows,bool pixel)1291 GMT_LOCAL void surface_suggest_sizes (struct GMT_CTRL *GMT, struct GMT_GRID *G, unsigned int factors[], unsigned int n_columns, unsigned int n_rows, bool pixel) {
1292 /* Calls gmt_optimal_dim_for_surface to determine if there are
1293 * better choices for n_columns, n_rows that might speed up calculations
1294 * by having many more common factors.
1295 *
1296 * W. H. F. Smith, 26 Feb 1992. */
1297
1298 unsigned int k;
1299 unsigned int n_sug = 0; /* Number of suggestions found */
1300 struct GMT_SURFACE_SUGGESTION *sug = NULL;
1301
1302 n_sug = gmt_optimal_dim_for_surface (GMT, factors, n_columns, n_rows, &sug);
1303
1304 if (n_sug) { /* We did find some suggestions, report them (up to the first 10 suggestions) */
1305 char region[GMT_LEN128] = {""}, buffer[GMT_LEN128] = {""};
1306 bool lat_bad = false;
1307 unsigned int m, save_range = GMT->current.io.geo.range;
1308 double w, e, s, n;
1309 GMT->current.io.geo.range = GMT_IS_GIVEN_RANGE; /* Override this setting explicitly */
1310 for (k = 0; k < n_sug && k < 10; k++) {
1311 m = sug[k].n_columns - (G->header->n_columns - 1); /* Additional nodes needed in x to give more factors */
1312 w = G->header->wesn[XLO] - (m/2)*G->header->inc[GMT_X]; /* Potential revised w/e extent */
1313 e = G->header->wesn[XHI] + (m/2)*G->header->inc[GMT_X];
1314 if (m%2) e += G->header->inc[GMT_X];
1315 m = sug[k].n_rows - (G->header->n_rows - 1); /* Additional nodes needed in y to give more factors */
1316 s = G->header->wesn[YLO] - (m/2)*G->header->inc[GMT_Y]; /* Potential revised s/n extent */
1317 n = G->header->wesn[YHI] + (m/2)*G->header->inc[GMT_Y];
1318 if (!lat_bad && gmt_M_y_is_lat (GMT, GMT_IN) && (s < -90.0 || n > 90.0))
1319 lat_bad = true;
1320 if (m%2) n += G->header->inc[GMT_Y];
1321 if (pixel) { /* Since we already added 1/2 pixel we need to undo that here so the report matches original phase */
1322 w -= G->header->inc[GMT_X] / 2.0; e -= G->header->inc[GMT_X] / 2.0;
1323 s -= G->header->inc[GMT_Y] / 2.0; n -= G->header->inc[GMT_Y] / 2.0;
1324 }
1325 gmt_ascii_format_col (GMT, buffer, w, GMT_OUT, GMT_X);
1326 sprintf (region, "-R%s/", buffer);
1327 gmt_ascii_format_col (GMT, buffer, e, GMT_OUT, GMT_X);
1328 strcat (region, buffer); strcat (region, "/");
1329 gmt_ascii_format_col (GMT, buffer, s, GMT_OUT, GMT_Y);
1330 strcat (region, buffer); strcat (region, "/");
1331 gmt_ascii_format_col (GMT, buffer, n, GMT_OUT, GMT_Y);
1332 strcat (region, buffer);
1333 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Hint: Choosing %s [n_columns = %d, n_rows = %d] might cut run time by a factor of %.8g\n",
1334 region, sug[k].n_columns, sug[k].n_rows, sug[k].factor);
1335 }
1336 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Hint: After completion you can recover the desired region via gmt grdcut\n");
1337 if (lat_bad) {
1338 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Note: One or more of the suggested south/north bounds exceed the allowable range [-90/90]\n");
1339 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "A workaround is to use -fx to only consider x as geographic longitudes\n");
1340 }
1341 gmt_M_free (GMT, sug);
1342 GMT->current.io.geo.range = save_range;
1343 }
1344 else
1345 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Cannot suggest any n_columns,n_rows better than your current -R -I settings.\n");
1346 return;
1347 }
1348
surface_init_parameters(struct SURFACE_INFO * C,struct SURFACE_CTRL * Ctrl)1349 GMT_LOCAL void surface_init_parameters (struct SURFACE_INFO *C, struct SURFACE_CTRL *Ctrl) {
1350 /* Place program options into the surface struct. This was done this way
1351 * since surface.c relied heavily on global variables which are a no-no
1352 * in GMT5. The simplest solution was to collect all those variables into
1353 * a single structure and pass a pointer to that structure to functions.
1354 */
1355 if (Ctrl->S.active) { /* Gave a search radius; adjust if minutes or seconds where specified */
1356 if (Ctrl->S.unit == 'm') Ctrl->S.radius /= 60.0;
1357 if (Ctrl->S.unit == 's') Ctrl->S.radius /= 3600.0;
1358 }
1359 C->radius = Ctrl->S.radius;
1360 C->relax_new = Ctrl->Z.value;
1361 C->relax_old = 1.0 - C->relax_new;
1362 C->max_iterations = Ctrl->N.value;
1363 C->radius = Ctrl->S.radius;
1364 C->limit_file[LO] = Ctrl->L.file[LO];
1365 C->limit_file[HI] = Ctrl->L.file[HI];
1366 C->set_limit[LO] = Ctrl->L.mode[LO];
1367 C->set_limit[HI] = Ctrl->L.mode[HI];
1368 C->limit[LO] = Ctrl->L.limit[LO];
1369 C->limit[HI] = Ctrl->L.limit[HI];
1370 C->boundary_tension = Ctrl->T.b_tension;
1371 C->interior_tension = Ctrl->T.i_tension;
1372 C->alpha = Ctrl->A.value;
1373 C->converge_limit = Ctrl->C.value;
1374 C->converge_mode = Ctrl->C.mode;
1375 C->n_alloc = GMT_INITIAL_MEM_ROW_ALLOC;
1376 C->z_rms = 1.0;
1377 C->r_z_rms = 1.0;
1378 C->mode_type[0] = 'I'; /* I means exclude data points when iterating */
1379 C->mode_type[1] = 'D'; /* D means include data points when iterating */
1380 C->n_columns = C->Grid->header->n_columns;
1381 C->n_rows = C->Grid->header->n_rows;
1382 C->nxny = C->Grid->header->nm;
1383 C->mx = C->Grid->header->mx;
1384 C->my = C->Grid->header->my;
1385 C->mxmy = C->Grid->header->size;
1386 gmt_M_memcpy (C->p, p, 20, unsigned int);
1387
1388 gmt_M_memcpy (C->info.wesn, C->Grid->header->wesn, 4, double);
1389 }
1390
surface_find_closest_point(double * x,double * y,double * z,uint64_t k,double x0,double y0,double half_dx,double half_dy,double * xx,double * yy,double * zz)1391 GMT_LOCAL double surface_find_closest_point (double *x, double *y, double *z, uint64_t k, double x0, double y0, double half_dx, double half_dy, double *xx, double *yy, double *zz) {
1392 /* Find the point (xx,yy) on the line from (x[k-1],y[k-1]) to (x[k], y[k]) that is closest to (x0,y0). If (xx,yy)
1393 * is outside the end of the line segment or outside the bin then we return r = DBL_MAX */
1394 double dx, dy, a, r= DBL_MAX; /* Initialize distance from (x0,y0) to nearest point (xx,yy) measured orthogonally onto break line */
1395 uint64_t km1 = k - 1;
1396 dx = x[k] - x[km1]; dy = y[k] - y[km1];
1397 if (gmt_M_is_zero (dx)) { /* Break line is vertical */
1398 if ((y[k] <= y0 && y[km1] > y0) || (y[km1] <= y0 && y[k] > y0)) { /* Nearest point is in same bin */
1399 *xx = x[k]; *yy = y0;
1400 r = fabs (*xx - x0);
1401 *zz = z[km1] + (z[k] - z[km1]) * (*yy - y[km1]) / dy;
1402 }
1403 }
1404 else if (gmt_M_is_zero (dy)) { /* Break line is horizontal */
1405 if ((x[k] <= x0 && x[km1] > x0) || (x[km1] <= x0 && x[k] > x0)) { /* Nearest point in same bin */
1406 *xx = x0; *yy = y[k];
1407 r = fabs (*yy - y0);
1408 *zz = z[km1] + (z[k] - z[km1]) * (*xx - x[km1]) / dx;
1409 }
1410 }
1411 else { /* General case. Nearest orthogonal point may or may not be in bin, in which case r > r_prev */
1412 a = dy / dx; /* Slope of line */
1413 *xx = (y0 - y[km1] + a * x[km1] + x0 / a) / (a + 1.0/a);
1414 *yy = a * (*xx - x[k]) + y[k];
1415 if ((x[k] <= *xx && x[km1] > *xx) || (x[km1] <= *xx && x[k] > *xx)) { /* Orthonormal point found between the end points of line */
1416 if (fabs (*xx-x0) < half_dx && fabs (*yy-y0) < half_dy) { /* Yes, within this bin */
1417 r = hypot (*xx - x0, *yy - y0);
1418 *zz = z[km1] + (z[k] - z[km1]) * (*xx - x[km1]) / dx;
1419 }
1420 }
1421 }
1422 return r;
1423 }
1424
surface_interpolate_add_breakline(struct GMT_CTRL * GMT,struct SURFACE_INFO * C,struct GMT_DATATABLE * T,char * file,bool fix_z,double z_level)1425 GMT_LOCAL void surface_interpolate_add_breakline (struct GMT_CTRL *GMT, struct SURFACE_INFO *C, struct GMT_DATATABLE *T, char *file, bool fix_z, double z_level) {
1426 int srow, scol;
1427 uint64_t new_n = 0, n_int = 0, nb = 0;
1428 uint64_t k = 0, n, kmax = 0, kmin = 0, row, seg, node_this, node_prev;
1429 size_t n_alloc, n_alloc_b;
1430 double dx, dy, dz, r, r_this, r_min, x0_prev, y0_prev, x0_this, y0_this;
1431 double xx, yy, zz, half_dx, half_dy, zmin = DBL_MAX, zmax = -DBL_MAX;
1432 double *xline = NULL, *yline = NULL, *zline = NULL;
1433 double *x = NULL, *y = NULL, *z = NULL, *xb = NULL, *yb = NULL, *zb = NULL;
1434 char fname1[GMT_LEN256] = {""}, fname2[GMT_LEN256] = {""};
1435 FILE *fp1 = NULL, *fp2 = NULL;
1436
1437 if (file) {
1438 sprintf (fname1, "%s.int", file);
1439 sprintf (fname2, "%s.final", file);
1440 if ((fp1 = fopen (fname1, "w")) == NULL) {
1441 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create file %s\n", fname1);
1442 return;
1443 }
1444 if ((fp2 = fopen (fname2, "w")) == NULL) {
1445 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create file %s\n", fname1);
1446 fclose (fp1);
1447 return;
1448 }
1449 }
1450 /* Add constraints from breaklines */
1451 /* Reduce breaklines to the nearest point per node of cells crossed */
1452
1453 n_alloc = n_alloc_b = GMT_INITIAL_MEM_ROW_ALLOC;
1454 xb = gmt_M_memory (GMT, NULL, n_alloc_b, double);
1455 yb = gmt_M_memory (GMT, NULL, n_alloc_b, double);
1456 zb = gmt_M_memory (GMT, NULL, n_alloc_b, double);
1457
1458 x = gmt_M_memory (GMT, NULL, n_alloc, double);
1459 y = gmt_M_memory (GMT, NULL, n_alloc, double);
1460 z = gmt_M_memory (GMT, NULL, n_alloc, double);
1461
1462 half_dx = 0.5 * C->inc[GMT_X]; half_dy = 0.5 * C->inc[GMT_Y];
1463 for (seg = 0; seg < T->n_segments; seg++) {
1464 xline = T->segment[seg]->data[GMT_X];
1465 yline = T->segment[seg]->data[GMT_Y];
1466 if (!fix_z) zline = T->segment[seg]->data[GMT_Z];
1467 /* 1. Interpolate the breakline to ensure there are points in every bin that it crosses */
1468 if (file) fprintf (fp1, "> Segment %d\n", (int)seg);
1469 for (row = k = 0, new_n = 1; row < T->segment[seg]->n_rows - 1; row++) {
1470 dx = xline[row+1] - xline[row];
1471 dy = yline[row+1] - yline[row];
1472 if (!fix_z) dz = zline[row+1] - zline[row];
1473 /* Given point spacing and grid spacing, how many points to interpolate? */
1474 n_int = lrint (hypot (dx, dy) * MAX (C->r_inc[GMT_X], C->r_inc[GMT_Y])) + 1;
1475 new_n += n_int;
1476 if (n_alloc <= new_n) {
1477 n_alloc += MAX (GMT_CHUNK, n_int);
1478 x = gmt_M_memory (GMT, x, n_alloc, double);
1479 y = gmt_M_memory (GMT, y, n_alloc, double);
1480 z = gmt_M_memory (GMT, z, n_alloc, double);
1481 }
1482
1483 dx /= n_int;
1484 dy /= n_int;
1485 if (!fix_z) dz /= n_int;
1486 for (n = 0; n < n_int; k++, n++) {
1487 x[k] = xline[row] + n * dx;
1488 y[k] = yline[row] + n * dy;
1489 z[k] = (fix_z) ? z_level : zline[row] + n * dz;
1490 if (file) fprintf (fp1, "%g\t%g\t%g\n", x[k], y[k], z[k]);
1491 }
1492 }
1493 x[k] = xline[row]; y[k] = yline[row]; z[k] = (fix_z) ? z_level : zline[row];
1494 if (file) fprintf (fp1, "%g\t%g\t%g\n", x[k], y[k], z[k]);
1495
1496 /* 2. Go along the (x,y,z), k = 1:new_n line and find the closest point to each bin node */
1497 if (file) fprintf (fp2, "> Segment %d\n", (int)seg);
1498 scol = (int)x_to_col (x[0], C->Grid->header->wesn[XLO], C->r_inc[GMT_X]);
1499 srow = (int)y_to_row (y[0], C->Grid->header->wesn[YLO], C->r_inc[GMT_Y], C->current_ny);
1500 node_this = row_col_to_node (srow, scol, C->current_mx); /* The bin we are in */
1501 x0_this = col_to_x (scol, C->Grid->header->wesn[XLO], C->Grid->header->wesn[XHI], C->inc[GMT_X], C->current_nx); /* Node center point */
1502 y0_this = row_to_y (srow, C->Grid->header->wesn[YLO], C->Grid->header->wesn[YHI], C->inc[GMT_Y], C->current_ny);
1503 r_min = hypot (x[0] - x0_this, y[0] - y0_this); /* Distance from node center to start of breakline */
1504 xb[nb] = x[0]; yb[nb] = y[0]; zb[nb] = z[0]; /* Add this as our "nearest" breakline point (so far) for this bin */
1505 for (k = 1; k < new_n; k++) {
1506 /* Reset what is the previous point */
1507 node_prev = node_this;
1508 x0_prev = x0_this; y0_prev = y0_this;
1509 scol = (int)x_to_col (x[k], C->Grid->header->wesn[XLO], C->r_inc[GMT_X]);
1510 srow = (int)y_to_row (y[k], C->Grid->header->wesn[YLO], C->r_inc[GMT_Y], C->current_ny);
1511 x0_this = col_to_x (scol, C->Grid->header->wesn[XLO], C->Grid->header->wesn[XHI], C->inc[GMT_X], C->current_nx); /* Node center point */
1512 y0_this = row_to_y (srow, C->Grid->header->wesn[YLO], C->Grid->header->wesn[YHI], C->inc[GMT_Y], C->current_ny);
1513 node_this = row_col_to_node (srow, scol, C->current_mx);
1514 r_this = hypot (x[k] - x0_this, y[k] - y0_this);
1515 if (node_this == node_prev) { /* Both points in same bin, see if 2nd point is closer */
1516 if (r_this < r_min) { /* This point is closer than previous point */
1517 xb[nb] = x[k]; yb[nb] = y[k]; zb[nb] = z[k];
1518 r_min = r_this;
1519 }
1520 }
1521 /* Find point on line closest to prev bin center */
1522 r = surface_find_closest_point (x, y, z, k, x0_prev, y0_prev, half_dx, half_dy, &xx, &yy, &zz);
1523 if (r < r_min) { /* Yes, closer than previous point */
1524 xb[nb] = xx; yb[nb] = yy; zb[nb] = zz;
1525 r_min = r;
1526 }
1527 if (node_this != node_prev) { /* Find point on line closest to this bin center */
1528 if (file) fprintf (fp2, "%g\t%g\t%g\n", xb[nb], yb[nb], zb[nb]);
1529 nb++; /* OK, moving on from this bin */
1530 if (nb > n_alloc_b) {
1531 n_alloc_b += GMT_CHUNK;
1532 xb = gmt_M_memory (GMT, xb, n_alloc_b, double);
1533 yb = gmt_M_memory (GMT, yb, n_alloc_b, double);
1534 zb = gmt_M_memory (GMT, zb, n_alloc_b, double);
1535 }
1536 xb[nb] = x[k]; yb[nb] = y[k]; zb[nb] = z[k];
1537 r_min = r_this;
1538 r = surface_find_closest_point (x, y, z, k, x0_this, y0_this, half_dx, half_dy, &xx, &yy, &zz);
1539 if (r < r_min) { /* Yes, closer than previous point */
1540 xb[nb] = xx; yb[nb] = yy; zb[nb] = zz;
1541 r_min = r;
1542 }
1543 }
1544 }
1545 if (file) fprintf (fp2, "%g\t%g\t%g\n", xb[nb], yb[nb], zb[nb]);
1546 nb++;
1547 }
1548 if (file) {
1549 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Reinterpolated breakline saved to file %s\n", fname1);
1550 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Final breakline constraints saved to file %s\n", fname2);
1551 fclose (fp1);
1552 fclose (fp2);
1553 }
1554
1555 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Found %d breakline points, reinterpolated to %d points, reduced to %d points\n", (int)T->n_records, (int)new_n, (int)nb);
1556
1557 /* Now append the interpolated breakline to the C data structure */
1558
1559 k = C->npoints;
1560 C->data = gmt_M_memory (GMT, C->data, k+nb, struct SURFACE_DATA);
1561 C->z_mean *= k; /* It was already computed, reset it to the sum so we can add more and recalculate the mean */
1562 if (C->set_limit[LO] == DATA) /* Lower limit should equal minimum data found. Start with what we have so far and change if we find lower values */
1563 zmin = C->limit[LO];
1564 if (C->set_limit[HI] == DATA) /* Upper limit should equal maximum data found. Start with what we have so far and change if we find higher values */
1565 zmax = C->limit[HI];
1566
1567 for (n = 0; n < nb; n++) {
1568
1569 if (gmt_M_is_dnan (zb[n])) continue;
1570
1571 scol = (int)x_to_col (xb[n], C->Grid->header->wesn[XLO], C->r_inc[GMT_X]);
1572 if (scol < 0 || scol >= C->current_nx) continue;
1573 srow = (int)y_to_row (yb[n], C->Grid->header->wesn[YLO], C->r_inc[GMT_Y], C->current_ny);
1574 if (srow < 0 || srow >= C->current_ny) continue;
1575
1576 C->data[k].index = row_col_to_index (srow, scol, C->current_nx);
1577 C->data[k].x = (gmt_grdfloat)xb[n];
1578 C->data[k].y = (gmt_grdfloat)yb[n];
1579 C->data[k].z = (gmt_grdfloat)zb[n];
1580 C->data[k].kind = SURFACE_BREAKLINE; /* Mark as breakline constraint */
1581 if (zmin > zb[n]) zmin = z[n], kmin = k;
1582 if (zmax < zb[n]) zmax = z[n], kmax = k;
1583 k++;
1584 C->z_mean += zb[n];
1585 }
1586
1587 if (k != (C->npoints + nb)) /* We had some NaNs */
1588 C->data = gmt_M_memory (GMT, C->data, k, struct SURFACE_DATA);
1589
1590 C->npoints = k;
1591 C->z_mean /= k;
1592
1593 if (C->set_limit[LO] == DATA) /* Update our lower data-driven limit to the new minimum found */
1594 C->limit[LO] = C->data[kmin].z;
1595 if (C->set_limit[HI] == DATA) /* Update our upper data-driven limit to the new maximum found */
1596 C->limit[HI] = C->data[kmax].z;
1597
1598 gmt_M_free (GMT, x);
1599 gmt_M_free (GMT, y);
1600 gmt_M_free (GMT, z);
1601 gmt_M_free (GMT, xb);
1602 gmt_M_free (GMT, yb);
1603 gmt_M_free (GMT, zb);
1604 }
1605
New_Ctrl(struct GMT_CTRL * GMT)1606 static void *New_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */
1607 struct SURFACE_CTRL *C;
1608
1609 C = gmt_M_memory (GMT, NULL, 1, struct SURFACE_CTRL);
1610
1611 /* Initialize values whose defaults are not 0/false/NULL */
1612 C->N.value = SURFACE_MAX_ITERATIONS;
1613 C->A.value = 1.0; /* Real xinc == yinc in terms of distances */
1614 C->W.file = strdup ("surface_log.txt");
1615 C->Z.value = SURFACE_OVERRELAXATION;
1616
1617 return (C);
1618 }
1619
Free_Ctrl(struct GMT_CTRL * GMT,struct SURFACE_CTRL * C)1620 static void Free_Ctrl (struct GMT_CTRL *GMT, struct SURFACE_CTRL *C) { /* Deallocate control structure */
1621 if (!C) return;
1622 gmt_M_str_free (C->G.file);
1623 if (C->D.file) gmt_M_str_free (C->D.file);
1624 if (C->L.file[LO]) gmt_M_str_free (C->L.file[LO]);
1625 if (C->L.file[HI]) gmt_M_str_free (C->L.file[HI]);
1626 if (C->M.arg) gmt_M_str_free (C->M.arg);
1627 if (C->W.file) gmt_M_str_free (C->W.file);
1628 gmt_M_free (GMT, C);
1629 }
1630
usage(struct GMTAPI_CTRL * API,int level)1631 static int usage (struct GMTAPI_CTRL *API, int level) {
1632 unsigned int ppm;
1633 const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
1634 if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
1635 GMT_Usage (API, 0, "usage: %s [<table>] -G%s %s %s [-A<aspect_ratio>|m] [-C<convergence_limit>] "
1636 "[-D<breakline>[+z[<zlevel>]]] [%s] [-Ll|u<limit>] [-M<radius>] [-N<n_iterations>] [-Q] "
1637 "[-S<search_radius>[m|s]] [-T[b|i]<tension>] [%s] [-W[<logfile>]] [-Z<over_relaxation>] "
1638 "[%s] [%s] [%s] [%s] [%s] [%s] [%s [%s] [%s] [%s] %s[%s] [%s]\n",
1639 name, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT, GMT_J_OPT, GMT_V_OPT, GMT_a_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_e_OPT, GMT_f_OPT,
1640 GMT_h_OPT, GMT_i_OPT, GMT_qi_OPT, GMT_r_OPT, GMT_w_OPT, GMT_x_OPT, GMT_colon_OPT, GMT_PAR_OPT);
1641
1642 if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);
1643 ppm = urint (SURFACE_CONV_LIMIT / 1e-6); /* Default convergence criteria */
1644
1645 GMT_Message (API, GMT_TIME_NONE, " REQUIRED ARGUMENTS:\n");
1646 GMT_Option (API, "<");
1647 gmt_outgrid_syntax (API, 'G', "Sets name of the output grid file");
1648 GMT_Option (API, "I,R");
1649 GMT_Message (API, GMT_TIME_NONE, "\n OPTIONAL ARGUMENTS:\n");
1650 GMT_Usage (API, 1, "\n-A<aspect_ratio>|m");
1651 GMT_Usage (API, -2, "Set <aspect-ratio> [Default = 1 gives an isotropic solution], "
1652 "i.e., <xinc> and <yinc> are assumed to give derivatives of equal weight; if not, specify "
1653 "<aspect_ratio> such that <yinc> = <xinc> / <aspect_ratio>. "
1654 "If gridding lon,lat use -Am to set <aspect_ratio> = cosine(middle of lat range).");
1655 GMT_Usage (API, 1, "\n-C<convergence_limit>");
1656 GMT_Usage (API, -2, "Set final convergence limit; iteration stops when max |change| < <convergence_limit>. "
1657 "Default will choose %g of the rms of your z data after removing L2 plane (%u ppm precision). "
1658 "Enter your own convergence limit in the same units as your z data.", SURFACE_CONV_LIMIT, ppm);
1659 GMT_Usage (API, 1, "\n-D<breakline>[+z[<zlevel>]]");
1660 GMT_Usage (API, -2, "Use xyz data in the <breakline> file as a 'soft breakline'. Optional modifier:");
1661 GMT_Usage (API, 3, "+z Override any z from the <breakline> file with the appended <z_level> [0].");
1662 GMT_Usage (API, 1, "\n%s", GMT_J_OPT);
1663 GMT_Usage (API, -2, "Select the data map projection. This projection is only used to add a CRS info to the "
1664 "grid formats that support it, i.e., netCDF, GeoTIFF, and others supported by GDAL.");
1665 GMT_Usage (API, 1, "\n-Ll|u<limit>");
1666 GMT_Usage (API, -2, "Constrain the range of output values; append directive and value, repeatable:");
1667 GMT_Usage (API, 3, "l: Set lower limit; forces solution to be >= <limit>.");
1668 GMT_Usage (API, 3, "u: Set upper limit; forces solution to be <= <limit>.");
1669 GMT_Usage (API, -2, "Note: <limit> can be any number, or the letter d for min (or max) input data value, "
1670 "or the filename of a grid with bounding values [Default solution is unconstrained]. "
1671 "Example: -Ll0 enforces a non-negative solution.");
1672 gmt_dist_syntax (API->GMT, "M<radius>", "Set maximum radius for masking the grid away from data points [no masking].");
1673 GMT_Usage (API, -2, "For Cartesian grids with different x and y units you may append <xlim>/<ylim>; "
1674 "this fills all nodes within the rectangular area of the given half-widths. "
1675 "One can also achieve the rectangular selection effect by using the -M<n_cells>c "
1676 "form. Here <n_cells> means the number of cells around the data point. As an example, "
1677 "-M0c means that only the cell where the point lies is retained, -M1c keeps one cell "
1678 "beyond that (i.e. makes a 3x3 neighborhood), and so on.");
1679 GMT_Usage (API, 1, "\n-N<n_iterations>");
1680 GMT_Usage (API, -2, "Set maximum number of iterations in the final cycle; default = %d.", SURFACE_MAX_ITERATIONS);
1681 GMT_Usage (API, 1, "\n-S<search_radius>[m|s]");
1682 GMT_Usage (API, -2, "Set <search_radius> to initialize grid; default = 0 will skip this step. "
1683 "This step is slow and not needed unless grid dimensions are pathological; "
1684 "i.e., have few or no common factors. "
1685 "Append m or s to give <search_radius> in minutes or seconds.");
1686 GMT_Usage (API, 1, "\n-T[b|i]<tension>");
1687 GMT_Usage (API, -2, "Add tension to the gridding equation; use a value between 0 and 1. "
1688 "Default = 0 gives minimum curvature (smoothest; bicubic) solution. "
1689 "1 gives a harmonic spline solution (local max/min occur only at data points). "
1690 "Typically, 0.25 or more is good for potential field (smooth) data; "
1691 "0.5-0.75 or so for topography. We encourage you to experiment. Optional directives:");
1692 GMT_Usage (API, 3, "b: Set tension in boundary conditions only.");
1693 GMT_Usage (API, 3, "i: Set tension in interior equations only.n");
1694 GMT_Usage (API, -2, "Note: Without a directive we set tension for both to same value.");
1695 GMT_Usage (API, 1, "\n-Q Query for grid sizes that might run faster than your selected -R -I, then exit.");
1696 GMT_Option (API, "V");
1697 GMT_Usage (API, 1, "\n-W[<logfile>]");
1698 GMT_Usage (API, -2, "Write convergence information to a log file [surface_log.txt].");
1699 GMT_Usage (API, 1, "\n-Z<over_relaxation>");
1700 GMT_Usage (API, -2, "Change over-relaxation parameter [Default = %g]. Use a value "
1701 "between 1 and 2. Larger number accelerates convergence but can be unstable. "
1702 "Use 1 if you want to be sure to have (slow) stable convergence.", SURFACE_OVERRELAXATION);
1703 GMT_Option (API, "a,bi3,di,e,f,h,i,qi,r,w,x,:,.");
1704 if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Note: Geographic data with 360-degree range use periodic boundary condition in longitude. "
1705 "For additional details, see Smith & Wessel, Geophysics, 55, 293-305, 1990.");
1706
1707 return (GMT_MODULE_USAGE);
1708 }
1709
parse(struct GMT_CTRL * GMT,struct SURFACE_CTRL * Ctrl,struct GMT_OPTION * options)1710 static int parse (struct GMT_CTRL *GMT, struct SURFACE_CTRL *Ctrl, struct GMT_OPTION *options) {
1711 /* Parse the options provided and set parameters in CTRL structure.
1712 * Any GMT common options will override values set previously by other commands.
1713 * It also replaces any file names specified as input or output with the data ID
1714 * returned when registering these sources/destinations with the API.
1715 */
1716
1717 unsigned int n_errors = 0, k, end;
1718 char modifier, *c = NULL, *d = NULL;
1719 struct GMT_OPTION *opt = NULL;
1720 struct GMTAPI_CTRL *API = GMT->parent;
1721
1722 for (opt = options; opt; opt = opt->next) {
1723 switch (opt->option) {
1724
1725 case '<': /* Skip input files */
1726 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(opt->arg))) n_errors++;;
1727 break;
1728
1729 /* Processes program-specific parameters */
1730
1731 case 'A':
1732 n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
1733 Ctrl->A.active = true;
1734 if (opt->arg[0] == 'm')
1735 Ctrl->A.mode = 1;
1736 else
1737 Ctrl->A.value = atof (opt->arg);
1738 break;
1739 case 'C':
1740 n_errors += gmt_M_repeated_module_option (API, Ctrl->C.active);
1741 Ctrl->C.active = true;
1742 Ctrl->C.value = atof (opt->arg);
1743 if (strchr (opt->arg, '%')) { /* Gave convergence in percent */
1744 Ctrl->C.mode = BY_PERCENT;
1745 Ctrl->C.value *= 0.01;
1746 }
1747 break;
1748 case 'D':
1749 n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
1750 if ((d = strstr (opt->arg, "+d"))) {
1751 d[0] = '\0'; /* Temporarily chop off +d part */
1752 Ctrl->D.debug = true;
1753 }
1754 if ((c = strstr (opt->arg, "+z"))) {
1755 c[0] = '\0'; /* Temporarily chop off +z part */
1756 if (c[2]) Ctrl->D.z = atof (&c[2]); /* Get the constant z-value [0] */
1757 Ctrl->D.fix_z = true;
1758 }
1759 Ctrl->D.active = true;
1760 if (opt->arg[0]) Ctrl->D.file = strdup (opt->arg);
1761 if (GMT_Get_FilePath (API, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->D.file))) n_errors++;
1762 if (c) c[0] = '+'; /* Restore original string */
1763 if (d) d[0] = '+'; /* Restore original string */
1764 break;
1765 case 'G':
1766 n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
1767 Ctrl->G.active = true;
1768 if (opt->arg[0]) Ctrl->G.file = strdup (opt->arg);
1769 if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file))) n_errors++;
1770 break;
1771 case 'I':
1772 n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
1773 Ctrl->I.active = true;
1774 n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
1775 break;
1776 case 'J': /* We have this gal here be cause it needs to be processed separately (after -R) */
1777 n_errors += gmt_M_repeated_module_option (API, Ctrl->J.active);
1778 Ctrl->J.active = true;
1779 Ctrl->J.projstring = strdup(opt->arg);
1780 break;
1781 case 'L': /* Set limits */
1782 switch (opt->arg[0]) {
1783 case 'l': case 'u': /* Lower or upper limits */
1784 end = (opt->arg[0] == 'l') ? LO : HI; /* Which one it is */
1785 n_errors += gmt_M_repeated_module_option (API, Ctrl->L.active[end]);
1786 Ctrl->L.active[end] = true;
1787 n_errors += gmt_M_check_condition (GMT, opt->arg[1] == 0, "Option -L%c: No argument given\n", opt->arg[0]);
1788 Ctrl->L.file[end] = strdup (&opt->arg[1]);
1789 if (!gmt_access (GMT, Ctrl->L.file[end], F_OK)) /* File exists */
1790 Ctrl->L.mode[end] = SURFACE;
1791 else if (Ctrl->L.file[end][0] == 'd') /* Use data minimum */
1792 Ctrl->L.mode[end] = DATA;
1793 else {
1794 Ctrl->L.mode[end] = VALUE; /* Use given value */
1795 Ctrl->L.limit[end] = atof (&opt->arg[1]);
1796 }
1797 break;
1798 default:
1799 n_errors++;
1800 break;
1801 }
1802 break;
1803 case 'M':
1804 n_errors += gmt_M_repeated_module_option (API, Ctrl->M.active);
1805 Ctrl->M.active = true;
1806 Ctrl->M.arg = strdup (opt->arg);
1807 break;
1808 case 'N':
1809 n_errors += gmt_M_repeated_module_option (API, Ctrl->N.active);
1810 Ctrl->N.active = true;
1811 Ctrl->N.value = atoi (opt->arg);
1812 break;
1813 case 'Q':
1814 n_errors += gmt_M_repeated_module_option (API, Ctrl->Q.active);
1815 Ctrl->Q.active = true;
1816 break;
1817 case 'S':
1818 n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
1819 Ctrl->S.active = true;
1820 Ctrl->S.radius = atof (opt->arg);
1821 Ctrl->S.unit = opt->arg[strlen(opt->arg)-1];
1822 if (Ctrl->S.unit == 'c' && gmt_M_compat_check (GMT, 4)) {
1823 GMT_Report (API, GMT_MSG_COMPAT, "Unit c for seconds is deprecated; use s instead.\n");
1824 Ctrl->S.unit = 's';
1825 }
1826 if (!strchr ("sm ", Ctrl->S.unit)) {
1827 GMT_Report (API, GMT_MSG_ERROR, "Option -S: Unrecognized unit %c\n", Ctrl->S.unit);
1828 n_errors++;
1829 }
1830 break;
1831 case 'T':
1832 n_errors += gmt_M_repeated_module_option (API, Ctrl->T.active);
1833 Ctrl->T.active = true;
1834 k = 0;
1835 if (gmt_M_compat_check (GMT, 4)) { /* GMT4 syntax allowed for upper case */
1836 modifier = opt->arg[strlen(opt->arg)-1];
1837 if (modifier == 'B') modifier = 'b';
1838 else if (modifier == 'I') modifier = 'i';
1839 if (!(modifier == 'b' || modifier == 'i'))
1840 modifier = opt->arg[0], k = 1;
1841 }
1842 else {
1843 modifier = opt->arg[0];
1844 k = 1;
1845 }
1846 if (modifier == 'b') {
1847 Ctrl->T.b_tension = atof (&opt->arg[k]);
1848 }
1849 else if (modifier == 'i') {
1850 Ctrl->T.i_tension = atof (&opt->arg[k]);
1851 }
1852 else if (modifier == '.' || (modifier >= '0' && modifier <= '9')) {
1853 Ctrl->T.i_tension = Ctrl->T.b_tension = atof (opt->arg);
1854 }
1855 else {
1856 GMT_Report (API, GMT_MSG_ERROR, "Option -T: Unrecognized modifier %c\n", modifier);
1857 n_errors++;
1858 }
1859 break;
1860 case 'W':
1861 n_errors += gmt_M_repeated_module_option (API, Ctrl->W.active);
1862 Ctrl->W.active = true;
1863 if (opt->arg[0]) { /* Specified named log file */
1864 gmt_M_str_free (Ctrl->W.file);
1865 Ctrl->W.file = strdup (opt->arg);
1866 }
1867 break;
1868 case 'Z':
1869 n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
1870 Ctrl->Z.active = true;
1871 Ctrl->Z.value = atof (opt->arg);
1872 break;
1873
1874 default: /* Report bad options */
1875 n_errors += gmt_default_error (GMT, opt->option);
1876 break;
1877 }
1878 }
1879
1880 n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
1881 n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0,
1882 "Option -I: Must specify positive increment(s)\n");
1883 n_errors += gmt_M_check_condition (GMT, Ctrl->N.value < 1, "Option -N: Max iterations must be nonzero\n");
1884 n_errors += gmt_M_check_condition (GMT, Ctrl->Z.value < 0.0 || Ctrl->Z.value > 2.0,
1885 "Option -Z: Relaxation value must be 1 <= z <= 2\n");
1886 n_errors += gmt_M_check_condition (GMT, !Ctrl->G.file && !Ctrl->Q.active, "Option -G: Must specify output grid file\n");
1887 n_errors += gmt_M_check_condition (GMT, Ctrl->A.mode && gmt_M_is_cartesian (GMT, GMT_IN),
1888 "Option -Am: Requires geographic input data\n");
1889 n_errors += gmt_check_binary_io (GMT, 3);
1890
1891 return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
1892 }
1893
1894 #define bailout(code) {gmt_M_free_options (mode); return (code);}
1895 #define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}
1896
GMT_surface(void * V_API,int mode,void * args)1897 EXTERN_MSC int GMT_surface (void *V_API, int mode, void *args) {
1898 int error = 0, key, one = 1, end;
1899 unsigned int old_verbose;
1900 char *limit[2] = {"lower", "upper"};
1901 double wesn[6];
1902
1903 struct SURFACE_INFO C;
1904 struct SURFACE_CTRL *Ctrl = NULL;
1905 struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
1906 struct GMT_OPTION *options = NULL;
1907 struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */
1908
1909 /*----------------------- Standard module initialization and parsing ----------------------*/
1910
1911 if (API == NULL) return (GMT_NOT_A_SESSION);
1912 if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */
1913 options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */
1914
1915 if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error); /* Give usage if requested */
1916
1917 /* Parse the command-line arguments */
1918
1919 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 */
1920 if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
1921 Ctrl = New_Ctrl (GMT); /* Allocate and initialize a new control structure */
1922 if ((error = parse (GMT, Ctrl, options)) != 0) Return (error);
1923
1924 /*---------------------------- This is the surface main code ----------------------------*/
1925
1926 gmt_M_tic(GMT);
1927 old_verbose = GMT->current.setting.verbose;
1928
1929 gmt_enable_threads (GMT); /* Set number of active threads, if supported */
1930 /* Some initializations and defaults setting */
1931 gmt_M_memset (&C, 1, struct SURFACE_INFO);
1932
1933 gmt_M_memcpy (C.wesn_orig, GMT->common.R.wesn, 4, double); /* Save original region in case user specified -r */
1934 gmt_M_memcpy (wesn, GMT->common.R.wesn, 4, double); /* Save specified region */
1935 C.periodic = (gmt_M_x_is_lon (GMT, GMT_IN) && gmt_M_360_range (wesn[XLO], wesn[XHI]));
1936 if (C.periodic && gmt_M_180_range (wesn[YLO], wesn[YHI])) {
1937 /* Trying to grid global geographic data - this is not something surface can do */
1938 GMT_Report (API, GMT_MSG_ERROR, "You are attempting to grid a global geographic data set, but surface cannot handle poles.\n");
1939 GMT_Report (API, GMT_MSG_ERROR, "It will do its best but it remains a Cartesian calculation which affects nodes near the poles.\n");
1940 GMT_Report (API, GMT_MSG_ERROR, "Because the grid is flaggged as geographic, the (repeated) pole values will be averaged upon writing to file.\n");
1941 GMT_Report (API, GMT_MSG_ERROR, "This may introduce a jump at either pole which will distort the grid near the poles.\n");
1942 GMT_Report (API, GMT_MSG_ERROR, "Consider spherical gridding instead with greenspline or sphinterpolate.\n");
1943 }
1944 if (GMT->common.R.registration == GMT_GRID_PIXEL_REG) { /* Pixel registration request. Use the trick of offsetting area by x_inc(y_inc) / 2 */
1945 /* Note that the grid remains node-registered and only gets tagged as pixel-registered upon writing the final grid to file */
1946 wesn[XLO] += GMT->common.R.inc[GMT_X] / 2.0; wesn[XHI] += GMT->common.R.inc[GMT_X] / 2.0;
1947 wesn[YLO] += GMT->common.R.inc[GMT_Y] / 2.0; wesn[YHI] += GMT->common.R.inc[GMT_Y] / 2.0;
1948 /* n_columns,n_rows remain the same for now but nodes are in "pixel" position. We reset to original wesn and reduce n_columns,n_rows by 1 when we write result */
1949 }
1950 if (Ctrl->A.mode) Ctrl->A.value = cosd (0.5 * (wesn[YLO] + wesn[YHI])); /* Set cos of middle latitude as aspect ratio */
1951
1952 if ((C.Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_ONLY, NULL, wesn, NULL,
1953 GMT_GRID_NODE_REG, GMT_NOTSET, NULL)) == NULL) Return (API->error);
1954
1955 if (C.Grid->header->n_columns < 4 || C.Grid->header->n_rows < 4) {
1956 GMT_Report (API, GMT_MSG_ERROR, "Grid must have at least 4 nodes in each direction (you have %d by %d) - abort.\n", C.Grid->header->n_columns, C.Grid->header->n_rows);
1957 Return (GMT_RUNTIME_ERROR);
1958 }
1959
1960 surface_init_parameters (&C, Ctrl); /* Pass parameters from parsing control to surface information structure C */
1961
1962 /* Determine the initial and intermediate grid dimensions */
1963 C.current_stride = gmt_gcd_euclid (C.n_columns-1, C.n_rows-1);
1964
1965 if (Ctrl->Q.active && old_verbose < GMT_MSG_INFORMATION) /* Temporarily escalate verbosity to INFORMATION */
1966 GMT->current.setting.verbose = GMT_MSG_INFORMATION;
1967 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION) || Ctrl->Q.active) {
1968 sprintf (C.format, "Grid domain: W: %s E: %s S: %s N: %s n_columns: %%d n_rows: %%d [", GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out, GMT->current.setting.format_float_out);
1969 (GMT->common.R.registration == GMT_GRID_PIXEL_REG) ? strcat (C.format, "pixel registration]\n") : strcat (C.format, "gridline registration]\n");
1970 GMT_Report (API, GMT_MSG_INFORMATION, C.format, C.wesn_orig[XLO], C.wesn_orig[XHI], C.wesn_orig[YLO], C.wesn_orig[YHI], C.n_columns-one, C.n_rows-one);
1971 }
1972 if (C.current_stride == 1) GMT_Report (API, GMT_MSG_WARNING, "Your grid dimensions are mutually prime. Convergence is very unlikely.\n");
1973 if ((C.current_stride == 1 && gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) || Ctrl->Q.active) surface_suggest_sizes (GMT, C.Grid, C.factors, C.n_columns-1, C.n_rows-1, GMT->common.R.registration == GMT_GRID_PIXEL_REG);
1974 if (Ctrl->Q.active) { /* Reset verbosity and bail */
1975 GMT->current.setting.verbose = old_verbose;
1976 Return (GMT_NOERROR);
1977 }
1978
1979 /* Set current_stride = 1, read data, setting indices. Then throw
1980 away data that can't be used in the end game, limiting the
1981 size of data arrays and Briggs->b[6] structure/array. */
1982
1983 C.current_stride = 1;
1984 surface_set_grid_parameters (&C);
1985 if (surface_read_data (GMT, &C, options))
1986 Return (GMT_RUNTIME_ERROR);
1987 if (Ctrl->D.active) { /* Append breakline dataset */
1988 struct GMT_DATASET *Lin = NULL;
1989 char *file = (Ctrl->D.debug) ? Ctrl->D.file : NULL;
1990 if (Ctrl->D.fix_z) { /* Either provide a fixed z value or override whatever input file may supply with this value */
1991 //GMT->common.b.ncol[GMT_IN] = 0; /* So Set_Columns will work */
1992 if ((error = GMT_Set_Columns (GMT->parent, GMT_IN, 2, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) /* Only read 2 columns */
1993 Return (GMT_RUNTIME_ERROR);
1994 }
1995 if ((Lin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, GMT_READ_NORMAL, NULL, Ctrl->D.file, NULL)) == NULL)
1996 Return (API->error);
1997 if (Lin->n_columns < 2) {
1998 GMT_Report (API, GMT_MSG_ERROR, "Input file %s has %d column(s) but at least 2 are needed\n", Ctrl->D.file, (int)Lin->n_columns);
1999 Return (GMT_DIM_TOO_SMALL);
2000 }
2001 surface_interpolate_add_breakline (GMT, &C, Lin->table[0], file, Ctrl->D.fix_z, Ctrl->D.z); /* Pass the single table since we read a single file */
2002 }
2003
2004 surface_throw_away_unusables (GMT, &C); /* Eliminate data points that will not serve as constraints */
2005 surface_remove_planar_trend (GMT, &C); /* Fit best-fitting plane and remove it from the data; plane will be restored at the end */
2006 key = surface_rescale_z_values (GMT, &C); /* Divide residual data by their rms value */
2007
2008 if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, C.Grid) != GMT_NOERROR) Return (API->error);
2009 if (key == 1) { /* Data lie exactly on a plane; write a grid with the plan and exit */
2010 gmt_M_free (GMT, C.data); /* The data set is no longer needed */
2011 if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL,
2012 0, 0, C.Grid) == NULL) Return (API->error); /* Get a grid of zeros... */
2013 surface_restore_planar_trend (&C); /* ...and restore the plane we found */
2014 if ((error = surface_write_grid (GMT, &C, Ctrl->G.file)) != 0) /* Write this grid */
2015 Return (error);
2016 Return (GMT_NOERROR); /* Clean up and return */
2017 }
2018
2019 surface_load_constraints (GMT, &C, true); /* Set lower and upper constraint grids, if requested */
2020
2021 /* Set up factors and reset current_stride to its initial (and largest) value */
2022
2023 C.current_stride = gmt_gcd_euclid (C.n_columns-1, C.n_rows-1);
2024 C.n_factors = gmt_get_prime_factors (GMT, C.current_stride, C.factors);
2025 surface_set_grid_parameters (&C);
2026 while (C.current_nx < 4 || C.current_ny < 4) { /* Must have at least a grid of 4x4 */
2027 surface_smart_divide (&C);
2028 surface_set_grid_parameters (&C);
2029 }
2030 surface_set_offset (&C); /* Initialize the node-jumps across rows for this grid size */
2031 surface_set_index (GMT, &C); /* Determine the nearest data constraint for this grid size */
2032
2033 if (Ctrl->W.active) { /* Want to log convergence information to file */
2034 if ((C.fp_log = gmt_fopen (GMT, Ctrl->W.file, "w")) == NULL) {
2035 GMT_Report (API, GMT_MSG_ERROR, "Unable to create log file %s.\n", Ctrl->W.file);
2036 Return (GMT_ERROR_ON_FOPEN);
2037 }
2038 C.logging = true;
2039 fprintf (C.fp_log, "#grid\tmode\tgrid_iteration\tchange\tlimit\ttotal_iteration\n");
2040 }
2041
2042 /* Now the data are ready to go for the first iteration. */
2043
2044 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) { /* Report on memory usage for this run */
2045 size_t mem_use, mem_total;
2046 mem_use = mem_total = C.npoints * sizeof (struct SURFACE_DATA);
2047 GMT_Report (API, GMT_MSG_INFORMATION, "------------------------------------------\n");
2048 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Memory for data array", gmt_memory_use (mem_use, 1));
2049 mem_use = sizeof (struct GMT_GRID) + C.mxmy * sizeof (gmt_grdfloat); mem_total += mem_use;
2050 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Memory for final grid", gmt_memory_use (mem_use, 1));
2051 for (end = LO; end <= HI; end++) if (C.set_limit[end]) { /* Will need to keep a lower|upper surface constrain grid */
2052 mem_total += mem_use;
2053 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Memory for constraint grid", gmt_memory_use (mem_use, 1));
2054 }
2055 mem_use = C.npoints * sizeof (struct SURFACE_BRIGGS) ; mem_total += mem_use;
2056 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Memory for Briggs coefficients", gmt_memory_use (mem_use, 1));
2057 mem_use = C.mxmy; mem_total += mem_use;
2058 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Memory for node status", gmt_memory_use (mem_use, 1));
2059 GMT_Report (API, GMT_MSG_INFORMATION, "------------------------------------------\n");
2060 GMT_Report (API, GMT_MSG_INFORMATION, "%-31s: %9s\n", "Total memory use", gmt_memory_use (mem_total, 1));
2061 GMT_Report (API, GMT_MSG_INFORMATION, "==========================================\n");
2062 }
2063
2064 /* Allocate the memory needed to perform the gridding */
2065
2066 C.Briggs = gmt_M_memory (GMT, NULL, C.npoints, struct SURFACE_BRIGGS);
2067 C.status = gmt_M_memory (GMT, NULL, C.mxmy, char);
2068 C.fraction = gmt_M_memory (GMT, NULL, C.current_stride, double);
2069 if (GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_DATA_ONLY, NULL, NULL, NULL, 0, 0, C.Grid) == NULL)
2070 Return (API->error);
2071 if (C.radius > 0) surface_initialize_grid (GMT, &C); /* Fill in nodes with a weighted average in a search radius */
2072 GMT_Report (API, GMT_MSG_INFORMATION, "Grid\tMode\tIteration\tMax Change\tConv Limit\tTotal Iterations\n");
2073
2074 surface_set_coefficients (GMT, &C); /* Initialize the coefficients needed in the finite-difference expressions */
2075
2076 /* Here is the main multigrid loop, were we first grid using a coarse grid and the
2077 * progressively refine the grid until we reach the final configuration. */
2078
2079 C.previous_stride = C.current_stride;
2080 surface_find_nearest_constraint (GMT, &C); /* Assign nearest data value to nodes and evaluate Briggs coefficients */
2081 surface_iterate (GMT, &C, GRID_DATA); /* Grid the data using the data constraints */
2082
2083 while (C.current_stride > 1) { /* More intermediate grids remain, go to next */
2084 surface_smart_divide (&C); /* Set the new current_stride */
2085 surface_set_grid_parameters (&C); /* Update node book-keeping constants */
2086 surface_set_offset (&C); /* Reset the node-jumps across rows for this grid size */
2087 surface_set_index (GMT, &C); /* Recompute the index values for the nearest data points */
2088 fill_in_forecast (GMT, &C); /* Expand the grid and fill it via bilinear interpolation */
2089 surface_iterate (GMT, &C, GRID_NODES); /* Grid again but only to improve on the bilinear guesses */
2090 surface_find_nearest_constraint (GMT, &C); /* Assign nearest data value to nodes and evaluate Briggs coefficients */
2091 surface_iterate (GMT, &C, GRID_DATA); /* Grid the data but now use the data constraints */
2092 C.previous_stride = C.current_stride; /* Remember previous stride before we smart-divide again */
2093 }
2094
2095 if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) surface_check_errors (GMT, &C); /* Report on mean misfit and curvature */
2096
2097 surface_restore_planar_trend (&C); /* Restore the least-square plane we removed earlier */
2098
2099 if (Ctrl->W.active) /* Close the log file */
2100 gmt_fclose (GMT, C.fp_log);
2101
2102 /* Clean up after ourselves */
2103
2104 gmt_M_free (GMT, C.Briggs);
2105 gmt_M_free (GMT, C.status);
2106 gmt_M_free (GMT, C.fraction);
2107 for (end = LO; end <= HI; end++) if (C.set_limit[end]) { /* Free lower|upper surface constrain grids */
2108 if (GMT_Destroy_Data (API, &C.Bound[end]) != GMT_NOERROR)
2109 GMT_Report (API, GMT_MSG_ERROR, "Failed to free grid with %s bounds\n", limit[end]);
2110 }
2111
2112 if (Ctrl->M.active) { /* Want to mask the grid first */
2113 char input[GMT_VF_LEN] = {""}, mask[GMT_VF_LEN] = {""}, cmd[GMT_LEN256] = {""};
2114 static char *V_level = GMT_VERBOSE_CODES;
2115 struct GMT_GRID *Gmask = NULL;
2116 struct GMT_VECTOR *V = NULL;
2117 double *data[2] = {NULL, NULL};
2118 uint64_t row, col, ij, dim[3] = {2, C.npoints, GMT_DOUBLE}; /* ncols, nrows, type */
2119
2120 if ((V = GMT_Create_Data (API, GMT_IS_VECTOR, GMT_IS_POINT, GMT_CONTAINER_ONLY, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
2121 Return (API->error);
2122 }
2123 for (col = 0; col < 2; col++)
2124 data[col] = gmt_M_memory (GMT, NULL, C.npoints, double);
2125 for (ij = 0; ij < C.npoints; ij++) {
2126 data[GMT_X][ij] = C.data[ij].x;
2127 data[GMT_Y][ij] = C.data[ij].y;
2128 }
2129 gmt_M_free (GMT, C.data);
2130 GMT_Put_Vector (API, V, GMT_X, GMT_DOUBLE, data[GMT_X]);
2131 GMT_Put_Vector (API, V, GMT_Y, GMT_DOUBLE, data[GMT_Y]);
2132 /* Create a virtual file for reading the input data grid */
2133 if (GMT_Open_VirtualFile (API, GMT_IS_DATASET|GMT_VIA_VECTOR, GMT_IS_POINT, GMT_IN|GMT_IS_REFERENCE, V, input) == GMT_NOTSET) {
2134 Return (API->error);
2135 }
2136 /* Create a virtual file to hold the mask grid */
2137 if (GMT_Open_VirtualFile (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_OUT|GMT_IS_REFERENCE, NULL, mask) == GMT_NOTSET) {
2138 Return (API->error);
2139 }
2140 gmt_disable_bghio_opts (GMT); /* Do not want any -b -g -h -i -o to affect the reading input file in grdmask */
2141 /* Hardwire -rg since internally, all grids are gridline registered (until output at least) */
2142 sprintf (cmd, "%s -G%s -R%g/%g/%g/%g -I%g/%g -NNaN/1/1 -S%s -V%c -rg --GMT_HISTORY=readonly",
2143 input, mask, wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI], GMT->common.R.inc[GMT_X],
2144 GMT->common.R.inc[GMT_Y], Ctrl->M.arg, V_level[GMT->current.setting.verbose]);
2145 GMT_Report (API, GMT_MSG_INFORMATION, "Masking grid nodes away from data points via grdmask\n");
2146 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Calling grdsample with args %s\n", cmd);
2147 if (GMT_Call_Module (API, "grdmask", GMT_MODULE_CMD, cmd) != GMT_NOERROR) { /* Resample the file */
2148 GMT_Report (API, GMT_MSG_ERROR, "Unable to mask the intermediate grid - exiting\n");
2149 Return (API->error);
2150 }
2151 if (GMT_Close_VirtualFile (API, input) == GMT_NOTSET) {
2152 Return (API->error);
2153 }
2154 if (GMT_Destroy_Data (API, &V) != GMT_NOERROR) { /* Done with the data set */
2155 Return (API->error);
2156 }
2157 gmt_M_free (GMT, data[GMT_X]); gmt_M_free (GMT, data[GMT_Y]);
2158 if ((Gmask = GMT_Read_VirtualFile (API, mask)) == NULL) { /* Load in the mask grid */
2159 Return (API->error);
2160 }
2161 /* Apply the mask */
2162 gmt_M_grd_loop (GMT, Gmask, row, col, ij) C.Grid->data[ij] *= Gmask->data[ij];
2163 if (GMT_Destroy_Data (API, &Gmask) != GMT_NOERROR) { /* Done with the mask */
2164 Return (API->error);
2165 }
2166 gmt_reenable_bghio_opts (GMT); /* Recover settings provided by user (if -b -g -h -i were used at all) */
2167 }
2168 else
2169 gmt_M_free (GMT, C.data);
2170
2171 /* Only gave -J to set the proj4 flag. Do it here only because otherwise would interfere with computations(!!) */
2172 if (Ctrl->J.active) {
2173 gmt_parse_common_options (GMT, "J", 'J', Ctrl->J.projstring); /* Has to be processed independently of -R */
2174 C.Grid->header->ProjRefPROJ4 = gmt_export2proj4 (GMT); /* Convert the GMT -J<...> into a proj4 string */;
2175 free (Ctrl->J.projstring);
2176 }
2177
2178 if ((error = surface_write_grid (GMT, &C, Ctrl->G.file)) != 0) /* Write the output grid */
2179 Return (error);
2180
2181 gmt_M_toc(GMT,""); /* Print total run time, but only if -Vt was set */
2182
2183 Return (GMT_NOERROR);
2184 }
2185