/*--------------------------------------------------------------------
 *	$Id: grdcut.c,v 1.3.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
 *--------------------------------------------------------------------*/
/*
 * grdcut.c reads a grdfile and writes a portion within it
 * to a new file.
 *
 * Author:	Walter Smith
 * Date:	5 august, 1988
 * Version:	3.4.1
 *
 */
 
#include "gmt.h"

float	*grd;

main (int argc, char **argv)
{
	BOOLEAN	error = FALSE, global = FALSE;

	int	i, nx_old, ny_old, nx_new, ny_new, one_or_zero;

	double	w_new = 0.0, e_new = 0.0, s_new = 0.0, n_new = 0.0;
	double	w_old, e_old, s_old, n_old;

	char *grd_in, *grd_out, format[BUFSIZ];

	struct GRD_HEADER header, test_header;

	argc = GMT_begin (argc, argv);
	
	grd_in = grd_out = CNULL;
	
	/* Check and interpret the command line arguments */
	
	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], &w_new, &e_new, &s_new, &n_new);
					break;
					
	 			case 'G':
	 				grd_out = &argv[i][2];
					break;
				default:		/* Options not recognized */
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
	 		grd_in = argv[i];
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr,"grdcut %s - Extract subsets from grdfiles\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdcut <input_grd> -G<output_grd> -R<west/east/south/north> [-V]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<input_grd> is file to extract a subset from.\n");
		fprintf (stderr, "\t-G specifies output grdfile\n");
		GMT_explain_option ('R');
		fprintf (stderr, "\t   Obviously, the WESN you specify must be within the WESN of the input file.\n");
		fprintf (stderr, "\t   If in doubt, run grdinfo first and check range of old file.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}

	/* Check that the options selected make sense */
	
	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (!grd_out) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G option:  Must specify output file\n", GMT_program);
		error++;
	}
	if (!grd_in) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

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


	if (GMT_read_grd_info (grd_in, &header)) {
		fprintf (stderr, "%s: Error opening file %s\n", GMT_program, grd_in);
		exit (EXIT_FAILURE);
	}
		
	if (s_new < header.y_min || s_new > header.y_max) error = TRUE;
	if (n_new < header.y_min || n_new > header.y_max) error = TRUE;
	global = (fabs (header.x_max - header.x_min) == 360.0);
		
	if ( !global && ((w_new < header.x_min) || (e_new > header.x_max)) ) error = TRUE;
	if (error) {
		fprintf (stderr, "%s: Subset exceeds data domain!\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	/* Make sure output grid is kosher */

	test_header.x_min = w_new;	test_header.x_max = e_new;	test_header.x_inc = header.x_inc;
	test_header.y_min = s_new;	test_header.y_max = n_new;	test_header.y_inc = header.y_inc;
	GMT_grd_RI_verify (&test_header, 1);

	GMT_grd_init (&header, argc, argv, TRUE);

	w_old = header.x_min;	e_old = header.x_max;
	s_old = header.y_min;	n_old = header.y_max;
	nx_old = header.nx;	ny_old = header.ny;
	one_or_zero = (header.node_offset) ? 0 : 1;
	nx_new = irint ((e_new - w_new) / header.x_inc) + one_or_zero;
	ny_new = irint ((n_new - s_new) / header.y_inc) + one_or_zero;
			
	grd = (float *) GMT_memory (VNULL, (size_t)(nx_new * ny_new), sizeof (float), GMT_program);
	if (GMT_read_grd (grd_in, &header, grd, w_new, e_new, s_new, n_new, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error reading file %s\n", GMT_program, grd_in);
		exit (EXIT_FAILURE);
	}

	if (gmtdefs.verbose) {
		sprintf (format, "\t%s\t%s\t%s\t%s\t%s\t%s\t%%d\t%%d\n", gmtdefs.d_format, gmtdefs.d_format,
			gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
		fprintf (stderr, "%s: File spec:\tW E S N dx dy nx ny:\n", GMT_program);
		fprintf (stderr, "%s: Old:", GMT_program);
		fprintf (stderr, format, w_old, e_old, s_old, n_old, header.x_inc, header.y_inc, nx_old, ny_old);
		fprintf (stderr, "%s: New:", GMT_program);
		fprintf (stderr, format, w_new, e_new, s_new, n_new, header.x_inc, header.y_inc, nx_new, ny_new);
	}

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