/*--------------------------------------------------------------------
 *	$Id: grdinfo.c,v 1.2.4.4 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
 *--------------------------------------------------------------------*/
/*
 * 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-OCT-1999
 * Version:	3.1	Based on 3.0
 * 		3.3.2	Added -D option
 * 		3.3.3	Added -F option via M. Odegard
 * 		3.3.5	12-JUL-2000
 * Version:	3.4.1
 */

#include "gmt.h"

void GMT_format_ddmmss (double val, char *format);

float *a;
char *type[2] = { "Normal", "Pixel"};

main (int argc, char **argv)
{
	int nfiles = 0, k, i, j, i_min, i_max, nm, n_nan = 0, n;
	
	BOOLEAN error = FALSE, l1 = FALSE, l2 = FALSE, quick = TRUE, find_max = FALSE;
	BOOLEAN world_form = FALSE, geographic = FALSE, columns = FALSE;
	
	double x_min, y_min, z_min, x_max, y_max, z_max;
	double mean, median, sum2, stdev, scale, rms, x, half;
	
	char file[BUFSIZ], format[BUFSIZ];
	
	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 '\0':
					error += GMT_get_common_args (argv[i], 0, 0, 0, 0);
					break;
				
				/* Supplemental parameters */
			
				case 'C':
					columns = TRUE;
					break;
				case 'D':
					geographic = TRUE;
					break;
				case 'F':
					world_form = 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;
				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] [-D] [-F] [-L1] [-L[2]] [-M]\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-D reports domain in dd:mm:ss format [Default is decimal]\n");
		fprintf (stderr, "\t-F reports domain in world mapping format [Default is generic]\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");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}
	
	if (nfiles == 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify one or more input files\n", GMT_program);
		error++;
	}
	if (geographic && world_form) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify only one of -D and -F\n", GMT_program);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

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

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

	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);

		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) {
			sprintf (format, "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%d\t%d", file, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, grd.nx, grd.ny);
			fprintf (GMT_stdout, format, grd.x_min, grd.x_max, grd.y_min, grd.y_max, grd.z_min, grd.z_max, grd.x_inc, grd.y_inc);
			if (find_max) {
				sprintf (format, "\t%s\t%s\t%s\t%s", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
				fprintf (GMT_stdout, format, x_min, y_min, x_max, y_max);
			}
			if (l1) {
				sprintf (format, "\t%s\t%s\n", gmtdefs.d_format, gmtdefs.d_format);
				fprintf (GMT_stdout, format, median, scale);
			}
			if (l2) {
				sprintf (format, "\t%s\t%s\t%s", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
				fprintf (GMT_stdout, format, mean, stdev, rms);
			}
			if (find_max) fprintf (GMT_stdout, "\t%d", n_nan);
			fprintf (GMT_stdout, "\n");
		}
		else {
			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 (geographic) {
				char item1[80], item2[80], item3[80];
				GMT_format_ddmmss (grd.x_min, item1);
				GMT_format_ddmmss (grd.x_max, item2);
				GMT_format_ddmmss (grd.x_inc, item3);
				fprintf (GMT_stdout, "%s: x_min: %s x_max: %s x_inc: %s units: %s nx: %d\n", file, item1, item2, item3, grd.x_units, grd.nx);
				GMT_format_ddmmss (grd.y_min, item1);
				GMT_format_ddmmss (grd.y_max, item2);
				GMT_format_ddmmss (grd.y_inc, item3);
				fprintf (GMT_stdout, "%s: y_min: %s y_max: %s y_inc: %s units: %s ny: %d\n", file, item1, item2, item3, grd.y_units, grd.ny);
			}
			else 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 {
				sprintf (format, "%s: x_min: %s x_max: %s x_inc: %s units: %s nx: %d\n", file, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, grd.x_units, grd.nx);
				fprintf (GMT_stdout, format, grd.x_min, grd.x_max, grd.x_inc);
				sprintf (format, "%s: y_min: %s y_max: %s y_inc: %s units: %s ny: %d\n", file, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, grd.y_units, grd.ny);
				fprintf (GMT_stdout, format, grd.y_min, grd.y_max, grd.y_inc);
			}
			
			if (find_max) {
				if (z_min == -DBL_MAX) z_min = GMT_d_NaN;
				if (z_max == +DBL_MAX) z_max = GMT_d_NaN;
				sprintf (format, "%s: z_min: %s at x = %s y = %s z_max: %s at x = %s y = %s units: %s\n", file, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format, grd.z_units);
				fprintf (GMT_stdout, format, z_min, x_min, y_min, z_max, x_max, y_max);
			}
			else if (world_form) {
				fprintf (GMT_stdout, "%s: zmin: %lg\n", file, grd.z_min);
				fprintf (GMT_stdout, "%s: zmax: %lg\n", file, grd.z_max);
				fprintf (GMT_stdout, "%s: units: %s\n", file, grd.z_units);
			}
			else {
				sprintf (format, "%s: z_min: %s z_max: %s units: %s\n", file, gmtdefs.d_format, gmtdefs.d_format, grd.z_units);
				fprintf (GMT_stdout, format, grd.z_min, grd.z_max);
			}

			sprintf (format, "%s: scale_factor: %s add_offset: %s\n", file, gmtdefs.d_format, gmtdefs.d_format);
			fprintf (GMT_stdout, format, grd.z_scale_factor, grd.z_add_offset);
			if (n_nan) fprintf (GMT_stdout, "%s: %d nodes set to NaN\n", file, n_nan);
			if (l1) {
				sprintf (format, "%s: median: %s scale: %s\n", file, gmtdefs.d_format, gmtdefs.d_format);
				fprintf (GMT_stdout, format, median, scale);
			}
			if (l2) {
				sprintf (format, "%s: mean: %s stdev: %s rms: %s\n", file, gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
				fprintf (GMT_stdout, format, mean, stdev, rms);
			}
		}
	}
	
	GMT_free ((void *)a);
	
	GMT_end (argc, argv);
}

void GMT_format_ddmmss (double val, char *format)
{
	int sign, minutes, seconds, ival;
	double frac_sec, f_seconds;

	/* Return a formatted dd:mm:ss or dd:mm:ss.fff string */
	
	sign = (val < 0.0) ? -1 : 1;
	val = fabs (val);
	ival = (int)val;	/* Truncate to integer in the direction toward 0 */
	minutes = seconds = 0;
	if ((val - (double) ival) > SMALL) {
		minutes = (int)floor (((val - ival) * 60.0) + SMALL);
		if (minutes == 60) {
			minutes = 0;
			ival = irint (val);
		}
		f_seconds = (val - (double)ival - (double)minutes / 60.0) * 3600.0;
		seconds = (int)floor (f_seconds + SMALL);
		frac_sec = f_seconds - (double)seconds;
		if (seconds == 60 && frac_sec < SMALL) {
			seconds = 0;
			minutes++;
			if (minutes == 60) {
				minutes = 0;
				ival = irint (val);
			}
		}
	}

	if (frac_sec > SMALL)	/* Must preserve second fraction */
		sprintf (format, "%d:%2.2d:%2.2d.%3.3d", sign * ival, minutes, seconds, irint (1000.0 * frac_sec));
	else
		sprintf (format, "%d:%2.2d:%2.2d", sign * ival, minutes, seconds);
}