/*--------------------------------------------------------------------
 *	$Id: grdproject.c,v 1.4.4.3 2002/02/27 17:41:10 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
 *--------------------------------------------------------------------*/
/*
 * grdproject reads a geographical grdfile and evaluates the grid at new grid positions
 * specified by the map projection and new dx/dy values using a weighted average of all
 * points within the search radius. Optionally, grdproject may perform the inverse
 * transformation, going from rectangular coordinates to geographical.
 *
 * Author:	Paul Wessel
 * Date:	15-JUL-2000
 * Ver:		3.4.1
 *
 */

#include "gmt.h"

float *geo, *rect;

main (int argc, char **argv)
{
	int i, dpi, nx, ny, unit = 0;
	
	BOOLEAN error = FALSE, inverse = FALSE, n_set = FALSE, set_n = FALSE, one_to_one = FALSE;
	BOOLEAN d_set = FALSE, e_set = FALSE, m_set = FALSE, map_center = FALSE, offset, toggle_offset = FALSE;
	
	double w, e, s, n, x_inc = 0.0, y_inc = 0.0, max_radius = 0.0, one_or_zero;
	double xmin, xmax, ymin, ymax, inch_to_unit, unit_to_inch, fwd_scale, inv_scale;
	
	char *infile, *outfile, format[BUFSIZ], unit_name[80], scale_unit_name[80];
	
	struct GRD_HEADER g_head, r_head;
	
	argc = GMT_begin (argc, argv);
	
	infile = outfile = CNULL;
	w = e = s = n = 0.0;
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */
			
				case 'J':
				case 'R':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], &w, &e, &s, &n);
					break;
				
				/* Supplemental parameters */
				
				case 'C':
					map_center = TRUE;
					break;
				case 'D':
					GMT_getinc (&argv[i][2], &x_inc, &y_inc);
					d_set = TRUE;
					break;
				case 'E':
					dpi = atoi (&argv[i][2]);
					e_set = TRUE;
					break;
				case 'A':
					one_to_one = TRUE;
					unit = GMT_check_scalingopt (argv[i], scale_unit_name);
					break;
				case 'F':
					toggle_offset = TRUE;
					break;
				case 'G':
					outfile = &argv[i][2];
					break;
				case 'I':
					inverse = TRUE;
					break;
				case 'M':	/* Directly specify units */
					GMT_set_measure_unit (argv[i]);
					m_set = TRUE;
					break;
				case 'N':
					sscanf (&argv[i][2], "%d/%d", &nx, &ny);
					if (ny == 0) ny = nx;
					n_set = TRUE;
					break;
				case 'S':
					max_radius = atof (&argv[i][2]);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else 
			infile = argv[i];
	}
	
	if ((d_set + e_set + n_set) == 0) n_set = set_n = TRUE;

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdproject %s - Project geographical grid to/from rectangular grid\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdproject <in_grdfile> -J<parameters> -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-A[k|m|n|i|c|p]] [-C] [-D<dx[m|c]>[/dy[m|c]]] [-E<dpi>] [-F] [-G<out_grdfile>] [-I] [-Mc|i|m]\n");
		fprintf (stderr, "\t[-N<nx/ny>] [-S<radius>] [-V]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<in_grdfile> is data set to be transformed\n");
		GMT_explain_option ('J');
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A force projected values to be in actual meters [Default uses the given map scale]\n");
		fprintf (stderr, "\t   Specify another unit by appending k (km), m (miles), n (nautical miles), i (inch), c (cm), or p (points)\n");
		fprintf (stderr, "\t-C coordinates relative to projection center [Default is relative to lower left corner]\n");
		fprintf (stderr, "\t-D sets the grid spacing for the new grid\n");
		fprintf (stderr, "\t-E sets dpi for output grid\n");
		fprintf (stderr, "\t-F toggle between pixel and grid registration  [Default is same as input]\n");
		fprintf (stderr, "\t-G name of output grid\n");
		fprintf (stderr, "\t-I Inverse transformation from rectangular to geographical\n");
		fprintf (stderr, "\t-M Temporarily reset MEASURE_UNIT to be c (cm), i (inch), m (meter), or p (point)\n");
		fprintf (stderr, "\t   Cannot be used if -A is set.\n");
		fprintf (stderr, "\t-N sets the number of nodes for the new grid\n");
		fprintf (stderr, "\t   Only one of -D, -E, and -N can be specified!\n");
		fprintf (stderr, "\t   If none are specified, nx,ny of the input file is used\n");
		fprintf (stderr, "\t-S sets the search radius in projected units [Default avoids aliasing]\n");
		GMT_explain_option ('V');
		
		exit (EXIT_FAILURE);
	}
	
	if (!infile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (!outfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", GMT_program);
		error++;
	}
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (max_radius < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S:  Max search radius, if specified, must be positive\n", GMT_program);
		error++;
	}
	if ((m_set + one_to_one) == 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Can specify only one of -A and -M\n", GMT_program);
		error++;
	}
	if ((d_set + e_set + n_set) != 1) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify only one of -D, -E, or -N\n", GMT_program);
		error++;
	}
	if (d_set && (x_inc <= 0.0 || y_inc <= 0.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -D option.  Must specify positive increment(s)\n", GMT_program);
		error++;
	}
	if (n_set && !set_n && (nx <= 0 || ny <= 0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  Must specify positive integers\n", GMT_program);
		error++;
	}
	if (e_set && dpi <= 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option.  Must specify positive dpi\n", GMT_program);
		error++;
	}
	
	if (error) exit (EXIT_FAILURE);

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

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

	GMT_map_setup (w, e, s, n);

	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 chosen units */
		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;
		if (unit) {	/* Change the 1:1 unit used */
			xmin *= fwd_scale;
			xmax *= fwd_scale;
			ymin *= fwd_scale;
			ymax *= fwd_scale;
		}
	}
	else {	/* Convert inches to chosen MEASURE */
		xmin *= inch_to_unit;
		xmax *= inch_to_unit;
		ymin *= inch_to_unit;
		ymax *= inch_to_unit;
	}

	GMT_grd_init (&r_head, argc, argv, FALSE);
	GMT_grd_init (&g_head, argc, argv, FALSE);

	sprintf (format, "(%s/%s/%s/%s)", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);

	if (inverse) {	/* Transforming from rectangular projection to geographical */
	
		/* if (!project_info.region) d_swap (s, e); */  /* Got w/s/e/n, make into w/e/s/n */

		g_head.x_min = w;	g_head.x_max = e;	g_head.y_min = s;	g_head.y_max = n;
	
		if (GMT_read_grd_info (infile, &r_head)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, infile);
			exit (EXIT_FAILURE);
		}
		rect = (float *) GMT_memory (VNULL, (size_t)(r_head.nx * r_head.ny), sizeof (float), GMT_program);
		if (GMT_read_grd (infile, &r_head, rect, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", infile, GMT_program);
			exit (EXIT_FAILURE);
		}
		offset = r_head.node_offset;		/* Same as input */
		if (toggle_offset) offset = !offset;	/* Toggle */
		one_or_zero = (offset) ? 0 : 1;
		if (set_n) {
			nx = r_head.nx;
			ny = r_head.ny;
		}
		GMT_grdproject_init (&g_head, x_inc, y_inc, nx, ny, dpi, offset);
		geo = (float *) GMT_memory (VNULL, (size_t)(g_head.nx * g_head.ny), sizeof (float), GMT_program);
		if (gmtdefs.verbose) {
			fprintf (stderr, "%s:  Transform ", GMT_program);
			fprintf (stderr, format, g_head.x_min, g_head.x_max, g_head.y_min, g_head.y_max);
			fprintf (stderr, " <-- ");
			fprintf (stderr, format, xmin, xmax, ymin, ymax);
			fprintf (stderr, " [%s]\n", unit_name);
		}

		/* Modify input rect header if -A, -C, -M have been set */
		
		if (one_to_one) {	/* Convert from 1:1 scale */
			if (unit) {	/* Undo the 1:1 unit used */
				r_head.x_min *= inv_scale;
				r_head.x_max *= inv_scale;
				r_head.y_min *= inv_scale;
				r_head.y_max *= inv_scale;
				max_radius *= inv_scale;
			}
			r_head.x_min *= project_info.x_scale;
			r_head.x_max *= project_info.x_scale;
			r_head.y_min *= project_info.y_scale;
			r_head.y_max *= project_info.y_scale;
			max_radius *= project_info.x_scale;
		}
		else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
			r_head.x_min *= unit_to_inch;
			r_head.x_max *= unit_to_inch;
			r_head.y_min *= unit_to_inch;
			r_head.y_max *= unit_to_inch;
			max_radius *= unit_to_inch;
		}
		if (map_center) {	/* Then correct so lower left corner is (0,0) */
			r_head.x_min += project_info.x0;
			r_head.x_max += project_info.x0;
			r_head.y_min += project_info.y0;
			r_head.y_max += project_info.y0;
		}
		r_head.x_inc = (r_head.x_max - r_head.x_min) / (r_head.nx - one_or_zero);
		r_head.y_inc = (r_head.y_max - r_head.y_min) / (r_head.ny - one_or_zero);
		
		/* rect xy values and max_radius are here in GMT projected inches */
		
		GMT_init_search_radius (&max_radius, &r_head, &g_head, TRUE);
		GMT_grd_inverse (geo, &g_head, rect, &r_head, max_radius);
			
		if (GMT_write_grd (outfile, &g_head, geo, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
			exit (EXIT_FAILURE);
		}
		
	}
	else {	/* Forward projection from geographical to rectangular grid */
	
		if (GMT_read_grd_info (infile, &g_head)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, infile);
			exit (EXIT_FAILURE);
		}

		geo = (float *) GMT_memory (VNULL, (size_t)(g_head.nx * g_head.ny), sizeof (float), GMT_program);
		if (GMT_read_grd (infile, &g_head, geo, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error reading file %s\n", GMT_program, infile);
			exit (EXIT_FAILURE);
		}
	
		r_head.x_min = project_info.xmin;	r_head.x_max = project_info.xmax;
		r_head.y_min = project_info.ymin;	r_head.y_max = project_info.ymax;
		if (one_to_one) {	/* Convert from 1:1 scale */
			if (unit) {	/* Undo the 1:1 unit used */
				x_inc *= inv_scale;
				y_inc *= inv_scale;
				max_radius *= inv_scale;
			}
			x_inc *= project_info.x_scale;
			y_inc *= project_info.y_scale;
			max_radius *= project_info.x_scale;
		}
		else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
			x_inc *= unit_to_inch;
			y_inc *= unit_to_inch;
			max_radius *= unit_to_inch;
		}
		if (set_n) {
			nx = g_head.nx;
			ny = g_head.ny;
		}
		
		if (gmtdefs.verbose) {
			fprintf (stderr, "%s:  Transform ", GMT_program);
			fprintf (stderr, format, g_head.x_min, g_head.x_max, g_head.y_min, g_head.y_max);
			fprintf (stderr, " --> ");
			fprintf (stderr, format, xmin, xmax, ymin, ymax);
			fprintf (stderr, " [%s]\n", unit_name);
		}

		offset = g_head.node_offset;		/* Same as input */
		if (toggle_offset) offset = !offset;	/* Toggle */
		one_or_zero = (offset) ? 0 : 1;

		GMT_grdproject_init (&r_head, x_inc, y_inc, nx, ny, dpi, offset);
		rect = (float *) GMT_memory (VNULL, (size_t)(r_head.nx * r_head.ny), sizeof (float), GMT_program);
		GMT_init_search_radius (&max_radius, &r_head, &g_head, FALSE);
		GMT_grd_forward (geo, &g_head, rect, &r_head, max_radius);
		
		/* Modify output rect header if -A, -C, -M have been set */

		if (map_center) {	/* Change origin from lower left to projection center */
			r_head.x_min -= project_info.x0;
			r_head.x_max -= project_info.x0;
			r_head.y_min -= project_info.y0;
			r_head.y_max -= project_info.y0;
		}
		if (one_to_one) {	/* Convert to 1:1 scale */
			r_head.x_min /= project_info.x_scale;
			r_head.x_max /= project_info.x_scale;
			r_head.y_min /= project_info.y_scale;
			r_head.y_max /= project_info.y_scale;
			if (unit) {	/* Change the 1:1 unit used */
				r_head.x_min *= fwd_scale;
				r_head.x_max *= fwd_scale;
				r_head.y_min *= fwd_scale;
				r_head.y_max *= fwd_scale;
			}
		}
		else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
			r_head.x_min *= unit_to_inch;
			r_head.x_max *= unit_to_inch;
			r_head.y_min *= unit_to_inch;
			r_head.y_max *= unit_to_inch;
		}
		r_head.x_inc = (r_head.x_max - r_head.x_min) / (r_head.nx - one_or_zero);
		r_head.y_inc = (r_head.y_max - r_head.y_min) / (r_head.ny - one_or_zero);

		/* rect xy values are here in GMT projected units chosen by user */

		if (GMT_write_grd (outfile, &r_head, rect, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error writing file %s\n", GMT_program, outfile);
			exit (EXIT_FAILURE);
		}
		
	}
	
	GMT_free ((void *)geo);
	GMT_free ((void *)rect);
	
	GMT_end (argc, argv);
}