/*-------------------------------------------------------------------- * $Id: xyz2grd.c,v 1.4.4.2 2002/02/27 17:41:10 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING file for copying and redistribution conditions. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; version 2 of the License. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * xyz2grd reads a xyz file from standard input and creates the * corresponding grd-file. The input file does not have to be * sorted. xyz2grd will report if some nodes are missing. nx and * ny are computed as: * nx = rint((east-west)/dx) + 1, ny = rint((north-south)/dy) + 1 * * Author: Paul Wessel * Date: 3-JAN-1991-2000 * Version: 2.0 Based on old v1.x * Version: 3.1 Based on old 3.0 * Version: 3.2 Added -S option for swapping only * 3.3.1 -S option need to set fopen modes to rb an wb * 3.3.4 -L option finally fixed * 3.3.5 -L option finally fixed, part II. Also added -A[z|n] * Version: 3.3.6 * 3.4 Now can handle multiple input files (except with -S) * Version: 3.4.1 */ #include "gmt.h" main (int argc, char **argv) { BOOLEAN error = FALSE, got_input = FALSE, pixel = FALSE, z_only = FALSE, count = FALSE, average = TRUE; BOOLEAN nodata_set = FALSE, b_only = FALSE, just_swap = FALSE, done = FALSE, nofile = TRUE; int i, ij, gmt_ij, ii, jj, z, nm, n_read = 0, n_filled = 0, n_used = 0, one_or_zero, fno, n_args, *flag; int n_empty = 0, n_stuffed = 0, n_bad = 0, n_expected_fields, n_fields, n_confused = 0, n_files = 0; double w, e, s, n, dx = 0.0, dy = 0.0, *in, d_value, off, idx, idy, no_data_d; float no_data_f, *a; char *grdfile, *zfile, line[BUFSIZ], input[BUFSIZ]; FILE *fp = NULL, *fps = NULL; struct GRD_HEADER grd; struct GMT_Z_IO r; argc = GMT_begin (argc, argv); grdfile = zfile = CNULL; input[0] = 0; w = e = s = n = 0.0; no_data_f = GMT_f_NaN; no_data_d = GMT_d_NaN; GMT_init_z_io (&r, TRUE); for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch (argv[i][1]) { /* Common parameters */ case 'H': case 'R': case 'V': case ':': case '\0': error += GMT_get_common_args (argv[i], &w, &e, &s, &n); break; case 'b': error += GMT_io_selection (&argv[i][2]); b_only = TRUE; break; case 'A': if (argv[i][2] == 'n') count = TRUE; else if (argv[i][2] == '\0' || argv[i][2] == 'z') average = FALSE; else { fprintf (stderr, "%s: GMT SYNTAX ERROR -A option: Select -An or -A[z]\n", GMT_program); error++; } break; case 'D': strcpy (input, &argv[i][2]); got_input = TRUE; break; case 'F': pixel = TRUE; break; case 'G': grdfile = &argv[i][2]; break; case 'I': GMT_getinc (&argv[i][2], &dx, &dy); break; case 'L': GMT_geographic_in = TRUE; break; case 'N': if (!argv[i][2]) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Must specify value or NaN\n", GMT_program); error++; } else { no_data_d = (argv[i][2] == 'N' || argv[i][2] == 'n') ? GMT_d_NaN : atof (&argv[i][2]); no_data_f = (float)no_data_d; nodata_set = TRUE; } break; case 'S': just_swap = TRUE; if (argv[i][2]) zfile = &argv[i][2]; break; case 'Z': error += GMT_parse_z_io (&argv[i][2], &r, TRUE); z_only = TRUE; break; default: error = TRUE; GMT_default_error (argv[i][1]); break; } } else n_files++; } if (argc == 1 || GMT_quick) { fprintf (stderr, "xyz2grd %s - Converting [xy]z data to a GMT grdfile\n\n", GMT_VERSION); fprintf (stderr, "usage: xyz2grd [<[xy]zfile(s)>] -G -I[/] -R\n"); fprintf (stderr, "\t[-A[n|z]] [-D] [-F] [-H[]] [-L]\n"); fprintf (stderr, "\t[-N] [-S[]] [-:] [-bi[s][]]\n"); if (GMT_quick) exit (EXIT_FAILURE); fprintf (stderr, "\txyzfile(s) is a 1- or 3-column file [or standard input] with [x,y,]z values\n"); fprintf (stderr, "\t-G to name the output grdfile.\n"); fprintf (stderr, "\t-I specifies grid size(s). Append m [or c] to and/or for minutes [or seconds]\n"); GMT_explain_option ('R'); fprintf (stderr, "\n\tOPTIONS:\n"); fprintf (stderr, "\t-A (or -Az): Add multiple entries at the same node.\n"); fprintf (stderr, "\t Append n (-An): Count number of multiple entries per node instead.\n"); fprintf (stderr, "\t [Default (no -A option) will compute mean values]\n"); fprintf (stderr, "\t-D to enter header information. Specify '=' to get default value\n"); fprintf (stderr, "\t-F will force pixel registration [Default is grid registration]\n"); GMT_explain_option ('H'); fprintf (stderr, "\t-L means that x is longitude, i.e. assumed to be periodic in 360\n"); fprintf (stderr, "\t-N set value for nodes without input xyz triplet [Default is NaN]\n"); fprintf (stderr, "\t Z-table entries that equal are replaced by NaN\n"); fprintf (stderr, "\t-S Requires -Z and will swap the byte-order of the input data and write\n"); fprintf (stderr, "\t the result to (or stdout if no file given). No grdfile created!\n"); fprintf (stderr, "\t For this option, only one input file (or stdin) is allowed.\n"); GMT_explain_option ('V'); fprintf (stderr, "\t-Z sets exact specification of incoming 1-column z-table\n"); fprintf (stderr, "\t If data is in row format, state if first row is at T(op) or B(ottom)\n"); fprintf (stderr, "\t Then, append L or R to indicate starting point in row\n"); fprintf (stderr, "\t If data is in column format, state if first columns is L(left) or R(ight)\n"); fprintf (stderr, "\t Then, append T or B to indicate starting point in column\n"); fprintf (stderr, "\t To skip a header of size bytes, append s [ = 0]\n"); fprintf (stderr, "\t To swap byte-order in 2-byte words, append w\n"); fprintf (stderr, "\t Append x if gridline-registered, periodic data in x without repeating column at xmax\n"); fprintf (stderr, "\t Append y if gridline-registered, periodic data in y without repeating row at ymax\n"); fprintf (stderr, "\t Specify one of the following data types (all binary except a):\n"); fprintf (stderr, "\t a Ascii\n"); fprintf (stderr, "\t c signed 1-byte character\n"); fprintf (stderr, "\t u unsigned 1-byte character\n"); fprintf (stderr, "\t h signed short 2-byte integer\n"); fprintf (stderr, "\t H unsigned short 2-byte integer\n"); fprintf (stderr, "\t i signed 4-byte integer\n"); fprintf (stderr, "\t I unsigned 4-byte integer\n"); fprintf (stderr, "\t l signed long (4- or 8-byte) integer\n"); fprintf (stderr, "\t f 4-byte floating point single precision\n"); fprintf (stderr, "\t d 8-byte floating point double precision\n"); fprintf (stderr, "\t [Default format is scanline orientation in ASCII representation: -ZTLa]\n"); fprintf (stderr, "\t This option assumes all nodes have data values.\n"); GMT_explain_option (':'); GMT_explain_option ('i'); GMT_explain_option ('n'); fprintf (stderr, "\t Default is 3 input columns\n"); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (just_swap) { /* Reading and writing binary file */ if (n_files > 1) { fprintf (stderr, "%s: GMT SYNTAX ERROR: -S can only handle one input file\n", GMT_program); error++; } else { strcpy (GMT_io.r_mode, "rb"); strcpy (GMT_io.w_mode, "wb"); } } if (b_only && z_only) GMT_io.binary[0] = FALSE; if (!just_swap) { if (!project_info.region_supplied) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (dx <= 0.0 || dy <= 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -I option. Must specify positive increment(s)\n", GMT_program); error++; } if (!grdfile) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -G: Must specify output file\n", GMT_program); error++; } } if (just_swap && !z_only) { fprintf (stderr, "%s: GMT SYNTAX ERROR option -S: Must also specify -Z\n", GMT_program); error++; } if ((GMT_io.binary[0] || r.binary) && gmtdefs.io_header) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data cannot have header -H\n", GMT_program); error++; } if (GMT_io.binary[0] && GMT_io.ncol[0] == 0) GMT_io.ncol[0] = 3; if (GMT_io.binary[0] && GMT_io.ncol[0] < 3) { fprintf (stderr, "%s: GMT SYNTAX ERROR. Binary input data (-bi) must have at least 3 columns\n", GMT_program); error++; } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ if (b_only && z_only) fprintf (stderr, "%s: GMT Warning. -Z overrides -bi\n", GMT_program); if (GMT_io.binary[0] && gmtdefs.verbose) { char *type[2] = {"double", "single"}; fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[0], type[GMT_io.single_precision[0]]); } GMT_grd_init (&grd, argc, argv, FALSE); /* if (!project_info.region) d_swap (s, e); */ /* Got w/s/e/n, make into w/e/s/n */ /* Decode grd information given, if any */ if (got_input) GMT_decode_grd_h_info (input, &grd); if (pixel) { grd.node_offset = 1; one_or_zero = 0; off = 0.0; } else { grd.node_offset = 0; one_or_zero = 1; off = 0.5; } if (just_swap) r.swab = TRUE; else { grd.nx = irint ((e-w)/dx) + one_or_zero; grd.ny = irint ((n-s)/dy) + one_or_zero; grd.x_min = w; grd.x_max = e; grd.y_min = s; grd.y_max = n; grd.x_inc = dx; grd.y_inc = dy; GMT_grd_RI_verify (&grd, 1); if (gmtdefs.verbose) fprintf (stderr, "%s: nx = %d ny = %d\n", GMT_program, grd.nx, grd.ny); nm = grd.nx * grd.ny; a = (float *) GMT_memory (VNULL, (size_t)nm, sizeof (float), GMT_program); flag = (int *) GMT_memory (VNULL, (size_t)nm, sizeof (int), GMT_program); idx = 1.0 / dx; idy = 1.0 / dy; n_expected_fields = (z_only) ? 1 : 3; z = (z_only) ? 0 : 2; ij = -1; /* Will be incremented to 0 or set first time around */ } GMT_set_z_io (&r, &grd); if (n_files > 0) nofile = FALSE; else n_files = 1; n_args = (argc > 1) ? argc : 2; for (fno = 1; !done && fno < n_args; fno++) { /* Loop over input files, if any */ if (!nofile && argv[fno][0] == '-') continue; if (nofile) { /* Just read standard input */ fp = GMT_stdin; done = TRUE; if (gmtdefs.verbose) fprintf (stderr, "%s: Reading from standard input\n", GMT_program); #ifdef SET_IO_MODE GMT_setmode (0); #endif } else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) { fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]); continue; } if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]); if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) GMT_fgets (line, BUFSIZ, fp); if (just_swap) { /* Just swap data and exit */ if (gmtdefs.verbose) fprintf (stderr, "%s: Swapping data bytes only\n", GMT_program); if (!zfile) { fps = GMT_stdout; #ifdef SET_IO_MODE GMT_setmode (1); #endif } else if ((fps = GMT_fopen (zfile, GMT_io.w_mode)) == NULL) { fprintf (stderr, "%s: Cannot create file %s\n", GMT_program, zfile); exit (EXIT_FAILURE); } if (r.skip) fseek (fp, (long)r.skip, SEEK_CUR); while ((r.read_item) (fp, &d_value)) (r.write_item) (fps, d_value); if (fp != GMT_stdin) GMT_fclose (fp); if (fps != GMT_stdout) GMT_fclose (fps); exit (EXIT_SUCCESS); /* We are done here */ } if (z_only) { /* Read separately because of all the possible formats */ if (nodata_set && GMT_is_dnan (no_data_d)) nodata_set = FALSE; /* No point testing */ if (r.skip) fseek (fp, (long)r.skip, SEEK_CUR); while ((r.read_item) (fp, &d_value)) { ij++; if (ij == r.n_expected) { fprintf (stderr, "%s: More than %d records, only %d was expected (aborting)!\n", GMT_program, ij, r.n_expected); exit (EXIT_FAILURE); } (r.get_gmt_ij) (&r, ij, &gmt_ij); if (nodata_set && d_value == no_data_d) a[gmt_ij] = GMT_f_NaN; else a[gmt_ij] = (float)d_value; } ij++; if (ij != r.n_expected) { /* Input amount doesnt match expecatations */ fprintf (stderr, "%s: Found %d records, but %d was expected (aborting)!\n", GMT_program, ij, r.n_expected); exit (EXIT_FAILURE); } GMT_check_z_io (&r, a); /* This fills in missing periodic row or column */ } else { /* Get x, y, z */ while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) { /* Not yet EOF */ if (GMT_io.status & GMT_IO_MISMATCH) { fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields, n_expected_fields, n_read); exit (EXIT_FAILURE); } n_read++; if (in[1] < s || in[1] > n) continue; /* Outside y-range */ if (in[0] < w || in[0] > e) continue; /* Outside y-range */ /* Ok, we are inside the region - process data */ ii = (int)floor ((in[0] - w) * idx + off); if (ii == grd.nx) ii--, n_confused++; jj = (int)floor ((n - in[1]) * idy + off); if (jj == grd.ny) jj--, n_confused++; ij = jj * grd.nx + ii; a[ij] += (float)in[z]; /* Add up incase we must average */ flag[ij]++; n_used++; } } if (fp != GMT_stdin) GMT_fclose (fp); } if (!z_only) { /* xyz data could have resulted in duplicates */ for (ij = 0; ij < nm; ij++) { /* Check if all nodes got one value only */ if (flag[ij] == 1) { if (count) a[ij] = 1.0; n_filled++; } else if (flag[ij] == 0) { n_empty++; a[ij] = no_data_f; } else { /* More than 1 value went to this node */ if (count) a[ij] = (float)flag[ij]; else if (average) a[ij] /= (float)flag[ij]; n_filled++; n_stuffed++; } } if (GMT_geographic_in && (fabs (fabs (grd.x_max - grd.x_min) - 360.0) < GMT_CONV_LIMIT) && !grd.node_offset) { /* Make sure longitudes got replicated */ int j, ij_left, ij_right; for (j = 0; j < grd.ny; j++) { ij_left = j * grd.nx; ij_right = ij_left + grd.nx - 1; if (flag[ij_left] && !flag[ij_right]) { a[ij_right] = a[ij_left]; n_empty--; } else if (flag[ij_right] && !flag[ij_left]) { a[ij_left] = a[ij_right]; n_empty--; } } } if (gmtdefs.verbose) { sprintf (line, "%s\n", gmtdefs.d_format); fprintf (stderr, "%s: n_read: %d n_used: %d n_filled: %d n_empty: %d set to ", GMT_program, n_read, n_used, n_filled, n_empty); (GMT_is_dnan (no_data_d)) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, no_data_d); if (n_bad) fprintf (stderr, "%s: %d records unreadable\n", GMT_program, n_bad); if (n_stuffed) fprintf (stderr, "%s: Warning - %d nodes had multiple entries that were averaged\n", GMT_program, n_stuffed); if (n_confused) fprintf (stderr, "%s: Warning - %d values gave bad indeces: Pixel vs gridline confusion?\n", GMT_program, n_confused); } } if (GMT_write_grd (grdfile, &grd, a, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) { fprintf (stderr, "%s: Error writing file %s\n", GMT_program, grdfile); exit (EXIT_FAILURE); } GMT_free ((void *)a); GMT_free ((void *)flag); GMT_end (argc, argv); }