/*--------------------------------------------------------------------
 *	$Id: mapproject.c,v 1.4.4.5 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
 *--------------------------------------------------------------------*/
/*
 * mapproject reads a pair of coordinates [+ optional data fields] from
 * standard input or file(s) and transforms the coordinates according to the
 * map projection selected. See man gmt-system for projections currently supported.
 *
 * The default is to expect longitude, latitude, [and optional datavalues],
 * and return x, y, [ and optional datavalues], but if the -I option is used,
 * the reverse is true.  Specifying -C means that the origin of projected coordinates
 * should be set to origin of projection  [Default origin is lower left corner of "map"].
 * If your data is lat lon instead of lon lat [Default], use -: to toggle x/y -> y/x.
 * Note that only unprojected values are affected by the -: switch.  True x,y values are
 * always printed out as x,y.
 *
 *
 * Author:	Paul Wessel
 * Date:	1-MAR-1990
 * Version:	3.4.1
 *
 */
 
#include "gmt.h"

main (int argc, char **argv)
{
	int i, k, fno, n = 0, n_read = 0, n_files = 0, x, y, n_args, unit = 0;
	int n_fields, n_expected_fields, *n_output;
	
	BOOLEAN error = FALSE, inverse = FALSE, suppress = FALSE, one_to_one = FALSE;
	BOOLEAN map_center = FALSE, nofile = TRUE, done = FALSE, first = TRUE;
	BOOLEAN pure_ascii = FALSE, d_set = FALSE;
	
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, x_in, y_in;
	double xmin, xmax, ymin, ymax, *in, *out, fwd_scale, inv_scale;
	double x_in_min, x_in_max, y_in_min, y_in_max, inch_to_unit, unit_to_inch;
	double x_out_min, x_out_max, y_out_min, y_out_max, u_scale;
	
	char stuff[BUFSIZ], line[BUFSIZ], format[BUFSIZ], format1[BUFSIZ], unit_name[80], scale_unit_name[80];
	
	FILE *fp = NULL;
	out = (double *)NULL;
	
	argc = GMT_begin (argc, argv);
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'H':
				case 'J':
				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 'C':
					map_center = TRUE;
					break;
				case 'D':
					d_set = TRUE;
					GMT_set_measure_unit (argv[i]);
					break;
				case 'F':
					one_to_one = TRUE;
					unit = GMT_check_scalingopt (argv[i], scale_unit_name);
					break;
				case 'I':
					inverse = TRUE;
					break;
				case 'M':               /* Multiple line segments */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'S':
					suppress = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "mapproject %s - Forward and Inverse map transformation of 2-D coordinates\n\n", GMT_VERSION);
		fprintf (stderr, "usage: mapproject <infiles> -J<parameters> -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-C] [-Dc|i|m|p] [-F[k|m|n|i|c|p]] [-H[<nrec>]] [-I] [-M[<flag>]] [-S] [-V] [-:] [-bi[s][<n>]] [-bo[s][<n>]]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tinfiles (in ASCII or binary) has 2 or more columns.  If no file(s) is given, standard input is read.\n");
		GMT_explain_option ('J');
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C returns x/y relative to projection center [Default is relative to lower left corner]\n");
		fprintf (stderr, "\t-D Temporarily reset MEASURE_UNIT to be c (cm), i (inch), m (meter), or p (point)\n");
		fprintf (stderr, "\t   Cannot be used if -F is set.\n");
		fprintf (stderr, "\t-F force projected values to be in actual meters [Default uses the given plot scale]\n");
		fprintf (stderr, "\t   Specify unit by appending k (km), m (miles), n (nautical miles), i (inch), c (cm), or p (points)\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I means Inverse, i.e., get lon/lat from x/y input. [Default is lon/lat -> x/y]\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-S means Suppress points outside region\n");
		GMT_explain_option ('V');
		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 (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\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 ((d_set + one_to_one) == 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Can specify only one of -D and -F\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]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (1);
#endif

	GMT_init_scales (unit, &fwd_scale, &inv_scale, &inch_to_unit, &unit_to_inch, unit_name);

	u_scale = (inverse) ? inv_scale : fwd_scale;

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

	GMT_map_setup (west, east, south, north);
		
	if (gmtdefs.verbose) {
		sprintf (format, "%s/%s/%s/%s", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		xmin = (map_center) ? project_info.xmin - project_info.x0 : project_info.xmin;
		xmax = (map_center) ? project_info.xmax - project_info.x0 : project_info.xmax;
		ymin = (map_center) ? project_info.ymin - project_info.y0 : project_info.ymin;
		ymax = (map_center) ? project_info.ymax - project_info.y0 : project_info.ymax;
		if (one_to_one) {	/* Convert to GMT inches */
			strncpy (unit_name, scale_unit_name, 80);
			xmin /= project_info.x_scale;
			xmax /= project_info.x_scale;
			ymin /= project_info.y_scale;
			ymax /= project_info.y_scale;
		}
		
		/* Convert inches to chosen MEASURE */
		xmin *= inch_to_unit;
		xmax *= inch_to_unit;
		ymin *= inch_to_unit;
		ymax *= inch_to_unit;

		fprintf (stderr, "%s:  Transform ", GMT_program);
		fprintf (stderr, format, project_info.w, project_info.e, project_info.s, project_info.n);
		(inverse) ? fprintf (stderr, " <- ") : fprintf (stderr, " -> ");
		fprintf (stderr, format, xmin, xmax, ymin, ymax);
		fprintf (stderr, " [%s]\n", unit_name);
	}
		
	/* Now we are ready to take on some input values */
	
	x = (gmtdefs.xy_toggle) ? 1 : 0;	y = 1 - x;		/* Set up which columns have x and y */
	if (MAPPING) {
		GMT_geographic_in = (inverse) ? FALSE : TRUE;
		GMT_geographic_out = !GMT_geographic_in;
	}
	else
		GMT_geographic_in = GMT_geographic_out = FALSE;

	n_expected_fields = (GMT_io.ncol[0]) ? GMT_io.ncol[0] : BUFSIZ;
	if (GMT_io.ncol[1])	/* Want binary output to be limitied to first ncol[1] columns */
		n_output = &GMT_io.ncol[1];
	else	/* Default to the number of input fields */
		n_output = &n_expected_fields;
	
	sprintf (format1, "%s\t%s\t%%s\n", gmtdefs.d_format, gmtdefs.d_format);
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
	
	x_in_min = y_in_min = x_out_min = y_out_min = DBL_MAX;
	x_in_max = y_in_max = x_out_max = y_out_max = -DBL_MAX;
	
	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 (0);
#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) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (line, BUFSIZ, fp);
				if (first && !GMT_io.binary[1]) fprintf (GMT_stdout, "%s", line);
			}
			first = FALSE;
		}
	
		if (inverse) {		/* Do inverse transformation */

			while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

				while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) {
					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;

				if (GMT_io.status & GMT_IO_MISMATCH && n_fields < 2) {
					fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d (skipped)\n", GMT_program, n_fields, n_expected_fields, n_read);
					continue;
				}
				if (!out) out = (double *) GMT_memory (VNULL, (size_t)n_expected_fields, sizeof (double), GMT_program);

				if (gmtdefs.verbose) {
					x_in = in[0];
					y_in = in[1];
				}
				if (one_to_one) {	/* Convert from 1:1 scale */
					if (unit) {
						in[0] *= u_scale;
						in[1] *= u_scale;
					}
					in[0] *= project_info.x_scale;
					in[1] *= project_info.y_scale;
				}
				else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from whatever to inch */
					in[0] *= unit_to_inch;
					in[1] *= unit_to_inch;
				}
				if (map_center) {	/* Then correct so lower left corner is (0,0) */
					in[0] += project_info.x0;
					in[1] += project_info.y0;
				}
				GMT_xy_to_geo (&out[0], &out[1], in[0], in[1]);
				n_read++;
				if (suppress && GMT_map_outside (out[0], out[1])) continue;
				if (gmtdefs.verbose) {
					x_in_min = MIN (x_in_min, x_in);
					x_in_max = MAX (x_in_max, x_in);
					y_in_min = MIN (y_in_min, y_in);
					y_in_max = MAX (y_in_max, y_in);
					x_out_min = MIN (x_out_min, out[0]);
					x_out_max = MAX (x_out_max, out[0]);
					y_out_min = MIN (y_out_min, out[1]);
					y_out_max = MAX (y_out_max, out[1]);
				}

				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, "%*s %*s %[^\n]", stuff);
					fprintf (GMT_stdout, format1, out[x], out[y], stuff);
				}
				else {	/* Simply copy other columns and output */
					for (k = 2; k < *n_output; k++) out[k] = in[k];
					GMT_output (GMT_stdout, *n_output, out);
				}
				n++;
			}
		}
		else {		/* Do forward transformation */

			while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && !(GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

				while ((GMT_io.status & GMT_IO_SEGMENT_HEADER) && !(GMT_io.status & GMT_IO_EOF)) {
					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;

				if (GMT_io.status & GMT_IO_MISMATCH) {
					fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d (skipped)\n", GMT_program, n_fields, n_expected_fields, n_read);
					continue;
				}
				if (!out) out = (double *) GMT_memory (VNULL, (size_t)n_expected_fields, sizeof (double), GMT_program);
				n_read++;

				/* Because -: is processed in GMT_input we use [0] for lon and [1] for y always */

				if (suppress && GMT_map_outside (in[0], in[1])) continue;
				GMT_geo_to_xy (in[0], in[1], &out[0], &out[1]);
				if (map_center) {	/* Change origin from lower left to projection center */
					out[0] -= project_info.x0;
					out[1] -= project_info.y0;
				}
				if (one_to_one) {	/* Convert to 1:1 scale */
					out[0] /= project_info.x_scale;
					out[1] /= project_info.y_scale;
					if (unit) {
						out[0] *= u_scale;
						out[1] *= u_scale;
					}
				}
				else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
					out[0] *= inch_to_unit;
					out[1] *= inch_to_unit;
				}
				if (gmtdefs.verbose) {
					x_in_min = MIN (x_in_min, in[0]);
					x_in_max = MAX (x_in_max, in[0]);
					y_in_min = MIN (y_in_min, in[1]);
					y_in_max = MAX (y_in_max, in[1]);
					x_out_min = MIN (x_out_min, out[0]);
					x_out_max = MAX (x_out_max, out[0]);
					y_out_min = MIN (y_out_min, out[1]);
					y_out_max = MAX (y_out_max, out[1]);
				}
				if (pure_ascii && n_expected_fields > 2) {
					/* Special case: Ascii input 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, "%*s %*s %[^\n]", stuff);
					fprintf (GMT_stdout, format1, out[0], out[1], stuff);
				}
				else {	/* Simply copy other columns and output */
					for (k = 2; k < *n_output; k++) out[k] = in[k];
					GMT_output (GMT_stdout, *n_output, out);
				}
				n++;
			}
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	if (gmtdefs.verbose && n_read > 0) {
		sprintf (format, "%%s: Input extreme values:  Xmin: %s Xmax: %s Ymin: %s Ymax %s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, x_in_min, x_in_max, y_in_min, y_in_max);
		sprintf (format, "%%s: Output extreme values:  Xmin: %s Xmax: %s Ymin: %s Ymax %s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, format, GMT_program, x_out_min, x_out_max, y_out_min, y_out_max);
		if (inverse)
			fprintf (stderr, "%s: Mapped %d x-y pairs [%s] to lon-lat\n", GMT_program, n, unit_name);
		else
			fprintf (stderr, "%s: Mapped %d lon-lat pairs to x-y [%s]\n", GMT_program, n, unit_name);
		if (suppress && n != n_read) fprintf (stderr, "%s: %d fell outside region\n", GMT_program, n_read - n);
	}

	GMT_free ((void *)out);

	GMT_end (argc, argv);
}