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 *
19 * G M T _ C D F . C R O U T I N E S
20 *
21 * Takes care of all grd input/output built on NCAR's netCDF routines (which is
22 * an XDR implementation)
23 * Most functions will return with error message if an internal error is returned.
24 * There functions are only called indirectly via the GMT_* grdio functions.
25 *
26 * Author: Paul Wessel
27 * Date: 1-JAN-2010
28 * Version: 5
29 *
30 * Public functions (5):
31 *
32 * gmt_cdf_read_grd_info : Read header from file
33 * gmt_cdf_read_grd : Read header and data set from file
34 * gmt_cdf_update_grd_info : Update header in existing file
35 * gmt_cdf_write_grd_info : Write header to new file
36 * gmt_cdf_write_grd : Write header and data set to new file
37 *
38 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
39
40 #include "gmt_dev.h"
41 #include "gmt_internals.h"
42
gmt_cdf_grd_info(struct GMT_CTRL * GMT,int ncid,struct GMT_GRID_HEADER * header,char job)43 int gmt_cdf_grd_info (struct GMT_CTRL *GMT, int ncid, struct GMT_GRID_HEADER *header, char job) {
44 int err; /* Implicitly by gmt_M_err_trap */
45 int nm[2];
46 double dummy[2];
47 char text[GMT_GRID_COMMAND_LEN320+GMT_GRID_REMARK_LEN160];
48 size_t limit = 2147483647U; /* 2^31 - 1 is the max length of a 1-D array in netCDF */
49 nc_type z_type;
50 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
51 /* Dimension ids, variable ids, etc. */
52 int side_dim, xysize_dim, x_range_id, y_range_id, z_range_id, inc_id, nm_id, z_id, dims[1];
53
54 /* Define and get dimensions and variables */
55
56 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Enter gmt_cdf_grd_info with argument %c\n", (int)job);
57
58 if (job == 'w') {
59 if (header->nm > limit) { /* Print error message and let things crash */
60 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Your grid contains more than 2^31 - 1 nodes (%" PRIu64 ") and cannot be stored with the deprecated GMT netCDF format.\n", header->nm);
61 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Please choose another grid format such as the default netCDF 4 COARDS-compliant grid format.\n");
62 return (GMT_DIM_TOO_LARGE);
63 }
64 gmt_M_err_trap (nc_def_dim (ncid, "side", 2U, &side_dim));
65 gmt_M_err_trap (nc_def_dim (ncid, "xysize", header->nm, &xysize_dim));
66
67 dims[0] = side_dim;
68 gmt_M_err_trap (nc_def_var (ncid, "x_range", NC_DOUBLE, 1, dims, &x_range_id));
69 gmt_M_err_trap (nc_def_var (ncid, "y_range", NC_DOUBLE, 1, dims, &y_range_id));
70 gmt_M_err_trap (nc_def_var (ncid, "z_range", NC_DOUBLE, 1, dims, &z_range_id));
71 gmt_M_err_trap (nc_def_var (ncid, "spacing", NC_DOUBLE, 1, dims, &inc_id));
72 gmt_M_err_trap (nc_def_var (ncid, "dimension", NC_LONG, 1, dims, &nm_id));
73
74 switch (header->type) {
75 case GMT_GRID_IS_CB: z_type = NC_BYTE; break;
76 case GMT_GRID_IS_CS: z_type = NC_SHORT; break;
77 case GMT_GRID_IS_CI: z_type = NC_INT; break;
78 case GMT_GRID_IS_CF: z_type = NC_FLOAT; break;
79 case GMT_GRID_IS_CD: z_type = NC_DOUBLE; break;
80 default: z_type = NC_NAT;
81 }
82
83 dims[0] = xysize_dim;
84 gmt_M_err_trap (nc_def_var (ncid, "z", z_type, 1, dims, &z_id));
85 }
86 else {
87 gmt_M_err_trap (nc_inq_varid (ncid, "x_range", &x_range_id));
88 gmt_M_err_trap (nc_inq_varid (ncid, "y_range", &y_range_id));
89 gmt_M_err_trap (nc_inq_varid (ncid, "z_range", &z_range_id));
90 gmt_M_err_trap (nc_inq_varid (ncid, "spacing", &inc_id));
91 gmt_M_err_trap (nc_inq_varid (ncid, "dimension", &nm_id));
92 gmt_M_err_trap (nc_inq_varid (ncid, "z", &z_id));
93 gmt_M_err_trap (nc_inq_vartype (ncid, z_id, &z_type));
94 switch (z_type) {
95 case NC_BYTE: header->type = GMT_GRID_IS_CB; HH->orig_datatype = GMT_CHAR; break;
96 case NC_SHORT: header->type = GMT_GRID_IS_CS; HH->orig_datatype = GMT_SHORT; break;
97 case NC_INT: header->type = GMT_GRID_IS_CI; HH->orig_datatype = GMT_INT; break;
98 case NC_FLOAT: header->type = GMT_GRID_IS_CF; HH->orig_datatype = GMT_FLOAT; break;
99 case NC_DOUBLE: header->type = GMT_GRID_IS_CD; HH->orig_datatype = GMT_DOUBLE; break;
100 default: header->type = k_grd_unknown_fmt; break;
101 }
102 }
103 HH->z_id = z_id;
104
105 /* Get or assign attributes */
106
107 gmt_M_memset (text, GMT_GRID_COMMAND_LEN320+GMT_GRID_REMARK_LEN160, char);
108
109 if (job == 'u') gmt_M_err_trap (nc_redef (ncid));
110
111 if (job == 'r') {
112 int reg;
113 gmt_M_err_trap (nc_get_att_text (ncid, x_range_id, "units", header->x_units));
114 gmt_M_err_trap (nc_get_att_text (ncid, y_range_id, "units", header->y_units));
115 gmt_M_err_trap (nc_get_att_text (ncid, z_range_id, "units", header->z_units));
116 gmt_M_err_trap (nc_get_att_double (ncid, z_id, "scale_factor", &header->z_scale_factor));
117 gmt_M_err_trap (nc_get_att_double (ncid, z_id, "add_offset", &header->z_add_offset));
118 gmt_M_err_trap (nc_get_att_int (ncid, z_id, "node_offset", ®));
119 header->registration = reg;
120 #ifdef DOUBLE_PRECISION_GRID
121 nc_get_att_double (ncid, z_id, "_FillValue", &header->nan_value);
122 #else
123 nc_get_att_float (ncid, z_id, "_FillValue", &header->nan_value);
124 #endif
125 gmt_M_err_trap (nc_get_att_text (ncid, NC_GLOBAL, "title", header->title));
126 gmt_M_err_trap (nc_get_att_text (ncid, NC_GLOBAL, "source", text));
127 strncpy (header->command, text, GMT_GRID_COMMAND_LEN320-1);
128 strncpy (header->remark, &text[GMT_GRID_COMMAND_LEN320], GMT_GRID_REMARK_LEN160-1);
129
130 gmt_M_err_trap (nc_get_var_double (ncid, x_range_id, dummy));
131 header->wesn[XLO] = dummy[0];
132 header->wesn[XHI] = dummy[1];
133 gmt_M_err_trap (nc_get_var_double (ncid, y_range_id, dummy));
134 header->wesn[YLO] = dummy[0];
135 header->wesn[YHI] = dummy[1];
136 gmt_M_err_trap (nc_get_var_double (ncid, inc_id, dummy));
137 header->inc[GMT_X] = dummy[0];
138 header->inc[GMT_Y] = dummy[1];
139 gmt_M_err_trap (nc_get_var_int (ncid, nm_id, nm));
140 header->n_columns = nm[0];
141 header->n_rows = nm[1];
142 gmt_M_err_trap (nc_get_var_double (ncid, z_range_id, dummy));
143 header->z_min = dummy[0];
144 header->z_max = dummy[1];
145 //HH->row_order = k_nc_start_north; /* N->S but that is just when read. We don't want to set this for writing later */
146 }
147 else {
148 int reg;
149 strncpy (text, header->command, GMT_GRID_COMMAND_LEN320-1);
150 strncpy (&text[GMT_GRID_COMMAND_LEN320], header->remark, GMT_GRID_REMARK_LEN160-1);
151 gmt_M_err_trap (nc_put_att_text (ncid, x_range_id, "units", GMT_GRID_UNIT_LEN80, header->x_units));
152 gmt_M_err_trap (nc_put_att_text (ncid, y_range_id, "units", GMT_GRID_UNIT_LEN80, header->y_units));
153 gmt_M_err_trap (nc_put_att_text (ncid, z_range_id, "units", GMT_GRID_UNIT_LEN80, header->z_units));
154 gmt_M_err_trap (nc_put_att_double (ncid, z_id, "scale_factor", NC_DOUBLE, 1U, &header->z_scale_factor));
155 gmt_M_err_trap (nc_put_att_double (ncid, z_id, "add_offset", NC_DOUBLE, 1U, &header->z_add_offset));
156 if (z_type != NC_FLOAT && z_type != NC_DOUBLE)
157 header->nan_value = rintf (header->nan_value); /* round to integer */
158 #ifdef DOUBLE_PRECISION_GRID
159 gmt_M_err_trap (nc_put_att_double (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value));
160 #else
161 gmt_M_err_trap (nc_put_att_float (ncid, z_id, "_FillValue", z_type, 1U, &header->nan_value));
162 #endif
163 reg = header->registration;
164 gmt_M_err_trap (nc_put_att_int (ncid, z_id, "node_offset", NC_LONG, 1U, ®));
165 gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "title", GMT_GRID_TITLE_LEN80, header->title));
166 gmt_M_err_trap (nc_put_att_text (ncid, NC_GLOBAL, "source", GMT_GRID_COMMAND_LEN320+GMT_GRID_REMARK_LEN160, text));
167
168 gmt_M_err_trap (nc_enddef (ncid));
169
170 dummy[0] = header->wesn[XLO]; dummy[1] = header->wesn[XHI];
171 gmt_M_err_trap (nc_put_var_double (ncid, x_range_id, dummy));
172 dummy[0] = header->wesn[YLO]; dummy[1] = header->wesn[YHI];
173 gmt_M_err_trap (nc_put_var_double (ncid, y_range_id, dummy));
174 dummy[0] = header->inc[GMT_X]; dummy[1] = header->inc[GMT_Y];
175 gmt_M_err_trap (nc_put_var_double (ncid, inc_id, dummy));
176 nm[0] = header->n_columns; nm[1] = header->n_rows;
177 gmt_M_err_trap (nc_put_var_int (ncid, nm_id, nm));
178 if (header->z_min <= header->z_max) {
179 dummy[0] = header->z_min; dummy[1] = header->z_max;
180 }
181 else {
182 dummy[0] = 0.0; dummy[1] = 0.0;
183 }
184 gmt_M_err_trap (nc_put_var_double (ncid, z_range_id, dummy));
185 }
186 return (GMT_NOERROR);
187 }
188
gmt_cdf_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)189 int gmt_cdf_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
190 int ncid, err;
191 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
192 if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE);
193 gmt_M_err_trap (gmt_nc_open (GMT, HH->name, NC_NOWRITE, &ncid));
194 gmt_M_err_trap (gmt_cdf_grd_info (GMT, ncid, header, 'r'));
195 gmt_M_err_trap (gmt_nc_close (GMT, ncid));
196 return (GMT_NOERROR);
197 }
198
gmt_cdf_update_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)199 int gmt_cdf_update_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
200 int ncid, old_fill_mode, err;
201 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
202 if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE);
203 gmt_M_err_trap (gmt_nc_open (GMT, HH->name, NC_WRITE, &ncid));
204 gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode));
205 gmt_M_err_trap (gmt_cdf_grd_info (GMT, ncid, header, 'u'));
206 gmt_M_err_trap (gmt_nc_close (GMT, ncid));
207 return (GMT_NOERROR);
208 }
209
gmt_cdf_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)210 int gmt_cdf_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
211 int ncid, old_fill_mode, err;
212 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
213 if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE);
214 gmt_M_err_trap (gmt_nc_create (GMT, HH->name, NC_CLOBBER, &ncid));
215 gmt_M_err_trap (nc_set_fill (ncid, NC_NOFILL, &old_fill_mode));
216 gmt_M_err_trap (gmt_cdf_grd_info (GMT, ncid, header, 'w'));
217 gmt_M_err_trap (gmt_nc_close (GMT, ncid));
218 return (GMT_NOERROR);
219 }
220
gmt_cdf_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)221 int gmt_cdf_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
222 /* header: grid structure header
223 * grid: array with final grid
224 * wesn: Sub-region to extract [Use entire file if 0,0,0,0]
225 * padding: # of empty rows/columns to add on w, e, s, n of grid, respectively
226 * complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only)
227 * Note: The file has only real values, we simply allow space in the complex array
228 * for real and imaginary parts when processed by grdfft etc.
229 *
230 * Reads a subset of a grid file and optionally pads the array with extra rows and columns
231 * header values for n_columns and n_rows are reset to reflect the dimensions of the logical array,
232 * not the physical size (i.e., the padding is not counted in n_columns and n_rows)
233 */
234
235 bool check;
236 int ncid, j, err, first_col, last_col, first_row, last_row;
237 unsigned int i, width_in, height_in;
238 unsigned int width_out, *actual_col = NULL;
239 uint64_t ij, kk, imag_offset;
240 size_t start[1], edge[1];
241 gmt_grdfloat *tmp = NULL;
242 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
243
244 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_in, &height_in, &first_col, &last_col, &first_row, &last_row, &actual_col), HH->name);
245 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
246
247 width_out = width_in; /* Width of output array */
248 if (pad[XLO] > 0) width_out += pad[XLO];
249 if (pad[XHI] > 0) width_out += pad[XHI];
250
251 /* Open the NetCDF file */
252
253 if (!strcmp (HH->name,"=")) {
254 gmt_M_free (GMT, actual_col);
255 return (GMT_GRDIO_NC_NO_PIPE);
256 }
257 if ((err = gmt_nc_open (GMT, HH->name, NC_NOWRITE, &ncid)) != 0) {
258 gmt_M_free (GMT, actual_col);
259 return err;
260 }
261 check = !isnan (header->nan_value);
262
263 /* Load data row by row. The data in the file is stored in the same
264 * "upside down" fashion as within GMT. The first row is the top row */
265
266 tmp = gmt_M_memory (GMT, NULL, header->n_columns, gmt_grdfloat);
267
268 edge[0] = header->n_columns;
269 ij = imag_offset + pad[YHI] * width_out + pad[XLO];
270 header->z_min = DBL_MAX;
271 header->z_max = -DBL_MAX;
272 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
273
274 for (j = first_row; j <= last_row; j++, ij += width_out) {
275 start[0] = j * header->n_columns;
276 if ((err = gmt_nc_get_vara_grdfloat (ncid, HH->z_id, start, edge, tmp))) { /* Get one row */
277 gmt_M_free (GMT, actual_col);
278 gmt_M_free (GMT, tmp);
279 gmt_nc_close (GMT, ncid);
280 return (err);
281 }
282 for (i = 0; i < width_in; i++) { /* Check for and handle NaN proxies */
283 kk = ij+i;
284 grid[kk] = tmp[actual_col[i]];
285 if (check && grid[kk] == header->nan_value)
286 grid[kk] = GMT->session.f_NaN;
287 if (gmt_M_is_fnan (grid[kk])) {
288 HH->has_NaNs = GMT_GRID_HAS_NANS;
289 continue;
290 }
291 header->z_min = MIN (header->z_min, (double)grid[kk]);
292 header->z_max = MAX (header->z_max, (double)grid[kk]);
293 }
294 }
295
296 header->n_columns = width_in;
297 header->n_rows = height_in;
298 gmt_M_memcpy (header->wesn, wesn, 4, double);
299
300 gmt_M_free (GMT, actual_col);
301 gmt_M_free (GMT, tmp);
302 gmt_M_err_trap (gmt_nc_close (GMT, ncid));
303
304 return (GMT_NOERROR);
305 }
306
gmt_cdf_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)307 int gmt_cdf_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
308 /* header: grid structure header
309 * grid: array with final grid
310 * wesn: Sub-region to write out [Use entire file if 0,0,0,0]
311 * padding: # of empty rows/columns to add on w, e, s, n of grid, respectively
312 * complex_mode: &1 | &2 if complex array is to hold real (1) and imaginary (2) parts (otherwise read as real only)
313 * Note: The file has only real values, we simply allow space in the complex array
314 * for real and imaginary parts when processed by grdfft etc.
315 */
316
317 size_t start[1], edge[1];
318 int ncid, old_fill_mode, err, first_col, last_col, first_row, last_row;
319 long *tmp_i = NULL;
320 unsigned int i, *actual_col = NULL;
321 unsigned int j, width_out, height_out, width_in;
322 uint64_t ij, nr_oor = 0, imag_offset;
323 gmt_grdfloat *tmp_f = NULL;
324 double limit[2] = {-FLT_MAX, FLT_MAX};
325 gmt_grdfloat value;
326 nc_type z_type;
327 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
328
329 if (!strcmp (HH->name,"=")) return (GMT_GRDIO_NC_NO_PIPE); /* Cannot do piping on netCDF files */
330
331 /* Determine the value to be assigned to missing data, if not already done so */
332
333 switch (header->type) {
334 case GMT_GRID_IS_CB:
335 if (isnan (header->nan_value)) header->nan_value = CHAR_MIN;
336 limit[0] = CHAR_MIN - 0.5; limit[1] = CHAR_MAX + 0.5;
337 z_type = NC_BYTE; break;
338 case GMT_GRID_IS_CS:
339 if (isnan (header->nan_value)) header->nan_value = SHRT_MIN;
340 limit[0] = SHRT_MIN - 0.5; limit[1] = SHRT_MAX + 0.5;
341 z_type = NC_SHORT; break;
342 case GMT_GRID_IS_CI:
343 if (isnan (header->nan_value)) header->nan_value = INT_MIN;
344 limit[0] = INT_MIN - 0.5; limit[1] = INT_MAX + 0.5;
345 z_type = NC_INT; break;
346 case GMT_GRID_IS_CF:
347 z_type = NC_FLOAT; break;
348 case GMT_GRID_IS_CD:
349 z_type = NC_DOUBLE; break;
350 default:
351 z_type = NC_NAT;
352 }
353
354 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_out, &height_out, &first_col, &last_col, &first_row, &last_row, &actual_col), HH->name);
355 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
356
357 width_in = width_out; /* Physical width of input array */
358 if (pad[XLO] > 0) width_in += pad[XLO];
359 if (pad[XHI] > 0) width_in += pad[XHI];
360
361 gmt_M_memcpy (header->wesn, wesn, 4, double);
362 header->n_columns = width_out;
363 header->n_rows = height_out;
364
365 /* Write grid header */
366
367 if ((err = gmt_nc_create (GMT, HH->name, NC_CLOBBER, &ncid))) {
368 gmt_M_free (GMT, actual_col);
369 return (err);
370 }
371 if ((err = nc_set_fill (ncid, NC_NOFILL, &old_fill_mode))) {
372 gmt_M_free (GMT, actual_col);
373 return (err);
374 }
375 if ((err = gmt_cdf_grd_info (GMT, ncid, header, 'w'))) {
376 gmt_M_free (GMT, actual_col);
377 return (err);
378 }
379
380 /* Set start position for writing grid */
381
382 edge[0] = width_out;
383 ij = first_col + pad[XLO] + (uint64_t)(first_row + pad[YHI]) * width_in;
384 header->z_min = DBL_MAX;
385 header->z_max = -DBL_MAX;
386
387 /* Store z-variable */
388
389 if (z_type == NC_FLOAT || z_type == NC_DOUBLE) {
390 tmp_f = gmt_M_memory (GMT, NULL, width_in, gmt_grdfloat);
391 for (j = 0; j < height_out; j++, ij += width_in) {
392 start[0] = j * width_out;
393 for (i = 0; i < width_out; i++) {
394 value = grid[ij+actual_col[i]+imag_offset];
395 if (!isfinite (value)) {
396 if (isinf(value))
397 nr_oor++; /* out of gmt_grdfloat range */
398 tmp_f[i] = header->nan_value;
399 }
400 else {
401 tmp_f[i] = value;
402 header->z_min = MIN (header->z_min, (double)tmp_f[i]);
403 header->z_max = MAX (header->z_max, (double)tmp_f[i]);
404 }
405 }
406 if ((err = gmt_nc_put_vara_grdfloat (ncid, HH->z_id, start, edge, tmp_f))) {
407 gmt_M_free (GMT, actual_col);
408 gmt_M_free (GMT, tmp_f);
409 return (err);
410 }
411 }
412 gmt_M_free (GMT, tmp_f);
413 }
414 else {
415 tmp_i = gmt_M_memory (GMT, NULL, width_in, long);
416 for (j = 0; j < height_out; j++, ij += width_in) {
417 start[0] = j * width_out;
418 for (i = 0; i < width_out; i++) {
419 value = grid[ij+actual_col[i]+imag_offset];
420 if (!isfinite (value)) {
421 if (isinf(value))
422 nr_oor++; /* out of gmt_grdfloat range */
423 tmp_i[i] = lrintf (header->nan_value);
424 }
425 else if (value <= limit[0] || value >= limit[1]) {
426 tmp_i[i] = lrintf (header->nan_value);
427 nr_oor++;
428 }
429 else {
430 tmp_i[i] = lrintf (value);
431 header->z_min = MIN (header->z_min, (double)tmp_i[i]);
432 header->z_max = MAX (header->z_max, (double)tmp_i[i]);
433 }
434 }
435 if ((err = nc_put_vara_long (ncid, HH->z_id, start, edge, tmp_i))) {
436 gmt_M_free (GMT, actual_col);
437 gmt_M_free (GMT, tmp_i);
438 return (err);
439 }
440 }
441 gmt_M_free (GMT, tmp_i);
442 }
443
444 if (nr_oor > 0) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%" PRIu64 " out-of-range grid values converted to _FillValue [%s]\n", nr_oor, HH->name);
445
446 gmt_M_free (GMT, actual_col);
447
448 if (header->z_min <= header->z_max) {
449 limit[0] = header->z_min; limit[1] = header->z_max;
450 }
451 else {
452 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "No valid values in grid [%s]\n", HH->name);
453 limit[0] = 0.0; limit[1] = 0.0;
454 }
455 gmt_M_err_trap (nc_put_var_double (ncid, HH->z_id - 3, limit));
456
457 /* Close grid */
458
459 gmt_M_err_trap (gmt_nc_close (GMT, ncid));
460
461 return (GMT_NOERROR);
462 }
463