/*--------------------------------------------------------------------
 *	$Id: grdtrack.c,v 1.19 2004/07/16 17:07:02 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
 *--------------------------------------------------------------------*/
/*
 * grdtrack reads a xyfile, opens the 2d binary gridded grdfile, 
 * and samples the dataset at the xy positions with a bilinear or bicubic
 * interpolant.  This new data is added to the input as an extra column
 * and printed to standard output.  In order to evaluate derivatives along
 * the edges of the grdfile region, we assume natural bicubic spline
 * boundary conditions (d2z/dn2 = 0, n being the normal to the edge;
 * d2z/dxdy = 0 in the corners).  Rectangles of size x_inc by y_inc are 
 * mapped to [0,1] x [0,1] by affine transformation, and the interpolation
 * done on the normalized rectangle.
 *
 * Author:	Walter H F Smith
 * Date:	23-SEP-1993
 * 
 * Based on the original grdtrack, which had this authorship/date/history:
 *
 * Author:	Paul Wessel
 * Date:	29-JUN-1988
 * Revised:	5-JAN-1990	PW: Updated to v.2.0
 *		4-AUG-1993	PW: Added -Q
 *		14-AUG-1998	PW: GMT 3.1
 *  Modified:	10 Jul 2000 3.3.5  by PW to allow plain -L to indicate geographic coordinates
 * Version:	4
 */

#include "gmt.h"

main (int argc, char **argv)
{
	int i, j, ix, iy, mx, my, nx, ny, n_read = 0, n_fields, n_points = 0, one_or_zero;
	int n_output = 0, n_expected_fields = 0;
	
	BOOLEAN error = FALSE, bilinear = FALSE, suppress = FALSE, node = FALSE, z_only = FALSE, pure_ascii = FALSE;
	
	double value, west, east, south, north, threshold = 1.0, i_dx, i_dy, half, *in, *out;
	
	float *f;

	char *grdfile, stuff[BUFSIZ], line[BUFSIZ];
	
	FILE *fp = NULL;
	
	struct GRD_HEADER grd;
	
	struct GMT_EDGEINFO edgeinfo;
	
	struct GMT_BCR bcr;

	grdfile = CNULL;
	west = east = south = north = 0.0;
	stuff[0] = 0;
	out = (double *)NULL;
	
	argc = GMT_begin (argc, argv);
	GMT_boundcond_init (&edgeinfo);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			
				/* Common parameters */
			
				case 'H':
				case 'R':
				case 'V':
				case ':':
				case 'b':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;

				/* Supplemental parameters */
				
				case 'G':
					grdfile = &argv[i][2];
					break;
				case 'L':
					if (argv[i][2]) {
						error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]);
						if (edgeinfo.gn) {
							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;
						}
					}
					else {
						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 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':
					node = TRUE;
					break;
				case 'Q':
					bilinear = TRUE;
					threshold = (argv[i][2]) ? atof (&argv[i][2]) : 1.0;
					break;
				case 'S':
					suppress = TRUE;
					break;
				case 'Z':
					z_only = TRUE;
					break;

				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else if ((fp = GMT_fopen (argv[i], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[i]);
			exit (EXIT_FAILURE);
		}
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdtrack %s - Sampling of a 2-D gridded netCDF grdfile along 1-D trackline\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdtrack <xyfile> -G<grdfile> [-H[<nrec>]] [-L<flag>] [-M[<flag>]] [-N]\n"); 
		fprintf (stderr, "\t[-Q[<value>]] [-R<west/east/south/north>[r]] [-S] [-V] [-Z] [-:] [-bi[s][<n>]]\n\t[-bo[s][<n>]] [-f[i|o]<colinfo>]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<xyfile> is an multicolumn ASCII file with (lon,lat) in the first two columns\n");
		fprintf (stderr, "\t-G <grdfile> is the name of the 2-D binary data set to sample\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-L sets boundary conditions.  <flag> can be either\n");
		fprintf (stderr, "\t   g for geographic boundary conditions\n");
		fprintf (stderr, "\t   or one or both of\n");
		fprintf (stderr, "\t   x for periodic boundary conditions on x\n");
		fprintf (stderr, "\t   y for periodic boundary conditions on y\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N Report value at nearest node instead of interpolating\n");
		fprintf (stderr, "\t-Q Quick mode, use bilinear rather than bicubic interpolation.\n");
		fprintf (stderr, "\t   Optionally, append <value> in the 0 < value <= 1 range.\n");
		fprintf (stderr, "\t   [Default = 1 requires all 4 nodes to be non-NaN.], <value> = 0.5\n");
		fprintf (stderr, "\t   will interpolate about 1/2 way from a non-NaN to a NaN node, while\n");
		fprintf (stderr, "\t   0.1 will go about 90%% of the way, etc.\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t-S Suppress output when result equals NaN\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t-Z only output z-values [Default gives all columns]\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (bilinear && (threshold <= 0.0 || threshold > 1.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -Q:  threshold must be in <0,1] range\n", GMT_program);
		error++;
	}
	if (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify input file\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && 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] = 2;
        if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < 2) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least 2 columns\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */

	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]]);
	}

	pure_ascii = !(GMT_io.binary[GMT_IN] || GMT_io.binary[GMT_OUT]);

	if (fp == NULL) {
		fp = GMT_stdin;
		if (gmtdefs.verbose) fprintf (stderr, "%s: Reads from standard input\n", GMT_program);
#ifdef SET_IO_MODE
		GMT_setmode (GMT_IN);
#endif
	}

	if (GMT_read_grd_info (grdfile, &grd)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}
	
	if (west == east) {	/* No subset asked for */
		west = grd.x_min;
		east = grd.x_max;
		south = grd.y_min;
		north = grd.y_max;
	}
	one_or_zero = (grd.node_offset) ? 0 : 1;
	half = (grd.node_offset) ? 0.5 : 0.0;
	nx = irint ( (east - west) / grd.x_inc) + one_or_zero;
	ny = irint ( (north - south) / grd.y_inc) + one_or_zero;
	mx = nx + 4;	my = ny + 4;
	i_dx = 1.0 / grd.x_inc;
	i_dy = 1.0 / grd.y_inc;
	
	f = (float *) GMT_memory (VNULL, (size_t)(mx * my), sizeof (float), GMT_program);

	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;
	if (GMT_read_grd (grdfile, &grd, f, west, east, south, north, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grdfile);
		exit (EXIT_FAILURE);
	}

	GMT_boundcond_param_prep (&grd, &edgeinfo);
	
	project_info.w = west;	project_info.e = east;
	project_info.s = south;	project_info.n = north;
	
	/* Initialize bcr structure:  */

	GMT_bcr_init (&grd, GMT_pad, bilinear, threshold, &bcr);

	/* Set boundary conditions  */
	
	GMT_boundcond_set (&grd, &edgeinfo, GMT_pad, f);
	
	if (gmtdefs.io_header[GMT_IN]) {	/* First echo headers, if any */
		for (i = 0; i < gmtdefs.n_header_recs - 1; i++) {
			GMT_fgets (line, BUFSIZ, fp);
			if (!GMT_io.binary[GMT_OUT] && gmtdefs.io_header[GMT_OUT]) fprintf (GMT_stdout, "%s", line);
		}
		GMT_fgets (line, BUFSIZ, fp);
		line[strlen(line)-1] = 0;
		if (!GMT_io.binary[GMT_OUT] && gmtdefs.io_header[GMT_OUT]) fprintf (GMT_stdout, "%s\tsample\n", line);
	}

	ix = (gmtdefs.xy_toggle[1]);	iy = 1 - ix;	/* These are used for output purposes only */
	n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;

	while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {

		while (GMT_io.status & GMT_IO_SEGMENT_HEADER) {
			GMT_write_segmentheader (GMT_stdout, n_expected_fields);
			n_fields = GMT_input (fp,  &n_expected_fields, &in);
		}
		if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */

		if ((GMT_io.status & GMT_IO_MISMATCH) && n_fields < 2) {
			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);
		}
		if (n_output == 0) n_output = n_expected_fields + 1;

		n_read++;

		/* If point is outside grd area, shift it using periodicity or skip if not periodic. */

		while ( (in[1] < grd.y_min) && (edgeinfo.nyp > 0) ) in[1] += (grd.y_inc * edgeinfo.nyp);
		if (in[1] < grd.y_min) continue;

		while ( (in[1] > grd.y_max) && (edgeinfo.nyp > 0) ) in[1] -= (grd.y_inc * edgeinfo.nyp);
		if (in[1] > grd.y_max) continue;

		while ( (in[0] < grd.x_min) && (edgeinfo.nxp > 0) ) in[0] += (grd.x_inc * edgeinfo.nxp);
		if (in[0] < grd.x_min) continue;

		while ( (in[0] > grd.x_max) && (edgeinfo.nxp > 0) ) in[0] -= (grd.x_inc * edgeinfo.nxp);
		if (in[0] > grd.x_max) continue;
		
		if (node) {
			i = irint ((in[0] - grd.x_min) * i_dx - half) + one_or_zero;
			j = irint ((grd.y_max - in[1]) * i_dy - half) + one_or_zero;
			value = f[(j+GMT_pad[3])*mx+i+GMT_pad[0]];
		}
		else
			value = GMT_get_bcr_z(&grd, in[0], in[1], f, &edgeinfo, &bcr);

		if (suppress && GMT_is_dnan (value)) continue;

		if (!out) out = (double *) GMT_memory (VNULL, (size_t)n_output, sizeof (double), GMT_program);

		if (pure_ascii && n_expected_fields > 2) {
			/* Special case: Ascii i/o and at least 3 columns:
			   Columns beyond first two could be text strings */

			/* First get rid of any commas that may cause grief */
			for (i = 0; GMT_io.current_record[i]; i++) if (GMT_io.current_record[i] == ',') GMT_io.current_record[i] = ' ';
			sscanf (GMT_io.current_record, "%*f %*f %[^\n]", stuff);
			GMT_ascii_output_one (GMT_stdout, in[ix], ix);	fprintf (GMT_stdout, "\t");
			GMT_ascii_output_one (GMT_stdout, in[iy], iy);	fprintf (GMT_stdout, "\t");
			fprintf (GMT_stdout, "%s\t", stuff);
			GMT_ascii_output_one (GMT_stdout, value, 2);	fprintf (GMT_stdout, "\n");
		}
		else if (z_only) {	/* Simply print out value */
			GMT_output (GMT_stdout, 1, &value);
		}
		else {	/* Simply copy other columns, append value, and output */
			for (i = 0; i < n_expected_fields; i++) out[i] = in[i];
			out[i] = value;
			GMT_output (GMT_stdout, n_output, out);
		}

		n_points++;
	}
	GMT_fclose (fp);
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Sampled %d points from grid %s (%d x %d)\n", GMT_program,
		n_points, grdfile, grd.nx, grd.ny);
	
	GMT_free ((void *)f);
	GMT_free ((void *)out);
	
	GMT_end (argc, argv);
}
