1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17 /* Contains the read/write functions needed to handle ERSI Arc/Info ASCII Exchange grids.
18 * Based on previous code in grd2xyz and xyz2grd. A few limitations of the format:
19 * 1) Grid spacing must be square (dx = dy) since only one cell_size record.
20 * 2) NaNs must be stored via a proxy value [Auto-defaults to -9999 if not set].
21 *
22 * Paul Wessel, June 2010.
23 *
24 * 3) Read also the so called ESRI .HDR format (a binary raw file plus a companion header)
25 * 4) Recognizes GTOPO30 and SRTM30 files from their names and use info coded in the
26 * name to fill the header struct.
27 *
28 * Joaquim Luis, Mars 2011.
29 */
30
31 #include "gmt_common_byteswap.h"
32
33 /* Private function visible only in this file */
34
gmtesriio_write_info(struct GMT_CTRL * GMT,FILE * fp,struct GMT_GRID_HEADER * header)35 GMT_LOCAL int gmtesriio_write_info (struct GMT_CTRL *GMT, FILE *fp, struct GMT_GRID_HEADER *header) {
36 /* Note: fp is an open file pointer passed in; it will be closed upstream and not here */
37 char record[GMT_BUFSIZ] = {""}, item[GMT_LEN64] = {""};
38
39 snprintf (record, GMT_BUFSIZ, "ncols %d\nnrows %d\n", header->n_columns, header->n_rows);
40 gmt_M_fputs (record, fp); /* Write a text record */
41 if (header->registration == GMT_GRID_PIXEL_REG) { /* Pixel format */
42 sprintf (record, "xllcorner ");
43 snprintf (item, GMT_LEN64, GMT->current.setting.format_float_out, header->wesn[XLO]);
44 strcat (record, item); strcat (record, "\n");
45 gmt_M_fputs (record, fp); /* Write a text record */
46 sprintf (record, "yllcorner ");
47 snprintf (item, GMT_LEN64, GMT->current.setting.format_float_out, header->wesn[YLO]);
48 strcat (record, item); strcat (record, "\n");
49 gmt_M_fputs (record, fp); /* Write a text record */
50 }
51 else { /* Gridline format */
52 sprintf (record, "xllcenter ");
53 snprintf (item, GMT_LEN64, GMT->current.setting.format_float_out, header->wesn[XLO]);
54 strcat (record, item); strcat (record, "\n");
55 gmt_M_fputs (record, fp); /* Write a text record */
56 sprintf (record, "yllcenter ");
57 snprintf (item, GMT_LEN64, GMT->current.setting.format_float_out, header->wesn[YLO]);
58 strcat (record, item); strcat (record, "\n");
59 gmt_M_fputs (record, fp); /* Write a text record */
60 }
61 sprintf (record, "cellsize ");
62 snprintf (item, GMT_LEN64, GMT->current.setting.format_float_out, header->inc[GMT_X]);
63 strcat (record, item); strcat (record, "\n");
64 gmt_M_fputs (record, fp); /* Write a text record */
65 if (isnan (header->nan_value)) {
66 GMT_Report (GMT->parent, GMT_MSG_WARNING, "ESRI Arc/Info ASCII Interchange file must use proxy for NaN; default to -9999\n");
67 header->nan_value = -9999.0f;
68 }
69 snprintf (record, GMT_BUFSIZ, "nodata_value %ld\n", lrintf (header->nan_value));
70 gmt_M_fputs (record, fp); /* Write a text record */
71
72 return (GMT_NOERROR);
73 }
74
gmtesriio_read_info_hdr(struct GMT_CTRL * GMT,struct GMT_GRID_HEADER * header)75 GMT_LOCAL int gmtesriio_read_info_hdr (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
76 /* Parse the contents of a .HDR file */
77 int nB, error;
78 char record[GMT_BUFSIZ];
79 FILE *fp = NULL;
80 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
81
82 if ((fp = gmt_fopen (GMT, header->title, "r")) == NULL) return (GMT_GRDIO_OPEN_FAILED);
83
84 header->registration = GMT_GRID_NODE_REG;
85 header->z_scale_factor = 1.0;
86 header->z_add_offset = 0.0;
87
88 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* BYTEORDER */
89 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* LAYOUT */
90 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
91 if (sscanf (record, "%*s %d", &header->n_rows) != 1) {
92 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding NROWS record\n");
93 gmt_fclose (GMT, fp);
94 return (GMT_GRDIO_READ_FAILED);
95 }
96 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
97 if (sscanf (record, "%*s %d", &header->n_columns) != 1) {
98 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding NCOLS record\n");
99 gmt_fclose (GMT, fp);
100 return (GMT_GRDIO_READ_FAILED);
101 }
102 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
103 if (sscanf (record, "%*s %d", &nB) != 1) {
104 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding NBANDS record\n");
105 gmt_fclose (GMT, fp);
106 return (GMT_GRDIO_READ_FAILED);
107 }
108 if (nB != 1) {
109 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Cannot read file with number of Bands != 1 \n");
110 gmt_fclose (GMT, fp);
111 return (GMT_GRDIO_READ_FAILED);
112 }
113 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
114 if (sscanf (record, "%*s %d", &header->bits) != 1) {
115 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding NBITS record\n");
116 gmt_fclose (GMT, fp);
117 return (GMT_GRDIO_READ_FAILED);
118 }
119 if ( header->bits != 16 && header->bits != 32 ) {
120 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: This data type (%d bits) is not supported\n", header->bits);
121 gmt_fclose (GMT, fp);
122 return (GMT_GRDIO_READ_FAILED);
123 }
124 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* BANDROWBYTES */
125 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* TOTALROWBYTES */
126 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* BANDGAPBYTES */
127 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
128 while (strncmp (record, "NODATA", 6) ) /* Keep reading till find this keyword */
129 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
130 #ifdef DOUBLE_PRECISION_GRID
131 if (sscanf (record, "%*s %lf", &header->nan_value) != 1) {
132 #else
133 if (sscanf (record, "%*s %f", &header->nan_value) != 1) {
134 #endif
135 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding nan_value_value record\n");
136 gmt_fclose (GMT, fp);
137 return (GMT_GRDIO_READ_FAILED);
138 }
139 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
140 if (sscanf (record, "%*s %lf", &header->wesn[XLO]) != 1) {
141 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding ULXMAP record\n");
142 gmt_fclose (GMT, fp);
143 return (GMT_GRDIO_READ_FAILED);
144 }
145 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
146 if (sscanf (record, "%*s %lf", &header->wesn[YHI]) != 1) {
147 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding ULYMAP record\n");
148 gmt_fclose (GMT, fp);
149 return (GMT_GRDIO_READ_FAILED);
150 }
151 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
152 if (sscanf (record, "%*s %lf", &header->inc[GMT_X]) != 1) {
153 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding XDIM record\n");
154 gmt_fclose (GMT, fp);
155 return (GMT_GRDIO_READ_FAILED);
156 }
157 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
158 if (sscanf (record, "%*s %lf", &header->inc[GMT_Y]) != 1) {
159 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding YDIM record\n");
160 gmt_fclose (GMT, fp);
161 return (GMT_GRDIO_READ_FAILED);
162 }
163
164 gmt_fclose (GMT, fp);
165 HH->orig_datatype = (header->bits == 16) ? GMT_SHORT : GMT_INT;
166
167 header->wesn[XHI] = header->wesn[XLO] + (header->n_columns - 1 + header->registration) * header->inc[GMT_X];
168 header->wesn[YLO] = header->wesn[YHI] - (header->n_rows - 1 + header->registration) * header->inc[GMT_Y];
169
170 if ((error = gmt_M_err_fail (GMT, gmt_grd_RI_verify (GMT, header, 1), HH->name)))
171 return error;
172
173 return (GMT_NOERROR);
174 }
175
176 GMT_LOCAL int gmtesriio_read_info (struct GMT_CTRL *GMT, FILE *fp, struct GMT_GRID_HEADER *header) {
177 /* Note: fp is an open file pointer passed in; it will be closed upstream and not here */
178 int c, error;
179 char record[GMT_BUFSIZ];
180 FILE *fp2 = NULL, *fpBAK = NULL;
181 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
182
183 header->registration = GMT_GRID_NODE_REG;
184 header->z_scale_factor = 1.0;
185 header->z_add_offset = 0.0;
186
187 if (HH->flags[0] == 'M' || HH->flags[0] == 'I') { /* We are dealing with a ESRI .hdr file */
188 int error;
189 if ((error = gmtesriio_read_info_hdr (GMT, header)) != 0) /* Continue the work someplace else */
190 return (error);
191 else
192 return (GMT_NOERROR);
193 }
194 else if (HH->flags[0] == 'B' && HH->flags[1] == '0') { /* A GTOPO30 or SRTM30 file */
195 size_t len = strlen (header->title);
196
197 header->inc[GMT_X] = header->inc[GMT_Y] = 30.0 * GMT_SEC2DEG; /* 30 arc seconds */
198 header->wesn[YHI] = atof (&header->title[len-2]);
199 if ( header->title[len-3] == 'S' || header->title[len-3] == 's' ) header->wesn[YHI] *= -1;
200 c = header->title[len-3];
201 header->title[len-3] = '\0';
202 header->wesn[XLO] = atof (&header->title[len-6]);
203 header->title[len-3] = (char)c; /* Reset because this function is called at least twice */
204 if ( header->title[len-7] == 'W' || header->title[len-7] == 'w' ) header->wesn[XLO] *= -1;
205 if (header->wesn[YHI] > -60) {
206 header->wesn[YLO] = header->wesn[YHI] - 50;
207 header->wesn[XHI] = header->wesn[XLO] + 40;
208 header->n_columns = 4800;
209 header->n_rows = 6000;
210 }
211 else { /* Antarctica tiles cover 30 degrees of latitude and 60 degrees of longitude each have 3,600 rows and 7,200 columns */
212 header->wesn[YLO] = -90;
213 header->wesn[XHI] = header->wesn[XLO] + 60;
214 header->n_columns = 7200;
215 header->n_rows = 3600;
216 }
217 header->registration = GMT_GRID_PIXEL_REG;
218 gmt_set_geographic (GMT, GMT_IN);
219 gmtlib_grd_set_units (GMT, header);
220
221 /* Different sign of NaN value between GTOPO30 and SRTM30 grids */
222 if (strstr (HH->name, ".DEM") || strstr (HH->name, ".dem"))
223 header->nan_value = -9999.0f;
224 else
225 header->nan_value = 9999.0f;
226 header->bits = 16; /* Temp pocket to store number of bits */
227 HH->orig_datatype = GMT_SHORT;
228 return (GMT_NOERROR);
229 }
230 else if (HH->flags[0] == 'B' && HH->flags[1] == '1') { /* A SRTM3 or SRTM1 file */
231 size_t len = strlen (header->title);
232 struct stat F;
233
234 header->wesn[XLO] = atof (&header->title[len-3]);
235 if ( header->title[len-4] == 'W' || header->title[len-4] == 'w' ) header->wesn[XLO] *= -1;
236 c = header->title[len-4];
237 header->title[len-4] = '\0';
238 header->wesn[YLO] = atof (&header->title[len-6]);
239 header->title[len-4] = (char)c; /* Reset because this function is called at least twice */
240 if ( header->title[len-7] == 'S' || header->title[len-7] == 's' ) header->wesn[YLO] *= -1;
241 header->wesn[YHI] = header->wesn[YLO] + 1;
242 header->wesn[XHI] = header->wesn[XLO] + 1;
243 header->nan_value = -32768.0f;
244 header->bits = 16; /* Temp pocket to store number of bits */
245 HH->orig_datatype = GMT_SHORT;
246 if (stat (HH->name, &F)) /* Must finally find out if it is a 1 or 3 arcseconds file */
247 return (GMT_GRDIO_STAT_FAILED); /* Inquiry about file failed somehow */
248 if (F.st_size < 3e6) { /* Actually the true size is 2884802 */
249 header->inc[GMT_X] = header->inc[GMT_Y] = 3.0 * GMT_SEC2DEG; /* 3 arc seconds */
250 strcpy (header->remark, "Assumed to be a SRTM3 tile");
251 }
252 else {
253 header->inc[GMT_X] = header->inc[GMT_Y] = 1.0 * GMT_SEC2DEG; /* 1 arc second */
254 strcpy (header->remark, "Assumed to be a SRTM1 tile");
255 }
256 gmt_set_geographic (GMT, GMT_IN);
257 gmtlib_grd_set_units (GMT, header);
258 return (GMT_NOERROR);
259 }
260 else if ((HH->flags[0] == 'L' || HH->flags[0] == 'B') && HH->flags[1] == '2') { /* A Arc/Info BINARY file */
261 if ((fp2 = gmt_fopen (GMT, header->title, "r")) == NULL) {
262 return (GMT_GRDIO_OPEN_FAILED);
263 }
264 /* To use the same parsing header code as in the ASCII file case where header and data are in the
265 same file, we will swap the file pointers and undo the swap at the end of this function */
266 fpBAK = fp; /* Copy of the input argument '*fp' */
267 fp = fp2;
268 }
269
270 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
271 if (sscanf (record, "%*s %d", &header->n_columns) != 1) {
272 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding ncols record\n");
273 if (fpBAK) gmt_fclose (GMT, fp2);
274 return (GMT_GRDIO_READ_FAILED);
275 }
276 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
277 if (sscanf (record, "%*s %d", &header->n_rows) != 1) {
278 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding nrows record\n");
279 if (fpBAK) gmt_fclose (GMT, fp2);
280 return (GMT_GRDIO_READ_FAILED);
281 }
282 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
283 if (sscanf (record, "%*s %lf", &header->wesn[XLO]) != 1) {
284 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding xll record\n");
285 if (fpBAK) gmt_fclose (GMT, fp2);
286 return (GMT_GRDIO_READ_FAILED);
287 }
288 gmt_str_tolower (record);
289 if (!strncmp (record, "xllcorner", 9U)) header->registration = GMT_GRID_PIXEL_REG; /* Pixel grid */
290 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
291 if (sscanf (record, "%*s %lf", &header->wesn[YLO]) != 1) {
292 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding yll record\n");
293 if (fpBAK) gmt_fclose (GMT, fp2);
294 return (GMT_GRDIO_READ_FAILED);
295 }
296 gmt_str_tolower (record);
297 if (!strncmp (record, "yllcorner", 9U)) header->registration = GMT_GRID_PIXEL_REG; /* Pixel grid */
298 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
299 if (sscanf (record, "%*s %lf", &header->inc[GMT_X]) != 1) {
300 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding cellsize record\n");
301 if (fpBAK) gmt_fclose (GMT, fp2);
302 return (GMT_GRDIO_READ_FAILED);
303 }
304 /* Handle the optional nodata_value record */
305 c = fgetc (fp); /* Get first char of next line... */
306 ungetc (c, fp); /* ...and put it back where it came from */
307 if (c == 'n' || c == 'N') { /* Assume this is a nodata_value record since we found an 'n|N' */
308 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
309 #ifdef DOUBLE_PRECISION_GRID
310 if (sscanf (record, "%*s %lf", &header->nan_value) != 1) {
311 #else
312 if (sscanf (record, "%*s %f", &header->nan_value) != 1) {
313 #endif
314 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info ASCII Grid: Error decoding nan_value_value record\n");
315 if (fpBAK) gmt_fclose (GMT, fp2);
316 gmt_fclose (GMT, fp);
317 return (GMT_GRDIO_READ_FAILED);
318 }
319 }
320 header->inc[GMT_Y] = header->inc[GMT_X];
321 header->wesn[XHI] = header->wesn[XLO] + (header->n_columns - 1 + header->registration) * header->inc[GMT_X];
322 header->wesn[YHI] = header->wesn[YLO] + (header->n_rows - 1 + header->registration) * header->inc[GMT_Y];
323
324 if ((error = gmt_M_err_fail (GMT, gmt_grd_RI_verify (GMT, header, 1), HH->name)))
325 return error;
326
327 if (fpBAK) { /* Case of Arc/Info binary file with a separate header file. We still have things to do. */
328 char tmp[16];
329 /* Read an extra record containing the endianness info */
330 gmt_fgets (GMT, record, GMT_BUFSIZ, fp);
331 if (sscanf (record, "%*s %s", tmp) != 1) {
332 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Arc/Info BINARY Grid: Error decoding endianness record\n");
333 gmt_fclose (GMT, fp2);
334 return (GMT_GRDIO_READ_FAILED);
335 }
336 HH->flags[0] = (tmp[0] == 'L') ? 'L' : 'B';
337
338 header->bits = 32; /* Those float binary files */
339 HH->orig_datatype = GMT_FLOAT;
340 /* Ok, now as mentioned above undo the file pointer swapping (point again to data file) */
341 gmt_fclose (GMT, fp);
342 fp = fpBAK;
343 gmt_fclose (GMT, fp2);
344 }
345
346 if (fp2) gmt_fclose(GMT, fp2);
347 return (GMT_NOERROR);
348 }
349
350 int gmtlib_is_esri_grid (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
351 /* Determine if file is an ESRI Interchange ASCII file */
352 FILE *fp = NULL;
353 char record[GMT_BUFSIZ] = {""}, *e = NULL;
354 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
355
356 if (!strcmp (HH->name, "="))
357 return (GMT_GRDIO_PIPE_CODECHECK); /* Cannot check on pipes */
358 if ((e = gmt_get_ext (HH->name)) && !strcmp (e, GMT_TILE_EXTENSION_REMOTE)) return (-1); /* Watch out for .jp2 tiles since they may contain W|E|S|N codes as well and the ESRI check comes before GDAL check*/
359 if ((fp = gmt_fopen (GMT, HH->name, "r")) == NULL)
360 return (GMT_GRDIO_OPEN_FAILED);
361 if (fgets (record, GMT_BUFSIZ, fp) == NULL) { /* Just get first line. Not using gmt_fgets since we may be reading a binary file */
362 gmt_fclose (GMT, fp);
363 return (GMT_GRDIO_OPEN_FAILED);
364 }
365 gmt_fclose (GMT, fp);
366 if (strncmp (record, "ncols ", 6) ) {
367 /* Failed to find "ncols"; probably a binary file */
368 char *file = NULL;
369 size_t name_len;
370
371 HH->orig_datatype = GMT_SHORT; /* May be overridden below */
372 /* If it got here, see if a companion .hdr file exists (must test upper & lower cases names) */
373 file = strdup (HH->name);
374 gmt_chop_ext (file);
375 name_len = strlen (HH->name);
376 if (name_len < strlen(file) + 4) {
377 /* The file extension had less than 3 chars, which means that 1) it's not an esri file.
378 2) would corrupt the heap with the later strcat (file, ".hdr");
379 On Win this would later cause a crash upon freeing 'file' */
380 gmt_M_str_free (file);
381 return (-1);
382 }
383 if (isupper ((unsigned char) HH->name[name_len - 1]))
384 strcat (file, ".HDR");
385 else
386 strcat (file, ".hdr");
387
388 if (!gmt_access (GMT, file, F_OK)) { /* Now, if first line has BYTEORDER or ncols keywords we are in the game */
389 if ((fp = gmt_fopen (GMT, file, "r")) == NULL)
390 return (GMT_GRDIO_OPEN_FAILED);
391 gmt_fgets (GMT, record, GMT_BUFSIZ, fp); /* Just get first line */
392 gmt_fclose (GMT, fp);
393
394 if (!strncmp (record, "BYTEORDER", 9) ) {
395 sscanf (record, "%*s %c", &HH->flags[0]); /* Store the endianness flag here */
396 strncpy (header->title, file, GMT_GRID_TITLE_LEN80-1);
397 }
398 else if (!strncmp (record, "ncols ", 6) ) { /* Ah. A Arc/Info float binary file with a separate .hdr */
399 strncpy (header->title, file, GMT_GRID_TITLE_LEN80-1);
400 HH->flags[0] = 'L'; /* If is truly 'L' or 'B' we'll find only when parsing the whole header */
401 HH->flags[1] = '2'; /* Flag to let us know the file type */
402 HH->orig_datatype = GMT_FLOAT;
403 }
404 else { /* Cannot do anything with this data */
405 gmt_M_str_free (file);
406 return (-1);
407 }
408
409 gmt_M_str_free (file);
410 }
411 else {
412 /* No header file; see if filename contains w/e/s/n information, as in W|ExxxN|Syy.dem
413 * for GTOPO30 (e.g W020N90.DEM) or N|SyyW|Exxx.hgt for SRTM1|3 (e.g. N00E006.hgt) */
414 size_t len;
415
416 while (gmt_chop_ext (file)); /* Remove all extensions so we know exactly where to look */
417 len = strlen (file);
418 if ((file[len-3] == 'N' || file[len-3] == 'n' || file[len-3] == 'S' || file[len-3] == 's') &&
419 (file[len-7] == 'W' || file[len-7] == 'w' || file[len-7] == 'E' || file[len-7] == 'e')) {
420 /* It is a GTOPO30 or SRTM30 source file without a .hdr companion. */
421 /* see http://dds.cr.usgs.gov/srtm/version1/SRTM30/GTOPO30_Documentation */
422 HH->flags[0] = 'B'; /* GTOPO30 & SRTM30 are Big Endians */
423 HH->flags[1] = '0'; /* Flag to let us know the file type */
424 /* Store the file name with all extensions removed.
425 * We'll use this to create header from file name info */
426 strncpy (header->title, file, GMT_GRID_TITLE_LEN80-1);
427 strcpy (header->remark, "Assumed to be a GTOPO30 or SRTM30 tile");
428 HH->orig_datatype = GMT_SHORT;
429 }
430 else if (name_len > 3 && !(strncmp (&HH->name[name_len-4], ".hgt", 4) && strncmp (&HH->name[name_len-4], ".HGT", 4))) {
431 /* Probably a SRTM1|3 file. In gmtesriio_read_info we'll check further if it is a 1 or 3 sec */
432 if ((file[len-4] == 'E' || file[len-4] == 'e' || file[len-4] == 'W' || file[len-4] == 'w') &&
433 (file[len-7] == 'N' || file[len-7] == 'n' || file[len-7] == 'S' || file[len-7] == 's')) {
434 HH->flags[0] = 'B'; /* SRTM1|3 are Big Endians */
435 HH->flags[1] = '1'; /* Flag to let us know the file type */
436 /* Store the file name with all extensions removed.
437 * We'll use this to create header from file name info */
438 strncpy (header->title, file, GMT_GRID_TITLE_LEN80-1);
439 HH->orig_datatype = GMT_SHORT;
440 }
441 }
442 else {
443 /* Cannot do anything with this data */
444 gmt_M_str_free (file);
445 return (-1); /* Not this kind of file */
446 }
447 }
448 }
449
450 header->type = GMT_GRID_IS_EI;
451 return GMT_NOERROR;
452 }
453
454 int gmt_esri_read_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
455 int error;
456 FILE *fp = NULL;
457 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
458
459 if (!strcmp (HH->name, "=")) /* Pipe in from stdin */
460 fp = GMT->session.std[GMT_IN];
461 else if ((fp = gmt_fopen (GMT, HH->name, "r")) == NULL)
462 return (GMT_GRDIO_OPEN_FAILED);
463
464 if ((error = gmtesriio_read_info (GMT, fp, header)) != 0) {
465 gmt_fclose (GMT, fp);
466 return (error);
467 }
468
469 gmt_fclose (GMT, fp);
470
471 return (GMT_NOERROR);
472 }
473
474 int gmt_esri_write_grd_info (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header) {
475 FILE *fp = NULL;
476 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
477
478 if (!strcmp (HH->name, "=")) /* Write to stdout */
479 fp = GMT->session.std[GMT_OUT];
480 else if ((fp = gmt_fopen (GMT, HH->name, "w")) == NULL)
481 return (GMT_GRDIO_CREATE_FAILED);
482
483 gmtesriio_write_info (GMT, fp, header);
484
485 gmt_fclose (GMT, fp);
486
487 return (GMT_NOERROR);
488 }
489
490 int gmt_esri_read_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
491 int error;
492 bool check, is_binary = false, swap = false;
493 unsigned int col, height_in, ii, in_nx;
494 int row, first_col, last_col, first_row, last_row;
495 unsigned int row2, width_in, *actual_col = NULL;
496 unsigned int nBits = 32U;
497 uint64_t ij, kk, width_out, imag_offset, n_left = 0;
498 size_t n_expected;
499 char *r_mode = NULL;
500 int16_t *tmp16 = NULL;
501 gmt_grdfloat value, *tmp = NULL;
502 FILE *fp = NULL;
503 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
504
505 if (HH->flags[0]) { /* We are dealing with a ESRI .hdr file or GTOPO30, SRTM30, SRTM1|3 */
506 r_mode = "rb";
507 if (((HH->flags[0] == 'M' || HH->flags[0] == 'B') && !GMT_BIGENDIAN) ||
508 (HH->flags[0] == 'L' && GMT_BIGENDIAN))
509 swap = true;
510 nBits = header->bits;
511 is_binary = true;
512 }
513 else
514 r_mode = GMT->current.io.r_mode;
515
516 if (!strcmp (HH->name, "=")) /* Read from pipe */
517 fp = GMT->session.std[GMT_IN];
518 else if ((fp = gmt_fopen (GMT, HH->name, r_mode)) != NULL) {
519 if ((error = gmtesriio_read_info (GMT, fp, header)) != 0) {
520 gmt_fclose (GMT, fp);
521 return (error);
522 }
523 }
524 else
525 return (GMT_GRDIO_OPEN_FAILED);
526
527 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);
528 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
529
530 width_out = width_in; /* Width of output array */
531 if (pad[XLO] > 0) width_out += pad[XLO];
532 if (pad[XHI] > 0) width_out += pad[XHI];
533 n_expected = header->n_columns;
534
535 if (nBits == 32) /* Either an ASCII file or ESRI .HDR with NBITS = 32, in which case we assume it's a file of floats */
536 tmp = gmt_M_memory (GMT, NULL, n_expected, float);
537 else
538 tmp16 = gmt_M_memory (GMT, NULL, n_expected, int16_t);
539
540 header->z_min = DBL_MAX; header->z_max = -DBL_MAX;
541 HH->has_NaNs = GMT_GRID_NO_NANS; /* We are about to check for NaNs and if none are found we retain 1, else 2 */
542 if (is_binary) {
543 int n_rows = header->n_rows;
544 if (last_row - first_row + 1 != n_rows) /* We have a sub-region */
545 if (fseek (fp, (off_t) (first_row * n_expected * 4UL * nBits / 32UL), SEEK_CUR)) {
546 gmt_fclose (GMT, fp);
547 gmt_M_free (GMT, actual_col);
548 if (nBits == 32) gmt_M_free (GMT, tmp); else gmt_M_free (GMT, tmp16);
549 return (GMT_GRDIO_SEEK_FAILED);
550 }
551
552 ij = imag_offset + pad[YHI] * width_out + pad[XLO];
553
554 for (row = first_row; row <= last_row; row++, ij += width_out) {
555 if (nBits == 32) { /* Get one row */
556 if (gmt_M_fread (tmp, 4, n_expected, fp) < n_expected) {
557 gmt_fclose (GMT, fp);
558 gmt_M_free (GMT, actual_col);
559 gmt_M_free (GMT, tmp);
560 return (GMT_GRDIO_READ_FAILED);
561 }
562 }
563 else {
564 if (gmt_M_fread (tmp16, 2, n_expected, fp) < n_expected) {
565 gmt_fclose (GMT, fp);
566 gmt_M_free (GMT, actual_col);
567 gmt_M_free (GMT, tmp16);
568 return (GMT_GRDIO_READ_FAILED);
569 }
570 }
571 for (col = 0, kk = ij; col < width_in; col++, kk++) {
572 if (nBits == 32) {
573 if (swap) {
574 /* need to memcpy because casting from float* to uint32_t*
575 * violates strict-aliasing rules. */
576 uint32_t u;
577 memcpy (&u, &tmp[actual_col[col]], sizeof (uint32_t));
578 u = bswap32 (u);
579 memcpy (&tmp[actual_col[col]], &u, sizeof (uint32_t));
580 }
581 grid[kk] = tmp[actual_col[col]];
582 }
583 else {
584 if (swap) {
585 uint16_t *p = (uint16_t *)&tmp16[actual_col[col]]; /* here casting the pointer is allowed */
586 *p = bswap16 (*p);
587 }
588 grid[kk] = tmp16[actual_col[col]];
589 }
590 if (grid[kk] == header->nan_value) {
591 HH->has_NaNs = GMT_GRID_HAS_NANS;
592 grid[kk] = GMT->session.f_NaN;
593 }
594 else { /* Update z_min, z_max */
595 header->z_min = MIN (header->z_min, (double)grid[kk]);
596 header->z_max = MAX (header->z_max, (double)grid[kk]);
597 }
598 }
599 }
600
601 }
602 else { /* ASCII */
603 n_left = header->nm;
604
605 /* ESRI grids are scanline oriented (top to bottom), as are the GMT grids.
606 * NaNs are not allowed; they are represented by a nodata_value instead. */
607 col = row = 0; /* For the entire file */
608 row2 = 0; /* For the inside region */
609 check = !isnan (header->nan_value);
610 in_nx = header->n_columns;
611 header->n_columns = width_in; /* Needed to be set here due to gmt_M_ijp below */
612 #ifdef DOUBLE_PRECISION_GRID
613 while (fscanf (fp, "%lf", &value) == 1 && n_left) { /* We read all values and skip those not inside our w/e/s/n */
614 #else
615 while (fscanf (fp, "%f", &value) == 1 && n_left) { /* We read all values and skip those not inside our w/e/s/n */
616 #endif
617 tmp[col] = value; /* Build up a single input row */
618 col++;
619 if (col == in_nx) { /* End of input row */
620 if (row >= first_row && row <= last_row) { /* We want a piece (or all) of this row */
621 ij = imag_offset + gmt_M_ijp (header, row2, 0); /* First out index for this row */
622 for (ii = 0; ii < width_in; ii++) {
623 kk = ij + ii;
624 grid[kk] = (check && tmp[actual_col[ii]] == header->nan_value) ? GMT->session.f_NaN : tmp[actual_col[ii]];
625 if (gmt_M_is_fnan (grid[kk])) {
626 HH->has_NaNs = GMT_GRID_HAS_NANS;
627 continue;
628 }
629 /* Update z_min, z_max */
630 header->z_min = MIN (header->z_min, (double)grid[kk]);
631 header->z_max = MAX (header->z_max, (double)grid[kk]);
632 }
633 row2++;
634 }
635 col = 0, row++;
636 }
637 n_left--;
638 }
639 }
640 if (header->z_min == DBL_MAX && header->z_max == -DBL_MAX) /* No valid data values in the grid */
641 header->z_min = header->z_max = NAN;
642
643 gmt_fclose (GMT, fp);
644 gmt_M_free (GMT, actual_col);
645 gmt_M_free (GMT, tmp);
646 if (nBits != 32) gmt_M_free (GMT, tmp16);
647
648 if (n_left) {
649 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Expected % "PRIu64 " points, found only % "PRIu64 "\n", header->nm, header->nm - n_left);
650 return (GMT_GRDIO_READ_FAILED);
651 }
652
653 header->n_columns = width_in;
654 header->n_rows = height_in;
655 gmt_M_memcpy (header->wesn, wesn, 4, double);
656
657 return (GMT_NOERROR);
658 }
659
660 int gmt_esri_write_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode, bool floating) {
661 unsigned int i2, j, j2, width_out, height_out, last;
662 int first_col, last_col, first_row, last_row;
663 unsigned int i, *actual_col = NULL;
664 uint64_t ij, width_in, kk, imag_offset;
665 char item[GMT_LEN64], c[2] = {0, 0};
666 FILE *fp = NULL;
667 struct GMT_GRID_HEADER_HIDDEN *HH = gmt_get_H_hidden (header);
668
669 if (fabs (header->inc[GMT_X] / header->inc[GMT_Y] - 1.0) > GMT_CONV8_LIMIT)
670 return (GMT_GRDIO_ESRI_NONSQUARE); /* Only square pixels allowed */
671 if (!strcmp (HH->name, "=")) /* Write to pipe */
672 fp = GMT->session.std[GMT_OUT];
673 else if ((fp = gmt_fopen (GMT, HH->name, GMT->current.io.w_mode)) == NULL)
674 return (GMT_GRDIO_CREATE_FAILED);
675 else
676 gmtesriio_write_info (GMT, fp, header);
677
678 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);
679 (void)gmtlib_init_complex (header, complex_mode, &imag_offset); /* Set offset for imaginary complex component */
680
681 width_in = width_out; /* Physical width of input array */
682 if (pad[XLO] > 0) width_in += pad[XLO];
683 if (pad[XHI] > 0) width_in += pad[XHI];
684
685 gmt_M_memcpy (header->wesn, wesn, 4, double);
686
687 /* Store header information and array */
688
689 i2 = first_col + pad[XLO];
690 last = width_out - 1;
691 for (j = 0, j2 = first_row + pad[YHI]; j < height_out; j++, j2++) {
692 ij = imag_offset + j2 * width_in + i2;
693 c[0] = '\t';
694 for (i = 0; i < width_out; i++) {
695 if (i == last) c[0] = '\n';
696 kk = ij+actual_col[i];
697 if (gmt_M_is_fnan (grid[kk]))
698 snprintf (item, GMT_LEN64, "%ld%c", lrintf (header->nan_value), c[0]);
699 else if (floating) {
700 snprintf (item, GMT_LEN64-1, GMT->current.setting.format_float_out, grid[kk]);
701 strcat (item, c);
702 }
703 else
704 snprintf (item, GMT_LEN64, "%ld%c", lrint ((double)grid[kk]), c[0]);
705 gmt_M_fputs (item, fp);
706 }
707 }
708
709 gmt_M_free (GMT, actual_col);
710 gmt_fclose (GMT, fp);
711
712 return (GMT_NOERROR);
713 }
714
715 int gmt_esri_writei_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
716 /* Standard integer values on output only */
717 return (gmt_esri_write_grd (GMT, header, grid, wesn, pad, complex_mode, false));
718 }
719
720 int gmt_esri_writef_grd (struct GMT_CTRL *GMT, struct GMT_GRID_HEADER *header, gmt_grdfloat *grid, double wesn[], unsigned int *pad, unsigned int complex_mode) {
721 /* Write floating point on output */
722 return (gmt_esri_write_grd (GMT, header, grid, wesn, pad, complex_mode, true));
723 }
724