/*--------------------------------------------------------------------
 *	$Id: xyz2grd.c,v 1.19 2004/04/25 09:10:46 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 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)
 *		4
 */
 
#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, high_low = 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 'b':
					b_only = TRUE;
				case 'H':
				case 'R':
				case 'V':
				case ':':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], &w, &e, &s, &n);
					break;

				case 'A':
					switch (argv[i][2]) {
						case 'n':	/* Count the number of values for each node */
							count = TRUE;
							break;
						case 'l':	/* Return the lowest (minimum) value at each node */
							high_low = -1;
							break;
						case 'u':	/* Return the upper (maximum)  value at each node */
							high_low = +1;
							break;
						case 'z':	/* Return the mean value at each node */
						case '\0':
							average = FALSE;
							break;
						default:
							fprintf (stderr, "%s: GMT SYNTAX ERROR -A option:  Select -An, Al, Au or -A[z]\n", GMT_program);
							error++;
						break;
					}
					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':	/* Obsolete, but backward compatibility prevails [use -f instead] */
					GMT_io.in_col_type[0] = GMT_io.out_col_type[0] = GMT_IS_LON;
					GMT_io.in_col_type[1] = GMT_io.out_col_type[1] = GMT_IS_LAT;
					fprintf (stderr, "%s: Option -L is obsolete (but is processed correctly).  Please use -f instead\n", GMT_program);
					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<grdfile> -I<dx[m|c]>[/<dy[m|c]>] -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-A[n|z|u|l]] [-D<xunit/yunit/zunit/scale/offset/title/remark>] [-F]\n");
		fprintf (stderr, "\t[-H[<nrec>]] [-N<nodata>] [-S[<zfile]] [-V] [-Z[<flags>]] [-:]\n\t[-bi[s][<n>]] [-f[i|o]<colinfo>]\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 <dx> and/or <dy> 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   Append u (-Au): Keep maximum value if multiple entries per node.\n");
		fprintf (stderr, "\t   Append l (-Al): Keep minimum value if multiple entries per node.\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-N set value for nodes without input xyz triplet [Default is NaN]\n");
		fprintf (stderr, "\t   Z-table entries that equal <nodata> 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 <zfile> (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 <n> bytes, append s<n> [<n> = 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 ('f');
		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[GMT_IN] = 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[GMT_IN] || r.binary) && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = 3;
	if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < 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 .gmtcommands4 */

	if (b_only && z_only) fprintf (stderr, "%s: GMT Warning.  -Z overrides -bi\n", GMT_program);

	if (GMT_io.binary[GMT_IN] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[GMT_IN], type[GMT_io.single_precision[GMT_IN]]);
	}

	if (GMT_io.binary[GMT_IN] && GMT_io.in_col_type[2] & GMT_IS_RATIME && GMT_io.single_precision[GMT_IN]) {
		fprintf (stderr, "%s: Warning: Your single precision binary input data are unlikely to hold absolute time coordinates without serious truncation.\n", GMT_program);
		fprintf (stderr, "%s: Warning: You must use double precision when storing absolute time coordinates in binary data tables.\n", GMT_program);
	}
	if (GMT_io.in_col_type[2] & GMT_IS_RATIME) {
		fprintf (stderr, "%s: Warning: Single precision grid format is unlikely to hold absolute time coordinates without serious truncation.\n", GMT_program);
		fprintf (stderr, "%s: Warning: You should select a grid format that allows for double precision elements.\n", GMT_program);
	}
		
	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 (GMT_IN);
#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[GMT_IN]) 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 (GMT_OUT);
#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 doesn't match expectations */
				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 = GMT_x_to_i (in[0], w, idx, off, grd.nx);
				if (ii == grd.nx) ii--, n_confused++;
				jj = GMT_y_to_j (in[1], s, idy, off, grd.ny);
				if (jj == -1) jj++, n_confused++;
				ij = jj * grd.nx + ii;
				if (high_low) {	/* Always come here if looking for extreme values */
					if (flag[ij]) {	/* Already assigned the first value */
						if ((high_low == -1 && (in[z] < (double)a[ij])) || (high_low == +1 && (in[z] > (double)a[ij]))) a[ij] = (float)in[z];
					}
					else {	/* First time, just assign the current value */
						a[ij] = (float)in[z];
						flag[ij] = 1;
					}
				}
				else { 	/* Add up incase we must average */
					a[ij] += (float)in[z];
					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) {	/* This catches nodes with one value or the -Al|u single values */
				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_io.in_col_type[0] == GMT_IS_LON && (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 indices: 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);
}
