/*--------------------------------------------------------------------
 *	$Id: grdinfo.c,v 1.15 2004/08/26 02:54:15 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
 *--------------------------------------------------------------------*/
/*
 * grdinfo reads one or more grd file and [optionally] prints out various
 * statistics like mean/standard deviation and median/scale
 *
 * Author:	Paul Wessel
 * Date:	20-SEP-2001
  * Version:	4
 */

#include "gmt.h"

main (int argc, char **argv)
{
	int nfiles = 0, k, i, j, i_min, i_max, nm, n_nan = 0, n, n_grds = 0;
	
	BOOLEAN error = FALSE, l1 = FALSE, l2 = FALSE, quick = TRUE, find_max = FALSE;
	BOOLEAN world_form = FALSE, columns = FALSE, step = FALSE, t_range = FALSE;
	
	float *a;
	
	double x_min, y_min, z_min, x_max, y_max, z_max, global_xmin, global_xmax, global_ymin, global_ymax;
	double global_zmin, global_zmax, delta[3] = {0.0, 0.0, 0.0};
	double mean, median, sum2, stdev, scale, rms, x, half;
	
	char file[BUFSIZ], format[BUFSIZ], text[64];
	char *type[2] = { "Normal", "Pixel"};
	
	struct GRD_HEADER grd;
	
	argc = GMT_begin (argc, argv);

	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
				/* Common parameters */
			
				case 'V':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;
				
				/* Supplemental parameters */
			
				case 'C':
					columns = TRUE;
					break;
				case 'D':	/* Left for backwards compatibility, use -f instead */
					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;
					break;
				case 'F':
					world_form = TRUE;
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &delta[0], &delta[1]);
					step = TRUE;
					break;
				case 'M':
					quick = FALSE;
					find_max = TRUE;
					break;
				case 'L':
					quick = FALSE;
					if (argv[i][2] == 0 || argv[i][2] == '2')
						l2 = TRUE;
					else if (argv[i][2] == '1')
						l1 = TRUE;
					else {
						error = TRUE;
						fprintf (stderr, "%s: GMT SYNTAX ERROR -L option:  Choose between -L1 or -L2\n", GMT_program);
					}
					break;
				case 'T':
					t_range = TRUE;
					delta[2] = atof (&argv[i][2]);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			nfiles ++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdinfo %s - Extract information from netCDF grdfiles\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdinfo <grdfiles> [-C] [-F] [-I<dx>[/<dy>]] [-L1] [-L[2]] [-M]\n");
		fprintf (stderr, "	[-T<dz>] [-f[i|o]<colinfo>]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t<grdfiles> may be one or more netCDF grdfiles\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-C formats report in fields on a single line using the format\n");
		fprintf (stderr, "\t   file w e s n z0 z1 dx dy nx ny [x0 y0 x1 y1] [med scale] [mean std rms] [n_nan]\n");
		fprintf (stderr, "\t   (-M gives [x0 y0 x1 y1] and [n_nan]; -L1 gives [med scale]; -L2 gives [mean std rms])\n");
		fprintf (stderr, "\t-F reports domain in world mapping format [Default is generic]\n");
		fprintf (stderr, "\t-I returns textstring -Rw/e/s/n to nearest multiple of dx/dy.\n");
		fprintf (stderr, "\t   If -C is set then rounded off will occurr but no -R string is issued\n");
		fprintf (stderr, "\t-L1 reports median and L1-scale of data set\n");
		fprintf (stderr, "\t-L[2] reports mean, standard deviation, and rms of data set\n");
		fprintf (stderr, "\t-M searches for the global min and max locations (x0,y0) and (x1,y1)\n");
		fprintf (stderr, "\t-T given increment dz, return global -Tzmin/zmax/dz in multiples of dz\n");
		GMT_explain_option ('V');
		GMT_explain_option ('f');
		exit (EXIT_FAILURE);
	}
	
	if (nfiles == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one or more input files\n", GMT_program);
		error++;
	}
	if (t_range && delta[2] <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -T: Must specify a positive increment\n", GMT_program);
		error++;
	}
	if (step && (delta[0] <= 0.0 || delta[1] <= 0.0)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I: Must specify a positive increment(s)\n", GMT_program);
		error++;
	}
	if ((step || t_range) && !quick) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -M: Not compatible with -I or -T\n", GMT_program);
		error++;
	}
	if (t_range && step) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Only one of -I -T can be specified\n", GMT_program);
		error++;
	}

	if (error) exit (EXIT_FAILURE);

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

	a = (float *) GMT_memory (VNULL, (size_t)1, sizeof (float), GMT_program);

	global_xmin = global_ymin = global_zmin = +DBL_MAX;
	global_xmax = global_ymax = global_zmax = -DBL_MAX;

	for (k = 1; k < argc; k++) {	/* Loop over arguments, skip options */
	
		if (argv[k][0] == '-') continue;
		
		strcpy (file, argv[k]);
		for (j = 0; file[j]; j++) if (file[j] == '=') file[j] = 0;
		if (strcmp (file, "=") && access (file, R_OK)) {
			fprintf (stderr, "%s:  File %s not found\n", GMT_program, file);
			continue;
		}

		GMT_grd_init (&grd, argc, argv, FALSE);

		if (GMT_read_grd_info (argv[k], &grd)) {
			fprintf (stderr, "%s: Error opening file %s\n", GMT_program, file);
			continue;
		}
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Processing file %s\n", GMT_program, file);

		n_grds++;

		half = (grd.node_offset) ? 0.5 : 0.0;
		
		if (grd.z_min == grd.z_max && GMT_grd_i_format) quick = FALSE, find_max = TRUE;
		
		if (!quick) {	/* Must determine the location of global min and max values */

			nm = grd.nx * grd.ny;
			a = (float *) GMT_memory ((void *)a, (size_t)nm, sizeof (float), GMT_program);
			if (GMT_read_grd (argv[k], &grd, a, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) continue;
		
			z_min = DBL_MAX;	z_max = -DBL_MAX;
			mean = median = sum2 = 0.0;
			i_min = i_max = 0;
			n_nan = 0;
			for (i = 0; i < nm; i++) {
				if (GMT_is_fnan (a[i])) {
					n_nan++;
					continue;
				}
				if (find_max) {
					if (a[i] < z_min) {
						z_min = a[i];
						i_min = i;
					}
					if (a[i] > z_max) {
						z_max = a[i];
						i_max = i;
					}
				}
				if (l2) {
					mean += a[i];
					sum2 += a[i]*a[i];
				}
			}
		
			x_min = grd.x_min + ((i_min % grd.nx) + half) * grd.x_inc;
			y_min = grd.y_max - ((i_min / grd.nx) + half) * grd.y_inc;
			x_max = grd.x_min + ((i_max % grd.nx) + half) * grd.x_inc;
			y_max = grd.y_max - ((i_max / grd.nx) + half) * grd.y_inc;
		}

		if (l1) {	/* Calculate the median and L1 scale */
			qsort ((void *)a, (size_t)nm, sizeof (float), GMT_comp_float_asc);
			n = nm - n_nan;
			median = (n%2) ? a[n/2] : 0.5*(a[n/2-1] + a[n/2]);
			for (i = 0; i < n; i++) a[i] = (float)fabs (a[i] - median);
			qsort ((void *)a, (size_t)n, sizeof (float), GMT_comp_float_asc);
			scale = (n%2) ? 1.4826 * a[n/2] : 0.7413 * (a[n/2-1] + a[n/2]);
		}
		if (l2) {	/* Calculate the mean, standard deviation, and rms */
			x = (double)(nm - n_nan);
			stdev = (rint (x) > 1.0) ? sqrt((x*sum2 - mean*mean)/(x*(x-1))) : GMT_d_NaN;
			rms = (rint (x) > 0.0) ? sqrt (sum2 / x) : GMT_d_NaN;
			mean = (rint (x) > 0.0) ? mean / x : GMT_d_NaN;
		}

		/* OK, time to report results */

		if (columns && !step) {
			fprintf (GMT_stdout, "%s\t", file);
			GMT_ascii_output_one (GMT_stdout, grd.x_min, 0);	fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, grd.x_max, 0);	fprintf (GMT_stdout, "\t");
			GMT_ascii_output_one (GMT_stdout, grd.y_min, 1);	fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, grd.y_max, 1);	fprintf (GMT_stdout, "\t");
			GMT_ascii_output_one (GMT_stdout, grd.z_min, 2);	fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, grd.z_max, 2);	fprintf (GMT_stdout, "\t");
			GMT_ascii_format_one (text, grd.x_inc, GMT_io.out_col_type[0]);
			if (isalpha ((int)text[strlen(text)-1])) text[strlen(text)-1] = '\0';	/* Chop of trailing WESN flag here */
			fprintf (GMT_stdout, "%s\t", text);
			GMT_ascii_format_one (text, grd.y_inc, GMT_io.out_col_type[1]);
			if (isalpha ((int)text[strlen(text)-1])) text[strlen(text)-1] = '\0';	/* Chop of trailing WESN flag here */
			fprintf (GMT_stdout, "%s\t", text);
			GMT_ascii_output_one (GMT_stdout, (double)grd.nx, 2);	fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, (double)grd.ny, 2);

			if (find_max) {
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, x_min, 0);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, y_min, 1);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, x_max, 0);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, y_max, 1);
			}
			if (l1) {
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, median, 2);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, scale, 2);
			}
			if (l2) {
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, mean, 2);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, stdev, 2);
				fprintf (GMT_stdout, "\t");	GMT_ascii_output_one (GMT_stdout, rms, 2);
			}
			if (find_max) fprintf (GMT_stdout, "\t%d", n_nan);
			fprintf (GMT_stdout, "\n");
		}
		else if (!(t_range || step)) {
			fprintf (GMT_stdout, "%s: Title: %s\n", file, grd.title);
			fprintf (GMT_stdout, "%s: Command: %s\n", file, grd.command);
			fprintf (GMT_stdout, "%s: Remark: %s\n", file, grd.remark);
			fprintf (GMT_stdout, "%s: %s node registration used\n", file, type[grd.node_offset]);
			fprintf (GMT_stdout, "%s: grdfile format # %d\n", file, GMT_grd_i_format);
			if (world_form) {
				if ((fabs (grd.x_min) < 500.0) && (fabs (grd.x_max) < 500.0) && (fabs (grd.y_min) < 500.0) && (fabs (grd.y_max) < 500.0)) {
					fprintf (GMT_stdout, "%s: x_min: %.7f\n", file, grd.x_min);
					fprintf (GMT_stdout, "%s: x_max: %.7f\n", file, grd.x_max);
					fprintf (GMT_stdout, "%s: x_inc: %.7f\n", file, grd.x_inc);
					fprintf (GMT_stdout, "%s: units: %s\n", file, grd.x_units);
					fprintf (GMT_stdout, "%s: nx: %d\n", file, grd.nx);
					fprintf (GMT_stdout, "%s: y_min: %.7f\n", file, grd.y_min);
					fprintf (GMT_stdout, "%s: y_max: %.7f\n", file, grd.y_max);
					fprintf (GMT_stdout, "%s: y_inc: %.7f\n", file, grd.y_inc);
					fprintf (GMT_stdout, "%s: units: %s\n", file, grd.y_units);
					fprintf (GMT_stdout, "%s: ny: %d\n", file, grd.ny);
				}
				else {
					fprintf (GMT_stdout, "%s: x_min: %.2f\n", file, grd.x_min);
					fprintf (GMT_stdout, "%s: x_max: %.2f\n", file, grd.x_max);
					fprintf (GMT_stdout, "%s: x_inc: %.2f\n", file, grd.x_inc);
					fprintf (GMT_stdout, "%s: units: %s\n", file, grd.x_units);
					fprintf (GMT_stdout, "%s: nx: %d\n", file, grd.nx);
					fprintf (GMT_stdout, "%s: y_min: %.2f\n", file, grd.y_min);
					fprintf (GMT_stdout, "%s: y_max: %.2f\n", file, grd.y_max);
					fprintf (GMT_stdout, "%s: y_inc: %.2f\n", file, grd.y_inc);
					fprintf (GMT_stdout, "%s: units: %s\n", file, grd.y_units);
					fprintf (GMT_stdout, "%s: ny: %d\n", file, grd.ny);
				}
			}
			else {
				GMT_ascii_format_one (text, grd.x_inc, GMT_io.out_col_type[0]);
				if (isalpha ((int)text[strlen(text)-1])) text[strlen(text)-1] = '\0';	/* Chop of trailing WESN flag here */
				fprintf (GMT_stdout, "%s: x_min: ", file);
				GMT_ascii_output_one (GMT_stdout, grd.x_min, 0);
				fprintf (GMT_stdout, " x_max: ");
				GMT_ascii_output_one (GMT_stdout, grd.x_max, 0);
				fprintf (GMT_stdout, " x_inc: %s", text);
				fprintf (GMT_stdout, " units: %s nx: %d\n", grd.x_units, grd.nx);
				fprintf (GMT_stdout, "%s: y_min: ", file);
				GMT_ascii_output_one (GMT_stdout, grd.y_min, 1);
				fprintf (GMT_stdout, " y_max: ");
				GMT_ascii_output_one (GMT_stdout, grd.y_max, 1);
				GMT_ascii_format_one (text, grd.y_inc, GMT_io.out_col_type[1]);
				if (isalpha ((int)text[strlen(text)-1])) text[strlen(text)-1] = '\0';	/* Chop of trailing WESN flag here */
				fprintf (GMT_stdout, " y_inc: %s", text);
				fprintf (GMT_stdout, " units: %s ny: %d\n", grd.y_units, grd.ny);
			}
			
			if (find_max) {
				if (z_min == -DBL_MAX) z_min = GMT_d_NaN;
				if (z_max == +DBL_MAX) z_max = GMT_d_NaN;
				fprintf (GMT_stdout, "%s: z_min: ", file);
				GMT_ascii_output_one (GMT_stdout, z_min, 2);
				fprintf (GMT_stdout, " at x = ");
				GMT_ascii_output_one (GMT_stdout, x_min, 0);
				fprintf (GMT_stdout, " y = ");
				GMT_ascii_output_one (GMT_stdout, y_min, 1);
				fprintf (GMT_stdout, " z_max: ");
				GMT_ascii_output_one (GMT_stdout, z_max, 2);
				fprintf (GMT_stdout, " at x = ");
				GMT_ascii_output_one (GMT_stdout, x_max, 0);
				fprintf (GMT_stdout, " y = ");
				GMT_ascii_output_one (GMT_stdout, y_max, 1);
				fprintf (GMT_stdout, "\n");
			}
			else if (world_form) {
				fprintf (GMT_stdout, "%s: zmin: %g\n", file, grd.z_min);
				fprintf (GMT_stdout, "%s: zmax: %g\n", file, grd.z_max);
				fprintf (GMT_stdout, "%s: units: %s\n", file, grd.z_units);
			}
			else {
				fprintf (GMT_stdout, "%s: z_min: ", file);
				GMT_ascii_output_one (GMT_stdout, grd.z_min, 2);
				fprintf (GMT_stdout, " z_max: ");
				GMT_ascii_output_one (GMT_stdout, grd.z_max, 2);
				fprintf (GMT_stdout, " units: %s\n", grd.z_units);
			}

			GMT_ascii_format_one (text, grd.z_add_offset, GMT_io.out_col_type[2]);
			if (isalpha ((int)text[strlen(text)-1])) text[strlen(text)-1] = '\0';	/* Chop of trailing WESN flag here */
			sprintf (format, "%s: scale_factor: %s add_offset: %%s\n", file, gmtdefs.d_format);
			fprintf (GMT_stdout, format, grd.z_scale_factor, text);
			if (n_nan) fprintf (GMT_stdout, "%s: %d nodes set to NaN\n", file, n_nan);
			if (l1) {
				fprintf (GMT_stdout, "%s: median: ", file);	GMT_ascii_output_one (GMT_stdout, median, 2);
				fprintf (GMT_stdout, " scale: ");		GMT_ascii_output_one (GMT_stdout, scale, 2);
				fprintf (GMT_stdout, "\n");
			}
			if (l2) {
				fprintf (GMT_stdout, "%s: mean: ", file);	GMT_ascii_output_one (GMT_stdout, mean, 2);
				fprintf (GMT_stdout, " stdev: ");		GMT_ascii_output_one (GMT_stdout, stdev, 2);
				fprintf (GMT_stdout, " rms: ");			GMT_ascii_output_one (GMT_stdout, rms, 2);
				fprintf (GMT_stdout, "\n");
			}
		}
		else {
			if (grd.z_min < global_zmin) global_zmin = grd.z_min;
			if (grd.z_max > global_zmax) global_zmax = grd.z_max;
			if (grd.x_min < global_xmin) global_xmin = grd.x_min;
			if (grd.x_max > global_xmax) global_xmax = grd.x_max;
			if (grd.y_min < global_ymin) global_ymin = grd.y_min;
			if (grd.y_max > global_ymax) global_ymax = grd.y_max;
		}
	}
	
	if (global_zmin == -DBL_MAX) global_zmin = GMT_d_NaN;	/* Never got set */
	if (global_zmax == +DBL_MAX) global_zmax = GMT_d_NaN;
	
	if (columns && step) {
		global_xmin = floor (global_xmin / delta[0]) * delta[0];
		global_xmax = ceil  (global_xmax / delta[0]) * delta[0];
		global_ymin = floor (global_ymin / delta[1]) * delta[1];
		global_ymax = ceil  (global_ymax / delta[1]) * delta[1];
		fprintf (GMT_stdout, "%d\t", n_grds);
		GMT_ascii_output_one (GMT_stdout, global_xmin, 0);	fputc ('\t', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_xmax, 0);	fputc ('\t', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_ymin, 1);	fputc ('\t', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_ymax, 1);	fputc ('\t', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_zmin, 1);	fputc ('\t', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_zmax, 1);	fputc ('\n', GMT_stdout);
	}
	else if (t_range) {
		global_zmin = floor (global_zmin / delta[2]) * delta[2];
		global_zmax = ceil  (global_zmax / delta[2]) * delta[2];
		fprintf (GMT_stdout, "-T");
		GMT_ascii_output_one (GMT_stdout, global_zmin, 2);
		fprintf (GMT_stdout, "/");
		GMT_ascii_output_one (GMT_stdout, global_zmax, 2);
		fprintf (GMT_stdout, "/");
		GMT_ascii_output_one (GMT_stdout, delta[2], 2);
		fprintf (GMT_stdout, "\n");
	}
	else if (step) {
		global_xmin = floor (global_xmin / delta[0]) * delta[0];
		global_xmax = ceil  (global_xmax / delta[0]) * delta[0];
		global_ymin = floor (global_ymin / delta[1]) * delta[1];
		global_ymax = ceil  (global_ymax / delta[1]) * delta[1];
		fprintf (GMT_stdout, "-R");
		GMT_ascii_output_one (GMT_stdout, global_xmin, 0);	fputc ('/', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_xmax, 0);	fputc ('/', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_ymin, 1);	fputc ('/', GMT_stdout);
		GMT_ascii_output_one (GMT_stdout, global_ymax, 1);	fputc ('\n', GMT_stdout);
	}
	
	GMT_free ((void *)a);
	
	GMT_end (argc, argv);
}
