/*--------------------------------------------------------------------
 *	$Id: mapproject.c,v 1.42 2004/09/23 00:55:50 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
 *--------------------------------------------------------------------*/
/*
 * 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 the man page 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.  Option -G allows calculation of distances along track or
 * to a fixed point, while -L calculates shortest distances to a line.
 * Finally, datum conversions can also be done, alone or in series with a
 * map projection.
 *
 *
 * Author:	Paul Wessel
 * Date:	1-MAR-1990
 * Version:	4
 *
 */
 
#include "gmt.h"

PFD azimuth_func;
	
main (int argc, char **argv)
{
	int i, j, k, fno, n = 0, n_read = 0, n_files = 0, x, y, n_args, unit = 0, n_slash, g_report = 0;
	int n_fields, n_expected_fields, distance = 0, proj_type = 0, *n_output, save[2], n_lines, two;
	
	BOOLEAN error = FALSE, inverse = FALSE, suppress = FALSE, one_to_one = FALSE, ECEF_conv = FALSE;
	BOOLEAN map_center = FALSE, nofile = TRUE, done = FALSE, first = TRUE, datum_conv = FALSE;
	BOOLEAN pure_ascii = FALSE, d_set = FALSE, line_start = TRUE, do_az = FALSE, long_verbose = FALSE;
	BOOLEAN geodetic_calc = FALSE, do_geocentric, do_line_dist = FALSE, greenwich = FALSE;
	BOOLEAN datum_conv_only = FALSE, back_az = FALSE, double_whammy = FALSE, shift_xy = FALSE;
	
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, x_in, y_in, d, s, x0, y0;
	double xmin, xmax, ymin, ymax, *in, *out, fwd_scale, inv_scale, xtmp, ytmp;
	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, d_scale;
	double false_easting = 0.0, false_northing = 0.0;
	
	char line[BUFSIZ], line_file[BUFSIZ], format[BUFSIZ], unit_name[80], scale_unit_name[80];
	char txt_a[32], txt_b[32], c, *p;
	
	FILE *fp = NULL;

	struct GMT_LINES *xyline;

	PFD distance_func;
	PFI near_a_line;

	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 'V':
					if (argv[i][2] == 'l') long_verbose = TRUE;
				case 'R':
				case 'H':
				case 'J':
				case ':':
				case 'b':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */
				
				case 'A':
					do_az = TRUE;
	 				n = sscanf (&argv[i][2], "%c%[^/]/%s", &c, txt_a, txt_b);
					if (n < 3) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR.  Expected -Ab|B|f|F<lon0>/<lat0>\n", GMT_program);
						error++;
					}
					else {
						do_geocentric = back_az = FALSE;
						switch (c) {
							case 'B':
								do_geocentric = TRUE;
							case 'b':
								back_az = TRUE;
								break;
							case 'F':
								do_geocentric = TRUE;
							case 'f':
								break;
							default:
								fprintf (stderr, "%s: GMT SYNTAX ERROR.  Expected -Ab|B|f|F<lon0>/<lat0>\n", GMT_program);
								error++;
								break;
						}
						error += GMT_verify_expectations (GMT_io.in_col_type[0], GMT_scanf_arg (txt_a, GMT_io.in_col_type[0], &x0), txt_a);
						error += GMT_verify_expectations (GMT_io.in_col_type[1], GMT_scanf_arg (txt_b, GMT_io.in_col_type[1], &y0), txt_b);
					}
					break;
				case 'C':
					map_center = TRUE;
					if (argv[i][2]) {	/* Also gave shifts */
	 					n = sscanf (&argv[i][2], "%lf/%lf", &false_easting, &false_northing);
						if (n != 2) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR.  Expected -C[<false_easting>/<false_northing>]\n", GMT_program);
							error++;
						}
						shift_xy = TRUE;
					}
					break;
				case 'D':
					d_set = TRUE;
					GMT_set_measure_unit (argv[i]);
					break;
				case 'E':
					ECEF_conv = TRUE;
					if (GMT_ECEF_init (&argv[i][2])) {
						error++;
						fprintf (stderr, "%s: GMT SYNTAX ERROR -E: Usage -E[<datum>]\n", GMT_program);
					}
					break;
				case 'F':
					one_to_one = TRUE;
					unit = GMT_check_scalingopt (argv[i], scale_unit_name);
					break;
				case 'G':
					for (n_slash = 0, k = 2; argv[i][k]; k++) if (argv[i][k] == '/') n_slash++;
					if (n_slash == 2 || n_slash == 1) {	/* Got -Glon0/lat0[/units] */
						distance = 1;
	 					n = sscanf (&argv[i][2], "%[^/]/%[^/]/%c", txt_a, txt_b, &c);
						if (n < 2) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR.  Expected -G<lon0>/<lat0>[/e|E|k|K|m|M|n|N|c|C|d|D]\n", GMT_program);
							error++;
						}
						error += GMT_verify_expectations (GMT_io.in_col_type[0], GMT_scanf_arg (txt_a, GMT_io.in_col_type[0], &x0), txt_a);
						error += GMT_verify_expectations (GMT_io.in_col_type[1], GMT_scanf_arg (txt_b, GMT_io.in_col_type[1], &y0), txt_b);
						if (n_slash == 2)
							error += GMT_get_dist_scale (c, &d_scale, &proj_type, &distance_func);
						else
							error += GMT_get_dist_scale ('\0', &d_scale, &proj_type, &distance_func);
					}
					else {				/* Got -G[units] */
						distance = 2;
						error += GMT_get_dist_scale (argv[i][2], &d_scale, &proj_type, &distance_func);
					}
					break;
				case 'I':
					inverse = TRUE;
					break;
				case 'L':
	 				strcpy (line_file, &argv[i][2]);
					k = strlen (line_file) - 2;
					if (line_file[k] == '/' && strchr ("ekmndcC", line_file[k+1])) {
						error += GMT_get_dist_scale (line_file[k+1], &d_scale, &proj_type, &distance_func);
						line_file[k] = 0;
					}
					else
						error += GMT_get_dist_scale ('\0', &d_scale, &proj_type, &distance_func);
					do_line_dist = TRUE;
					break;
				case 'Q':
					if (argv[i][2] == 'e') g_report += 1;
					if (argv[i][2] == 'd') g_report += 2;
					if (argv[i][2] == '\0') g_report = 3;
					break;
				case 'M':               /* Multiple line segments */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'S':
					suppress = TRUE;
					break;
				case 'T':
					datum_conv = TRUE;
					if (GMT_datum_init (&argv[i][2])) {
						error++;
						fprintf (stderr, "%s: GMT SYNTAX ERROR -T: Usage -T[h]<from>[/<to>]\n", GMT_program);
					}
					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 transformations and geodesy\n\n", GMT_VERSION);
		fprintf (stderr, "usage: mapproject <infiles> -J<parameters> -R<west/east/south/north> [-C[<dx/dy>]]\n");
		fprintf (stderr, "\t[-Ab|B|f|F<lon0/lat0>] [-Dc|i|m|p] [-E[<datum>]] [-F[k|m|n|i|c|p]] [-G[<lon0/lat0>/][<unit>]\n");
		fprintf (stderr, "\t[-H[<nrec>]] [-I] [-L<line.xy>[/<unit>] [-M[<flag>]] [-S] [-T[h]<from>[/<to>]\n");
		fprintf (stderr, "\t[-V[l]] [-:] [-bi[s][<n>]] [-bo[s][<n>]] [-f[i|o]<colinfo>]\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-A Calculate azimuths from specified point to each input data point with -Af.\n");
		fprintf (stderr, "\t   Use -Ab to calculate backazimuths from data to specified point.\n");
		fprintf (stderr, "\t   Upper case B or F gives azimuths of geodesics using current ellipsoid.\n");
		fprintf (stderr, "\t-C returns x/y relative to projection center [Default is relative to lower left corner]\n");
		fprintf (stderr, "\t   Optionally append dx/dy to add (or subtract if -I) (i.e., false easting & northing) [0/0]\n");
		fprintf (stderr, "\t   Units are plot units unless -F is set in which case the unit is meters.\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-E Convert (lon, lat, h) to Earth Centered Earth Fixed (ECEF) coordinates [-I for inverse].\n");
		fprintf (stderr, "\t   Specify <datum> using datum ID (see -Qd or man page) or as <ellipsoid>:<dx,dy,dz>\n");
		fprintf (stderr, "\t   where <ellipsoid> may be ellipsoid ID (see -Qe or man page) or <major,inv_flattening>.\n");
		fprintf (stderr, "\t   If <datum> = - or not given we assume WGS-84.\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");
		fprintf (stderr, "\t-G Calculate distances to specified point OR culumlative distances along track (if point not given).\n");
		fprintf (stderr, "\t   Specify unit as m(e)ter, (k)m, (m)ile, (n)autical mile, (d)egree, or (c)artesian in user units.\n");
		fprintf (stderr, "\t   Unit C means Cartesian distances after first projecting the input coordinates (-R, -J).\n");
		fprintf (stderr, "\t   Units E, K, M, N, D mean geodesic distance using current ellipsoid [lower case is spherical].\n");
		fprintf (stderr, "\t   Default is meters on spherical earth.\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");
		fprintf (stderr, "\t-L Calculate minimum distances to specified line(s) in the file <line.xy>.\n");
		fprintf (stderr, "\t   Specify unit as m(e)ter, (k)m, (m)ile, (n)autical mile, (d)egree, or (c)artesian in user units.\n");
		fprintf (stderr, "\t   Unit C means Cartesian distances after first projecting the input coordinates (-R, -J).\n");
		fprintf (stderr, "\t   Calculations uses spherical approximations.  Default unit is meters.\n");
		GMT_explain_option ('M');
		fprintf (stderr, "\t-Q list projection parameters and exit.  For subsets [Default is all] use\n");
		fprintf (stderr, "\t   -Qe shows ellipsoid parameters\n");
		fprintf (stderr, "\t   -Qd shows datum parameters\n");
		fprintf (stderr, "\t-S means Suppress points outside region\n");
		fprintf (stderr, "\t-T means coordinate transformation from datum <from> to datum <to>.\n");
		fprintf (stderr, "\t   Prepend h if input data are lon, lat, height [Default sets height = 0].\n");
		fprintf (stderr, "\t   Specify datums using datum ID (see -Qd or man page) or as <ellipsoid>:<dx,dy,dz>\n");
		fprintf (stderr, "\t   where <ellipsoid> may be ellipsoid ID (see -Qe or man page) or <major,inv_flattening>.\n");
		fprintf (stderr, "\t   <from> = - means WGS-84.  If /<to> is not given we assume WGS-84.\n");
		fprintf (stderr, "\t   -T can be used as pre- or post- (-I) processing for -J -R.\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t   Append l for long verbose, reporting every 1000 points.\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 (g_report & 1) {	/* List ellipsoid parameters */
		fprintf (stderr, "GMT supports %d ellipsoids, given below (-> indicates default setting)\n", N_ELLIPSOIDS);
		fprintf (stderr, "  ID                      Date        a         1/f\n");
		fprintf (stderr, "--------------------------------------------------------\n");
		for (i = 0; i < N_ELLIPSOIDS; i++) {
			(i == gmtdefs.ellipsoid) ? fprintf (stderr, "->") : fprintf (stderr, "  ");
			fprintf (stderr, "%-23s %4d %.3lf %.9lf\n", gmtdefs.ref_ellipsoid[i].name, gmtdefs.ref_ellipsoid[i].date, gmtdefs.ref_ellipsoid[i].eq_radius, 1.0/gmtdefs.ref_ellipsoid[i].flattening);
		}
		fprintf (stderr, "--------------------------------------------------------\n");
	}
	if (g_report & 2) {	/* List datum parameters */
		fprintf (stderr, "GMT supports %d datums, given below (-> indicates default setting)\n", N_DATUMS);
		fprintf (stderr, "  ID  Name                               Ellipsoid                 x     y     z   Region\n");
		fprintf (stderr, "-----------------------------------------------------------------------------------------\n");
		for (i = 0; i < N_DATUMS; i++) {
			(!strcmp (gmtdefs.datum[i].name, "WGS 1984")) ? fprintf (stderr, "->") : fprintf (stderr, "  ");
			fprintf (stderr, "%3d %-34s %-23s %5.0lf %5.0lf %5.0lf %s\n", i, gmtdefs.datum[i].name, gmtdefs.datum[i].ellipsoid, gmtdefs.datum[i].xyz[0], gmtdefs.datum[i].xyz[1], gmtdefs.datum[i].xyz[2], gmtdefs.datum[i].region);
		}
		fprintf (stderr, "-----------------------------------------------------------------------------------------\n");
	}
	if (g_report) exit (EXIT_FAILURE);

	geodetic_calc = (distance || do_az || do_line_dist);
	
	if (datum_conv && (distance + ECEF_conv + do_line_dist) > 0) {	/* No good... */
		fprintf (stderr, "%s: GMT SYNTAX ERROR: -T cannot work with -E, -G or -L\n", GMT_program);
		error++;
	}
	if (geodetic_calc && inverse) {	/* No good... */
		fprintf (stderr, "%s: GMT SYNTAX ERROR: -A, -G, and -L cannot work with -I\n", GMT_program);
		error++;
	}
	if (project_info.projection < 0 && (distance || do_line_dist) && proj_type == 2) {	/* Must have -J */
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -J option with selected form of -G or -L\n", GMT_program);
		error++;
	}
	if (!project_info.region_supplied && !(geodetic_calc || datum_conv || ECEF_conv)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\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 ((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[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 (((datum_conv && GMT_datum.h_given) || ECEF_conv) && GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < 3) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  For -E or -T, binary input data (-bi) must have at least 3 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]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (GMT_OUT);
#endif
	if (datum_conv && project_info.projection && project_info.region_supplied) {	/* Do datum shift & project coordinates */
		double_whammy = TRUE;
		if (inverse) {	/* Need to set the ellipsoid to that of the old datum */
			if (GMT_datum.from.ellipsoid_id < 0) {
				gmtdefs.ellipsoid = N_ELLIPSOIDS - 1;
				gmtdefs.ref_ellipsoid[i].eq_radius = GMT_datum.from.a;
				gmtdefs.ref_ellipsoid[i].flattening = GMT_datum.from.f;
			}
			else
				gmtdefs.ellipsoid = GMT_datum.from.ellipsoid_id;
		}
		else {	/* Need to set the ellipsoid to that of the new datum */
			if (GMT_datum.to.ellipsoid_id < 0) {
				gmtdefs.ellipsoid = N_ELLIPSOIDS - 1;
				gmtdefs.ref_ellipsoid[i].eq_radius = GMT_datum.to.a;
				gmtdefs.ref_ellipsoid[i].flattening = GMT_datum.to.f;
			}
			else
				gmtdefs.ellipsoid = GMT_datum.to.ellipsoid_id;
		}
	}
	else
		datum_conv_only = datum_conv;

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

	if (distance) {	/* save output format in case -J changes it */
		save[0] = GMT_io.out_col_type[0];
		save[1] = GMT_io.out_col_type[1];
	}
	u_scale = (inverse) ? inv_scale : fwd_scale;

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

	if (project_info.projection < 0) {	/* Supply dummy linear proj */
		GMT_map_getproject ("x1d");
		if (!project_info.region_supplied) {
			west = 0.0;	east = 360.0;
			south = -90.0;	north = 90.0;
		}
	}
	GMT_map_setup (west, east, south, north);
	
	azimuth_func = (SPHERICAL || !do_geocentric) ? GMT_az_backaz_sphere : GMT_az_backaz_geodesic;
	 
	if (distance && proj_type < 2) {	/* Ensure we use the selected output coordinates */
		GMT_io.out_col_type[0] = save[0];
		GMT_io.out_col_type[1] = save[1];
	}

	if (gmtdefs.verbose && !(geodetic_calc || datum_conv)) {
		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;
		
		if (shift_xy) {
			xmin += false_easting;
			xmax += false_easting;
			ymin += false_northing;
			ymax += false_northing;
		}

		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);
	}
		
	if (GMT_io.in_col_type[0] & GMT_IS_GEO && proj_type == 0) {	/* Geographic data */
		GMT_distance_func = (PFD) GMT_great_circle_dist;
		near_a_line = (PFI) GMT_near_a_line_spherical;
		greenwich = (west < 0.0 && east > 0.0);
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
		near_a_line = (PFI) GMT_near_a_line_cartesian;
	}
	if (do_line_dist) {
		n_lines  = GMT_lines_init (line_file, &xyline, 0.0, greenwich);
		if (proj_type == 2) {	/* Must convert the line points first */
			for (i = 0; i < n_lines; i++) {
				for (j = 0; j < xyline[i].np; j++) {
					GMT_geo_to_xy (xyline[i].lon[j], xyline[i].lat[j], &xtmp, &ytmp);
					xyline[i].lon[j] = xtmp;
					xyline[i].lat[j] = ytmp;
				}
			}
		}
	}

	/* Now we are ready to take on some input values */
	
	x = (gmtdefs.xy_toggle[1]) ? 1 : 0;	y = 1 - x;		/* Set up which columns have x and y for output only*/
	if ((MAPPING || ECEF_conv) && inverse) {
		GMT_io.out_col_type[0] = GMT_IS_LON;	GMT_io.out_col_type[1] = GMT_IS_LAT;	/* Inverse projection expects x,y and gives lon, lat */
		GMT_io.in_col_type[0] = GMT_io.in_col_type[1] = GMT_IS_FLOAT;
	}
	if (datum_conv_only) {	/* Both in and out are geographic */
		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;
		GMT_io.in_col_type[2] = GMT_io.out_col_type[2] = GMT_IS_FLOAT;
	}

	n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;
	if (GMT_io.ncol[GMT_OUT])	/* Want binary output to be limitied to first ncol[1] columns */
		n_output = &GMT_io.ncol[GMT_OUT];
	else	/* Default to the number of input fields */
		n_output = &n_expected_fields;
		
	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;
	
	two = (ECEF_conv || (datum_conv && GMT_datum.h_given)) ? 3 : 2;	/* # of output points from conversion */
	
	if (shift_xy && one_to_one) {	/* Use same units in -C and -F */
		false_easting *= u_scale;
		false_northing *= u_scale;		
	}
	if (distance == 2 && proj_type == 2) {	/* Must project the fixed point here */
		GMT_geo_to_xy (x0, y0, &xtmp, &ytmp);
		if (map_center) {	/* Change origin from lower left to projection center */
			xtmp -= project_info.x0;
			ytmp -= project_info.y0;
		}
		if (one_to_one) {	/* Convert to 1:1 scale */
			xtmp /= project_info.x_scale;
			ytmp /= project_info.y_scale;
			if (unit) {
				xtmp *= u_scale;
				ytmp *= u_scale;
			}
		}
		else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
			xtmp *= inch_to_unit;
			ytmp *= inch_to_unit;
		}
		if (shift_xy) {
			xtmp += false_easting;
			ytmp += false_northing;
		}
		x0 = xtmp;
		y0 = ytmp;
	}
	
	n = 0;
	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 (GMT_IN);
#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[GMT_IN]) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (line, BUFSIZ, fp);
				if (first && !GMT_io.binary[GMT_OUT] && gmtdefs.io_header[GMT_OUT]) 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 (shift_xy) {
					in[0] -= false_easting;
					in[1] -= false_northing;
				}
				if (ECEF_conv) {
					GMT_ECEF_inverse (in, out);
				}
				else {
					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]);
				}
				if (double_whammy) {	/* Now apply datum shift */
					in[0] = out[0];
					in[1] = out[1];
					GMT_conv_datum (in, out);
				}
					
				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 */

					/* We will use strtok to step past the first 2 [or 3] columns.  The remainder
					 * will then be the user text that we want to preserve.  Since strtok places
					 * 0 to indicate start of next token we count our way to the start of the text. */
					
					strcpy (line, GMT_io.current_record);
					GMT_chop (line);
					p = strtok (line,  " \t,");	/* Returns xstring */
					p = strtok (CNULL, " \t,");	/* Returns ystring */
					if (ECEF_conv) {
						p = strtok (CNULL, " \t,");	/* Returns zstring */
						for (i = k = 0; k < 3; i++) if (line[i] == '\0') k++;	/* Find start of user text */
						GMT_ascii_output_one (GMT_stdout, out[x], 0);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						GMT_ascii_output_one (GMT_stdout, out[y], 1);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						GMT_ascii_output_one (GMT_stdout, out[2], 2);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						fprintf (GMT_stdout, "%s\n", &line[i]);
					}
					else {
						for (i = k = 0; k < 2; i++) if (line[i] == '\0') k++;	/* Find start of user text */
						GMT_ascii_output_one (GMT_stdout, out[x], 0);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						GMT_ascii_output_one (GMT_stdout, out[y], 1);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						fprintf (GMT_stdout, "%s\n", &line[i]);
					}
				}
				else {	/* Simply copy other columns and output */
					for (k = two; k < *n_output; k++) out[k] = in[k];
					GMT_output (GMT_stdout, *n_output, out);
				}
				n++;
				if (long_verbose && (n%1000) == 0) fprintf (stderr, "%s: Projected %d points\r", GMT_program, 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);
					line_start = TRUE;
				}
				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+1, 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;
				
				if (datum_conv_only) {
					GMT_conv_datum (in, out);
				}
				else if (ECEF_conv) {
					GMT_ECEF_forward (in, out);
				}
				else {
					if (double_whammy) {	/* Apply datum shift first */
						GMT_conv_datum (in, out);
						in[0] = out[0];
						in[1] = out[1];
					}
					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 (shift_xy) {
						out[0] += false_easting;
						out[1] += false_northing;
					}
				}
				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]);
				}
				if (geodetic_calc) {	/* Get either distances or azimuths */
					if (distance) {	/* Cumulative distances along track */
						if (distance == 2 && line_start)
							s = d = 0.0;
						else if (proj_type == 2)	/* Calculate Cartesian distances using projected units */
							s = hypot (x0 - out[0], y0 - out[1]);
						else if (proj_type == 1)	/* Plain Cartesian distances using input points */
							s = hypot (x0 - in[0], y0 - in[1]);
						else				/* Great circle distances */
							s = d_scale * (*distance_func) (x0, y0, in[0], in[1]);
						if (distance == 2) {
							line_start = FALSE;
							d += s;
							if (proj_type == 2) {	/* Calculate distances using projected units */
								x0 = out[0];
								y0 = out[1];
							}
							else {
								x0 = in[0];
								y0 = in[1];
							}
						}
						else
							d = s;
					}
					else if (do_line_dist) {	/* Compute closest distance to line */
						if (proj_type == 2)	/* Using projected coordinates */
							(void *) near_a_line (out[0], out[1], xyline, n_lines, TRUE, &d);
						else			/* Using input coordinates */
							(void *) near_a_line (in[0], in[1], xyline, n_lines, TRUE, &d);
						d *= d_scale;
					}
					else {	/* Azimuths */
						d = (*azimuth_func) (x0, y0, in[0], in[1], back_az);
					}
					if (pure_ascii && n_expected_fields > 2) {
						/* Special case: Ascii input and at least 3 columns:
					 	  Columns beyond first two could be text strings */

						/* We will use strtok to step past the first 2 [or 3] columns.  The remainder
						 * will then be the user text that we want to preserve.  Since strtok places
						 * 0 to indicate start of next token we count our way to the start of the text. */
					
						strcpy (line, GMT_io.current_record);
						GMT_chop (line);
						p = strtok (line,  " \t,");	/* Returns xstring */
						p = strtok (CNULL, " \t,");	/* Returns ystring */
						for (i = k = 0; k < 2; i++) if (line[i] == '\0') k++;	/* Find start of user text */

						GMT_ascii_output_one (GMT_stdout, in[x], 0);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						GMT_ascii_output_one (GMT_stdout, in[y], 1);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						fprintf (GMT_stdout, "%s", &line[i]);
						fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
						GMT_ascii_output_one (GMT_stdout, d, 2);
						fprintf (GMT_stdout, "\n");
					}
					else {	/* Simply copy other columns and output */
						for (k = 0; k < *n_output; k++) out[k] = in[k];
						out[*n_output] = d;
						
						GMT_output (GMT_stdout, *n_output + 1, out);
					}
				}
				else {
					if (gmtdefs.verbose) {
						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 */

						/* We will use strtok to step past the first 2 [or 3] columns.  The remainder
						 * will then be the user text that we want to preserve.  Since strtok places
						 * 0 to indicate start of next token we count our way to the start of the text. */
					
						strcpy (line, GMT_io.current_record);
						GMT_chop (line);
						p = strtok (line,  " \t,");	/* Returns xstring */
						p = strtok (CNULL, " \t,");	/* Returns ystring */
						if (ECEF_conv || (datum_conv && GMT_datum.h_given)) {
							p = strtok (CNULL, " \t,");	/* Returns zstring */
							for (i = k = 0; k < 3; i++) if (line[i] == '\0') k++;	/* Find start of user text */
							GMT_ascii_output_one (GMT_stdout, out[x], 0);
							fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							GMT_ascii_output_one (GMT_stdout, out[y], 1);
							fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							GMT_ascii_output_one (GMT_stdout, out[2], 2);
							fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							fprintf (GMT_stdout, "%s\n", &line[i]);
						}
						else {
							for (i = k = 0; k < 2; i++) if (line[i] == '\0') k++;	/* Find start of user text */
							GMT_ascii_output_one (GMT_stdout, out[x], 0);
							fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							GMT_ascii_output_one (GMT_stdout, out[y], 1);
							fprintf (GMT_stdout, "%s", gmtdefs.field_delimiter);
							fprintf (GMT_stdout, "%s\n", &line[i]);
						}
					}
					else {	/* Simply copy other columns and output */
						for (k = two; k < *n_output; k++) out[k] = in[k];
						GMT_output (GMT_stdout, *n_output, out);
					}
				}
				n++;
				if (long_verbose && (n%1000) == 0) fprintf (stderr, "%s: Projected %d points\r", GMT_program, n);
			}
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	if (gmtdefs.verbose && n_read > 0) {
		fprintf (stderr, "%s: Projected %d points\n", GMT_program, n);
		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);
		if (!geodetic_calc) {
			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) {
				if (ECEF_conv)
					fprintf (stderr, "%s: Mapped %d ECEF coordinates [m] to (lon,lat,h)\n", GMT_program, n);
				else
					fprintf (stderr, "%s: Mapped %d x-y pairs [%s] to lon-lat\n", GMT_program, n, unit_name);
			}
			else if (datum_conv && GMT_datum.h_given)
				fprintf (stderr, "%s: Datum-converted %d (lon,lat,h) triplets\n", GMT_program, n);
			else if (datum_conv)
				fprintf (stderr, "%s: Datum-converted %d (lon,lat) pairs\n", GMT_program, n);
			else if (ECEF_conv)
				fprintf (stderr, "%s: Mapped %d (lon,lat,h) triplets to ECEF coordinates [m]\n", GMT_program, n);
			else if (MAPPING)
				fprintf (stderr, "%s: Mapped %d lon-lat pairs to x-y [%s]\n", GMT_program, n, unit_name);
			else
				fprintf (stderr, "%s: Mapped %d data 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);
}
