/*--------------------------------------------------------------------
 *	$Id: grd2xyz.c,v 1.14 2004/04/25 09:10:45 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
 *--------------------------------------------------------------------*/
/*
 * grd2xyz.c reads a grd file and prints out the x,y,z values to
 * standard output.
 *
 * Author:	Paul Wessel
 * Date:	3-JAN-1991
 * Version:	4
 */

#include "gmt.h"

float *z;

main (int argc, char **argv)
{
	BOOLEAN error = FALSE, global = FALSE, z_only = FALSE, b_only = FALSE, suppress = FALSE;
	BOOLEAN first = TRUE, reverse = FALSE, weight = FALSE;

	int i, j, k, ij, nm, nx, ny, one_or_zero, gmt_ij, n_suppressed = 0, n_out, n_files = 0;
	int n_total = 0;

	double w, e, s, n, d, *x, *y, out[4], d_value;

	struct GRD_HEADER grd;

	struct GMT_Z_IO r;
	
	argc = GMT_begin (argc, argv);

	w = e = s = n = 0.0;
	GMT_init_z_io (&r, FALSE);
	
	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;
					
				/* Supplemental options */
				
				case 'L':
					GMT_geographic_out = TRUE;
					break;
				case 'Z':
					error += GMT_parse_z_io (&argv[i][2], &r, FALSE);
					z_only = TRUE;
					break;
				case 'S':
					suppress = TRUE;
					if (argv[i][2] == 'r') reverse = TRUE;
					break;
				case 'W':
					weight = TRUE;
					out[3] = (argv[i][2]) ? atof (&argv[i][2]) : 1.0;
					break;
					
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grd2xyz %s - Converting netCDF grdfile(s) to ASCII xyz data\n\n", GMT_VERSION);
		fprintf( stderr, "usage: grd2xyz grdfile(s) [-H[<nrec>]] [-L] [-Rw/e/s/n] [-S[r]] [-V] [-W[<weight>]]\n");
		fprintf( stderr, "\t[-Z[<flags>]] [-:] [-bo[s][<n>]] [-f[i|o]<colinfo>] > xyzfile\n");

		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tgrdfile is one or more grd files to convert\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-H Write 1 ASCII header record [Default is no header]\n");
		fprintf (stderr, "\t-L means that x is longitude, i.e. assumed to be periodic in 360\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t-S Suppress output for nodes whose z equals NaN [Default prints all nodes]\n");
		fprintf (stderr, "\t   Append r to reverse the suppression (only output NaN nodes)\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-W Write xyzw using supplied weight (or 1 if not given) [Default is xyz]\n");
		fprintf (stderr, "\t-Z sets exact specification of resulting 1-column output 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   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  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");
		GMT_explain_option (':');
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (n_files == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one input file\n", GMT_program);
		error++;
	}

	if (b_only && z_only) GMT_io.binary[GMT_OUT] = FALSE;

	if ((GMT_io.binary[GMT_OUT] || r.binary) && gmtdefs.io_header[GMT_OUT]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary output data cannot have header -H\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 -bo\n", GMT_program);

#ifdef SET_IO_MODE
		GMT_setmode (GMT_OUT);
#endif

	n_out = (weight) ? 4 : 3;
	for (k = 1; k < argc; k++) {
		if (argv[k][0] == '-') continue;	/* Skip the options */
		
		if (GMT_read_grd_info (argv[k], &grd)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, argv[k]);
			exit (EXIT_FAILURE);
		}

		if (gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[k]);

		nm = grd.nx * grd.ny;
		n_total += nm;
	
		if (e > w && n > s) {	/* Subset */
			global = (fabs (grd.x_max - grd.x_min) == 360.0);
			if (!global && (w < grd.x_min || e > grd.x_max)) error = TRUE;
			if (s < grd.y_min || n > grd.y_max) error = TRUE;
			if (error) {
				fprintf (stderr, "%s: GMT ERROR: Subset exceeds data domain!\n", GMT_program);
				exit (EXIT_FAILURE);
			}
			one_or_zero = (grd.node_offset) ? 0 : 1;
			nx = irint ((e - w) / grd.x_inc) + one_or_zero;
			ny = irint ((n - s) / grd.y_inc) + one_or_zero;
		
			z = (float *) GMT_memory (VNULL, (size_t) (nx * ny), sizeof (float), GMT_program);
		
			if (GMT_read_grd (argv[k], &grd, z, w, e, s, n, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error reading file %s\n", GMT_program, argv[k]);
				exit (EXIT_FAILURE);
			}
		}
		else {
			z = (float *) GMT_memory (VNULL, (size_t) nm, sizeof (float), GMT_program);

			if (GMT_read_grd (argv[k], &grd, z, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
				fprintf (stderr, "%s: Error reading file %s\n", GMT_program, argv[k]);
				exit (EXIT_FAILURE);
			}
		}

		GMT_set_z_io (&r, &grd);

		if (z_only) {
			if (gmtdefs.io_header[GMT_OUT] && !r.binary) fprintf (GMT_stdout, "%s\n", grd.z_units);

			for (ij = 0; ij < r.n_expected; ij++) {
				(r.get_gmt_ij) (&r, ij, &gmt_ij);
				d_value = z[gmt_ij];
				if (suppress && (GMT_is_dnan (d_value) + reverse) == 1) {
					n_suppressed++;
					continue;
				}
				if ((r.x_missing && r.gmt_i == r.x_period) || (r.y_missing && r.gmt_j == 0)) continue;
				(r.write_item) (GMT_stdout, d_value);
			}
		}
		else {
	
			x = (double *) GMT_memory (VNULL, (size_t) grd.nx, sizeof (double), GMT_program);
			y = (double *) GMT_memory (VNULL, (size_t) grd.ny, sizeof (double), GMT_program);

			/* Compute grid node positions once only */

			d = (grd.node_offset) ? 0.5 : 0.0;

			for (j = 0; j < grd.ny; j++) y[j] = GMT_j_to_y (j, grd.y_min, grd.y_max, grd.y_inc, d, grd.ny);
			for (i = 0; i < grd.nx; i++) x[i] = GMT_i_to_x (i, grd.x_min, grd.x_max, grd.x_inc, d, grd.nx);
		
			if (gmtdefs.io_header[GMT_OUT] && first) {
				if (!grd.x_units[0]) strcpy (grd.x_units, "x");
				if (!grd.y_units[0]) strcpy (grd.y_units, "y");
				if (!grd.z_units[0]) strcpy (grd.z_units, "z");
				if (gmtdefs.xy_toggle[0])
					fprintf (GMT_stdout, "%s\t%s\t%s", grd.y_units, grd.x_units, grd.z_units);
				else
					fprintf (GMT_stdout, "%s\t%s\t%s", grd.x_units, grd.y_units, grd.z_units);
				if (weight)
					fprintf (GMT_stdout, "\tweight\n");
				else
					fprintf (GMT_stdout, "\n");
				first = FALSE;
			}
	
			for (j = ij = 0; j < grd.ny; j++) for (i = 0; i < grd.nx; i++, ij++) {
				out[2] = z[ij];
				if (suppress && (GMT_is_dnan (out[2]) + reverse) == 1) {
					n_suppressed++;
					continue;
				}
				out[0] = x[i];	out[1] = y[j];
				GMT_output (GMT_stdout, n_out, out);
			}
			GMT_free ((void *)x);
			GMT_free ((void *)y);
		}
	
		GMT_free ((void *)z);
	}
		
	if (gmtdefs.verbose) fprintf (stderr, "%s: %d values extracted\n", GMT_program, n_total - n_suppressed);
	if (n_suppressed && gmtdefs.verbose) {
		if (reverse)
			fprintf (stderr, "%s: %d finite values suppressed\n", GMT_program, n_suppressed);
		else
			fprintf (stderr, "%s: %d NaN values suppressed\n", GMT_program, n_suppressed);
	}

	GMT_end (argc, argv);
}
