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 wi1552ll 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 U S T O M I O . C R O U T I N E S
20 *
21 * GMT has built-in support for several grid formats, which is extended
22 * to a broad selection of formats via GDAL. If there are very unusual
23 * and custom formats you wish to see supported by GMT then please make
24 * your case on the GMT trac website and we can start a dialog of whether
25 * or not we should support your suggested format.
26 *
27 * For gurus: To add another data format:
28 *
29 * 1. Write the five required routines (see below).
30 * 2. Edit gmt_grdio.h and add one entry to the Gmt_grid_id enum
31 * 3. Increment GMT_N_GRD_FORMATS (in gmt_grdio.h too)
32 *
33 * Author: Paul Wessel
34 * Date: 1-JAN-2010
35 * Version: 5
36 *
37 * Functions include:
38 *
39 * GMT_*_read_grd_info : Read header from file
40 * GMT_*_read_grd : Read header and data set from file
41 * GMT_*_update_grd_info : Update header in existing file
42 * GMT_*_write_grd_info : Write header to new file
43 * GMT_*_write_grd : Write header and data set to new file
44 *
45 * where * is a tag specific to a particular data format
46 *
47 * NOTE: 1. GMT assumes that gmtlib_read_grd_info has been called before calls
48 * to gmtlib_read_grd. This normally is done via GMT_Read_Data.
49 * 2. Some formats may permit pipes to be used. In that case GMT
50 * expects the filename to be "=" (the equal sign). It is the
51 * responsibility of the custom routines to test for "=" and
52 * give error messages if piping is not supported for this format
53 * (e.g., netcdf uses fseek and can therefore not use pipes; other
54 * formats may have similar limitations).
55 * 3. For most formats the write_grd_info and update_grd_info is the
56 * same function (but netCDF is one exception)
57 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
58
59 #include "gmt_dev.h"
60 #include "gmt_internals.h"
61
62 /* Defined in gmt_cdf.c */
63 int gmt_cdf_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
64 int gmt_cdf_update_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
65 int gmt_cdf_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
66 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);
67 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);
68
69 /* Defined in gmt_nc.c */
70 int gmt_nc_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
71 int gmt_nc_update_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
72 int gmt_nc_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header);
73 int gmt_nc_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode);
74 int gmt_nc_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode);
75
76 /* CUSTOM I/O FUNCTIONS FOR GRIDDED DATA FILES */
77
78 /*-----------------------------------------------------------
79 * Format : dummy
80 * Purpose :
81 * Use this function to direct all unsupported formats to.
82 * Functions : gmt_dummy_grd_info, gmt_dummy_grd_read
83 *-----------------------------------------------------------*/
84
gmt_dummy_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)85 int gmt_dummy_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
86 if (header) GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unknown grid format.\n");
87 return (GMT_GRDIO_UNKNOWN_FORMAT);
88 }
89
gmt_dummy_grd_read(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)90 int gmt_dummy_grd_read (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
91 if (header && grid && wesn && pad && complex_mode < 1024) GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unknown grid format.\n");
92 return (GMT_GRDIO_UNKNOWN_FORMAT);
93 }
94
95 /*-----------------------------------------------------------
96 * Format : rb
97 * Type : Standard 8-bit Sun rasterfiles
98 * Prefix : GMT_ras_
99 * Author : Paul Wessel, SOEST
100 * Date : 17-JUN-1999
101 * Purpose: Rasterfiles may often be used to store
102 * datasets of limited dynamic range where
103 * 8-bits is all that is needed. Since the
104 * usual information like w,e,s,n is not part
105 * of a rasterfile's header, we assign default
106 * values as follows:
107 * w = s = 0.
108 * e = ras_width;
109 * n = ras_height;
110 * dx = dy = 1
111 * Such files are always pixel registered
112 *
113 * Functions : gmt_ras_read_grd_info,
114 * gmt_ras_write_grd_info, gmt_ras_read_grd, gmt_ras_write_grd
115 *-----------------------------------------------------------*/
116
117 #define RAS_MAGIC 0x59a66a95
118
gmtcustomio_read_rasheader(FILE * fp,struct rasterfile * h)119 GMT_LOCAL int gmtcustomio_read_rasheader (FILE *fp, struct rasterfile *h) {
120 /* Reads the header of a Sun rasterfile byte by byte
121 since the format is defined as the byte order on the
122 PDP-11.
123 */
124
125 uint8_t byte[4];
126 int32_t i, value;
127
128 for (i = 0; i < 8; i++) {
129
130 if (gmt_M_fread (byte, sizeof (uint8_t), 4, fp) != 4)
131 return (GMT_GRDIO_READ_FAILED);
132
133 value = (byte[0] << 24) + (byte[1] << 16) + (byte[2] << 8) + byte[3];
134
135 switch (i) {
136 case 0: h->magic = value; break;
137 case 1: h->width = value; break;
138 case 2: h->height = value; break;
139 case 3: h->depth = value; break;
140 case 4: h->length = value; break;
141 case 5: h->type = value; break;
142 case 6: h->maptype = value; break;
143 case 7: h->maplength = value; break;
144 }
145 }
146
147 if (h->type == RT_OLD && h->length == 0) h->length = 2 * irint (ceil (h->width * h->depth / 16.0)) * h->height;
148
149 return (GMT_NOERROR);
150 }
151
gmtcustomio_write_rasheader(FILE * fp,struct rasterfile * h)152 GMT_LOCAL int gmtcustomio_write_rasheader (FILE *fp, struct rasterfile *h) {
153 /* Writes the header of a Sun rasterfile byte by byte
154 since the format is defined as the byte order on the
155 PDP-11.
156 */
157
158 int i;
159 uint8_t byte[4];
160 int32_t value;
161
162 if (h->type == RT_OLD && h->length == 0) {
163 h->length = 2 * irint (ceil (h->width * h->depth / 16.0)) * h->height;
164 h->type = RT_STANDARD;
165 }
166
167 for (i = 0; i < 8; i++) {
168
169 switch (i) {
170 case 0: value = h->magic; break;
171 case 1: value = h->width; break;
172 case 2: value = h->height; break;
173 case 3: value = h->depth; break;
174 case 4: value = h->length; break;
175 case 5: value = h->type; break;
176 case 6: value = h->maptype; break;
177 default: value = h->maplength; break;
178 }
179 byte[0] = (uint8_t)((value >> 24) & 0xFF);
180 byte[1] = (uint8_t)((value >> 16) & 0xFF);
181 byte[2] = (uint8_t)((value >> 8) & 0xFF);
182 byte[3] = (uint8_t)(value & 0xFF);
183
184 if (gmt_M_fwrite (byte, sizeof (uint8_t), 4, fp) != 4)
185 return (GMT_GRDIO_WRITE_FAILED);
186 }
187
188 return (GMT_NOERROR);
189 }
190
gmtlib_is_ras_grid(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)191 int gmtlib_is_ras_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
192 /* Determine if file is a Sun rasterfile */
193 FILE *fp = NULL;
194 struct rasterfile h;
195 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
196 if (!strcmp (HH->name, "="))
197 return (GMT_GRDIO_PIPE_CODECHECK); /* Cannot check on pipes */
198 if ((fp = gmt_fopen (GMT, HH->name, "rb")) == NULL)
199 return (GMT_GRDIO_OPEN_FAILED);
200 gmt_M_memset (&h, 1, struct rasterfile);
201 if (gmtcustomio_read_rasheader (fp, &h)) {
202 gmt_fclose (GMT, fp);
203 return (GMT_GRDIO_READ_FAILED);
204 }
205 if (h.magic != RAS_MAGIC) {
206 gmt_fclose (GMT, fp);
207 return (GMT_GRDIO_NOT_RAS);
208 }
209 if (h.type != 1 || h.depth != 8) {
210 gmt_fclose (GMT, fp);
211 return (GMT_GRDIO_NOT_8BIT_RAS);
212 }
213 header->type = GMT_GRID_IS_RB;
214 gmt_fclose (GMT, fp);
215 return GMT_NOERROR;
216 }
217
gmt_ras_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)218 int gmt_ras_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
219 FILE *fp = NULL;
220 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
221 struct rasterfile h;
222 unsigned char u;
223 int i;
224
225 if (!strcmp (HH->name, "=")) { /* Read from pipe */
226 #ifdef SET_IO_MODE
227 gmt_setmode (GMT, GMT_IN);
228 #endif
229 fp = GMT->session.std[GMT_IN];
230 }
231 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) == NULL)
232 return (GMT_GRDIO_OPEN_FAILED);
233
234 gmt_M_memset (&h, 1, struct rasterfile);
235 if (gmtcustomio_read_rasheader (fp, &h)) {
236 gmt_fclose (GMT, fp);
237 return (GMT_GRDIO_READ_FAILED);
238 }
239 if (h.magic != RAS_MAGIC) {
240 gmt_fclose (GMT, fp);
241 return (GMT_GRDIO_NOT_RAS);
242 }
243 if (h.type != 1 || h.depth != 8) {
244 gmt_fclose (GMT, fp);
245 return (GMT_GRDIO_NOT_8BIT_RAS);
246 }
247
248 for (i = 0; i < h.maplength; i++) {
249 if (gmt_M_fread (&u, sizeof (unsigned char), 1U, fp) < 1U) {
250 gmt_fclose (GMT, fp);
251 return (GMT_GRDIO_READ_FAILED); /* Skip colormap by reading since fp could be stdin */
252 }
253 }
254 gmt_fclose (GMT, fp);
255
256 /* Since we have no info on boundary values, just use integer size and steps = 1 */
257
258 header->wesn[XLO] = header->wesn[YLO] = 0.0;
259 header->wesn[XHI] = header->n_columns = h.width;
260 header->wesn[YHI] = header->n_rows = h.height;
261 header->inc[GMT_X] = header->inc[GMT_Y] = 1.0;
262 header->registration = GMT_GRID_PIXEL_REG; /* Always pixel format */
263 header->z_scale_factor = 1.0; header->z_add_offset = 0.0;
264 HH->orig_datatype = GMT_CHAR;
265
266 return (GMT_NOERROR);
267 }
268
gmt_ras_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)269 int gmt_ras_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
270 FILE *fp = NULL;
271 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
272 struct rasterfile h;
273
274 if (!strcmp (HH->name, "=")) { /* Write to pipe */
275 #ifdef SET_IO_MODE
276 gmt_setmode (GMT, GMT_OUT);
277 #endif
278 fp = GMT->session.std[GMT_OUT];
279 }
280 else if ((fp = gmt_fopen (GMT, HH->name, "rb+")) == NULL && (fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
281 return (GMT_GRDIO_CREATE_FAILED);
282
283 h.magic = RAS_MAGIC;
284 h.width = header->n_columns;
285 h.height = header->n_rows;
286 h.depth = 8;
287 h.length = header->n_rows * irint (ceil (header->n_columns/2.0)) * 2;
288 h.type = 1;
289 h.maptype = h.maplength = 0;
290
291 if (gmtcustomio_write_rasheader (fp, &h)) {
292 gmt_fclose (GMT, fp);
293 return (GMT_GRDIO_WRITE_FAILED);
294 }
295
296 gmt_fclose (GMT, fp);
297
298 return (GMT_NOERROR);
299 }
300
gmt_ras_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)301 int gmt_ras_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
302 /* header: grid structure header */
303 /* grid: array with final grid */
304 /* wesn: Sub-region to extract [Use entire file if 0,0,0,0] */
305 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
306 /* complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) */
307 /* Note: The file has only real values, we simply allow space in the complex array */
308 /* for real and imaginary parts when processed by grdfft etc. */
309
310 bool piping = false, check;
311 int j, first_col, last_col, first_row, last_row;
312 unsigned int i, width_in, height_in, *actual_row = NULL;
313 size_t n2;
314 uint64_t kk, ij, j2, width_out, imag_offset;
315 FILE *fp = NULL;
316 unsigned char *tmp = NULL;
317 struct rasterfile h;
318 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
319
320 gmt_M_memset (&h, 1, struct rasterfile);
321 if (!strcmp (HH->name, "=")) { /* Read from pipe */
322 #ifdef SET_IO_MODE
323 gmt_setmode (GMT, GMT_IN);
324 #endif
325 fp = GMT->session.std[GMT_IN];
326 piping = true;
327 }
328 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) != NULL) { /* Skip header */
329 if (gmtcustomio_read_rasheader (fp, &h)) {
330 gmt_fclose (GMT, fp);
331 return (GMT_GRDIO_READ_FAILED);
332 }
333 if (h.maplength && fseek (fp, (off_t) h.maplength, SEEK_CUR)) {
334 gmt_fclose (GMT, fp);
335 return (GMT_GRDIO_SEEK_FAILED);
336 }
337 }
338 else
339 return (GMT_GRDIO_OPEN_FAILED);
340
341 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
342
343 n2 = lrint (ceil (header->n_columns / 2.0)) * 2; /* Sun 8-bit rasters are stored using 16-bit words */
344 tmp = gmt_M_memory (GMT, NULL, n2, unsigned char);
345
346 check = !isnan (header->nan_value);
347
348 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_row), HH->name);
349
350 width_out = width_in; /* Width of output array */
351 if (pad[XLO] > 0) width_out += pad[XLO];
352 if (pad[XHI] > 0) width_out += pad[XHI];
353
354 if (piping) { /* Skip data by reading it */
355 for (j = 0; j < first_row; j++) {
356 if (gmt_M_fread (tmp, sizeof (unsigned char), n2, fp) < n2) {
357 gmt_M_free (GMT, actual_row);
358 gmt_M_free (GMT, tmp);
359 return (GMT_GRDIO_READ_FAILED);
360 }
361 }
362 }
363 else {/* Simply seek by it */
364 if (first_row && fseek (fp, (off_t) (first_row * n2 * sizeof (unsigned char)), SEEK_CUR)) {
365 gmt_fclose (GMT, fp);
366 gmt_M_free (GMT, actual_row);
367 gmt_M_free (GMT, tmp);
368 return (GMT_GRDIO_SEEK_FAILED);
369 }
370 }
371
372 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
373 HH->has_NaNs = GMT_GRID_NO_NANS; /* Cannot have nans in a raster */
374
375 for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
376 ij = imag_offset + (j2 + pad[3]) * width_out + pad[XLO];
377 if (gmt_M_fread ( tmp, sizeof (unsigned char), n2, fp) < n2) {
378 if (!piping) gmt_fclose (GMT, fp);
379 gmt_M_free (GMT, actual_row);
380 gmt_M_free (GMT, tmp);
381 return (GMT_GRDIO_READ_FAILED);
382 }
383 for (i = 0; i < width_in; i++) {
384 kk = ij + i;
385 grid[kk] = (gmt_grdfloat) tmp[actual_row[i]];
386 if (check && grid[kk] == header->nan_value)
387 grid[kk] = GMT->session.f_NaN;
388 if (gmt_M_is_fnan (grid[kk])) continue;
389 /* Update z min/max */
390 header->z_min = MIN (header->z_min, (double)grid[kk]);
391 header->z_max = MAX (header->z_max, (double)grid[kk]);
392 }
393 }
394 if (piping) { /* Skip data by reading it */
395 int n_rows = header->n_rows;
396 for (j = last_row + 1; j < n_rows; j++) if (gmt_M_fread ( tmp, sizeof (unsigned char), n2, fp) < n2) {
397 gmt_M_free (GMT, actual_row);
398 gmt_M_free (GMT, tmp);
399 return (GMT_GRDIO_READ_FAILED);
400 }
401 }
402 header->n_columns = width_in;
403 header->n_rows = height_in;
404 gmt_M_memcpy (header->wesn, wesn, 4, double);
405
406 if (!piping) gmt_fclose (GMT, fp);
407
408 gmt_M_free (GMT, actual_row);
409 gmt_M_free (GMT, tmp);
410 return (GMT_NOERROR);
411 }
412
gmt_ras_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)413 int gmt_ras_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
414 /* header: grid structure header */
415 /* grid: array with final grid */
416 /* wesn: Sub-region to write [Use entire file if 0,0,0,0] */
417 /* padding: # of empty rows/columns to remove on w, e, s, n of grid, respectively */
418 /* complex_mode: &1 | &2 if complex array is to hold real (1) and imaginary (2) parts (otherwise read as real only) */
419 /* Note: The file has only real values, we simply allow space in the complex array */
420 /* for real and imaginary parts when processed by grdfft etc. */
421 /* If 64 is added we write no header */
422
423 bool check;
424 unsigned int i, i2, j, width_out, height_out, n2, *actual_col = NULL;
425 int first_col, last_col, first_row, last_row;
426 uint64_t kk, ij, j2, width_in, imag_offset;
427 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
428
429 unsigned char *tmp = NULL;
430
431 FILE *fp = NULL;
432
433 struct rasterfile h;
434
435 if (!strcmp (HH->name, "=")) { /* Write to pipe */
436 #ifdef SET_IO_MODE
437 gmt_setmode (GMT, GMT_OUT);
438 #endif
439 fp = GMT->session.std[GMT_OUT];
440 }
441 else if ((fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
442 return (GMT_GRDIO_CREATE_FAILED);
443
444 h.magic = RAS_MAGIC;
445 h.width = header->n_columns;
446 h.height = header->n_rows;
447 h.depth = 8;
448 h.length = header->n_rows * irint (ceil (header->n_columns/2.0)) * 2;
449 h.type = 1;
450 h.maptype = h.maplength = 0;
451
452 n2 = irint (ceil (header->n_columns / 2.0)) * 2;
453 tmp = gmt_M_memory (GMT, NULL, n2, unsigned char);
454
455 check = !isnan (header->nan_value);
456
457 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);
458
459 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
460
461 width_in = width_out; /* Physical width of input array */
462 if (pad[XLO] > 0) width_in += pad[XLO];
463 if (pad[XHI] > 0) width_in += pad[XHI];
464
465 gmt_M_memcpy (header->wesn, wesn, 4, double);
466
467 h.width = header->n_columns;
468 h.height = header->n_rows;
469 h.length = header->n_rows * irint (ceil (header->n_columns/2.0)) * 2;
470
471 /* store header information and array */
472
473 if (gmtcustomio_write_rasheader (fp, &h)) {
474 gmt_fclose (GMT, fp);
475 gmt_M_free (GMT, actual_col);
476 gmt_M_free (GMT, tmp);
477 return (GMT_GRDIO_WRITE_FAILED);
478 }
479
480 i2 = first_col + pad[XLO];
481 for (j = 0, j2 = first_row + pad[3]; j < height_out; j++, j2++) {
482 ij = imag_offset + j2 * width_in + i2;
483 for (i = 0; i < width_out; i++) {
484 kk = ij + actual_col[i];
485 if (check && gmt_M_is_fnan (grid[kk])) grid[kk] = header->nan_value;
486 tmp[i] = (unsigned char) grid[kk];
487 }
488 if (gmt_M_fwrite (tmp, sizeof (unsigned char), n2, fp) < n2) {
489 gmt_fclose (GMT, fp);
490 gmt_M_free (GMT, actual_col);
491 gmt_M_free (GMT, tmp);
492 return (GMT_GRDIO_WRITE_FAILED);
493 }
494 }
495 gmt_fclose (GMT, fp);
496
497 gmt_M_free (GMT, actual_col);
498 gmt_M_free (GMT, tmp);
499
500 return (GMT_NOERROR);
501
502 }
503
504 /*-----------------------------------------------------------
505 * Format : bm
506 * Type : Native binary (bit) C file
507 * Prefix : GMT_bit_
508 * Author : Paul Wessel, SOEST
509 * Date : 27-JUN-1994
510 * Purpose: The native binary bit format is used
511 * primarily for piped i/o between programs
512 * that otherwise would use temporary, large
513 * intermediate grdfiles. Note that not all
514 * programs can support piping (they may have
515 * to re-read the file or accept more than one
516 * grdfile). Datasets containing ON/OFF information
517 * like bitmasks can be stored using bits
518 * We use 4-byte integers to store 32 bits at the time
519 * Functions : gmt_bit_read_grd, gmt_bit_write_grd
520 *-----------------------------------------------------------*/
521
522 /* sizeof the first part of the native header */
523 #define SIZEOF_NATIVE_GRD_HDR1 12 /* 3 * sizeof (unsigned) */
524 /* sizeof the last part of the native header */
525 #define SIZEOF_NATIVE_GRD_HDR2 880 /* sizeof (struct GMT_GRID_HEADER) - (size_t)header->wesn */
526 /* sizeof the complete native header */
527 #define SIZEOF_NATIVE_GRD_HDR 892 /* SIZEOF_NATIVE_GRD_HDR1 + SIZEOF_NATIVE_GRD_HD */
528
gmtcustomio_native_read_grd_header(FILE * fp,struct GMT_GRID_HEADER * header)529 GMT_LOCAL int gmtcustomio_native_read_grd_header (FILE *fp, struct GMT_GRID_HEADER *header) {
530 int err = GMT_NOERROR;
531 /* Because GMT_GRID_HEADER is not 64-bit aligned we must read it in parts */
532 if (gmt_M_fread (&header->n_columns, SIZEOF_NATIVE_GRD_HDR1, 1U, fp) != 1 ||
533 gmt_M_fread (&header->wesn[0], SIZEOF_NATIVE_GRD_HDR2, 1U, fp) != 1)
534 err = GMT_GRDIO_READ_FAILED;
535 return (err);
536 }
537
gmt_native_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)538 int gmt_native_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
539 /* Read GRD header structure from native binary file. This is used by
540 * all the native binary formats in GMT */
541
542 int status;
543 FILE *fp = NULL;
544 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
545
546 if (!strcmp (HH->name, "=")) { /* Read from pipe */
547 #ifdef SET_IO_MODE
548 gmt_setmode (GMT, GMT_IN);
549 #endif
550 fp = GMT->session.std[GMT_IN];
551 }
552 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) == NULL)
553 return (GMT_GRDIO_OPEN_FAILED);
554
555 status = gmtcustomio_native_read_grd_header (fp, header);
556 gmt_fclose (GMT, fp);
557 return status;
558 }
559
560 #define free_and_return(code) { gmt_free_header(GMT, &t_head); return (code); }
561
gmtlib_is_native_grid(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)562 int gmtlib_is_native_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
563 uint64_t mx, status, size;
564 off_t nm;
565 double item_size;
566 struct stat buf;
567 struct GMT_GRID_HEADER *t_head = gmt_get_header (GMT);
568 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
569 struct GMT_GRID_HEADER_HIDDEN *HHt = gmt_get_H_hidden (t_head);
570 if (!strcmp (HH->name, "=")) {
571 free_and_return (GMT_GRDIO_PIPE_CODECHECK); /* Cannot check on pipes */
572 }
573 if (stat (HH->name, &buf)) {
574 free_and_return (GMT_GRDIO_STAT_FAILED); /* Inquiry about file failed somehow */
575 }
576 strncpy (HHt->name, HH->name, GMT_LEN256);
577 if ((status = gmt_native_read_grd_info (GMT, t_head)) != 0) {
578 free_and_return (GMT_GRDIO_READ_FAILED); /* Failed to read header */
579 }
580 if (t_head->n_columns <= 0 || t_head->n_rows <= 0 || !(t_head->registration == GMT_GRID_NODE_REG || t_head->registration == GMT_GRID_PIXEL_REG)) {
581 free_and_return (GMT_GRDIO_BAD_VAL); /* Garbage for n_columns or n_rows */
582 }
583 if (t_head->wesn[XLO] >= t_head->wesn[XHI] || t_head->wesn[YLO] >= t_head->wesn[YHI]) {
584 free_and_return (GMT_GRDIO_BAD_VAL); /* Garbage for wesn */
585 }
586 nm = gmt_M_get_nm (GMT, t_head->n_columns, t_head->n_rows);
587 if (nm <= 0) {
588 free_and_return (GMT_GRDIO_BAD_VAL); /* Overflow for n_columns * n_rows? */
589 }
590 item_size = (double)((buf.st_size - GMT_GRID_HEADER_SIZE) / nm); /* Estimate size of elements */
591 size = lrint (item_size);
592 if (!doubleAlmostEqualZero (item_size, (double)size)) {
593 free_and_return (GMT_GRDIO_BAD_VAL); /* Size not an integer */
594 }
595
596 switch (size) {
597 case 0: /* Possibly bit map; check some more */
598 mx = lrint (ceil (t_head->n_columns / 32.0)); /* Number of 4-byte integers it would take per row to store the bits */
599 nm = 4 * mx * ((uint64_t)t_head->n_rows); /* Number of bytes to store all the rows */
600 if ((buf.st_size - GMT_GRID_HEADER_SIZE) == nm) /* Yes, it was a bit mask file */
601 header->type = GMT_GRID_IS_BM;
602 else { /* No, junk data */
603 free_and_return (GMT_GRDIO_BAD_VAL);
604 }
605 HH->orig_datatype = GMT_INT;
606 break;
607 case 1: /* 1-byte elements */
608 header->type = GMT_GRID_IS_BB;
609 HH->orig_datatype = GMT_CHAR;
610 break;
611 case 2: /* 2-byte short int elements */
612 header->type = GMT_GRID_IS_BS;
613 HH->orig_datatype = GMT_SHORT;
614 break;
615 case 4: /* 4-byte elements - could be int or float */
616 /* See if we can decide it is a float grid */
617 if (gmt_M_compat_check (GMT, 4)) {
618 GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Will try to determine if a native 4-byte grid is float or int but may be wrong.\n");
619 GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Please append =bf (float) or =bi (integer) to avoid this situation.\n");
620 /* Naive test to see if we can decide it is a float grid */
621 if ((t_head->z_scale_factor == 1.0 && t_head->z_add_offset == 0.0) || fabs((t_head->z_min/t_head->z_scale_factor) - rint(t_head->z_min/t_head->z_scale_factor)) > GMT_CONV8_LIMIT || fabs((t_head->z_max/t_head->z_scale_factor) - rint(t_head->z_max/t_head->z_scale_factor)) > GMT_CONV8_LIMIT) {
622 GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Based on header values we guessed the grid is 4-byte float. If wrong you must add =bi.\n");
623 header->type = GMT_GRID_IS_BF;
624 HH->orig_datatype = GMT_FLOAT;
625 }
626 else {
627 GMT_Report (GMT->parent, GMT_MSG_COMPAT, "Based on header values we guessed the grid is 4-byte int. If wrong you must add =bf.\n");
628 header->type = GMT_GRID_IS_BI;
629 HH->orig_datatype = GMT_INT;
630 }
631 }
632 else {
633 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot determine if a native 4-byte grid is float or int without more information.\n");
634 GMT_Report (GMT->parent, GMT_MSG_ERROR, "You must append =bf (float) or =bi (integer) to avoid this situation.\n");
635 free_and_return (GMT_GRDIO_NONUNIQUE_FORMAT);
636 }
637 break;
638 case 8: /* 8-byte elements */
639 header->type = GMT_GRID_IS_BD;
640 HH->orig_datatype = GMT_DOUBLE;
641 break;
642 default: /* Garbage */
643 free_and_return (GMT_GRDIO_BAD_VAL);
644 break;
645 }
646 free_and_return (GMT_NOERROR);
647 }
648
gmtcustomio_native_write_grd_header(FILE * fp,struct GMT_GRID_HEADER * header)649 GMT_LOCAL int gmtcustomio_native_write_grd_header (FILE *fp, struct GMT_GRID_HEADER *header) {
650 int err = GMT_NOERROR;
651 /* Because GMT_GRID_HEADER is not 64-bit aligned we must write it in parts */
652
653 if (gmt_M_fwrite (&header->n_columns, SIZEOF_NATIVE_GRD_HDR1, 1U, fp) != 1 ||
654 gmt_M_fwrite (&header->wesn[0], SIZEOF_NATIVE_GRD_HDR2, 1U, fp) != 1)
655 err = GMT_GRDIO_WRITE_FAILED;
656 return (err);
657 }
658
gmtcustomio_native_skip_grd_header(FILE * fp,struct GMT_GRID_HEADER * header)659 GMT_LOCAL int gmtcustomio_native_skip_grd_header (FILE *fp, struct GMT_GRID_HEADER *header) {
660 int err = GMT_NOERROR;
661 gmt_M_unused(header);
662 /* Because GMT_GRID_HEADER is not 64-bit aligned we must estimate the # of bytes in parts */
663
664 if (fseek (fp, SIZEOF_NATIVE_GRD_HDR, SEEK_SET))
665 err = GMT_GRDIO_SEEK_FAILED;
666 return (err);
667 }
668
gmt_bit_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)669 int gmt_bit_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
670 /* header: grid structure header */
671 /* grid: array with final grid */
672 /* wesn: Sub-region to extract [Use entire file if 0,0,0,0] */
673 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
674 /* complex_mode: &1 | &2 if complex array is to hold real (1) and imaginary (2) parts (otherwise read as real only) */
675 /* Note: The file has only real values, we simply allow space in the complex array */
676 /* for real and imaginary parts when processed by grdfft etc. */
677
678 int j, err, bit;
679 bool piping = false, check = false;
680 int first_col, last_col, first_row, last_row;
681 unsigned int i, width_in, height_in, mx, word;
682 unsigned int *actual_col = NULL;
683 uint64_t kk, ij, j2, width_out, imag_offset;
684 FILE *fp = NULL;
685 unsigned int *tmp = NULL, ival;
686 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
687
688 if (!strcmp (HH->name, "=")) { /* Read from pipe */
689 #ifdef SET_IO_MODE
690 gmt_setmode (GMT, GMT_IN);
691 #endif
692 fp = GMT->session.std[GMT_IN];
693 piping = true;
694 }
695 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) != NULL) { /* Skip header */
696 gmt_M_err_trap (gmtcustomio_native_skip_grd_header (fp, header));
697 }
698 else
699 return (GMT_GRDIO_OPEN_FAILED);
700
701 check = !isnan (header->nan_value);
702 mx = urint (ceil (header->n_columns / 32.0)); /* Whole multiple of 32-bit integers */
703
704 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);
705 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
706
707 width_out = width_in; /* Width of output array */
708 if (pad[XLO] > 0) width_out += pad[XLO];
709 if (pad[XHI] > 0) width_out += pad[XHI];
710
711 tmp = gmt_M_memory (GMT, NULL, mx, unsigned int);
712
713 if (piping) { /* Skip data by reading it */
714 for (j = 0; j < first_row; j++) if (gmt_M_fread (tmp, sizeof (unsigned int), mx, fp) < mx) {
715 gmt_M_free (GMT, actual_col);
716 gmt_M_free (GMT, tmp);
717 return (GMT_GRDIO_READ_FAILED);
718 }
719 }
720 else { /* Simply seek by it */
721 if (fseek (fp, (off_t) (first_row * mx * sizeof (unsigned int)), SEEK_CUR)) {
722 gmt_M_free (GMT, actual_col);
723 gmt_M_free (GMT, tmp);
724 gmt_fclose (GMT, fp);
725 return (GMT_GRDIO_SEEK_FAILED);
726 }
727 }
728
729 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
730 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
731 for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
732 if (gmt_M_fread ( tmp, sizeof (unsigned int), mx, fp) < mx) {
733 if (!piping) gmt_fclose (GMT, fp);
734 gmt_M_free (GMT, actual_col);
735 gmt_M_free (GMT, tmp);
736 return (GMT_GRDIO_READ_FAILED); /* Failed to get one row */
737 }
738 ij = imag_offset + (j2 + pad[YHI]) * width_out + pad[XLO];
739 for (i = 0; i < width_in; i++) {
740 kk = ij + i;
741 word = actual_col[i] / 32;
742 bit = actual_col[i] % 32;
743 ival = (tmp[word] >> bit) & 1;
744 grid[kk] = (gmt_grdfloat) ival;
745 if (check && grid[kk] == header->nan_value)
746 grid[kk] = GMT->session.f_NaN;
747 if (gmt_M_is_fnan (grid[kk])) {
748 HH->has_NaNs = GMT_GRID_HAS_NANS;
749 continue;
750 }
751 /* Update z min/max */
752 header->z_min = MIN (header->z_min, (double)grid[kk]);
753 header->z_max = MAX (header->z_max, (double)grid[kk]);
754 }
755 }
756 if (piping) { /* Skip data by reading it */
757 int n_rows = header->n_rows;
758 for (j = last_row + 1; j < n_rows; j++) {
759 if (gmt_M_fread ( tmp, sizeof (unsigned int), mx, fp) < mx) {
760 gmt_M_free (GMT, actual_col);
761 gmt_M_free (GMT, tmp);
762 return (GMT_GRDIO_READ_FAILED);
763 }
764 }
765 }
766
767 header->n_columns = width_in;
768 header->n_rows = height_in;
769 gmt_M_memcpy (header->wesn, wesn, 4, double);
770
771 if (!piping) gmt_fclose (GMT, fp);
772
773 gmt_M_free (GMT, actual_col);
774 gmt_M_free (GMT, tmp);
775 return (GMT_NOERROR);
776 }
777
gmt_bit_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)778 int gmt_bit_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
779 /* header: grid structure header */
780 /* grid: array with final grid */
781 /* wesn: Sub-region to extract [Use entire file if 0,0,0,0] */
782 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
783 /* complex_mode: &1 | &2 if complex array is to hold real (1) and imaginary (2) parts (otherwise read as real only) */
784 /* Note: The file has only real values, we simply allow space in the complex array */
785 /* If 64 is added we write no header*/
786
787 unsigned int i2, iu, ju, width_out, height_out, mx, word, *actual_col = NULL;
788 int first_col, last_col, first_row, last_row;
789 int i, j, bit, err;
790 bool check = false, do_header = true;
791 uint64_t kk, ij, j2, width_in, imag_offset;
792 unsigned int *tmp = NULL, ival;
793 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
794
795 FILE *fp = NULL;
796
797 if (!strcmp (HH->name, "=")) { /* Write to pipe */
798 #ifdef SET_IO_MODE
799 gmt_setmode (GMT, GMT_OUT);
800 #endif
801 fp = GMT->session.std[GMT_OUT];
802 }
803 else if ((fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
804 return (GMT_GRDIO_CREATE_FAILED);
805
806 check = !isnan (header->nan_value);
807
808 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);
809 do_header = gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
810
811 width_in = width_out; /* Physical width of input array */
812 if (pad[XLO] > 0) width_in += pad[XLO];
813 if (pad[XHI] > 0) width_in += pad[XHI];
814
815 gmt_M_memcpy (header->wesn, wesn, 4, double);
816
817 /* Find z_min/z_max */
818
819 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
820 for (j = first_row, j2 = pad[YHI]; j <= last_row; j++, j2++) {
821 for (i = first_col, i2 = pad[XLO]; i <= last_col; i++, i2++) {
822 ij = j2 * width_in + i2 + imag_offset;
823 if (gmt_M_is_fnan (grid[ij])) {
824 if (check)
825 grid[ij] = header->nan_value;
826 }
827 else {
828 ival = (unsigned)lrintf (grid[ij]);
829 if (ival > 1) ival = 1; /* Truncate to 1 */
830 header->z_min = MIN (header->z_min, (double)ival);
831 header->z_max = MAX (header->z_max, (double)ival);
832 }
833 }
834 }
835 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
836 header->z_min = header->z_max = NAN;
837
838 /* Store header information and array */
839
840 if (do_header) {
841 if ((err = gmtcustomio_native_write_grd_header (fp, header))) {
842 gmt_fclose (GMT, fp);
843 gmt_M_free (GMT, actual_col);
844 return err;
845 }
846 }
847
848 mx = urint (ceil (width_out / 32.0));
849 tmp = gmt_M_memory (GMT, NULL, mx, unsigned int);
850
851 i2 = first_col + pad[XLO];
852 for (ju = 0, j2 = first_row + pad[YHI]; ju < height_out; ju++, j2++) {
853 gmt_M_memset (tmp, mx, unsigned int);
854 ij = j2 * width_in + i2 + imag_offset;
855 for (iu = 0; iu < width_out; iu++) {
856 kk = ij + actual_col[iu];
857 word = iu / 32;
858 bit = iu % 32;
859 ival = (unsigned)lrintf (grid[kk]);
860 if (ival > 1) ival = 1; /* Truncate to 1 */
861 tmp[word] |= (ival << bit);
862 }
863 if (gmt_M_fwrite (tmp, sizeof (unsigned int), mx, fp) < mx) {
864 gmt_fclose (GMT, fp);
865 gmt_M_free (GMT, actual_col);
866 gmt_M_free (GMT, tmp);
867 return (GMT_GRDIO_WRITE_FAILED);
868 }
869 }
870
871 gmt_fclose (GMT, fp);
872
873 gmt_M_free (GMT, actual_col);
874 gmt_M_free (GMT, tmp);
875
876 return (GMT_NOERROR);
877 }
878
879 /*-----------------------------------------------------------
880 * Format : bb, bs, bi, bf, bd
881 * Type : Native binary C file
882 * Prefix : GMT_native_
883 * Author : Paul Wessel, SOEST
884 * Date : 17-JUN-1999
885 * Purpose: The native binary output format is used
886 * primarily for piped i/o between programs
887 * that otherwise would use temporary, large
888 * intermediate grdfiles. Note that not all
889 * programs can support piping (they may have
890 * to re-read the file or accept more than one
891 * grdfile). Datasets with limited range may
892 * be stored using 1-, 2-, or 4-byte integers
893 * which will reduce storage space and improve
894 * throughput.
895 *-----------------------------------------------------------*/
896
897 /* GMT 64-bit Modification:
898 *
899 * Read/write GRD header structure from native binary file. This is
900 * used by all the native binary formats in GMT. We isolate the I/O of
901 * the header structure here because of 32/64 bit issues of alignment.
902 * The GRD header is 892 bytes long, three 4-byte integers followed
903 * by ten 8-byte doubles and six character strings. This created a
904 * problem on 64-bit systems, where the GMT_GRID_HEADER structure was
905 * automatically padded with 4-bytes before the doubles. Taking
906 * advantage of the code developed to deal with the 32/64-bit anomaly
907 * (below), we have permanently added a 4-byte integer to the
908 * GMT_GRID_HEADER structure, but skip it when reading or writing the
909 * header.
910 */
911
gmt_native_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)912 int gmt_native_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
913 /* Write GRD header structure to native binary file. This is used by
914 * all the native binary formats in GMT */
915
916 int err;
917 FILE *fp = NULL;
918 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
919
920 if (!strcmp (HH->name, "=")) { /* Write to pipe */
921 #ifdef SET_IO_MODE
922 gmt_setmode (GMT, GMT_OUT);
923 #endif
924 fp = GMT->session.std[GMT_OUT];
925 }
926 else if ((fp = gmt_fopen (GMT, HH->name, "rb+")) == NULL && (fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
927 return (GMT_GRDIO_CREATE_FAILED);
928
929 err = gmtcustomio_native_write_grd_header (fp, header);
930 gmt_fclose (GMT, fp);
931
932 if (err)
933 return err;
934 else
935 return GMT_NOERROR;
936 }
937
gmt_native_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)938 int gmt_native_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
939 /* header: grid structure header */
940 /* grid: array with final grid */
941 /* wesn: Sub-region to extract [Use entire file if 0,0,0,0] */
942 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
943 /* complex_mode: &4 | &8 if complex array is to hold real (1) and imaginary (2) parts (otherwise read as real only) */
944 /* Note: The file has only real values, we simply allow space in the complex array */
945 /* for real and imaginary parts when processed by grdfft etc. */
946
947 int j, type; /* Data type */
948 bool piping = false; /* true if we read input pipe instead of from file */
949 bool check = false; /* true if nan-proxies are used to signify NaN (for non-floating point types) */
950 unsigned int err; /* Error code returned */
951 int first_col, last_col; /* First and last column to deal with */
952 int first_row, last_row; /* First and last row to deal with */
953 unsigned int height_in; /* Number of columns in subregion */
954 unsigned int i; /* Misc. counters */
955 unsigned int *k = NULL; /* Array with indices */
956 unsigned int width_in; /* Number of items in one row of the subregion */
957 uint64_t kk, ij, j2, width_out, imag_offset; /* Width of row as return (may include padding) */
958 size_t size; /* Length of data type */
959 size_t n_expected; /* Length of row to read */
960 FILE *fp = NULL; /* File pointer to data or pipe */
961 void *tmp = NULL; /* Array pointer for reading in rows of data */
962 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
963
964 if (!strcmp (HH->name, "=")) { /* Read from pipe */
965 #ifdef SET_IO_MODE
966 gmt_setmode (GMT, GMT_IN);
967 #endif
968 fp = GMT->session.std[GMT_IN];
969 piping = true;
970 }
971 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) != NULL) { /* Skip header */
972 gmt_M_err_trap (gmtcustomio_native_skip_grd_header (fp, header));
973 }
974 else
975 return (GMT_GRDIO_OPEN_FAILED);
976
977 type = GMT->session.grdformat[header->type][1];
978 size = gmtlib_grd_data_size (GMT, header->type, &header->nan_value);
979 check = !isnan (header->nan_value);
980
981 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
982
983 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_in, &height_in, &first_col, &last_col, &first_row, &last_row, &k), HH->name);
984
985 width_out = width_in; /* Width of output array */
986 if (pad[XLO] > 0) width_out += pad[XLO];
987 if (pad[XHI] > 0) width_out += pad[XHI];
988
989 /* Allocate memory for one row of data (for reading purposes) */
990
991 n_expected = header->n_columns;
992 tmp = gmt_M_memory (GMT, NULL, n_expected * size, char);
993
994 /* Now deal with skipping */
995
996 if (piping) { /* Skip data by reading it */
997 for (j = 0; j < first_row; j++) if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
998 gmt_M_free (GMT, k);
999 gmt_M_free (GMT, tmp);
1000 return (GMT_GRDIO_READ_FAILED);
1001 }
1002 }
1003 else { /* Simply seek over it */
1004 if (fseek (fp, (off_t) (first_row * n_expected * size), SEEK_CUR)) {
1005 gmt_fclose (GMT, fp);
1006 gmt_M_free (GMT, k);
1007 gmt_M_free (GMT, tmp);
1008 return (GMT_GRDIO_SEEK_FAILED);
1009 }
1010 }
1011
1012 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
1013 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
1014 for (j = first_row, j2 = 0; j <= last_row; j++, j2++) {
1015 if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
1016 if (!piping) gmt_fclose (GMT, fp);
1017 gmt_M_free (GMT, k);
1018 gmt_M_free (GMT, tmp);
1019 return (GMT_GRDIO_READ_FAILED); /* Get one row */
1020 }
1021 ij = imag_offset + (j2 + pad[YHI]) * width_out + pad[XLO];
1022 for (i = 0, kk = ij; i < width_in; i++, kk++) {
1023 grid[kk] = gmtlib_decode (GMT, tmp, k[i], type); /* Convert whatever to gmt_grdfloat */
1024 if (check && grid[kk] == header->nan_value)
1025 grid[kk] = GMT->session.f_NaN;
1026 if (gmt_M_is_fnan (grid[kk])) {
1027 HH->has_NaNs = GMT_GRID_HAS_NANS;
1028 continue;
1029 }
1030 /* Update z_min, z_max */
1031 header->z_min = MIN (header->z_min, (double)grid[kk]);
1032 header->z_max = MAX (header->z_max, (double)grid[kk]);
1033 }
1034 }
1035 if (piping) { /* Skip remaining data by reading them */
1036 int n_rows = header->n_rows;
1037 for (j = last_row + 1; j < n_rows; j++) {
1038 if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
1039 gmt_M_free (GMT, k);
1040 gmt_M_free (GMT, tmp);
1041 return (GMT_GRDIO_READ_FAILED);
1042 }
1043 }
1044 }
1045
1046 header->n_columns = width_in;
1047 header->n_rows = height_in;
1048 gmt_M_memcpy (header->wesn, wesn, 4, double);
1049
1050 if (!piping) gmt_fclose (GMT, fp);
1051
1052 gmt_M_free (GMT, k);
1053 gmt_M_free (GMT, tmp);
1054
1055 return (GMT_NOERROR);
1056 }
1057
gmt_native_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)1058 int gmt_native_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
1059 /* header: grid structure header */
1060 /* grid: array with final grid */
1061 /* wesn: Sub-region to write out [Use entire file if 0,0,0,0] */
1062 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
1063 /* complex_mode: &4 | &8 if complex array holds real (4) and imaginary (8) parts (otherwise write real only) */
1064 /* Note: The file has only real values, we simply have allowed for space in the complex array */
1065 /* for real and imaginary parts when processed by grdfft etc. */
1066 /* If 64 is added we write no header */
1067
1068 int err; /* Error code*/
1069 int i, j, type; /* Data type */
1070 bool check = false; /* true if nan-proxies are used to signify NaN (for non-floating point types) */
1071 bool do_header = true; /* true if we should write the header first */
1072 int first_col, last_col; /* First and last column to deal with */
1073 int first_row, last_row; /* First and last row to deal with */
1074 unsigned int width_out; /* Width of row as return (may include padding) */
1075 unsigned int height_out; /* Number of columns in subregion */
1076 unsigned int i2, ju, iu; /* Misc. counters */
1077 unsigned int *k = NULL; /* Array with indices */
1078 uint64_t ij, width_in, imag_offset, j2;
1079 size_t size; /* Length of data type */
1080 size_t n_expected; /* Length of row to read */
1081 FILE *fp = NULL; /* File pointer to data or pipe */
1082 void *tmp = NULL; /* Array pointer for writing in rows of data */
1083 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1084
1085 if (!strcmp (HH->name, "=")) { /* Write to pipe */
1086 #ifdef SET_IO_MODE
1087 gmt_setmode (GMT, GMT_OUT);
1088 #endif
1089 fp = GMT->session.std[GMT_OUT];
1090 }
1091 else if ((fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
1092 return (GMT_GRDIO_CREATE_FAILED);
1093
1094 type = GMT->session.grdformat[header->type][1];
1095 size = gmtlib_grd_data_size (GMT, header->type, &header->nan_value);
1096 check = !isnan (header->nan_value);
1097
1098 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_out, &height_out, &first_col, &last_col, &first_row, &last_row, &k), HH->name);
1099 do_header = gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
1100
1101 width_in = width_out; /* Physical width of input array */
1102 if (pad[XLO] > 0) width_in += pad[XLO];
1103 if (pad[XHI] > 0) width_in += pad[XHI];
1104
1105 gmt_M_memcpy (header->wesn, wesn, 4, double);
1106
1107 /* Find z_min/z_max */
1108
1109 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
1110 for (j = first_row, j2 = pad[YHI]; j <= last_row; j++, j2++) {
1111 for (i = first_col, i2 = pad[XLO]; i <= last_col; i++, i2++) {
1112 ij = imag_offset + (j2 * width_in + i2);
1113 if (gmt_M_is_fnan (grid[ij])) {
1114 if (check) grid[ij] = header->nan_value;
1115 }
1116 else {
1117 header->z_min = MIN (header->z_min, (double)grid[ij]);
1118 header->z_max = MAX (header->z_max, (double)grid[ij]);
1119 }
1120 }
1121 }
1122
1123 /* Round off to chosen type */
1124
1125 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
1126 header->z_min = header->z_max = NAN;
1127 else if (type != 'f' && type != 'd') {
1128 header->z_min = rint (header->z_min);
1129 header->z_max = rint (header->z_max);
1130 }
1131
1132 /* Store header information and array */
1133
1134 if (do_header) {
1135 if ((err = gmtcustomio_native_write_grd_header (fp, header)) != 0) {
1136 gmt_M_free (GMT, k); gmt_fclose (GMT, fp);
1137 return err;
1138 }
1139 }
1140
1141 /* Allocate memory for one row of data (for writing purposes) */
1142
1143 n_expected = header->n_columns;
1144 tmp = gmt_M_memory (GMT, NULL, n_expected * size, char);
1145
1146 i2 = first_col + pad[XLO];
1147 for (ju = 0, j2 = first_row + pad[YHI]; ju < height_out; ju++, j2++) {
1148 ij = imag_offset + j2 * width_in + i2;
1149 for (iu = 0; iu < width_out; iu++) gmtlib_encode (GMT, tmp, iu, grid[ij+k[iu]], type);
1150 if (gmt_M_fwrite (tmp, size, n_expected, fp) < n_expected) {
1151 gmt_M_free (GMT, k);
1152 gmt_M_free (GMT, tmp);
1153 gmt_fclose (GMT, fp);
1154 return (GMT_GRDIO_WRITE_FAILED);
1155 }
1156 }
1157
1158 gmt_M_free (GMT, k);
1159 gmt_M_free (GMT, tmp);
1160
1161 gmt_fclose (GMT, fp);
1162
1163 return (GMT_NOERROR);
1164 }
1165
gmtlib_encode(struct GMT_CTRL * GMT,void * vptr,uint64_t k,gmt_grdfloat z,unsigned int type)1166 void gmtlib_encode (struct GMT_CTRL *GMT, void *vptr, uint64_t k, gmt_grdfloat z, unsigned int type) {
1167 /* Place the z value in the array location of the (type) pointer */
1168 switch (type) {
1169 case 'b':
1170 ((char *)vptr)[k] = (char)lrintf (z);
1171 break;
1172 case 's':
1173 ((short int *)vptr)[k] = (short int)lrintf (z);
1174 break;
1175 case 'i':
1176 case 'm':
1177 ((int *)vptr)[k] = (int)lrintf (z);
1178 break;
1179 case 'f':
1180 ((float *)vptr)[k] = (float)z;
1181 break;
1182 case 'd':
1183 ((double *)vptr)[k] = (double)z;
1184 break;
1185 default:
1186 GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT: Bad call to gmtlib_encode\n");
1187 break;
1188 }
1189 }
1190
gmtlib_decode(struct GMT_CTRL * GMT,void * vptr,uint64_t k,unsigned int type)1191 gmt_grdfloat gmtlib_decode (struct GMT_CTRL *GMT, void *vptr, uint64_t k, unsigned int type) {
1192 /* Retrieve the z value from the array location of the (type) pointer */
1193 gmt_grdfloat fval;
1194
1195 switch (type) {
1196 case 'b':
1197 fval = (gmt_grdfloat)(((char *)vptr)[k]);
1198 break;
1199 case 's':
1200 fval = (gmt_grdfloat)(((short int *)vptr)[k]);
1201 break;
1202 case 'i':
1203 case 'm':
1204 fval = (gmt_grdfloat)(((int *)vptr)[k]);
1205 break;
1206 case 'f':
1207 fval = ((gmt_grdfloat *)vptr)[k];
1208 break;
1209 case 'd':
1210 fval = (gmt_grdfloat)(((double *)vptr)[k]);
1211 break;
1212 default:
1213 GMT_Report (GMT->parent, GMT_MSG_ERROR, "GMT: Bad call to gmtlib_decode\n");
1214 fval = GMT->session.f_NaN;
1215 break;
1216 }
1217
1218 return (fval);
1219 }
1220
1221 /*-----------------------------------------------------------
1222 * Format : sf, sd
1223 * Type : Surfer 6 and 7 (float) file
1224 * Prefix : GMT_srf_
1225 * Author : Joaquim Luis
1226 * Date : 09-SEP-1999
1227 * 27-AUG-2005 Added minimalist support to GS format 7
1228 * Type : Native binary (double) C file
1229 * Purpose: to transform to/from Surfer grid file format
1230 * Functions : gmt_srf_read_grd_info, gmt_srf_write_grd_info,
1231 * gmt_srf_write_grd_info, gmt_srf_read_grd, gmt_srf_write_grd
1232 *-----------------------------------------------------------*/
1233
1234 struct srf_header6 { /* Surfer 6 file header structure */
1235 char id[4]; /* ASCII Binary identifier (DSBB) */
1236 unsigned short int n_columns; /* Number of columns -- NOTE: original definition by GoldenSoft is "short int" */
1237 unsigned short int n_rows; /* Number of rows */
1238 double wesn[4]; /* Min/maximum x/y coordinates */
1239 double z_min; /* Minimum z value */
1240 double z_max; /* Maximum z value */
1241 };
1242
1243 /* The format 7 is rather more complicated. It has now headers that point to "sections"
1244 that may either be also headers or have real data. Besides that, is follows also the
1245 stupidity of storing the grid using doubles (I would bat that there are no more than 0-1
1246 Surfer users that really need to save their grids in doubles). The following header
1247 does not strictly follow the GS format description, but its enough for reading grids
1248 that do not contain break-lines (again my estimate is that it covers (100 - 1e4)% users).
1249
1250 Note: I had significant troubles to be able to read correctly the Surfer 7 format.
1251 In its basic and mostly used form (that is, without break-lines info) what we normally
1252 call a header, can be described by the srf_header7 structure below (but including
1253 the three commented lines). This would make that the header is composed of 2 char[4] and
1254 and 5 ints followed by doubles. The problem was that after the ints the doubles were not
1255 read correctly. It looked like everything was displaced by 4 bytes.
1256 I than found the note about the GMT 64-bit Modification and tried the same trick.
1257 While that worked with gcc compiled code, it crashed with the VC6 compiler.
1258 Since of the first 3 variables, only the first is important to find out which Surfer
1259 grid format we are dealing with, I removed them from the header definition (and jump
1260 12 bytes before reading the header). As a result the header has now one 4 byte char +
1261 trhee 4-bytes ints followed by 8-bytes doubles. With this organization the header is
1262 read correctly by gmtcustomio_read_srfheader7. Needless to say that I don't understand why the
1263 even number of 4-bytes variables before the 8-bytes caused that the doubles we incorrectly read.
1264
1265 Joaquim Luis 08-2005. */
1266
1267 struct srf_header7 { /* Surfer 7 file header structure */
1268 /*char id[4]; * ASCII Binary identifier (DSRB) */
1269 /*int idumb1; * Size of Header in bytes (is == 1) */
1270 /*int idumb2; * Version number of the file format. Currently must be set to 1*/
1271 char id2[4]; /* Tag ID indicating a grid section (GRID) */
1272 int len_g; /* Length in bytes of the grid section (72) */
1273 int n_rows; /* Number of rows */
1274 int n_columns; /* Number of columns */
1275 double x_min; /* Minimum x coordinate */
1276 double y_min; /* Minimum y coordinate */
1277 double x_inc; /* Spacing between columns */
1278 double y_inc; /* Spacing between rows */
1279 double z_min; /* Minimum z value */
1280 double z_max; /* Maximum z value */
1281 double rotation; /* not currently used */
1282 double no_value; /* If GS were cleverer this would be NaN */
1283 char id3[4]; /* Tag ID indicating a data section (DATA) */
1284 int len_d; /* Length in bytes of the DATA section */
1285 };
1286
gmtlib_is_srf_grid(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)1287 int gmtlib_is_srf_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
1288 FILE *fp = NULL;
1289 char id[5];
1290 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1291 if (!strcmp (HH->name, "="))
1292 return (GMT_GRDIO_PIPE_CODECHECK); /* Cannot check on pipes */
1293 if ((fp = gmt_fopen (GMT, HH->name, "rb")) == NULL)
1294 return (GMT_GRDIO_OPEN_FAILED);
1295 if (gmt_M_fread (id, sizeof (char), 4U, fp) < 4U) {
1296 gmt_fclose (GMT, fp);
1297 return (GMT_GRDIO_READ_FAILED);
1298 }
1299 gmt_fclose (GMT, fp);
1300 if (!strncmp (id, "DSBB", 4U))
1301 header->type = GMT_GRID_IS_SF;
1302 else if (!strncmp (id, "DSRB", 4U))
1303 header->type = GMT_GRID_IS_SD;
1304 else
1305 return (GMT_GRDIO_BAD_VAL); /* Neither */
1306 return GMT_NOERROR;
1307 }
1308
gmtcustomio_read_srfheader6(FILE * fp,struct srf_header6 * h)1309 GMT_LOCAL int gmtcustomio_read_srfheader6 (FILE *fp, struct srf_header6 *h) {
1310 /* Reads the header of a Surfer 6 gridfile */
1311 /* if (gmt_M_fread (h, sizeof (struct srf_header6), 1U, fp) < 1U) return (GMT_GRDIO_READ_FAILED); */
1312
1313 /* UPDATE: Because srf_header6 is not 64-bit aligned we must read it in parts */
1314 if (gmt_M_fread (h->id, 4*sizeof (char), 1U, fp) != 1U)
1315 return (GMT_GRDIO_READ_FAILED);
1316 if (gmt_M_fread (&h->n_columns, 2*sizeof (short int), 1U, fp) != 1U)
1317 return (GMT_GRDIO_READ_FAILED);
1318 if (gmt_M_fread (h->wesn, sizeof (double), 4U, fp) != 4U)
1319 return (GMT_GRDIO_READ_FAILED);
1320 if (gmt_M_fread(&h->z_min, 2*sizeof(double), 1U, fp) != 1U)
1321 return (GMT_GRDIO_READ_FAILED);
1322
1323 return GMT_NOERROR;
1324 }
1325
gmtcustomio_read_srfheader7(FILE * fp,struct srf_header7 * h)1326 GMT_LOCAL int gmtcustomio_read_srfheader7 (FILE *fp, struct srf_header7 *h) {
1327 /* Reads the header of a Surfer 7 gridfile */
1328
1329 if (fseek (fp, (off_t)(3*sizeof(int)), SEEK_SET))
1330 return (GMT_GRDIO_SEEK_FAILED); /* skip the first 12 bytes */
1331 /* if (gmt_M_fread (h, sizeof (struct srf_header7), 1U, fp) < 1U) return (GMT_GRDIO_READ_FAILED); */
1332
1333 /* UPDATE: Because srf_header6 is not 64-bit aligned we must read it in parts */
1334 if (gmt_M_fread (h->id2, 4*sizeof (char), 1U, fp) != 1)
1335 return (GMT_GRDIO_READ_FAILED);
1336 if (gmt_M_fread (&h->len_g, 3*sizeof (int), 1U, fp) != 1)
1337 return (GMT_GRDIO_READ_FAILED);
1338 if (gmt_M_fread (&h->x_min, 8*sizeof (double), 1U, fp) != 1)
1339 return (GMT_GRDIO_READ_FAILED);
1340 if (gmt_M_fread (h->id3, 4*sizeof (char), 1U, fp) != 1)
1341 return (GMT_GRDIO_READ_FAILED);
1342 if (gmt_M_fread (&h->len_d, sizeof (int), 1U, fp) != 1)
1343 return (GMT_GRDIO_READ_FAILED);
1344
1345 return GMT_NOERROR;
1346 }
1347
gmtcustomio_write_srfheader(FILE * fp,struct srf_header6 * h)1348 GMT_LOCAL int gmtcustomio_write_srfheader (FILE *fp, struct srf_header6 *h) {
1349 /* if (gmt_M_fwrite (h, sizeof (struct srf_header6), 1U, fp) < 1U) return (GMT_GRDIO_WRITE_FAILED); */
1350 /* UPDATE: Because srf_header6 is not 64-bit aligned we must write it in parts */
1351 if (gmt_M_fwrite (h->id, 4*sizeof (char), 1U, fp) != 1)
1352 return (GMT_GRDIO_WRITE_FAILED);
1353 if (gmt_M_fwrite (&h->n_columns, 2*sizeof (short int), 1U, fp) != 1)
1354 return (GMT_GRDIO_WRITE_FAILED);
1355 if (gmt_M_fwrite (h->wesn, sizeof (struct srf_header6) - ((size_t)h->wesn - (size_t)h->id), 1U, fp) != 1)
1356 return (GMT_GRDIO_WRITE_FAILED);
1357 return GMT_NOERROR;
1358 }
1359
gmt_srf_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)1360 int gmt_srf_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
1361 FILE *fp = NULL;
1362 struct srf_header6 h6;
1363 struct srf_header7 h7;
1364 char id[5];
1365 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1366
1367 if (!strcmp (HH->name, "=")) { /* Read from pipe */
1368 #ifdef SET_IO_MODE
1369 gmt_setmode (GMT, GMT_IN);
1370 #endif
1371 fp = GMT->session.std[GMT_IN];
1372 }
1373 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) == NULL) {
1374 gmt_fclose (GMT, fp);
1375 return (GMT_GRDIO_OPEN_FAILED);
1376 }
1377
1378 if (gmt_M_fread (id, sizeof (char), 4U, fp) < 4U) {
1379 gmt_fclose (GMT, fp);
1380 return (GMT_GRDIO_READ_FAILED);
1381 }
1382 if (fseek(fp, (off_t)0, SEEK_SET)) {
1383 gmt_fclose (GMT, fp);
1384 return (GMT_GRDIO_SEEK_FAILED);
1385 }
1386 if (strncmp (id, "DSBB", 4U) && strncmp (id, "DSRB", 4U)) {
1387 gmt_fclose (GMT, fp);
1388 return (GMT_GRDIO_NOT_SURFER);
1389 }
1390
1391 gmt_M_memset (&h6, 1, struct srf_header6);
1392 gmt_M_memset (&h7, 1, struct srf_header7);
1393 if (!strncmp (id, "DSBB", 4U)) { /* Version 6 format */
1394 if (gmtcustomio_read_srfheader6 (fp, &h6)) return (GMT_GRDIO_READ_FAILED);
1395 header->type = GMT_GRID_IS_SF;
1396 HH->orig_datatype = GMT_FLOAT;
1397 }
1398 else { /* Version 7 format */
1399 if (gmtcustomio_read_srfheader7 (fp, &h7)) return (GMT_GRDIO_READ_FAILED);
1400 if (h7.len_d != h7.n_columns * h7.n_rows * 8 || !strcmp (h7.id2, "GRID")) return (GMT_GRDIO_SURF7_UNSUPPORTED);
1401 header->type = GMT_GRID_IS_SD;
1402 HH->orig_datatype = GMT_DOUBLE;
1403 }
1404
1405 gmt_fclose (GMT, fp);
1406
1407 header->registration = GMT_GRID_NODE_REG; /* Grid node registration */
1408 if (header->type == GMT_GRID_IS_SF) {
1409 strcpy (header->title, "Grid originally in Surfer 6 format");
1410 header->n_columns = h6.n_columns; header->n_rows = h6.n_rows;
1411 header->wesn[XLO] = h6.wesn[XLO]; header->wesn[XHI] = h6.wesn[XHI];
1412 header->wesn[YLO] = h6.wesn[YLO]; header->wesn[YHI] = h6.wesn[YHI];
1413 header->z_min = h6.z_min; header->z_max = h6.z_max;
1414 header->inc[GMT_X] = gmt_M_get_inc (GMT, h6.wesn[XLO], h6.wesn[XHI], h6.n_columns, header->registration);
1415 header->inc[GMT_Y] = gmt_M_get_inc (GMT, h6.wesn[YLO], h6.wesn[YHI], h6.n_rows, header->registration);
1416 }
1417 else { /* Format GMT_GRID_IS_SD */
1418 strcpy (header->title, "Grid originally in Surfer 7 format");
1419 header->n_columns = h7.n_columns; header->n_rows = h7.n_rows;
1420 header->wesn[XLO] = h7.x_min; header->wesn[YLO] = h7.y_min;
1421 header->wesn[XHI] = h7.x_min + h7.x_inc * (h7.n_columns - 1);
1422 header->wesn[YHI] = h7.y_min + h7.y_inc * (h7.n_rows - 1);
1423 header->z_min = h7.z_min; header->z_max = h7.z_max;
1424 header->inc[GMT_X] = h7.x_inc; header->inc[GMT_Y] = h7.y_inc;
1425 }
1426 header->z_scale_factor = 1; header->z_add_offset = 0;
1427
1428 return (GMT_NOERROR);
1429 }
1430
gmt_srf_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)1431 int gmt_srf_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
1432 FILE *fp = NULL;
1433 struct srf_header6 h;
1434 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1435
1436 if (!strcmp (HH->name, "=")) /* Write to pipe */
1437 {
1438 #ifdef SET_IO_MODE
1439 gmt_setmode (GMT, GMT_OUT);
1440 #endif
1441 fp = GMT->session.std[GMT_OUT];
1442 }
1443 else if ((fp = gmt_fopen (GMT, HH->name, "rb+")) == NULL && (fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
1444 return (GMT_GRDIO_CREATE_FAILED);
1445
1446 /* coverity[buffer_size] */ /* For Coverity analysis. Do not remove this comment */
1447 gmt_strncpy (h.id, "DSBB", 4U);
1448 h.n_columns = (short int)header->n_columns; h.n_rows = (short int)header->n_rows;
1449 if (header->registration == GMT_GRID_PIXEL_REG) {
1450 h.wesn[XLO] = header->wesn[XLO] + header->inc[GMT_X]/2.0; h.wesn[XHI] = header->wesn[XHI] - header->inc[GMT_X]/2.0;
1451 h.wesn[YLO] = header->wesn[YLO] + header->inc[GMT_Y]/2.0; h.wesn[YHI] = header->wesn[YHI] - header->inc[GMT_Y]/2.0;
1452 }
1453 else
1454 gmt_M_memcpy (h.wesn, header->wesn, 4, double);
1455 h.z_min = header->z_min; h.z_max = header->z_max;
1456
1457 if (gmtcustomio_write_srfheader (fp, &h)) {
1458 gmt_fclose (GMT, fp);
1459 return (GMT_GRDIO_WRITE_FAILED);
1460 }
1461
1462 gmt_fclose (GMT, fp);
1463
1464 return (GMT_NOERROR);
1465 }
1466
gmt_srf_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)1467 int gmt_srf_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
1468 /* header: grid structure header */
1469 /* grid: array with final grid */
1470 /* w,e,s,n: Sub-region to extract [Use entire file if 0,0,0,0] */
1471 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
1472 /* complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) */
1473 /* Note: The file has only real values, we simply allow space in the complex array */
1474 /* for real and imaginary parts when processed by grdfft etc. */
1475
1476 int j, type, n_rows; /* Data type */
1477 bool piping = false; /* true if we read input pipe instead of from file */
1478 int first_col, last_col; /* First and last column to deal with */
1479 int first_row, last_row; /* First and last row to deal with */
1480 unsigned int width_in; /* Number of items in one row of the subregion */
1481 unsigned int height_in; /* Number of columns in subregion */
1482 unsigned int i; /* Misc. counters */
1483 unsigned int *k = NULL; /* Array with indices */
1484 uint64_t kk, ij, j2, width_out, imag_offset;
1485 size_t size; /* Length of data type */
1486 size_t n_expected; /* Length of a row */
1487 FILE *fp = NULL; /* File pointer to data or pipe */
1488 void *tmp = NULL; /* Array pointer for reading in rows of data */
1489 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1490
1491 /* Fixed nan_value in Surfer grids */
1492 header->nan_value = 1.70141e38f;
1493
1494 if (!strcmp (HH->name, "=")) { /* Read from pipe */
1495 #ifdef SET_IO_MODE
1496 gmt_setmode (GMT, GMT_IN);
1497 #endif
1498 fp = GMT->session.std[GMT_IN];
1499 piping = true;
1500 }
1501 else if ((fp = gmt_fopen (GMT, HH->name, "rb")) != NULL) { /* Skip header */
1502 if (header->type == GMT_GRID_IS_SF) { /* Surfer Version 6 */
1503 if (fseek (fp, (off_t) sizeof (struct srf_header6), SEEK_SET)) return (GMT_GRDIO_SEEK_FAILED);
1504 }
1505 else { /* Version 7 (skip also the first 12 bytes) */
1506 if (fseek (fp, (off_t) (3*sizeof(int) + sizeof (struct srf_header7)), SEEK_SET)) return (GMT_GRDIO_SEEK_FAILED);
1507 }
1508 }
1509 else
1510 return (GMT_GRDIO_OPEN_FAILED);
1511
1512 type = GMT->session.grdformat[header->type][1];
1513 size = gmtlib_grd_data_size (GMT, header->type, &header->nan_value);
1514
1515 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_in, &height_in, &first_col, &last_col, &first_row, &last_row, &k), HH->name);
1516 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
1517
1518 width_out = width_in; /* Width of output array */
1519 if (pad[XLO] > 0) width_out += pad[XLO];
1520 if (pad[XHI] > 0) width_out += pad[XHI];
1521
1522 n_rows = header->n_rows;
1523 if ( (last_row - first_row + 1) != n_rows) { /* We have a sub-region */
1524 /* Surfer grids are stored starting from Lower Left, which is contrary to
1525 the rest of GMT grids that start at Top Left. So we must do a flip here */
1526 first_row = n_rows - height_in - first_row;
1527 last_row = first_row + height_in - 1;
1528 }
1529
1530 /* Allocate memory for one row of data (for reading purposes) */
1531
1532 n_expected = header->n_columns;
1533 tmp = gmt_M_memory (GMT, NULL, n_expected * size, char);
1534
1535 /* Now deal with skipping */
1536
1537 if (piping) { /* Skip data by reading it */
1538 for (j = 0; j < first_row; j++)
1539 if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
1540 gmt_M_free (GMT, k);
1541 gmt_M_free (GMT, tmp);
1542 return (GMT_GRDIO_READ_FAILED);
1543 }
1544 }
1545 else { /* Simply seek over it */
1546 if (first_row && fseek (fp, (off_t) (first_row * n_expected * size), SEEK_CUR)) {
1547 gmt_fclose (GMT, fp);
1548 gmt_M_free (GMT, k);
1549 gmt_M_free (GMT, tmp);
1550 return (GMT_GRDIO_SEEK_FAILED);
1551 }
1552 }
1553
1554 header->z_min = DBL_MAX;
1555 header->z_max = -DBL_MAX;
1556 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
1557
1558 for (j = first_row, j2 = height_in-1; j <= last_row; j++, j2--) {
1559 if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
1560 if (!piping) gmt_fclose (GMT, fp);
1561 gmt_M_free (GMT, k);
1562 gmt_M_free (GMT, tmp);
1563 return (GMT_GRDIO_READ_FAILED); /* Failed to get one row */
1564 }
1565 ij = imag_offset + (j2 + pad[YHI]) * width_out + pad[XLO];
1566 for (i = 0; i < width_in; i++) {
1567 kk = ij + i;
1568 grid[kk] = gmtlib_decode (GMT, tmp, k[i], type); /* Convert whatever to gmt_grdfloat */
1569 if (grid[kk] >= header->nan_value) {
1570 HH->has_NaNs = GMT_GRID_HAS_NANS;
1571 grid[kk] = GMT->session.f_NaN;
1572 }
1573 else { /* Update z_min, z_max */
1574 header->z_min = MIN (header->z_min, (double)grid[kk]);
1575 header->z_max = MAX (header->z_max, (double)grid[kk]);
1576 }
1577 }
1578 }
1579 if (piping) { /* Skip remaining data by reading it */
1580 for (j = last_row + 1; j < n_rows; j++) {
1581 if (gmt_M_fread (tmp, size, n_expected, fp) < n_expected) {
1582 gmt_M_free (GMT, k);
1583 gmt_M_free (GMT, tmp);
1584 return (GMT_GRDIO_READ_FAILED);
1585 }
1586 }
1587 }
1588
1589 header->n_columns = width_in;
1590 header->n_rows = height_in;
1591 gmt_M_memcpy (header->wesn, wesn, 4, double);
1592
1593 if (!piping) gmt_fclose (GMT, fp);
1594
1595 gmt_M_free (GMT, k);
1596 gmt_M_free (GMT, tmp);
1597
1598 return (GMT_NOERROR);
1599 }
1600
gmt_srf_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)1601 int gmt_srf_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
1602 /* header: grid structure header */
1603 /* grid: array with final grid */
1604 /* wesnn: Sub-region to write out [Use entire file if 0,0,0,0] */
1605 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
1606 /* complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) */
1607 /* Note: The file has only real values, we simply allow space in the complex array */
1608 /* for real and imaginary parts when processed by grdfft etc. */
1609
1610 int i, j, type; /* Data type */
1611 int first_col, last_col; /* First and last column to deal with */
1612 int first_row, last_row; /* First and last row to deal with */
1613 unsigned int width_out; /* Width of row as return (may include padding) */
1614 unsigned int height_out; /* Number of columns in subregion */
1615 unsigned int i2, ju, iu; /* Misc. counters */
1616 unsigned int *k = NULL; /* Array with indices */
1617 uint64_t ij, kk, j2, width_in, imag_offset; /* Number of items in one row of the subregion */
1618 size_t size; /* Length of data type */
1619 size_t n_expected; /* Length of a row */
1620 FILE *fp = NULL; /* File pointer to data or pipe */
1621 void *tmp = NULL; /* Array pointer for writing in rows of data */
1622 struct srf_header6 h;
1623 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1624
1625 if (GMT->session.grdformat[header->type][1] == 'd') {
1626 #ifdef HAVE_GDAL
1627 GMT_Report(GMT->parent, GMT_MSG_INFORMATION,
1628 "Surfer 7 format in GMT is read-only but you can do it via GDAL by appending '=gd:GS7BG' to the file name\n");
1629 #else
1630 GMT_Report(GMT->parent, GMT_MSG_INFORMATION,
1631 "As mentioned in the manual, Surfer 7 format in GMT is read-only\n");
1632 #endif
1633 return (GMT_NOERROR);
1634 }
1635
1636 header->nan_value = 1.70141e38f; /* Fixed nan_value in Surfer grids */
1637
1638 if (!strcmp (HH->name, "=")) { /* Write to pipe */
1639 #ifdef SET_IO_MODE
1640 gmt_setmode (GMT, GMT_OUT);
1641 #endif
1642 fp = GMT->session.std[GMT_OUT];
1643 }
1644 else if ((fp = gmt_fopen (GMT, HH->name, "wb")) == NULL)
1645 return (GMT_GRDIO_CREATE_FAILED);
1646
1647 type = GMT->session.grdformat[header->type][1];
1648 size = gmtlib_grd_data_size (GMT, header->type, &header->nan_value);
1649
1650 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_out, &height_out, &first_col, &last_col, &first_row, &last_row, &k), HH->name);
1651 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
1652
1653 width_in = width_out; /* Physical width of input array */
1654 if (pad[XLO] > 0) width_in += pad[XLO];
1655 if (pad[XHI] > 0) width_in += pad[XHI];
1656
1657 gmt_M_memcpy (header->wesn, wesn, 4, double);
1658
1659 /* Find z_min/z_max */
1660
1661 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
1662 for (j = first_row, j2 = pad[YHI]; j <= last_row; j++, j2++) {
1663 ij = imag_offset + j2 * width_in;
1664 for (i = first_col, i2 = pad[XLO]; i <= last_col; i++, i2++) {
1665 kk = ij + i2;
1666 if (isnan (grid[kk]))
1667 grid[kk] = header->nan_value;
1668 else {
1669 header->z_min = MIN (header->z_min, (double)grid[kk]);
1670 header->z_max = MAX (header->z_max, (double)grid[kk]);
1671 }
1672 }
1673 }
1674 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
1675 header->z_min = header->z_max = NAN;
1676
1677 /* store header information and array */
1678
1679 /* coverity[buffer_size] */ /* For Coverity analysis. Do not remove this comment */
1680 gmt_strncpy (h.id, "DSBB", 4U);
1681 h.n_columns = (short int)header->n_columns; h.n_rows = (short int)header->n_rows;
1682 if (header->registration == GMT_GRID_PIXEL_REG) {
1683 h.wesn[XLO] = header->wesn[XLO] + header->inc[GMT_X]/2; h.wesn[XHI] = header->wesn[XHI] - header->inc[GMT_X]/2;
1684 h.wesn[YLO] = header->wesn[YLO] + header->inc[GMT_Y]/2; h.wesn[YHI] = header->wesn[YHI] - header->inc[GMT_Y]/2;
1685 }
1686 else
1687 gmt_M_memcpy (h.wesn, header->wesn, 4, double);
1688
1689 h.z_min = header->z_min; h.z_max = header->z_max;
1690
1691 if (gmt_M_fwrite (&h, sizeof (struct srf_header6), 1U, fp) != 1) {
1692 gmt_M_free (GMT, k);
1693 gmt_fclose (GMT, fp);
1694 return (GMT_GRDIO_WRITE_FAILED);
1695 }
1696
1697 /* Allocate memory for one row of data (for writing purposes) */
1698
1699 n_expected = header->n_columns;
1700 tmp = gmt_M_memory (GMT, NULL, n_expected * size, char);
1701
1702 i2 = first_col + pad[XLO];
1703 for (ju = 0, j2 = last_row + pad[YHI]; ju < height_out; ju++, j2--) {
1704 ij = imag_offset + j2 * width_in + i2;
1705 for (iu = 0; iu < width_out; iu++) gmtlib_encode (GMT, tmp, iu, grid[ij+k[iu]], type);
1706 if (gmt_M_fwrite (tmp, size, n_expected, fp) < n_expected) {
1707 gmt_fclose (GMT, fp);
1708 gmt_M_free (GMT, k);
1709 gmt_M_free (GMT, tmp);
1710 return (GMT_GRDIO_WRITE_FAILED);
1711 }
1712 }
1713
1714 gmt_M_free (GMT, k);
1715 gmt_M_free (GMT, tmp);
1716
1717 gmt_fclose (GMT, fp);
1718
1719 return (GMT_NOERROR);
1720 }
1721
1722 #ifdef HAVE_GDAL
1723 #include "gmt_gdalread.c"
1724 #include "gmt_gdalwrite.c"
1725 #include "gmt_ogrproj.c" /* For coordinate conversions but can "enter" here too */
1726 #if GDAL_VERSION_MAJOR >= 2
1727 #include "gmt_ogrread.c"
1728 #endif
1729 /* GDAL support */
1730 /*-----------------------------------------------------------
1731 * Format : gd
1732 * Type : GDAL compatible format
1733 * Prefix : GMT_gdal_
1734 * Author : Joaquim Luis
1735 * Date : 06-SEP-2009
1736 *
1737 * Purpose: to access data read trough the GDAL interface
1738 * Functions : gmt_gdal_read_grd_info, gmt_gdal_write_grd_info,
1739 * gmt_gdal_write_grd_info, gmt_gdal_read_grd, gmt_gdal_write_grd
1740 *-----------------------------------------------------------*/
1741
gmtcustomio_free_from_gdalread(struct GMT_CTRL * GMT,struct GMT_GDALREAD_OUT_CTRL * from_gdalread)1742 static inline void gmtcustomio_free_from_gdalread (struct GMT_CTRL *GMT, struct GMT_GDALREAD_OUT_CTRL *from_gdalread) {
1743 int i;
1744 gmt_M_free (GMT, from_gdalread->ColorMap);
1745 for (i = 0; i < from_gdalread->RasterCount; i++)
1746 gmt_M_str_free (from_gdalread->band_field_names[i].DataType); /* Those were allocated with strdup */
1747 gmt_M_str_free (from_gdalread->ProjRefPROJ4);
1748 gmt_M_str_free (from_gdalread->ProjRefWKT);
1749 gmt_M_free (GMT, from_gdalread->band_field_names);
1750 }
1751
gmt_gdal_read_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)1752 int gmt_gdal_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
1753 struct GMT_GDALREAD_IN_CTRL *to_gdalread = NULL;
1754 struct GMT_GDALREAD_OUT_CTRL *from_gdalread = NULL;
1755 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1756
1757 if (!strcmp (HH->name, "=")) {
1758 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Pipes cannot be used within the GDAL interface.\n");
1759 return (GMT_GRDIO_OPEN_FAILED);
1760 }
1761
1762 /* Allocate new control structures */
1763 to_gdalread = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALREAD_IN_CTRL);
1764 from_gdalread = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALREAD_OUT_CTRL);
1765
1766 to_gdalread->M.active = true; /* Metadata only */
1767 if (HH->pocket) { /* Have a band request. */
1768 to_gdalread->B.active = true;
1769 to_gdalread->B.bands = HH->pocket; /* Band parsing and error testing is done in gmt_gdalread */
1770 }
1771
1772 if (gmt_gdalread (GMT, HH->name, to_gdalread, from_gdalread)) {
1773 GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR reading file (metadata) with gdalread.\n");
1774 gmt_M_free (GMT, to_gdalread);
1775 gmtcustomio_free_from_gdalread (GMT, from_gdalread);
1776 gmt_M_free (GMT, from_gdalread);
1777 return (GMT_GRDIO_OPEN_FAILED);
1778 }
1779
1780 if (from_gdalread->UInt8.active)
1781 HH->orig_datatype = GMT_UCHAR;
1782 else if (from_gdalread->UInt16.active)
1783 HH->orig_datatype = GMT_USHORT;
1784 else if (from_gdalread->Int16.active)
1785 HH->orig_datatype = GMT_SHORT;
1786 else if (from_gdalread->UInt32.active)
1787 HH->orig_datatype = GMT_UINT;
1788 else if (from_gdalread->Int32.active)
1789 HH->orig_datatype = GMT_INT;
1790 else if (from_gdalread->Float.active)
1791 HH->orig_datatype = GMT_FLOAT;
1792 else if (from_gdalread->Double.active)
1793 HH->orig_datatype = GMT_DOUBLE;
1794 else if (!strcmp(from_gdalread->band_field_names->DataType, "Byte"))
1795 HH->orig_datatype = GMT_UCHAR;
1796
1797 HH->grdtype = gmtlib_get_grdtype (GMT, GMT_IN, header);
1798 header->type = GMT_GRID_IS_GD;
1799 header->registration = (int)from_gdalread->hdr[6]; /* Which registration? */
1800 strcpy (header->title, "Grid imported via GDAL");
1801 header->n_columns = from_gdalread->RasterXsize, header->n_rows = from_gdalread->RasterYsize;
1802 gmt_M_memcpy (header->wesn, from_gdalread->hdr, 4, double);
1803 header->inc[GMT_X] = from_gdalread->hdr[7];
1804 header->inc[GMT_Y] = from_gdalread->hdr[8];
1805 header->z_min = from_gdalread->hdr[4];
1806 header->z_max = from_gdalread->hdr[5];
1807 if (from_gdalread->band_field_names) {
1808 header->z_scale_factor = from_gdalread->band_field_names[0].ScaleOffset[0];
1809 header->z_add_offset = from_gdalread->band_field_names[0].ScaleOffset[1];
1810 header->nan_value = (gmt_grdfloat)from_gdalread->band_field_names[0].nodata;
1811 }
1812 else {
1813 header->z_scale_factor = 1.0;
1814 header->z_add_offset = 0.0;
1815 }
1816
1817 /* Make sure we don't leak due to a previous copy */
1818 gmt_M_str_free (header->ProjRefPROJ4);
1819 gmt_M_str_free (header->ProjRefWKT);
1820 if (from_gdalread->ProjRefPROJ4)
1821 /* Need to strdup because from_gdalread is freed later */
1822 header->ProjRefPROJ4 = strdup (from_gdalread->ProjRefPROJ4);
1823 if (from_gdalread->ProjRefWKT)
1824 header->ProjRefWKT = strdup (from_gdalread->ProjRefWKT);
1825
1826 gmt_M_free (GMT, to_gdalread);
1827 gmtcustomio_free_from_gdalread (GMT, from_gdalread);
1828 gmt_M_free (GMT, from_gdalread);
1829
1830 return (GMT_NOERROR);
1831 }
1832
gmt_gdal_write_grd_info(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)1833 int gmt_gdal_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
1834 gmt_M_unused(GMT); gmt_M_unused(header);
1835 return (GMT_NOERROR);
1836 }
1837
gmt_gdal_read_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)1838 int gmt_gdal_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
1839 /* header: grid structure header */
1840 /* grid: array with final grid */
1841 /* wesn: Sub-region to extract [Use entire file if 0,0,0,0] */
1842 /* padding: # of empty rows/columns to add on w, e, s, n of grid, respectively */
1843 /* complex_mode: &4 | &8 if complex array is to hold real (4) and imaginary (8) parts (otherwise read as real only) */
1844 /* Note: The file has only real values, we simply allow space in the complex array */
1845 /* for real and imaginary parts when processed by grdfft etc. */
1846
1847 struct GMT_GDALREAD_IN_CTRL *to_gdalread = NULL;
1848 struct GMT_GDALREAD_OUT_CTRL *from_gdalread = NULL;
1849 int nBand, subset;
1850 uint64_t i, j, row, col;
1851 char strR[GMT_LEN128] = {""};
1852 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
1853
1854 /* Allocate new control structures */
1855 to_gdalread = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALREAD_IN_CTRL);
1856 from_gdalread = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALREAD_OUT_CTRL);
1857
1858 if (complex_mode & GMT_GRID_IS_COMPLEX_MASK) {
1859 to_gdalread->Z.active = true; /* Force reading into a complex array */
1860 to_gdalread->Z.complex_mode = (int)((complex_mode & GMT_GRID_IS_COMPLEX_MASK) >> 2); /* Gives 0, 1, or 2 */
1861 }
1862
1863 subset = gmt_M_is_subset (GMT, header, wesn); /* We have a Sub-region demand */
1864 if (subset) { /* We have a Sub-region demand */
1865 to_gdalread->R.active = true;
1866 snprintf (strR, GMT_LEN128, "%.10f/%.10f/%.10f/%.10f", wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
1867 to_gdalread->R.region = strR;
1868 to_gdalread->registration.val = header->registration; /* Due to pix-reg only by GDAL we need to inform it about our reg type */
1869 to_gdalread->registration.x_inc = header->inc[GMT_X];
1870 to_gdalread->registration.y_inc = header->inc[GMT_Y];
1871 }
1872
1873 /* This code chunk fixed #173 (r10638) grdpaste of via-GDAL grids but later caused #403 */
1874 if (pad[XLO] > 0 || pad[XHI] > 0 || pad[YLO] > 0 || pad[YHI] > 0) {
1875 to_gdalread->mini_hdr.active = true;
1876 if (pad[XLO] >= header->n_columns - 1) { /* With -1 we account for both grid & pixel registration cases */
1877 to_gdalread->mini_hdr.offset = pad[XLO]; to_gdalread->mini_hdr.side[0] = 'r';
1878 to_gdalread->mini_hdr.mx = header->mx;
1879 if (gmt_M_check_condition (GMT, !header->mx, "Programming error, header.mx not set\n")) {
1880 gmt_M_free (GMT, to_gdalread); gmt_M_free (GMT, from_gdalread);
1881 return (GMT_N_COLS_NOT_SET);
1882 }
1883 }
1884 else if (pad[XHI] >= header->n_columns - 1) {
1885 to_gdalread->mini_hdr.offset = pad[XHI]; to_gdalread->mini_hdr.side[0] = 'l';
1886 to_gdalread->mini_hdr.mx = header->mx;
1887 if (gmt_M_check_condition (GMT, !header->mx, "Programming error, header.mx not set\n")) {
1888 gmt_M_free (GMT, to_gdalread); gmt_M_free (GMT, from_gdalread);
1889 return (GMT_N_COLS_NOT_SET);
1890 }
1891 }
1892 else if (pad[YLO] >= header->n_rows - 1) {
1893 to_gdalread->mini_hdr.offset = pad[YLO]; to_gdalread->mini_hdr.side[0] = 't';
1894 to_gdalread->mini_hdr.my = header->my;
1895 if (gmt_M_check_condition (GMT, !header->my, "Programming error, header.my not set\n")) {
1896 gmt_M_free (GMT, to_gdalread); gmt_M_free (GMT, from_gdalread);
1897 return (GMT_N_ROWS_NOT_SET);
1898 }
1899 }
1900 else if (pad[YHI] >= header->n_rows - 1) {
1901 to_gdalread->mini_hdr.offset = pad[YHI]; to_gdalread->mini_hdr.side[0] = 'b';
1902 to_gdalread->mini_hdr.my = header->my;
1903 if (gmt_M_check_condition (GMT, !header->my, "Programming error, header.my not set\n")) {
1904 gmt_M_free (GMT, to_gdalread); gmt_M_free (GMT, from_gdalread);
1905 return (GMT_N_ROWS_NOT_SET);
1906 }
1907 }
1908 else {
1909 /* Here we assume that all pad[0] ... pad[3] are equal. Otherwise ... */
1910 to_gdalread->mini_hdr.active = false; /* Undo above setting */
1911 to_gdalread->p.active = true;
1912 gmt_M_memcpy (to_gdalread->p.pad, pad, 4U, unsigned int);
1913 }
1914
1915 /* OK, now test if we are under the condition of #403 (very small grids).
1916 If yes, undo the mini_hdr solution ... and hope no more troubles arise. */
1917 if (to_gdalread->mini_hdr.active && (header->n_columns <= 4 || header->n_rows <= 4)) {
1918 to_gdalread->mini_hdr.active = false;
1919 to_gdalread->p.active = true;
1920 gmt_M_memcpy (to_gdalread->p.pad, pad, 4U, unsigned int);
1921 }
1922 }
1923 if (HH->pocket) { /* Have a band request. */
1924 to_gdalread->B.active = true;
1925 to_gdalread->B.bands = HH->pocket; /* Band parsing and error testing is done in gmt_gdalread */
1926 }
1927
1928 /* Tell gmt_gdalread that we already have the memory allocated and send in the *grid pointer */
1929 to_gdalread->f_ptr.active = true;
1930 to_gdalread->f_ptr.grd = grid;
1931
1932 /* If header->nan_value != NaN tell gdalread to replace nan_value by NaN (in floats only) */
1933 to_gdalread->N.nan_value = header->nan_value;
1934
1935 if (gmt_gdalread (GMT, HH->name, to_gdalread, from_gdalread)) {
1936 GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR reading file with gdalread.\n");
1937 gmt_M_free (GMT, to_gdalread); gmt_M_free (GMT, from_gdalread);
1938 return (GMT_GRDIO_OPEN_FAILED);
1939 }
1940
1941 if (to_gdalread->B.active) gmt_M_str_free (HH->pocket); /* It was allocated by strdup. Free it for an eventual reuse. */
1942
1943 if (subset) { /* We had a Sub-region demand */
1944 header->n_columns = from_gdalread->RasterXsize;
1945 header->n_rows = from_gdalread->RasterYsize;
1946 header->nm = gmt_M_grd_get_nm (header); /* Sets the number of actual data items */
1947 gmt_M_memcpy (header->wesn, from_gdalread->hdr, 4, double);
1948 header->z_min = from_gdalread->hdr[4];
1949 header->z_max = from_gdalread->hdr[5];
1950 }
1951
1952 header->registration = (int)from_gdalread->hdr[6]; /* Confirm registration. It may not be the same as reported by read_grd_info */
1953
1954 if (from_gdalread->Float.active) {
1955 if (!to_gdalread->f_ptr.active) /* If the float array was allocated inside _gdalread */
1956 gmt_M_memcpy (grid, from_gdalread->Float.data, header->size, gmt_grdfloat);
1957 else if (!to_gdalread->c_ptr.active && from_gdalread->UInt8.active) { /* If the char array was allocated inside _gdalread */
1958 for (j = 0; j < header->size; j++)
1959 grid[j] = (gmt_grdfloat)from_gdalread->UInt8.data[j];
1960 }
1961 }
1962 else {
1963 /* Convert everything else do float */
1964 nBand = 0; /* Need a solution to RGB or multiband files */
1965 i = nBand * header->size;
1966 if (from_gdalread->UInt8.active)
1967 for (j = 0; j < header->size; j++)
1968 grid[j] = (gmt_grdfloat)from_gdalread->UInt8.data[j+i];
1969 else if (from_gdalread->UInt16.active)
1970 for (j = 0; j < header->size; j++)
1971 grid[j] = (gmt_grdfloat)from_gdalread->UInt16.data[j+i];
1972 else if (from_gdalread->Int16.active)
1973 for (j = 0; j < header->size; j++)
1974 grid[j] = (gmt_grdfloat)from_gdalread->Int16.data[j+i];
1975 else if (from_gdalread->Int32.active)
1976 for (j = 0; j < header->size; j++)
1977 grid[j] = (gmt_grdfloat)from_gdalread->Int32.data[j+i];
1978 else {
1979 GMT_Report (GMT->parent, GMT_MSG_ERROR, "ERROR data type not supported with gdalread in gmt_customio.\n");
1980 gmt_M_free (GMT, to_gdalread);
1981 gmtcustomio_free_from_gdalread (GMT, from_gdalread);
1982 gmt_M_free (GMT, from_gdalread);
1983 return (GMT_GRDIO_OPEN_FAILED);
1984 }
1985 }
1986
1987 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
1988 if (from_gdalread->nodata && !gmt_M_is_dnan (from_gdalread->nodata)) { /* Data has a nodata value */
1989 /* Since all originally integer types were actually converted to float above, we can use
1990 the same test to search for nodata values other than NaN. Note that gmt_galread sets
1991 unknown nodata return by GDAL as NaN, so in those cases this block of code is not executed */
1992
1993 /* Pointer arithmetic solution that should be parallelizable (but those IFs ...) */
1994 if (subset) { /* We had a Sub-region demand so n_rows * n_columns == grid's allocated size */
1995 for (row = 0; row < header->n_rows; row++) {
1996 for (col = 0; col < header->n_columns; col++, grid++) {
1997 if (*grid == (gmt_grdfloat)from_gdalread->nodata) {
1998 *grid = GMT->session.f_NaN;
1999 HH->has_NaNs = GMT_GRID_HAS_NANS;
2000 }
2001 }
2002 }
2003 }
2004 else { /* Full region. We are scanning also the padding zone which is known to have only 0's but never mind */
2005 for (row = 0; row < header->my; row++) {
2006 for (col = 0; col < header->mx; col++, grid++) {
2007 if (*grid == (gmt_grdfloat)from_gdalread->nodata) {
2008 *grid = GMT->session.f_NaN;
2009 HH->has_NaNs = GMT_GRID_HAS_NANS;
2010 }
2011 }
2012 }
2013 }
2014 grid = &grid[0]; /* Put the pointer pointing back to first element in array */
2015 }
2016 header->nan_value = GMT->session.f_NaN;
2017
2018 if (from_gdalread->UInt8.active)
2019 gmt_M_free (GMT, from_gdalread->UInt8.data);
2020 else if (from_gdalread->Float.active && !to_gdalread->f_ptr.active) /* Do not release the *grid pointer */
2021 gmt_M_free (GMT, from_gdalread->Float.data);
2022 else if (from_gdalread->UInt16.active)
2023 gmt_M_free (GMT, from_gdalread->UInt16.data);
2024 else if (from_gdalread->Int16.active)
2025 gmt_M_free (GMT, from_gdalread->Int16.data);
2026 else if (from_gdalread->Int32.active)
2027 gmt_M_free (GMT, from_gdalread->Int32.data);
2028
2029 gmt_M_free (GMT, to_gdalread);
2030 gmtcustomio_free_from_gdalread (GMT, from_gdalread);
2031 gmt_M_free (GMT, from_gdalread);
2032
2033 return (GMT_NOERROR);
2034 }
2035
gmt_gdal_write_grd(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header,gmt_grdfloat * grid,double wesn[],unsigned int * pad,unsigned int complex_mode)2036 int gmt_gdal_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[],
2037 unsigned int *pad, unsigned int complex_mode) {
2038 uint64_t node = 0, ij, imag_offset, imsize;
2039 int first_col, last_col; /* First and last column to deal with */
2040 int first_row, last_row; /* First and last row to deal with */
2041 unsigned int width_out; /* Width of row as return (may include padding) */
2042 unsigned int height_out; /* Number of columns in subregion */
2043 unsigned int *k = NULL; /* Array with indices */
2044 unsigned int row, col;
2045 char driver[16], type[16], *pch;
2046 unsigned char *zu8 = NULL;
2047 short int *zi16 = NULL;
2048 unsigned short int *zu16 = NULL;
2049 int *zi32 = NULL;
2050 unsigned int *zu32 = NULL;
2051 struct GMT_GDALWRITE_CTRL *to_GDALW = NULL;
2052 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
2053 gmt_M_unused(pad);
2054 type[0] = '\0';
2055
2056 if (HH->pocket == NULL) {
2057 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot write with GDAL without knowing which driver to use.\n");
2058 return (GMT_NOERROR);
2059 }
2060
2061 gmt_M_err_pass (GMT, gmt_grd_prep_io (GMT, header, wesn, &width_out, &height_out, &first_col, &last_col, &first_row, &last_row, &k), HH->name);
2062 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
2063
2064 to_GDALW = gmt_M_memory (GMT, NULL, 1, struct GMT_GDALWRITE_CTRL);
2065 if ((pch = strstr(HH->pocket, "+c")) != NULL) { /* If we have a list of +c<options>, trim and save it */
2066 to_GDALW->co_options = strdup (pch);
2067 pch[0] = '\0';
2068 }
2069 sscanf (HH->pocket, "%[^/]/%s", driver, type);
2070 to_GDALW->driver = strdup(driver);
2071 if (header->ProjRefPROJ4) {to_GDALW->P.ProjRefPROJ4 = header->ProjRefPROJ4; to_GDALW->P.active = true;}
2072 if (header->ProjRefWKT) {to_GDALW->P.ProjRefWKT = header->ProjRefWKT; to_GDALW->P.active = true;}
2073 if (header->ProjRefEPSG) to_GDALW->P.ProjRefEPSG = header->ProjRefEPSG; // Not yet used
2074 to_GDALW->flipud = 0;
2075 if (gmt_M_is_geographic (GMT, GMT_IN))
2076 to_GDALW->geog = 1;
2077 else
2078 to_GDALW->geog = 0;
2079 to_GDALW->n_columns = width_out;
2080 to_GDALW->n_rows = height_out;
2081 to_GDALW->nXSizeFull = header->mx;
2082 to_GDALW->n_bands = header->n_bands;
2083 to_GDALW->registration = header->registration;
2084 to_GDALW->pad[0] = header->pad[XLO]; to_GDALW->pad[1] = header->pad[XHI];
2085 to_GDALW->pad[2] = header->pad[YLO]; to_GDALW->pad[3] = header->pad[YHI];
2086 to_GDALW->ULx = wesn[XLO];
2087 to_GDALW->ULy = wesn[YHI];
2088 to_GDALW->x_inc = gmt_M_get_inc (GMT, header->wesn[XLO], header->wesn[XHI], header->n_columns, header->registration);
2089 to_GDALW->y_inc = gmt_M_get_inc (GMT, header->wesn[YLO], header->wesn[YHI], header->n_rows, header->registration);
2090 to_GDALW->nan_value = header->nan_value;
2091 to_GDALW->command = header->command;
2092 to_GDALW->orig_type = HH->orig_datatype;
2093
2094 /* Lazy implementation of nodata value update as it doesn't check and apply on an eventual sub-region on output only */
2095 if (!isnan (header->nan_value)) {
2096 for (ij = 0; ij < header->size; ij++)
2097 if (isnan (grid[ij]))
2098 grid[ij] = header->nan_value;
2099 }
2100
2101 imsize = gmt_M_get_nm (GMT, width_out, height_out);
2102 if (!type[0] || gmt_strlcmp(type, "float32")) {
2103 /* We have to shift the grid pointer in order to use the GDALRasterIO ability to extract a subregion. */
2104 /* See: osgeo-org.1560.n6.nabble.com/gdal-dev-writing-a-subregion-with-GDALRasterIO-td4960500.html */
2105 to_GDALW->data = &grid[2 * header->mx + (header->pad[XLO] + first_col)+imag_offset];
2106 to_GDALW->type = strdup("float32");
2107 gmt_gdalwrite(GMT, HH->name, to_GDALW);
2108 gmt_M_str_free (to_GDALW->driver);
2109 gmt_M_str_free (to_GDALW->type);
2110 gmt_M_free (GMT, to_GDALW);
2111 gmt_M_free (GMT, k);
2112 return (GMT_NOERROR);
2113 }
2114 else if (gmt_strlcmp(type,"u8") || gmt_strlcmp(type,"u08")) {
2115 zu8 = gmt_M_memory(GMT, NULL, imsize, unsigned char);
2116 for (row = first_row; row < height_out; row++)
2117 for (col = first_col, ij = gmt_M_ijp (header, row, 0)+imag_offset; col < width_out; col++, ij++)
2118 zu8[node++] = (unsigned char)grid[ij];
2119
2120 to_GDALW->data = zu8;
2121 to_GDALW->type = strdup("uint8");
2122 }
2123 else if (gmt_strlcmp(type,"i16")) {
2124 zi16 = gmt_M_memory(GMT, NULL, imsize, short int);
2125 for (row = first_row; row < height_out; row++)
2126 for (col = first_col, ij = gmt_M_ijp (header, row, 0)+imag_offset; col < width_out; col++, ij++)
2127 zi16[node++] = (short int)grid[ij];
2128
2129 to_GDALW->data = zi16;
2130 to_GDALW->type = strdup("int16");
2131 }
2132 else if (gmt_strlcmp(type,"u16")) {
2133 zu16 = gmt_M_memory(GMT, NULL, imsize, unsigned short int);
2134 for (row = first_row; row < height_out; row++)
2135 for (col = first_col, ij = gmt_M_ijp (header, row, 0)+imag_offset; col < width_out; col++, ij++)
2136 zu16[node++] = (unsigned short int)grid[ij];
2137
2138 to_GDALW->data = zu16;
2139 to_GDALW->type = strdup("uint16");
2140 }
2141 else if (gmt_strlcmp(type,"i32")) {
2142 zi32 = gmt_M_memory(GMT, NULL, imsize, int);
2143 for (row = first_row; row < height_out; row++)
2144 for (col = first_col, ij = gmt_M_ijp (header, row, 0)+imag_offset; col < width_out; col++, ij++)
2145 zi32[node++] = (int)grid[ij];
2146
2147 to_GDALW->data = zi32;
2148 to_GDALW->type = strdup("int32");
2149 }
2150 else if (gmt_strlcmp(type,"u32")) {
2151 zu32 = gmt_M_memory(GMT, NULL, imsize, unsigned int);
2152 for (row = first_row; row < height_out; row++)
2153 for (col = first_col, ij = gmt_M_ijp (header, row, 0)+imag_offset; col < width_out; col++, ij++)
2154 zu32[node++] = (unsigned int)grid[ij];
2155
2156 to_GDALW->data = zu32;
2157 to_GDALW->type = strdup("uint32");
2158 }
2159 else {
2160 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unknown or unsupported data type code in gmt_customio for writing file with GDAL.\n");
2161 gmt_M_free (GMT, k); gmt_M_free (GMT, to_GDALW->data); gmt_M_str_free (to_GDALW->driver);
2162 gmt_M_str_free (to_GDALW->type); gmt_M_free (GMT, to_GDALW);
2163 return (GMT_GRDIO_OPEN_FAILED);
2164 }
2165
2166 gmt_gdalwrite(GMT, HH->name, to_GDALW);
2167
2168 gmt_M_free (GMT, k);
2169 gmt_M_free (GMT, to_GDALW->data);
2170 gmt_M_str_free (to_GDALW->driver);
2171 gmt_M_str_free (to_GDALW->type);
2172 gmt_M_free (GMT, to_GDALW);
2173 return (GMT_NOERROR);
2174 }
2175
2176 #endif
2177
2178 /* Add custom code here */
2179
2180 /* 12: NOAA NGDC MGG Format */
2181 #include "gmt_mgg_header2.c"
2182
2183 /* 21: Atlantic Geoscience Center format */
2184 #include "gmt_agc_io.c"
2185
2186 /* 23: ESRI Arc/Info ASCII interchange format */
2187 #include "gmt_esri_io.c"
2188
gmtlib_grdio_init(struct GMT_CTRL * GMT)2189 void gmtlib_grdio_init (struct GMT_CTRL *GMT) {
2190 unsigned int id;
2191
2192 /* First element is empty */
2193
2194 id = k_grd_unknown_fmt;
2195 GMT->session.grdformat[id] = "Unknown grid format";
2196 GMT->session.readinfo[id] = &gmt_dummy_grd_info;
2197 GMT->session.updateinfo[id] = &gmt_dummy_grd_info;
2198 GMT->session.writeinfo[id] = &gmt_dummy_grd_info;
2199 GMT->session.readgrd[id] = &gmt_dummy_grd_read;
2200 GMT->session.writegrd[id] = &gmt_dummy_grd_read;
2201
2202 /* FORMAT: GMT netCDF-based (byte) grdio (COARDS compliant) */
2203
2204 id = GMT_GRID_IS_NB;
2205 GMT->session.grdformat[id] = "nb = GMT netCDF format (8-bit integer), " GMT_NC_CONVENTION;
2206 GMT->session.readinfo[id] = &gmt_nc_read_grd_info;
2207 GMT->session.updateinfo[id] = &gmt_nc_update_grd_info;
2208 GMT->session.writeinfo[id] = &gmt_nc_write_grd_info;
2209 GMT->session.readgrd[id] = &gmt_nc_read_grd;
2210 GMT->session.writegrd[id] = &gmt_nc_write_grd;
2211
2212 /* FORMAT: GMT netCDF-based (short) grdio (COARDS compliant) */
2213
2214 id = GMT_GRID_IS_NS;
2215 GMT->session.grdformat[id] = "ns = GMT netCDF format (16-bit integer), " GMT_NC_CONVENTION;
2216 GMT->session.readinfo[id] = &gmt_nc_read_grd_info;
2217 GMT->session.updateinfo[id] = &gmt_nc_update_grd_info;
2218 GMT->session.writeinfo[id] = &gmt_nc_write_grd_info;
2219 GMT->session.readgrd[id] = &gmt_nc_read_grd;
2220 GMT->session.writegrd[id] = &gmt_nc_write_grd;
2221
2222 /* FORMAT: GMT netCDF-based (int) grdio (COARDS compliant) */
2223
2224 id = GMT_GRID_IS_NI;
2225 GMT->session.grdformat[id] = "ni = GMT netCDF format (32-bit integer), " GMT_NC_CONVENTION;
2226 GMT->session.readinfo[id] = &gmt_nc_read_grd_info;
2227 GMT->session.updateinfo[id] = &gmt_nc_update_grd_info;
2228 GMT->session.writeinfo[id] = &gmt_nc_write_grd_info;
2229 GMT->session.readgrd[id] = &gmt_nc_read_grd;
2230 GMT->session.writegrd[id] = &gmt_nc_write_grd;
2231
2232 /* FORMAT: GMT netCDF-based (float) grdio (COARDS compliant) */
2233
2234 id = GMT_GRID_IS_NF;
2235 GMT->session.grdformat[id] = "nf = GMT netCDF format (32-bit float), " GMT_NC_CONVENTION;
2236 GMT->session.readinfo[id] = &gmt_nc_read_grd_info;
2237 GMT->session.updateinfo[id] = &gmt_nc_update_grd_info;
2238 GMT->session.writeinfo[id] = &gmt_nc_write_grd_info;
2239 GMT->session.readgrd[id] = &gmt_nc_read_grd;
2240 GMT->session.writegrd[id] = &gmt_nc_write_grd;
2241
2242 /* FORMAT: GMT netCDF-based (double) grdio (COARDS compliant) */
2243
2244 id = GMT_GRID_IS_ND;
2245 GMT->session.grdformat[id] = "nd = GMT netCDF format (64-bit float), " GMT_NC_CONVENTION;
2246 GMT->session.readinfo[id] = &gmt_nc_read_grd_info;
2247 GMT->session.updateinfo[id] = &gmt_nc_update_grd_info;
2248 GMT->session.writeinfo[id] = &gmt_nc_write_grd_info;
2249 GMT->session.readgrd[id] = &gmt_nc_read_grd;
2250 GMT->session.writegrd[id] = &gmt_nc_write_grd;
2251
2252 /* FORMAT: GMT netCDF-based (byte) grdio */
2253
2254 id = GMT_GRID_IS_CB;
2255 GMT->session.grdformat[id] = "cb = GMT netCDF format (8-bit integer, deprecated)";
2256 GMT->session.readinfo[id] = &gmt_cdf_read_grd_info;
2257 GMT->session.updateinfo[id] = &gmt_cdf_update_grd_info;
2258 GMT->session.writeinfo[id] = &gmt_cdf_write_grd_info;
2259 GMT->session.readgrd[id] = &gmt_cdf_read_grd;
2260 GMT->session.writegrd[id] = &gmt_cdf_write_grd;
2261
2262 /* FORMAT: GMT netCDF-based (short) grdio */
2263
2264 id = GMT_GRID_IS_CS;
2265 GMT->session.grdformat[id] = "cs = GMT netCDF format (16-bit integer, deprecated)";
2266 GMT->session.readinfo[id] = &gmt_cdf_read_grd_info;
2267 GMT->session.updateinfo[id] = &gmt_cdf_update_grd_info;
2268 GMT->session.writeinfo[id] = &gmt_cdf_write_grd_info;
2269 GMT->session.readgrd[id] = &gmt_cdf_read_grd;
2270 GMT->session.writegrd[id] = &gmt_cdf_write_grd;
2271
2272 /* FORMAT: GMT netCDF-based (int) grdio */
2273
2274 id = GMT_GRID_IS_CI;
2275 GMT->session.grdformat[id] = "ci = GMT netCDF format (32-bit integer, deprecated)";
2276 GMT->session.readinfo[id] = &gmt_cdf_read_grd_info;
2277 GMT->session.updateinfo[id] = &gmt_cdf_update_grd_info;
2278 GMT->session.writeinfo[id] = &gmt_cdf_write_grd_info;
2279 GMT->session.readgrd[id] = &gmt_cdf_read_grd;
2280 GMT->session.writegrd[id] = &gmt_cdf_write_grd;
2281
2282 /* FORMAT: GMT netCDF-based (float) grdio */
2283
2284 id = GMT_GRID_IS_CF;
2285 GMT->session.grdformat[id] = "cf = GMT netCDF format (32-bit float, deprecated)";
2286 GMT->session.readinfo[id] = &gmt_cdf_read_grd_info;
2287 GMT->session.updateinfo[id] = &gmt_cdf_update_grd_info;
2288 GMT->session.writeinfo[id] = &gmt_cdf_write_grd_info;
2289 GMT->session.readgrd[id] = &gmt_cdf_read_grd;
2290 GMT->session.writegrd[id] = &gmt_cdf_write_grd;
2291
2292 /* FORMAT: GMT netCDF-based (double) grdio */
2293
2294 id = GMT_GRID_IS_CD;
2295 GMT->session.grdformat[id] = "cd = GMT netCDF format (64-bit float, deprecated)";
2296 GMT->session.readinfo[id] = &gmt_cdf_read_grd_info;
2297 GMT->session.updateinfo[id] = &gmt_cdf_update_grd_info;
2298 GMT->session.writeinfo[id] = &gmt_cdf_write_grd_info;
2299 GMT->session.readgrd[id] = &gmt_cdf_read_grd;
2300 GMT->session.writegrd[id] = &gmt_cdf_write_grd;
2301
2302 /* FORMAT: GMT native binary (bit) grdio */
2303
2304 id = GMT_GRID_IS_BM;
2305 GMT->session.grdformat[id] = "bm = GMT native, C-binary format (bit-mask)";
2306 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2307 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2308 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2309 GMT->session.readgrd[id] = &gmt_bit_read_grd;
2310 GMT->session.writegrd[id] = &gmt_bit_write_grd;
2311
2312 /* FORMAT: GMT native binary (byte) grdio */
2313
2314 id = GMT_GRID_IS_BB;
2315 GMT->session.grdformat[id] = "bb = GMT native, C-binary format (8-bit integer)";
2316 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2317 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2318 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2319 GMT->session.readgrd[id] = &gmt_native_read_grd;
2320 GMT->session.writegrd[id] = &gmt_native_write_grd;
2321
2322 /* FORMAT: GMT native binary (short) grdio */
2323
2324 id = GMT_GRID_IS_BS;
2325 GMT->session.grdformat[id] = "bs = GMT native, C-binary format (16-bit integer)";
2326 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2327 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2328 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2329 GMT->session.readgrd[id] = &gmt_native_read_grd;
2330 GMT->session.writegrd[id] = &gmt_native_write_grd;
2331
2332 /* FORMAT: GMT native binary (int) grdio */
2333
2334 id = GMT_GRID_IS_BI;
2335 GMT->session.grdformat[id] = "bi = GMT native, C-binary format (32-bit integer)";
2336 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2337 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2338 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2339 GMT->session.readgrd[id] = &gmt_native_read_grd;
2340 GMT->session.writegrd[id] = &gmt_native_write_grd;
2341
2342 /* FORMAT: GMT native binary (float) grdio */
2343
2344 id = GMT_GRID_IS_BF;
2345 GMT->session.grdformat[id] = "bf = GMT native, C-binary format (32-bit float)";
2346 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2347 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2348 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2349 GMT->session.readgrd[id] = &gmt_native_read_grd;
2350 GMT->session.writegrd[id] = &gmt_native_write_grd;
2351
2352 /* FORMAT: GMT native binary (double) grdio */
2353
2354 id = GMT_GRID_IS_BD;
2355 GMT->session.grdformat[id] = "bd = GMT native, C-binary format (64-bit float)";
2356 GMT->session.readinfo[id] = &gmt_native_read_grd_info;
2357 GMT->session.updateinfo[id] = &gmt_native_write_grd_info;
2358 GMT->session.writeinfo[id] = &gmt_native_write_grd_info;
2359 GMT->session.readgrd[id] = &gmt_native_read_grd;
2360 GMT->session.writegrd[id] = &gmt_native_write_grd;
2361
2362 /* FORMAT: SUN 8-bit standard rasterfile grdio */
2363
2364 id = GMT_GRID_IS_RB;
2365 GMT->session.grdformat[id] = "rb = SUN rasterfile format (8-bit standard)";
2366 GMT->session.readinfo[id] = &gmt_ras_read_grd_info;
2367 GMT->session.updateinfo[id] = &gmt_ras_write_grd_info;
2368 GMT->session.writeinfo[id] = &gmt_ras_write_grd_info;
2369 GMT->session.readgrd[id] = &gmt_ras_read_grd;
2370 GMT->session.writegrd[id] = &gmt_ras_write_grd;
2371
2372 /* FORMAT: NOAA NGDC MGG grid format */
2373
2374 id = GMT_GRID_IS_RF;
2375 GMT->session.grdformat[id] = "rf = GEODAS grid format GRD98 (NGDC)";
2376 GMT->session.readinfo[id] = &gmt_mgg2_read_grd_info;
2377 GMT->session.updateinfo[id] = &gmt_mgg2_write_grd_info;
2378 GMT->session.writeinfo[id] = &gmt_mgg2_write_grd_info;
2379 GMT->session.readgrd[id] = &gmt_mgg2_read_grd;
2380 GMT->session.writegrd[id] = &gmt_mgg2_write_grd;
2381
2382 /* FORMAT: GMT native binary (float) grdio (Surfer format) */
2383
2384 id = GMT_GRID_IS_SF;
2385 GMT->session.grdformat[id] = "sf = Golden Software Surfer format 6 (32-bit float)";
2386 GMT->session.readinfo[id] = &gmt_srf_read_grd_info;
2387 GMT->session.updateinfo[id] = &gmt_srf_write_grd_info;
2388 GMT->session.writeinfo[id] = &gmt_srf_write_grd_info;
2389 GMT->session.readgrd[id] = &gmt_srf_read_grd;
2390 GMT->session.writegrd[id] = &gmt_srf_write_grd;
2391
2392 /* FORMAT: GMT native binary (double) grdio (Surfer format) */
2393
2394 id = GMT_GRID_IS_SD;
2395 GMT->session.grdformat[id] = "sd = Golden Software Surfer format 7 (64-bit float, read-only)";
2396 GMT->session.readinfo[id] = &gmt_srf_read_grd_info;
2397 GMT->session.updateinfo[id] = &gmt_srf_write_grd_info;
2398 GMT->session.writeinfo[id] = &gmt_srf_write_grd_info;
2399 GMT->session.readgrd[id] = &gmt_srf_read_grd;
2400 GMT->session.writegrd[id] = &gmt_srf_write_grd;
2401
2402 /* FORMAT: GMT native binary (float) grdio (AGC format) */
2403
2404 id = GMT_GRID_IS_AF;
2405 GMT->session.grdformat[id] = "af = Atlantic Geoscience Center format AGC (32-bit float)";
2406 GMT->session.readinfo[id] = &gmt_agc_read_grd_info;
2407 GMT->session.updateinfo[id] = &gmt_agc_write_grd_info;
2408 GMT->session.writeinfo[id] = &gmt_agc_write_grd_info;
2409 GMT->session.readgrd[id] = &gmt_agc_read_grd;
2410 GMT->session.writegrd[id] = &gmt_agc_write_grd;
2411
2412 /* FORMAT: ESRI Arc/Info ASCII Interchange Grid format (integer) */
2413
2414 id = GMT_GRID_IS_EI;
2415 GMT->session.grdformat[id] = "ei = ESRI Arc/Info ASCII Grid Interchange format (ASCII integer)";
2416 GMT->session.readinfo[id] = &gmt_esri_read_grd_info;
2417 GMT->session.updateinfo[id] = &gmt_esri_write_grd_info;
2418 GMT->session.writeinfo[id] = &gmt_esri_write_grd_info;
2419 GMT->session.readgrd[id] = &gmt_esri_read_grd;
2420 GMT->session.writegrd[id] = &gmt_esri_writei_grd;
2421
2422 /* FORMAT: ESRI Arc/Info ASCII Interchange Grid format (float) */
2423
2424 id = GMT_GRID_IS_EF;
2425 GMT->session.grdformat[id] = "ef = ESRI Arc/Info ASCII Grid Interchange format (ASCII float)";
2426 GMT->session.readinfo[id] = &gmt_esri_read_grd_info;
2427 GMT->session.updateinfo[id] = &gmt_esri_write_grd_info;
2428 GMT->session.writeinfo[id] = &gmt_esri_write_grd_info;
2429 GMT->session.readgrd[id] = &gmt_esri_read_grd;
2430 GMT->session.writegrd[id] = &gmt_esri_writef_grd;
2431
2432 /* FORMAT: Import via the GDAL interface */
2433
2434 id = GMT_GRID_IS_GD;
2435 #ifdef HAVE_GDAL
2436 GMT->session.grdformat[id] = "gd = Import/export through GDAL";
2437 GMT->session.readinfo[id] = &gmt_gdal_read_grd_info;
2438 GMT->session.updateinfo[id] = &gmt_gdal_write_grd_info;
2439 GMT->session.writeinfo[id] = &gmt_gdal_write_grd_info;
2440 GMT->session.readgrd[id] = &gmt_gdal_read_grd;
2441 GMT->session.writegrd[id] = &gmt_gdal_write_grd;
2442 #else
2443 GMT->session.grdformat[id] = "gd = Import/export through GDAL (not supported)";
2444 GMT->session.readinfo[id] = &gmt_dummy_grd_info;
2445 GMT->session.updateinfo[id] = &gmt_dummy_grd_info;
2446 GMT->session.writeinfo[id] = &gmt_dummy_grd_info;
2447 GMT->session.readgrd[id] = &gmt_dummy_grd_read;
2448 GMT->session.writegrd[id] = &gmt_dummy_grd_read;
2449 #endif
2450
2451 /* ----------------------------------------------
2452 * ADD CUSTOM FORMATS BELOW AS THEY ARE NEEDED */
2453 }
2454
2455 /* Comparator for qsort in gmt_grdformats_sorted() */
gmtcustomio_compare_grd_fmt_strings(const void * a,const void * b)2456 GMT_LOCAL int gmtcustomio_compare_grd_fmt_strings (const void* a, const void* b) {
2457 const char **ia = (const char **)a;
2458 const char **ib = (const char **)b;
2459 return strncmp(*ia, *ib, 2);
2460 }
2461
2462 /* This function sorts the grid formats alphabetically by their type id */
gmt_grdformats_sorted(struct GMT_CTRL * Ctrl)2463 char ** gmt_grdformats_sorted (struct GMT_CTRL *Ctrl) {
2464 static char *formats_sorted[GMT_N_GRD_FORMATS];
2465 static bool sorted = false;
2466
2467 if (sorted)
2468 return formats_sorted;
2469
2470 /* copy array with char pointers to type id strings: */
2471 memcpy (formats_sorted, Ctrl->session.grdformat, GMT_N_GRD_FORMATS * sizeof(char*));
2472 /* sort pointers beginning from the 2nd element: */
2473 qsort (formats_sorted + 1, GMT_N_GRD_FORMATS - 1, sizeof(char*), &gmtcustomio_compare_grd_fmt_strings);
2474 sorted = true;
2475
2476 return formats_sorted;
2477 }
2478