/*--------------------------------------------------------------------
 *	$Id: grdproject.c,v 1.28 2006/04/13 08:34:41 pwessel Exp $
 *
 *	Copyright (c) 1991-2006 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:		4
 *
 */

#include "gmt.h"

struct GRDPROJECT_CTRL {
	struct A {	/* -A[k|m|n|i|c|p] */
		BOOLEAN active;
		char unit;
	} A;
	struct C {	/* -C[<dx/dy>] */
		BOOLEAN active;
		double easting, northing;
	} C;
	struct D {	/* -Ddx[/dy] */
		BOOLEAN active;
		double xinc, yinc;
	} D;
	struct E {	/* -E<dpi> */
		BOOLEAN active;
		int dpi;
	} E;
	struct F {	/* -F */
		BOOLEAN active;
	} F;
	struct G {	/* -G */
		BOOLEAN active;
		char *file;
	} G;
	struct I {	/* -I */
		BOOLEAN active;
	} I;
	struct M {	/* -Mc|i|m */
		BOOLEAN active;
		char unit;
	} M;
	struct N {	/* -N<nx/ny> */
		BOOLEAN active;
		int nx, ny;
	} N;
	struct S {	/* -S<radius> */
		BOOLEAN active;
		BOOLEAN bilinear;
		BOOLEAN careful;
		double radius;
	} S;
};

int main (int argc, char **argv)
{
	BOOLEAN error = FALSE, set_n = FALSE, shift_xy = FALSE, offset;

	char *infile, format[BUFSIZ], unit_name[80], scale_unit_name[80];

	int i, unit = 0;
	
	size_t nm;
	
	float *geo, *rect;
	
	double w, e, s, n;
	double xmin, xmax, ymin, ymax, inch_to_unit, unit_to_inch, fwd_scale, inv_scale;

	struct GRD_HEADER g_head, r_head;
	struct GMT_EDGEINFO edgeinfo;
	struct GRDPROJECT_CTRL *Ctrl;

	void *New_Grdproject_Ctrl (), Free_Grdproject_Ctrl (struct GRDPROJECT_CTRL *C);

	argc = GMT_begin (argc, argv);

	Ctrl = (struct GRDPROJECT_CTRL *)New_Grdproject_Ctrl ();	/* Allocate and initialize a new control structure */
	
	infile = 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_parse_common_options (argv[i], &w, &e, &s, &n);
					break;

				/* Supplemental parameters */

				case 'C':
					Ctrl->C.active = TRUE;
					if (argv[i][2]) {	/* Also gave shifts */
	 					n = sscanf (&argv[i][2], "%lf/%lf", &Ctrl->C.easting, &Ctrl->C.northing);
						if (n != 2) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR.  Expected -C[<false_easting>/<false_northing>]\n", GMT_program);
							error++;
						}
					}
					break;
				case 'D':
					Ctrl->D.active = TRUE;
					if (GMT_getinc (&argv[i][2], &Ctrl->D.xinc, &Ctrl->D.yinc)) {
						GMT_inc_syntax ('D', 1);
						error = TRUE;
					}
					break;
				case 'E':
					Ctrl->E.active = TRUE;
					Ctrl->E.dpi = atoi (&argv[i][2]);
					break;
				case 'A':
					Ctrl->A.active = TRUE;
					Ctrl->A.unit = argv[i][2];
					break;
				case 'F':
					Ctrl->F.active = TRUE;
					break;
				case 'G':
					Ctrl->G.file = strdup (&argv[i][2]);
					break;
				case 'I':
					Ctrl->I.active = TRUE;
					break;
				case 'M':	/* Directly specify units */
					Ctrl->M.active = TRUE;
					Ctrl->M.unit = argv[i][2];
					break;
				case 'N':
					sscanf (&argv[i][2], "%d/%d", &Ctrl->N.nx, &Ctrl->N.ny);
					if (Ctrl->N.ny == 0) Ctrl->N.ny = Ctrl->N.nx;
					Ctrl->N.active = TRUE;
					break;
				case 'S':
					Ctrl->S.active = TRUE;
					if (argv[i][2] == '+')
						Ctrl->S.careful = TRUE;
					else if (argv[i][2] == '-') {
						Ctrl->S.careful = TRUE;
						Ctrl->S.bilinear = TRUE;
					}
					else
						Ctrl->S.radius = atof (&argv[i][2]);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else 
			infile = argv[i];
	}

	if ((Ctrl->D.active + Ctrl->E.active + Ctrl->N.active) == 0) Ctrl->N.active = set_n = TRUE;

	if (argc == 1 || GMT_give_synopsis_and_exit) {
		fprintf (stderr, "grdproject %s - Project geographical grid to/from rectangular grid\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdproject <in_grdfile> %s %s\n", GMT_J_OPT, GMT_Rgeo_OPT);
		fprintf (stderr, "\t[-A[k|m|n|i|c|p]] [-C[<dx/dy>]] [-D%s] [-E<dpi>] [-F]\n", GMT_inc_OPT);
		fprintf (stderr, "\t[-G<out_grdfile>] [-I] [-Mc|i|m] [-N<nx/ny>] [-S<radius>] [-V]\n\n");

		if (GMT_give_synopsis_and_exit) 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   Optionally append dx/dy to add (or subtract if -I) (i.e., false easting & northing) [0/0]\n");
		GMT_inc_syntax ('D', 0);
		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 (!Ctrl->G.file) {
		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 (Ctrl->S.radius < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -S:  Max search radius, if specified, must be positive\n", GMT_program);
		error++;
	}
	if ((Ctrl->M.active + Ctrl->A.active) == 2) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Can specify only one of -A and -M\n", GMT_program);
		error++;
	}
	if ((Ctrl->D.active + Ctrl->E.active + Ctrl->N.active) != 1) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify only one of -D, -E, or -N\n", GMT_program);
		error++;
	}
	if (Ctrl->D.active && (Ctrl->D.xinc <= 0.0 || Ctrl->D.yinc < 0.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -D option.  Must specify positive increment(s)\n", GMT_program);
		error++;
	}
	if (Ctrl->N.active && !set_n && (Ctrl->N.nx <= 0 || Ctrl->N.ny <= 0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -N option.  Must specify positive integers\n", GMT_program);
		error++;
	}
	if (Ctrl->E.active && Ctrl->E.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 .gmtcommands4 */

	if (Ctrl->M.active) GMT_set_measure_unit ('M', Ctrl->M.unit);
	shift_xy = !(Ctrl->C.easting == 0.0 && Ctrl->C.northing == 0.0);
	
	unit = GMT_check_scalingopt ('A', Ctrl->A.unit, scale_unit_name);
	GMT_init_scales (unit, &fwd_scale, &inv_scale, &inch_to_unit, &unit_to_inch, unit_name);

	GMT_map_setup (w, e, s, n);

	xmin = (Ctrl->C.active) ? project_info.xmin - project_info.x0 : project_info.xmin;
	xmax = (Ctrl->C.active) ? project_info.xmax - project_info.x0 : project_info.xmax;
	ymin = (Ctrl->C.active) ? project_info.ymin - project_info.y0 : project_info.ymin;
	ymax = (Ctrl->C.active) ? project_info.ymax - project_info.y0 : project_info.ymax;
	if (Ctrl->A.active) {	/* 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;
	}
	if (shift_xy) {
		xmin += Ctrl->C.easting;
		xmax += Ctrl->C.easting;
		ymin += Ctrl->C.northing;
		ymax += Ctrl->C.northing;
	}

	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 (Ctrl->I.active) {	/* 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);
		}

		if (Ctrl->S.careful) {
			GMT_boundcond_init (&edgeinfo);
			nm = ((size_t)(4 + r_head.nx)) * ((size_t)(4 + r_head.ny));
			GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;
		}
		else
			nm = ((size_t)r_head.nx) * ((size_t)r_head.ny);

		rect = (float *) GMT_memory (VNULL, nm, 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 (Ctrl->F.active) offset = !offset;	/* Toggle */
		if (set_n) {
			Ctrl->N.nx = r_head.nx;
			Ctrl->N.ny = r_head.ny;
		}
		GMT_grdproject_init (&g_head, Ctrl->D.xinc, Ctrl->D.yinc, Ctrl->N.nx, Ctrl->N.ny, Ctrl->E.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 (shift_xy) {
			r_head.x_min -= Ctrl->C.easting;
			r_head.x_max -= Ctrl->C.easting;
			r_head.y_min -= Ctrl->C.northing;
			r_head.y_max -= Ctrl->C.northing;

		}
		if (Ctrl->A.active) {	/* 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;
				Ctrl->S.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;
			Ctrl->S.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;
			Ctrl->S.radius *= unit_to_inch;
		}
		if (Ctrl->C.active) {	/* 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 = GMT_get_inc (r_head.x_min, r_head.x_max, r_head.nx, r_head.node_offset);
		r_head.y_inc = GMT_get_inc (r_head.y_min, r_head.y_max, r_head.ny, r_head.node_offset);

		/* rect xy values and Ctrl->S.radius are here in GMT projected inches */

		GMT_init_search_radius (&Ctrl->S.radius, &r_head, &g_head, TRUE);
		if (Ctrl->S.careful)
			GMT_grd_project (rect, &r_head, geo, &g_head, &edgeinfo, Ctrl->S.bilinear, TRUE);
		else
			GMT_grd_inverse (geo, &g_head, rect, &r_head, Ctrl->S.radius);

		if (Ctrl->S.careful) GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0;

		if (GMT_write_grd (Ctrl->G.file, &g_head, geo, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error writing file %s\n", GMT_program, Ctrl->G.file);
			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);
		}

		if (Ctrl->S.careful) {
			GMT_boundcond_init (&edgeinfo);
			nm = ((size_t)(4 + g_head.nx)) * ((size_t)(4 + g_head.ny));
			GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;
		}
		else
			nm = ((size_t)g_head.nx) * ((size_t)g_head.ny);

		geo = (float *) GMT_memory (VNULL, nm, 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 (Ctrl->A.active) {	/* Convert from 1:1 scale */
			if (unit) {	/* Undo the 1:1 unit used */
				Ctrl->D.xinc *= inv_scale;
				Ctrl->D.yinc *= inv_scale;
				Ctrl->S.radius *= inv_scale;
			}
			Ctrl->D.xinc *= project_info.x_scale;
			Ctrl->D.yinc *= project_info.y_scale;
			Ctrl->S.radius *= project_info.x_scale;
		}
		else if (gmtdefs.measure_unit != GMT_INCH) {	/* Convert from inch to whatever */
			Ctrl->D.xinc *= unit_to_inch;
			Ctrl->D.yinc *= unit_to_inch;
			Ctrl->S.radius *= unit_to_inch;
		}
		if (set_n) {
			Ctrl->N.nx = g_head.nx;
			Ctrl->N.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 (Ctrl->F.active) offset = !offset;	/* Toggle */

		GMT_grdproject_init (&r_head, Ctrl->D.xinc, Ctrl->D.yinc, Ctrl->N.nx, Ctrl->N.ny, Ctrl->E.dpi, offset);
		rect = (float *) GMT_memory (VNULL, (size_t)(r_head.nx * r_head.ny), sizeof (float), GMT_program);
		GMT_init_search_radius (&Ctrl->S.radius, &r_head, &g_head, FALSE);
		if (Ctrl->S.careful)
			GMT_grd_project (geo, &g_head, rect, &r_head, &edgeinfo, Ctrl->S.bilinear, FALSE);
		else
			GMT_grd_forward (geo, &g_head, rect, &r_head, Ctrl->S.radius);

		/* Modify output rect header if -A, -C, -M have been set */

		if (Ctrl->C.active) {	/* 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 (Ctrl->A.active) {	/* 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;
		}
		if (shift_xy) {
			r_head.x_min += Ctrl->C.easting;
			r_head.x_max += Ctrl->C.easting;
			r_head.y_min += Ctrl->C.northing;
			r_head.y_max += Ctrl->C.northing;

		}
		r_head.x_inc = GMT_get_inc (r_head.x_min, r_head.x_max, r_head.nx, r_head.node_offset);
		r_head.y_inc = GMT_get_inc (r_head.y_min, r_head.y_max, r_head.ny, r_head.node_offset);

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

		if (Ctrl->S.careful) GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0;

		if (GMT_write_grd (Ctrl->G.file, &r_head, rect, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
			fprintf (stderr, "%s: Error writing file %s\n", GMT_program, Ctrl->G.file);
			exit (EXIT_FAILURE);
		}

	}

	GMT_free ((void *)geo);
	GMT_free ((void *)rect);

	Free_Grdproject_Ctrl (Ctrl);	/* Deallocate control structure */

	GMT_end (argc, argv);

	exit (EXIT_SUCCESS);
}

void *New_Grdproject_Ctrl () {	/* Allocate and initialize a new control structure */
	struct GRDPROJECT_CTRL *C;
	
	C = (struct GRDPROJECT_CTRL *) GMT_memory (VNULL, 1, sizeof (struct GRDPROJECT_CTRL), "New_Grdproject_Ctrl");
	
	/* Initialize values whose defaults are not 0/FALSE/NULL */
		
	return ((void *)C);
}

void Free_Grdproject_Ctrl (struct GRDPROJECT_CTRL *C) {	/* Deallocate control structure */
	if (C->G.file) GMT_free ((void *)C->G.file);	
	GMT_free ((void *)C);	
}
