/*--------------------------------------------------------------------
 *	$Id: grd2xyz.c,v 1.3.4.2 2002/02/27 17:58:55 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
 *--------------------------------------------------------------------*/
/*
 * 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:	3.4.1
 */

#include "gmt.h"

float *z;

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

	int i, j, ij, nm, nx, ny, one_or_zero, gmt_ij, n_suppressed = 0;

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

	char *grdfile = CNULL;

	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 '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 '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;
					
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			grdfile = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grd2xyz %s - Converting a netCDF grdfile to an ASCII xyz-file\n\n", GMT_VERSION);
		fprintf( stderr, "usage: grd2xyz grdfile [-H[<nrec>]] [-L] [-Rw/e/s/n] [-S[r]] [-V] [-Z[<flags>]] [-:] [-bo[s][<n>]] > xyzfile\n");

		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tgrdfile is the grd file 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-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 ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}

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

	if ((GMT_io.binary[1] || r.binary) && gmtdefs.io_header) {
		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 .gmtcommands */

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

	if (GMT_read_grd_info (grdfile, &grd)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}

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

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	nm = grd.nx * grd.ny;
	
	if (e > w && n > s) {
		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 (grdfile, &grd, z, w, e, s, n, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grdfile);
			exit (EXIT_FAILURE);
		}
	}
	else {
		z = (float *) GMT_memory (VNULL, (size_t) nm, sizeof (float), GMT_program);

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

	GMT_set_z_io (&r, &grd);

	if (z_only) {
		if (gmtdefs.io_header && !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] = (j == (grd.ny-1)) ? grd.y_min + d * grd.y_inc : grd.y_max - (j + d) * grd.y_inc;

		for (i = 0; i < grd.nx; i++) x[i] = (i == (grd.nx-1)) ? grd.x_max - d * grd.x_inc: grd.x_min + (i + d) * grd.x_inc;
		if (gmtdefs.io_header) {
			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)
				fprintf (GMT_stdout, "%s\t%s\t%s\n", grd.y_units, grd.x_units, grd.z_units);
			else
				fprintf (GMT_stdout, "%s\t%s\t%s\n", grd.x_units, grd.y_units, grd.z_units);
		}
	
		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, 3, 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, nm - 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);
}