/*--------------------------------------------------------------------
 *	$Id: grdtrack.c,v 1.5.4.3 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
 *--------------------------------------------------------------------*/
/*
 * 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:	3.4.1
 */

#include "gmt.h"

main (int argc, char **argv)
{
	int i, 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, z_only = FALSE, pure_ascii = FALSE;
	
	double value, west, east, south, north, *in, *out;
	
	float *f;

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

	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 '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;

				/* Supplemental parameters */
				
				case 'b':
					error += GMT_io_selection (&argv[i][2]);
					break;
				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_geographic_in = TRUE;
					}
					else
						GMT_geographic_in = TRUE;
					break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					break;
				case 'Q':
					bilinear = TRUE;
					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>]] [-Q]\n");
		fprintf (stderr, "[-Rw/s/e/n[r]] [-S] [-V] [-Z] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\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");
		fprintf (stderr, "\t   If no flag is given, x is assumed to be longitude and periodic in 360.\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-Q Quick mode, use bilinear rather than bicubic interpolation\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 ('.');
		exit (EXIT_FAILURE);
	}
	
	if (!grdfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify input file\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[0] && 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] = 2;
        if (GMT_io.binary[0] && GMT_io.ncol[0] < 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 .gmtcommands */

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

	pure_ascii = !(GMT_io.binary[0] || GMT_io.binary[1]);

	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 (0);
#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;
	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;
	
	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);

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

	ix = (gmtdefs.xy_toggle);	iy = 1 - ix;
	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : 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;
		
		value = GMT_get_bcr_z(&grd, in[0], in[1], f, &edgeinfo);

		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);
			fprintf (GMT_stdout, format1, in[ix], in[iy], stuff, value);
		}
		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);
}