/*--------------------------------------------------------------------
 *	$Id: grdsample.c,v 1.4.4.2 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
 *--------------------------------------------------------------------*/
/*
 * grdsample reads a grdfile and evaluates the grid at new grid positions
 * specified by new dx/dy values using a 2-D Taylor expansion of order 3.
 * In order to evaluate derivatives along the edges of the surface, I assume 
 * natural bi-cubic spline conditions, i.e. both the second and third normal 
 * derivatives are zero, and that the dxdy derivative in the corners are zero, too.
 *
 * Author:	Paul Wessel
 * Date:	19-JUL-1989
 * Revised:	6-JAN-1990	PW: Updated to v.2.0
 * Revised:	16-JUN-1998	PW: Updated to v.3.1
 * Version:	3.4.1
 */

#include "gmt.h"

main (int argc, char **argv)
{
	int i, j, ij, one;
	
	BOOLEAN error = FALSE, greenwich = FALSE, offset = FALSE, bilinear = FALSE;
	BOOLEAN area_set = FALSE, n_set = FALSE, inc_set = FALSE, toggle = FALSE;
	
	double *lon, lat, dx2, dy2;
	
	float *a, *b;

	char *infile, *outfile, format[BUFSIZ];
	
	struct GRD_HEADER grd_a, grd_b;

	struct GMT_EDGEINFO edgeinfo;

	argc = GMT_begin (argc, argv);
	
	infile = outfile = CNULL;
	
	GMT_grd_init (&grd_b, argc, argv, FALSE);

	GMT_boundcond_init (&edgeinfo);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'R':
				case 'V':
				case '\0':
					error += GMT_get_common_args (argv[i], &grd_b.x_min, &grd_b.x_max, &grd_b.y_min, &grd_b.y_max);
					break;
				
				/* Supplemental parameters */
			
				case 'F':
					offset = TRUE;
					break;
				case 'G':
					outfile = &argv[i][2];
					break;
				case 'N':
					sscanf (&argv[i][2], "%d/%d", &grd_b.nx, &grd_b.ny);
					if (grd_b.ny == 0) grd_b.ny = grd_b.nx;
					n_set = TRUE;
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &grd_b.x_inc, &grd_b.y_inc);
					inc_set = TRUE;
					break;
				case 'L':
					error += GMT_boundcond_parse (&edgeinfo, &argv[i][2]);
					break;
				case 'Q':
					bilinear = TRUE;
					break;
				case 'T':	/* Convert from pixel file <-> gridfile */
					toggle = TRUE;
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else 
			infile = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdsample %s - Resample a gridded file onto a new grid\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdsample <old_grdfile> -G<new_grdfile> [-F] [-I<dx>[m|c][/<dy>[m|c]]]\n");
		fprintf (stderr, "\t[-L<flag>] [-N<nx/ny>] [-Q] [-R<west/east/south/north>] [-T] [-V]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<old_grdfile> is data set to be resampled\n");
		fprintf (stderr, "\t-G sets the name of the interpolated output grdfile\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-F force pixel node registration  [Default is centered]\n");
		fprintf (stderr, "\t-I sets the grid spacing (dx, dy) for the new grid\n");
		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-N specifies number of columns (nx) and rows (ny) of new grid\n");
		fprintf (stderr, "\t-Q do quick, bilinear interpolation [Default is bicubic]\n");
		fprintf (stderr, "\t-R specifies a subregion [Default is old region]\n");
		fprintf (stderr, "\t-T Toggles between grid registration and pixel registration\n");
		GMT_explain_option ('V');
		fprintf (stderr, "\t   One only of -N -I must be specified\n");
		
		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:  Must specify output file\n", GMT_program);
		error++;
	}
	if (!toggle) {
		if (grd_b.x_min != grd_b.x_max && grd_b.y_min != grd_b.y_max) area_set = TRUE;
		if (inc_set && n_set) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  Only one of -I, -N may be specified\n", GMT_program);
			error++;
		}
		if (n_set && (grd_b.nx <= 0 || grd_b.ny <= 0)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -N:  Must specify positive integers\n", GMT_program);
			error++;
		}
		if (inc_set && (grd_b.x_inc <= 0.0 || grd_b.y_inc <= 0.0)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR -I:  Must specify positive increments\n", GMT_program);
			error++;
		}
		if (!(inc_set || n_set)) {
			fprintf (stderr, "%s: GMT SYNTAX ERROR:  One of -I, -N must be specified\n", GMT_program);
			error++;
		}
	}
	
	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands */
	
	if (GMT_read_grd_info (infile, &grd_a)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, infile);
		exit (EXIT_FAILURE);
	}
	
	GMT_boundcond_param_prep (&grd_a, &edgeinfo);

	if (toggle) {
		offset = !grd_a.node_offset;	/* Change to the opposite of what it is */
		grd_b.nx = (offset) ? grd_a.nx - 1 : grd_a.nx + 1;
		grd_b.ny = (offset) ? grd_a.ny - 1 : grd_a.ny + 1;
		area_set = inc_set = FALSE;
	}
	
	a = (float *) GMT_memory (VNULL, (size_t)((grd_a.nx + 4) * (grd_a.ny + 4)), sizeof(float), GMT_program);
	
	if (area_set) {
		if (grd_b.y_min < grd_a.y_min || grd_b.y_max > grd_a.y_max) {
			fprintf (stderr, "%s:  Selected region exceeds the boundaries of the grdfile!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
		else if (!edgeinfo.nxp && (grd_b.x_min < grd_a.x_min || grd_b.x_max > grd_a.x_max)) {
			fprintf (stderr, "%s:  Selected region exceeds the boundaries of the grdfile!\n", GMT_program);
			exit (EXIT_FAILURE);
		}
	}
	
	if (!offset && !toggle) offset = grd_a.node_offset;
	one = !offset;
	grd_b.node_offset = offset;
	
	if (!area_set) {
		grd_b.x_min = grd_a.x_min;
		grd_b.x_max = grd_a.x_max;
		grd_b.y_min = grd_a.y_min;
		grd_b.y_max = grd_a.y_max;
	}
	
	if (edgeinfo.nxp && grd_b.x_min < 0.0 && grd_b.x_max > 0.0)
		greenwich = TRUE;
	else if (edgeinfo.nxp && grd_b.x_max > 360.0) {
		greenwich = TRUE;
		grd_b.x_min -= 360.0;
		grd_b.x_max -= 360.0;
	}
	if (inc_set) {
		grd_b.nx = irint ((grd_b.x_max - grd_b.x_min) / grd_b.x_inc) + one;
		grd_b.ny = irint ((grd_b.y_max - grd_b.y_min) / grd_b.y_inc) + one;
		grd_b.x_inc = (grd_b.x_max - grd_b.x_min) / (grd_b.nx - one);
		grd_b.y_inc = (grd_b.y_max - grd_b.y_min) / (grd_b.ny - one);
	}
	else {
		grd_b.x_inc = (grd_b.x_max - grd_b.x_min) / (grd_b.nx - one);
		grd_b.y_inc = (grd_b.y_max - grd_b.y_min) / (grd_b.ny - one);
	}
	
	GMT_grd_RI_verify (&grd_b, 1);

	b = (float *) GMT_memory (VNULL, (size_t)(grd_b.nx * grd_b.ny), sizeof(float), GMT_program);
	
	sprintf (format, "%%s: New grid (%s/%s/%s/%s) nx = %%d ny = %%d dx = %s dy = %s\n", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
	if (gmtdefs.verbose) fprintf (stderr, format, GMT_program, grd_b.x_min, grd_b.x_max, grd_b.y_min, grd_b.y_max, grd_b.nx, grd_b.ny, grd_b.x_inc, grd_b.y_inc);
	
	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 2;	/* Leave room for 2 empty boundary rows/cols */
	
	if (GMT_read_grd (infile, &grd_a, a, grd_a.x_min, grd_a.x_max, grd_a.y_min, grd_a.y_max, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, infile);
		exit (EXIT_FAILURE);
	}

	/* Initialize bcr structure:  */

	GMT_bcr_init (&grd_a, GMT_pad, bilinear);

	/* Set boundary conditions  */
	
	GMT_boundcond_set (&grd_a, &edgeinfo, GMT_pad, a);

	/* Precalculate longitudes */
	
	dx2 = 0.5 * grd_b.x_inc;
	dy2 = 0.5 * grd_b.y_inc;
	lon = (double *) GMT_memory (VNULL, (size_t)grd_b.nx, sizeof (double), GMT_program);
	for (i = 0; i < grd_b.nx; i++) {
		lon[i] = grd_b.x_min + (i * grd_b.x_inc) + ((offset) ? dx2 : 0.0);
		if (edgeinfo.nxp && greenwich && lon[i] > 180.0) lon[i] -= 360.0;
	}

	for (j = ij = 0; j < grd_b.ny; j++) {
		lat = grd_b.y_max - (j * grd_b.y_inc);
		if (offset) lat -= dy2;
		for (i = 0; i < grd_b.nx; i++, ij++) b[ij] = (float)GMT_get_bcr_z (&grd_a, lon[i], lat, a, &edgeinfo);
	}
	
	GMT_pad[0] = GMT_pad[1] = GMT_pad[2] = GMT_pad[3] = 0;	/* No boundary rows/cols on output */
	if (GMT_write_grd (outfile, &grd_b, b, 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 *)a);
	GMT_free ((void *)b);
	GMT_free ((void *)lon);
	
	GMT_end (argc, argv);
}