1 /*
2 * Code donated by David Divens, NOAA/NGDC
3 * Distributed under the GNU Public License (see LICENSE.TXT for details)
4 *--------------------------------------------------------------------*/
5
6 /* Comments added by P. Wessel:
7 *
8 * 1) GRD98 can now support pixel grids since 11/24/2006.
9 * 2) Learned that 1/x_inc and 1/y_inc are stored as integers. This means
10 * GRD98 imposes restrictions on x_inc & y_inc.
11 * 3) Added full support for padding and complex_mode on 3/17/2011
12 *
13 * Public functions (5):
14 *
15 * gmtlib_is_mgg2_grid :
16 * gmt_mgg2_read_grd_info :
17 * gmt_mgg2_write_grd_info :
18 * gmt_mgg2_read_grd :
19 * gmt_mgg2_write_grd :
20 */
21
22 #include "gmt_mgg_header2.h"
23
24 #define MGG_BYTE_SIZE
25
26 /* Local functions */
27
gmtmggheader2_swap_word(void * ptr)28 GMT_LOCAL void gmtmggheader2_swap_word (void* ptr) {
29 unsigned char *tmp = ptr;
30 unsigned char a = tmp[0];
31 tmp[0] = tmp[1];
32 tmp[1] = a;
33 }
34
gmtmggheader2_swap_long(void * ptr)35 GMT_LOCAL void gmtmggheader2_swap_long (void *ptr) {
36 unsigned char *tmp = ptr;
37 unsigned char a = tmp[0];
38 tmp[0] = tmp[3];
39 tmp[3] = a;
40 a = tmp[1];
41 tmp[1] = tmp[2];
42 tmp[2] = a;
43 }
44
gmtmggheader2_swap_mgg_header(MGG_GRID_HEADER_2 * header)45 GMT_LOCAL int gmtmggheader2_swap_mgg_header (MGG_GRID_HEADER_2 *header) {
46 int i, version;
47 /* Determine if swapping is needed */
48 if (header->version == (GRD98_MAGIC_NUM + GRD98_VERSION)) return (0); /* Version matches, No need to swap */
49 version = header->version;
50 gmtmggheader2_swap_long (&version);
51 if (version != (GRD98_MAGIC_NUM + GRD98_VERSION)) return (GMT_NOTSET); /* Cannot make sense of header */
52 /* Here we come when we do need to swap */
53 gmtmggheader2_swap_long (&header->version);
54 gmtmggheader2_swap_long (&header->length);
55 gmtmggheader2_swap_long (&header->dataType);
56 gmtmggheader2_swap_long (&header->latDeg);
57 gmtmggheader2_swap_long (&header->latMin);
58 gmtmggheader2_swap_long (&header->latSec);
59 gmtmggheader2_swap_long (&header->latSpacing);
60 gmtmggheader2_swap_long (&header->latNumCells);
61 gmtmggheader2_swap_long (&header->lonDeg);
62 gmtmggheader2_swap_long (&header->lonMin);
63 gmtmggheader2_swap_long (&header->lonSec);
64 gmtmggheader2_swap_long (&header->lonSpacing);
65 gmtmggheader2_swap_long (&header->lonNumCells);
66 gmtmggheader2_swap_long (&header->minValue);
67 gmtmggheader2_swap_long (&header->maxValue);
68 gmtmggheader2_swap_long (&header->gridRadius);
69 gmtmggheader2_swap_long (&header->precision);
70 gmtmggheader2_swap_long (&header->nanValue);
71 gmtmggheader2_swap_long (&header->numType);
72 gmtmggheader2_swap_long (&header->waterDatum);
73 gmtmggheader2_swap_long (&header->dataLimit);
74 gmtmggheader2_swap_long (&header->cellRegistration);
75 for (i = 0; i < GRD98_N_UNUSED; i++) gmtmggheader2_swap_long (&header->unused[i]);
76 return (1); /* Signal we need to swap the data also */
77 }
78
gmtmggheader2_dms2degrees(int deg,int min,int sec)79 GMT_LOCAL double gmtmggheader2_dms2degrees (int deg, int min, int sec) {
80 double decDeg = (double)deg;
81
82 decDeg += (double)min * GMT_MIN2DEG;
83 decDeg += (double)sec * GMT_SEC2DEG;
84
85 return decDeg;
86 }
87
gmtmggheader2_degrees2dms(double degrees,int * deg,int * min,int * sec)88 GMT_LOCAL void gmtmggheader2_degrees2dms (double degrees, int *deg, int *min, int *sec) {
89 /* Round off to the nearest half second */
90 if (degrees < 0) degrees -= (0.5 * GMT_SEC2DEG);
91
92 *deg = (int)degrees;
93 degrees -= *deg;
94
95 degrees *= GMT_DEG2MIN_F;
96 *min = (int)(degrees);
97 degrees -= *min;
98
99 *sec = (int)(degrees * GMT_MIN2SEC_F);
100 }
101
gmtmggheader2_GMTtoMGG2(struct GMT_GRID_HEADER * gmt,MGG_GRID_HEADER_2 * mgg)102 GMT_LOCAL int gmtmggheader2_GMTtoMGG2 (struct GMT_GRID_HEADER *gmt, MGG_GRID_HEADER_2 *mgg) {
103 double f;
104 gmt_M_memset (mgg, 1, MGG_GRID_HEADER_2);
105
106 mgg->version = GRD98_MAGIC_NUM + GRD98_VERSION;
107 mgg->length = sizeof (MGG_GRID_HEADER_2);
108 mgg->dataType = 1;
109
110 mgg->cellRegistration = gmt->registration;
111 mgg->lonNumCells = gmt->n_columns;
112 f = gmt->inc[GMT_X] * GMT_DEG2SEC_F;
113 mgg->lonSpacing = irint(f);
114 if (fabs (f - (double)mgg->lonSpacing) > GMT_CONV8_LIMIT) return (GMT_GRDIO_GRD98_XINC);
115 gmtmggheader2_degrees2dms(gmt->wesn[XLO], &mgg->lonDeg, &mgg->lonMin, &mgg->lonSec);
116
117 mgg->latNumCells = gmt->n_rows;
118 f = gmt->inc[GMT_Y] * GMT_DEG2SEC_F;
119 mgg->latSpacing = irint(gmt->inc[GMT_Y] * GMT_DEG2SEC_F);
120 if (fabs (f - (double)mgg->latSpacing) > GMT_CONV8_LIMIT) return (GMT_GRDIO_GRD98_YINC);
121 gmtmggheader2_degrees2dms(gmt->wesn[YHI], &mgg->latDeg, &mgg->latMin, &mgg->latSec);
122
123 /* Default values */
124 mgg->gridRadius = -1;
125 mgg->precision = GRD98_DEFAULT_PREC;
126 mgg->nanValue = GRD98_NAN_VALUE;
127 mgg->numType = sizeof (int);
128 mgg->minValue = irint(gmt->z_min * mgg->precision);
129 mgg->maxValue = irint(gmt->z_max * mgg->precision);
130
131 /* Data fits in two byte boundary */
132 if ((-SHRT_MAX <= mgg->minValue) && (mgg->maxValue <= SHRT_MAX)) {
133 mgg->numType = sizeof (short);
134 mgg->nanValue = (short)SHRT_MIN;
135 }
136 #ifdef MGG_BYTE_SIZE
137 /* Data fits in one byte boundary */
138 if ((gmt->z_min >= 0) && (gmt->z_max <= 127)) {
139 mgg->numType = sizeof (char);
140 mgg->nanValue = (char)255;
141 mgg->precision = 1;
142 mgg->minValue = irint (gmt->z_min);
143 mgg->maxValue = irint (gmt->z_max);
144 }
145 #endif
146 return (GMT_NOERROR);
147 }
148
gmtmggheader2_MGG2toGMT(MGG_GRID_HEADER_2 * mgg,struct GMT_GRID_HEADER * gmt)149 GMT_LOCAL void gmtmggheader2_MGG2toGMT (MGG_GRID_HEADER_2 *mgg, struct GMT_GRID_HEADER *gmt) {
150 int one_or_zero;
151 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (gmt);
152
153 /* Do not memset the gmt header since it has the file name set */
154
155 gmt->type = GMT_GRID_IS_RF;
156 gmt->registration = mgg->cellRegistration;
157 one_or_zero = 1 - gmt->registration;
158 gmt->n_columns = mgg->lonNumCells;
159 gmt->wesn[XLO] = gmtmggheader2_dms2degrees(mgg->lonDeg, mgg->lonMin, mgg->lonSec);
160 gmt->inc[GMT_X] = gmtmggheader2_dms2degrees(0, 0, mgg->lonSpacing);
161 gmt->wesn[XHI] = gmt->wesn[XLO] + (gmt->inc[GMT_X] * (gmt->n_columns - one_or_zero));
162
163 gmt->n_rows = mgg->latNumCells;
164 gmt->wesn[YHI] = gmtmggheader2_dms2degrees(mgg->latDeg, mgg->latMin, mgg->latSec);
165 gmt->inc[GMT_Y] = gmtmggheader2_dms2degrees(0, 0, mgg->latSpacing);
166 gmt->wesn[YLO] = gmt->wesn[YHI] - (gmt->inc[GMT_Y] * (gmt->n_rows - one_or_zero));
167
168 gmt->z_min = (double)mgg->minValue / (double)mgg->precision;
169 gmt->z_max = (double)mgg->maxValue / (double)mgg->precision;
170 gmt->z_scale_factor = 1.0;
171 gmt->z_add_offset = 0.0;
172 switch (mgg->numType) {
173 case 1: HH->orig_datatype = GMT_CHAR; break;
174 case 2: HH->orig_datatype = GMT_SHORT; break;
175 case 4: HH->orig_datatype = GMT_INT; break;
176 case -4: HH->orig_datatype = GMT_FLOAT; break;
177 default: HH->orig_datatype = GMT_INT; break;
178 }
179 }
180
181 /*----------------------------------------------------------|
182 * Public functions that are part of the GMT Devel library |
183 *----------------------------------------------------------|
184 */
185
gmtlib_is_mgg2_grid(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)186 int gmtlib_is_mgg2_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
187 /* Determine if file is a GRD98 file */
188 FILE *fp = NULL;
189 MGG_GRID_HEADER_2 mggHeader;
190 int ok;
191 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
192
193 if (!strcmp (HH->name, "="))
194 return (GMT_GRDIO_PIPE_CODECHECK); /* Cannot check on pipes */
195 if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.r_mode)) == NULL)
196 return (GMT_GRDIO_OPEN_FAILED);
197
198 gmt_M_memset (&mggHeader, 1, MGG_GRID_HEADER_2);
199 if (gmt_M_fread (&mggHeader, sizeof (MGG_GRID_HEADER_2), 1U, fp) != 1) {
200 gmt_fclose (GMT, fp);
201 return (GMT_GRDIO_READ_FAILED);
202 }
203
204 /* Swap header bytes if necessary; ok is 0|1 if successful and GMT_NOTSET if bad file */
205 ok = gmtmggheader2_swap_mgg_header (&mggHeader);
206
207 /* Check the magic number and size of header */
208 if (ok == GMT_NOTSET) {
209 gmt_fclose (GMT, fp);
210 return (GMT_NOTSET); /* Not this kind of file */
211 }
212 header->type = GMT_GRID_IS_RF;
213 gmt_fclose (GMT, fp);
214 return GMT_NOERROR;
215 }
216
gmt_mgg2_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)217 int gmt_mgg2_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
218 FILE *fp = NULL;
219 MGG_GRID_HEADER_2 mggHeader;
220 int ok;
221 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
222
223 if (!strcmp (HH->name, "="))
224 fp = GMT->session.std[GMT_IN];
225 else if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.r_mode)) == NULL)
226 return (GMT_GRDIO_OPEN_FAILED);
227
228 gmt_M_memset (&mggHeader, 1, MGG_GRID_HEADER_2);
229 if (gmt_M_fread (&mggHeader, sizeof (MGG_GRID_HEADER_2), 1U, fp) != 1) {
230 gmt_fclose (GMT, fp);
231 return (GMT_GRDIO_READ_FAILED);
232 }
233
234 /* Swap header bytes if necessary; ok is 0|1 if successful and GMT_NOTSET if bad file */
235 ok = gmtmggheader2_swap_mgg_header (&mggHeader);
236
237 /* Check the magic number and size of header */
238 if (ok == GMT_NOTSET) {
239 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unrecognized header, expected 0x%04X saw 0x%04X\n",
240 GRD98_MAGIC_NUM + GRD98_VERSION, mggHeader.version);
241 gmt_fclose (GMT, fp);
242 return (GMT_GRDIO_GRD98_BADMAGIC);
243 }
244
245 if (mggHeader.length != sizeof (MGG_GRID_HEADER_2)) {
246 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Invalid grid header size, expected %d, found %d\n",
247 (int)sizeof (MGG_GRID_HEADER_2), mggHeader.length);
248 gmt_fclose (GMT, fp);
249 return (GMT_GRDIO_GRD98_BADLENGTH);
250 }
251
252 gmt_fclose (GMT, fp);
253
254 gmtmggheader2_MGG2toGMT (&mggHeader, header);
255
256 return (GMT_NOERROR);
257 }
258
gmt_mgg2_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)259 int gmt_mgg2_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
260 FILE *fp = NULL;
261 MGG_GRID_HEADER_2 mggHeader;
262 int err;
263 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
264
265 if (!strcmp (HH->name, "="))
266 fp = GMT->session.std[GMT_OUT];
267 else if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.w_mode)) == NULL)
268 return (GMT_GRDIO_CREATE_FAILED);
269
270 if ((err = gmtmggheader2_GMTtoMGG2 (header, &mggHeader)) != 0) {
271 gmt_fclose (GMT, fp);
272 return (err);
273 }
274
275 if (gmt_M_fwrite (&mggHeader, sizeof (MGG_GRID_HEADER_2), 1U, fp) != 1) {
276 gmt_fclose (GMT, fp);
277 return (GMT_GRDIO_WRITE_FAILED);
278 }
279
280 gmt_fclose (GMT, fp);
281
282 return (GMT_NOERROR);
283 }
284
gmt_mgg2_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)285 int gmt_mgg2_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
286 MGG_GRID_HEADER_2 mggHeader;
287 FILE *fp = NULL;
288 int *tLong = NULL;
289 short *tShort = NULL;
290 char *tChar = NULL;
291 gmt_grdfloat *tFloat = NULL;
292 bool piping = false, is_float = false;
293 int j, first_col, last_col, first_row, last_row, swap_all = 0;
294 unsigned int width_in, height_in;
295 unsigned int i, width_out, *actual_col = NULL;
296 uint64_t kk, ij, j2, imag_offset;
297 off_t long_offset; /* For fseek only */
298 size_t n_expected; /* Items in one row */
299 size_t size; /* Item size */
300 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
301
302 gmt_M_memset (&mggHeader, 1, MGG_GRID_HEADER_2);
303 if (!strcmp (HH->name, "=")) {
304 fp = GMT->session.std[GMT_IN];
305 piping = true;
306 }
307 else if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.r_mode)) != NULL) {
308 if (gmt_M_fread (&mggHeader, sizeof (MGG_GRID_HEADER_2), 1U, fp) != 1) {
309 gmt_fclose (GMT, fp);
310 return (GMT_GRDIO_READ_FAILED);
311 }
312 swap_all = gmtmggheader2_swap_mgg_header (&mggHeader);
313 if (swap_all == GMT_NOTSET) {
314 gmt_fclose (GMT, fp);
315 return (GMT_GRDIO_GRD98_BADMAGIC);
316 }
317 if (mggHeader.numType == 0) mggHeader.numType = sizeof (int);
318 }
319 else
320 return (GMT_GRDIO_OPEN_FAILED);
321
322 is_float = (mggHeader.numType < 0 && abs (mggHeader.numType) == (int)sizeof (gmt_grdfloat)); /* Float file */
323
324 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);
325 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
326
327 width_out = width_in; /* Width of output array */
328 if (pad[XLO] > 0) width_out += pad[XLO];
329 if (pad[XHI] > 0) width_out += pad[XHI];
330
331 n_expected = header->n_columns;
332 tLong = gmt_M_memory (GMT, NULL, n_expected, int);
333 tShort = (short *)tLong; tChar = (char *)tLong; tFloat = (gmt_grdfloat *)tLong;
334 size = abs (mggHeader.numType);
335
336 if (piping) { /* Skip data by reading it */
337 for (j = 0; j < first_row; j++) if (gmt_M_fread ( tLong, size, n_expected, fp) != n_expected) {
338 gmt_fclose (GMT, fp);
339 gmt_M_free (GMT, actual_col); gmt_M_free (GMT, tLong);
340 return (GMT_GRDIO_READ_FAILED);
341 }
342 }
343 else if (first_row) { /* Simply seek by it */
344 long_offset = (off_t)first_row * (off_t)n_expected * size;
345 if (fseek (fp, long_offset, 1)) {
346 gmt_fclose (GMT, fp);
347 gmt_M_free (GMT, actual_col); gmt_M_free (GMT, tLong);
348 return (GMT_GRDIO_SEEK_FAILED);
349 }
350 }
351
352 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
353 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
354 for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
355 if (gmt_M_fread (tLong, size, n_expected, fp) != n_expected) {
356 gmt_fclose (GMT, fp);
357 gmt_M_free (GMT, actual_col); gmt_M_free (GMT, tLong);
358 return (GMT_GRDIO_READ_FAILED);
359 }
360 ij = imag_offset + (j2 + pad[YHI]) * width_out + pad[XLO];
361 for (i = 0; i < width_in; i++) {
362 kk = ij + i;
363 if (mggHeader.numType == sizeof (int)) {
364 if (swap_all) gmtmggheader2_swap_long (&tLong[actual_col[i]]);
365 if (tLong[actual_col[i]] == mggHeader.nanValue) grid[kk] = GMT->session.f_NaN;
366 else grid[kk] = (gmt_grdfloat) tLong[actual_col[i]] / (gmt_grdfloat) mggHeader.precision;
367 }
368 else if (is_float) {
369 if (swap_all) gmtmggheader2_swap_long (&tLong[actual_col[i]]);
370 if (tLong[actual_col[i]] == mggHeader.nanValue) grid[kk] = GMT->session.f_NaN;
371 else grid[kk] = tFloat[actual_col[i]];
372 }
373
374 else if (mggHeader.numType == sizeof (short)) {
375 if (swap_all) gmtmggheader2_swap_word(&tShort[actual_col[i]]);
376 if (tShort[actual_col[i]] == mggHeader.nanValue) grid[kk] = GMT->session.f_NaN;
377 else grid[kk] = (gmt_grdfloat) tShort[actual_col[i]] / (gmt_grdfloat) mggHeader.precision;
378 }
379
380 else if (mggHeader.numType == sizeof (char)) {
381 if (tChar[actual_col[i]] == mggHeader.nanValue) grid[kk] = GMT->session.f_NaN;
382 else grid[kk] = (gmt_grdfloat) tChar[actual_col[i]] / (gmt_grdfloat) mggHeader.precision;
383 }
384 else {
385 gmt_fclose (GMT, fp);
386 gmt_M_free (GMT, actual_col);
387 gmt_M_free (GMT, tLong);
388 return (GMT_GRDIO_UNKNOWN_TYPE);
389 }
390 if (gmt_M_is_fnan (grid[kk])) {
391 HH->has_NaNs = GMT_GRID_HAS_NANS;
392 continue;
393 }
394 /* Update z_min, z_max */
395 header->z_min = MIN (header->z_min, (double)grid[kk]);
396 header->z_max = MAX (header->z_max, (double)grid[kk]);
397 }
398 }
399 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
400 header->z_min = header->z_max = NAN;
401 if (piping) { /* Skip data by reading it */
402 int n_rows = header->n_rows;
403 for (j = last_row + 1; j < n_rows; j++) {
404 if (gmt_M_fread ( tLong, size, n_expected, fp) != n_expected) {
405 gmt_M_free (GMT, tLong);
406 gmt_M_free (GMT, actual_col);
407 gmt_fclose (GMT, fp);
408 return (GMT_GRDIO_READ_FAILED);
409 }
410 }
411 }
412
413 gmt_M_free (GMT, tLong);
414 gmt_M_free (GMT, actual_col);
415
416 header->n_columns = width_in;
417 header->n_rows = height_in;
418 gmt_M_memcpy (header->wesn, wesn, 4, double);
419
420 gmt_fclose (GMT, fp);
421
422 return (GMT_NOERROR);
423 }
424
gmt_mgg2_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)425 int gmt_mgg2_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
426 MGG_GRID_HEADER_2 mggHeader;
427 bool is_float = false, check;
428 int i, j, err;
429 unsigned int i2, ju, iu, width_out, height_out, width_in, *actual_col = NULL;
430 int first_col, last_col, first_row, last_row;
431 uint64_t ij, kk, j2, imag_offset;
432 size_t size;
433
434 int *tLong = NULL;
435 short *tShort = NULL;
436 char *tChar = NULL;
437 gmt_grdfloat *tFloat = NULL;
438 FILE *fp = NULL;
439 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
440
441 if (!strcmp (HH->name, "="))
442 fp = GMT->session.std[GMT_OUT];
443 else if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.w_mode)) == NULL)
444 return (GMT_GRDIO_CREATE_FAILED);
445
446 check = !isnan (header->nan_value);
447
448 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);
449 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
450
451 width_in = width_out; /* Physical width of input array */
452 if (pad[XLO] > 0) width_in += pad[XLO];
453 if (pad[XHI] > 0) width_in += pad[XHI];
454
455 gmt_M_memcpy (header->wesn, wesn, 4, double);
456
457 /* Find xmin/zmax */
458
459 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
460 for (j = first_row, j2 = pad[YHI]; j <= last_row; j++, j2++) {
461 ij = imag_offset + j2 * width_in;
462 for (i = first_col, i2 = pad[XLO]; i <= last_col; i++, i2++) {
463 kk = ij + i2;
464 if (isnan (grid[kk])) {
465 if (check) grid[kk] = header->nan_value;
466 }
467 else {
468 header->z_min = MIN (header->z_min, (double)grid[kk]);
469 header->z_max = MAX (header->z_max, (double)grid[kk]);
470 }
471 }
472 }
473 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
474 header->z_min = header->z_max = NAN;
475
476 /* store header information and array */
477 if ((err = gmtmggheader2_GMTtoMGG2(header, &mggHeader)) != 0) {
478 gmt_fclose (GMT, fp);
479 gmt_M_free (GMT, actual_col);
480 return (err);
481 }
482 if (gmt_M_fwrite (&mggHeader, sizeof (MGG_GRID_HEADER_2), 1U, fp) != 1) {
483 gmt_M_free (GMT, actual_col);
484 gmt_fclose (GMT, fp);
485 return (GMT_GRDIO_WRITE_FAILED);
486 }
487 is_float = (mggHeader.numType < 0 && abs (mggHeader.numType) == (int)sizeof (gmt_grdfloat)); /* Float file */
488
489 tLong = gmt_M_memory (GMT, NULL, width_in, int);
490 tShort = (short *) tLong; tChar = (char *)tLong; tFloat = (gmt_grdfloat *) tLong;
491
492 i2 = first_col + pad[XLO];
493 size = abs (mggHeader.numType);
494 for (ju = 0, j2 = first_row + pad[YHI]; ju < height_out; ju++, j2++) {
495 ij = imag_offset + j2 * width_in + i2;
496 for (iu = 0; iu < width_out; iu++) {
497 kk = ij+actual_col[iu];
498 if (gmt_M_is_fnan (grid[kk])) {
499 if (mggHeader.numType == sizeof (int)) tLong[iu] = mggHeader.nanValue;
500 else if (is_float) tFloat[iu] = (gmt_grdfloat)mggHeader.nanValue;
501 else if (mggHeader.numType == sizeof (short)) tShort[iu] = (short)mggHeader.nanValue;
502 else if (mggHeader.numType == sizeof (char)) tChar[iu] = (char)mggHeader.nanValue;
503 else {
504 gmt_M_free (GMT, tLong);
505 gmt_M_free (GMT, actual_col);
506 gmt_fclose (GMT, fp);
507 return (GMT_GRDIO_UNKNOWN_TYPE);
508 }
509 } else {
510 if (grid[kk] > -0.1 && grid[kk] < 0) grid[kk] = (gmt_grdfloat)(-0.1);
511
512 if (mggHeader.numType == sizeof (int))
513 tLong[iu] = (int)rint ((double)grid[kk] * mggHeader.precision);
514 else if (is_float)
515 tFloat[iu] = grid[kk];
516 else if (mggHeader.numType == sizeof (short))
517 tShort[iu] = (short) rint((double)grid[kk] * mggHeader.precision);
518 else if (mggHeader.numType == sizeof (char))
519 tChar[iu] = (char) rint((double)grid[kk] * mggHeader.precision);
520 else {
521 gmt_M_free (GMT, tLong);
522 gmt_M_free (GMT, actual_col);
523 gmt_fclose (GMT, fp);
524 return (GMT_GRDIO_UNKNOWN_TYPE);
525 }
526 }
527 }
528 if (gmt_M_fwrite (tLong, size, width_out, fp) != width_out) {
529 gmt_M_free (GMT, tLong);
530 gmt_M_free (GMT, actual_col);
531 gmt_fclose (GMT, fp);
532 return (GMT_GRDIO_WRITE_FAILED);
533 }
534 }
535 gmt_M_free (GMT, tLong);
536 gmt_M_free (GMT, actual_col);
537
538 gmt_fclose (GMT, fp);
539
540 return (GMT_NOERROR);
541 }
542