/*--------------------------------------------------------------------
 *
 *	Copyright (c) 1991-2024 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
 *	See LICENSE.TXT 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 Lesser General Public License as published by
 *	the Free Software Foundation; version 3 or any later version.
 *
 *	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 Lesser General Public License for more details.
 *
 *	Contact info: www.generic-mapping-tools.org
 *--------------------------------------------------------------------*/
/*
 * Brief synopsis: xyz2grd reads a xyz file from standard input and creates the
 * corresponding grd-file. The input file does not have to be
 * sorted. xyz2grd will report if some nodes are missing.
 *
 * Author:	Paul Wessel
 * Date:	1-JAN-2010
 * Version:	6 API
 */

#include "gmt_dev.h"
#include "longopt/xyz2grd_inc.h"

#define THIS_MODULE_CLASSIC_NAME	"xyz2grd"
#define THIS_MODULE_MODERN_NAME	"xyz2grd"
#define THIS_MODULE_LIB		"core"
#define THIS_MODULE_PURPOSE	"Convert data table to a grid"
#define THIS_MODULE_KEYS	"<D{,SD),GG}"
#define THIS_MODULE_NEEDS	"R"
#define THIS_MODULE_OPTIONS "-:JRVbdefhiqrw" GMT_OPT("FH")

struct XYZ2GRD_CTRL {
	struct XYZ2GRD_In {
		char *file;	/* Deprecated use with -E only */
	} In;
	struct XYZ2GRD_A {	/* -A[f|l|n|m|r|s|u|z] */
		bool active;
		char mode;
	} A;
	struct XYZ2GRD_D {	/* -D[+x<xname>][+yyname>][+z<zname>][+s<scale>][+ooffset>][+n<invalid>][+t<title>][+r<remark>] */
		bool active;
		char *information;
	} D;
	struct XYZ2GRD_E {	/* -E[<nodata>] */
		bool active;
		bool set;
		double nodata;
	} E;
	struct XYZ2GRD_G {	/* -G<output_grdfile> */
		bool active;
		char *file;
	} G;
	struct XYZ2GRD_I {	/* -I (for checking only) */
		bool active;
	} I;
	struct XYZ2GRD_S {	/* -S */
		bool active;
		char *file;
	} S;
	struct GMT_PARSE_Z_IO Z;
};

static void *New_Ctrl (struct GMT_CTRL *GMT) {	/* Allocate and initialize a new control structure */
	struct XYZ2GRD_CTRL *C = NULL;

	C = gmt_M_memory (GMT, NULL, 1, struct XYZ2GRD_CTRL);

	/* Initialize values whose defaults are not 0/false/NULL */
	C->Z.type = 'a';
	C->Z.format[0] = 'T';	C->Z.format[1] = 'L';
	return (C);
}

static void Free_Ctrl (struct GMT_CTRL *GMT, struct XYZ2GRD_CTRL *C) {	/* Deallocate control structure */
	if (!C) return;
	gmt_M_str_free (C->In.file);
	gmt_M_str_free (C->D.information);
	gmt_M_str_free (C->G.file);
	gmt_M_str_free (C->S.file);
	gmt_M_free (GMT, C);
}

static int usage (struct GMTAPI_CTRL *API, int level) {
	const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
	if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
	GMT_Usage (API, 0, "usage: %s [<table>] -G%s %s %s [-A[d|f|l|m|n|r|S|s|u|z]] [%s] [%s] [-S[<zfile>]] [%s] "
		"[-Z[<flags>]] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s] [%s]\n", name, GMT_OUTGRID, GMT_I_OPT, GMT_Rgeo_OPT,
		GMT_GRDEDIT2D, GMT_J_OPT, GMT_V_OPT, GMT_bi_OPT, GMT_di_OPT, GMT_e_OPT, GMT_f_OPT, GMT_h_OPT, GMT_i_OPT,
		GMT_qi_OPT, GMT_r_OPT, GMT_w_OPT, GMT_colon_OPT, GMT_PAR_OPT);

	if (level == GMT_SYNOPSIS) return (GMT_MODULE_SYNOPSIS);

	GMT_Message (API, GMT_TIME_NONE, "\n  REQUIRED ARGUMENTS:\n");
	GMT_Option (API, "<");
	gmt_outgrid_syntax (API, 'G', "Set name of the output grid file");
	GMT_Option (API, "IR");
	GMT_Message (API, GMT_TIME_NONE, "\n  OPTIONAL ARGUMENTS:\n");
	GMT_Usage (API, 1, "\n-A[d|f|l|m|n|r|S|s|u|z]");
	GMT_Usage (API, -2, "Determine what to do if multiple entries are found for a node [Default will compute mean values]:");
	GMT_Usage (API, 3, "d: Compute the range (between min and max) of multiple entries per node.");
	GMT_Usage (API, 3, "f: Keep first value if multiple entries per node.");
	GMT_Usage (API, 3, "l: Keep lower (minimum) value if multiple entries per node.");
	GMT_Usage (API, 3, "m: Compute mean of multiple entries per node.");
	GMT_Usage (API, 3, "n: Count number of multiple entries per node.");
	GMT_Usage (API, 3, "r: Compute RMS of multiple entries per node.");
	GMT_Usage (API, 3, "S: Compute standard deviation of multiple entries per node.");
	GMT_Usage (API, 3, "s: Keep last value if multiple entries per node.");
	GMT_Usage (API, 3, "u: Keep upper (maximum) value if multiple entries per node.");
	GMT_Usage (API, 3, "z: Sum multiple entries at the same node.");
	gmt_grd_info_syntax (API->GMT, 'D');
	GMT_Option (API, "J");
	GMT_Usage (API, 1, "\n-S[<zfile>]");
	GMT_Usage (API, -2, "Swap the byte-order of the input data and write result to <zfile> (or standard output if no "
		"file given).  Requires -Z, and no grid file created! For this option, only one input file (or standard input) "
		"is allowed. Note: Cannot handle swapping of 64-bit integers.");
	GMT_Option (API, "V");
	GMT_Usage (API, 1, "\n-Z[<flags>]");
	GMT_Usage (API, -2, "Set exact specification of incoming 1-column z-table. Unless -S is used, append two chars that "
		"indicate file organization:");
	GMT_Usage (API, 3, "%s If data is in row format, state if first row is at T(op) or B(ottom). Then, append L or R to "
		"indicate starting point in row.", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "%s If data is in column format, state if first columns is L(left) or R(ight). Then, append T or "
		"B to indicate starting point in column.", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "%s To skip a header of size <n> bytes, append s<n> [<n> = 0].", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "%s To swap the byte-order of each word, append w.", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "%s Append x if gridline-registered, periodic data in x without repeating column at xmax.", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "%s Append y if gridline-registered, periodic data in y without repeating row at ymax.", GMT_LINE_BULLET);
	GMT_Usage (API, -2, "Specify one of the following data types (all binary except a):");
	GMT_Usage (API, 3, "A:  ASCII (multiple floating point values per record).", GMT_LINE_BULLET);
	GMT_Usage (API, 3, "a:  ASCII (one value per record).");
	GMT_Usage (API, 3, "c:  int8_t, signed 1-byte character.");
	GMT_Usage (API, 3, "u:  uint8_t, unsigned 1-byte character.");
	GMT_Usage (API, 3, "h:  int16_t, signed 2-byte integer.");
	GMT_Usage (API, 3, "H:  uint16_t, unsigned 2-byte integer.");
	GMT_Usage (API, 3, "i:  int32_t, signed 4-byte integer.");
	GMT_Usage (API, 3, "I:  uint32_t, unsigned 4-byte integer.");
	GMT_Usage (API, 3, "l:  int64_t, signed long (8-byte) integer.");
	GMT_Usage (API, 3, "L:  uint64_t, unsigned long (8-byte) integer.");
	GMT_Usage (API, 3, "f:  4-byte floating point single precision.");
	GMT_Usage (API, 3, "d:  8-byte floating point double precision.");
	GMT_Usage (API, -2, "[Default format is scanline orientation in ASCII representation: -ZTLa]. This option assumes "
		"all nodes have data values.");
	GMT_Option (API, "bi3,di");
	if (gmt_M_showusage (API)) GMT_Usage (API, -2, "Also sets value for nodes without input xyz triplet [Default is NaN].");
	GMT_Option (API, "e,f,h,i,qi,r,w,:,.");

	return (GMT_MODULE_USAGE);
}

static int parse (struct GMT_CTRL *GMT, struct XYZ2GRD_CTRL *Ctrl, struct GMT_Z_IO *io, struct GMT_OPTION *options) {
	/* This parses the options provided to xyz2grd and sets parameters in CTRL.
	 * Any GMT common options will override values set previously by other commands.
	 * It also replaces any file names specified as input or output with the data ID
	 * returned when registering these sources/destinations with the API.
	 */

	unsigned int n_errors = 0, n_files = 0;
	uint64_t n_req;
	bool do_grid, b_only = false, got_E = false;
	char *ptr_to_arg = NULL;
	struct GMT_OPTION *opt = NULL;
	struct GMTAPI_CTRL *API = GMT->parent;

	memset (io, 0, sizeof(struct GMT_Z_IO)); /* Initialize with zero */

	for (opt = options; opt; opt = opt->next) {	/* Look for deprecated -E first */
		if (opt->option == 'E') got_E = true;
	}

	for (opt = options; opt; opt = opt->next) {
		switch (opt->option) {

			case '<':	/* Input files */
				if (n_files++ > 0) break;
				if (got_E)
					n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_DATASET, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->In.file));
				break;

			/* Processes program-specific parameters */

			case 'A':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->A.active);
				if (!opt->arg[0] && gmt_M_compat_check (GMT, 4)) {	/* In GMT4, just -A implied -Az */
					GMT_Report (API, GMT_MSG_COMPAT, "Option -A is deprecated; use -Az instead.\n");
					Ctrl->A.mode = 'z';
				}
				else if (!strchr ("dflmnrSsuz", opt->arg[0])) {
					GMT_Report (API, GMT_MSG_ERROR, "Select -Ad, -Af, -Al, -Am, -An, -Ar, -AS, -As, -Au, or -Az\n");
					n_errors++;
				}
				else
					Ctrl->A.mode = opt->arg[0];
				break;
			case 'D':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->D.active);
				n_errors += gmt_get_required_string (GMT, opt->arg, opt->option, 0, &Ctrl->D.information);
				break;
			case 'E':
				if (gmt_M_compat_check (GMT, 4)) {
					GMT_Report (API, GMT_MSG_COMPAT, "Option -E is deprecated; use grdconvert instead.\n");
					n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active);
					if (opt->arg[0]) {
						Ctrl->E.nodata = atof (opt->arg);
						Ctrl->E.set = true;
					}
				}
				else
					n_errors += gmt_default_option_error (GMT, opt);
				break;
			case 'G':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->G.active);
				n_errors += gmt_get_required_file (GMT, opt->arg, opt->option, 0, GMT_IS_GRID, GMT_OUT, GMT_FILE_LOCAL, &(Ctrl->G.file));
				break;
			case 'I':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->I.active);
				n_errors += gmt_parse_inc_option (GMT, 'I', opt->arg);
				break;
			case 'N':
				if (gmt_M_compat_check (GMT, 4)) {	/* Honor old -N<value> option */
					GMT_Report (API, GMT_MSG_COMPAT, "Option -N is deprecated; use GMT common option -di<nodata> instead.\n");
					if (opt->arg[0]) {
						char arg[GMT_LEN64] = {""};
						sprintf (arg, "i%s", opt->arg);
						n_errors += gmt_parse_d_option (GMT, arg);
					}
					else {
						GMT_Report (API, GMT_MSG_ERROR, "Option -N: Must specify value or NaN\n");
						n_errors++;
					}
				}
				else
					n_errors += gmt_default_option_error (GMT, opt);
				break;
			case 'S':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->S.active);
				if (opt->arg[0]) Ctrl->S.file = strdup (opt->arg);
				break;
			case 'Z':
				n_errors += gmt_M_repeated_module_option (API, Ctrl->Z.active);
				ptr_to_arg = opt->arg;
				break;

			default:	/* Report bad options */
				n_errors += gmt_default_option_error (GMT, opt);
				if (opt->option == 'b') b_only = true;
				break;
		}
	}

	if (Ctrl->Z.active) {
		if (Ctrl->S.active) Ctrl->Z.not_grid = true;	/* The row/col organization does not apply */
		n_errors += gmt_parse_z_io (GMT, ptr_to_arg, &Ctrl->Z);
	}

	n_errors += gmt_M_check_condition (GMT, Ctrl->S.active && !Ctrl->Z.active, "Option -S: Must also specify -Z\n");
	if (Ctrl->S.active) {	/* Reading and writing binary file */
		if (n_files > 1) {
			GMT_Report (API, GMT_MSG_ERROR, "-S can only handle one input file\n");
			n_errors++;
		}
		else {
			strcpy (GMT->current.io.r_mode, "rb");
			strcpy (GMT->current.io.w_mode, "wb");
		}
		Ctrl->Z.swab = 1;	/* Only swap on input */
	}

	if (Ctrl->Z.active) {
		n_errors += gmt_init_z_io (GMT, Ctrl->Z.format, Ctrl->Z.repeat, Ctrl->Z.swab, Ctrl->Z.skip, Ctrl->Z.type, io);
		GMT->common.b.type[GMT_IN] = Ctrl->Z.type;
		if (b_only) {
			GMT->common.b.active[GMT_IN] = false;
			GMT_Report (API, GMT_MSG_ERROR, "-Z overrides -bi\n");
		}
	}

	do_grid = !(Ctrl->S.active || Ctrl->E.active);
	if (do_grid) {
		n_errors += gmt_M_check_condition (GMT, !GMT->common.R.active[RSET], "Must specify -R option\n");
		n_errors += gmt_M_check_condition (GMT, GMT->common.R.inc[GMT_X] <= 0.0 || GMT->common.R.inc[GMT_Y] <= 0.0, "Option -I: Must specify positive increment(s)\n");
	}
	n_errors += gmt_M_check_condition (GMT, !(Ctrl->S.active || Ctrl->G.active), "Must specify either -G or -S\n");
	n_req = (Ctrl->Z.active) ? 1 : ((Ctrl->A.mode == 'n') ? 2 : 3);
	n_errors += gmt_check_binary_io (GMT, n_req);

	return (n_errors ? GMT_PARSE_ERROR : GMT_NOERROR);
}

#define bailout(code) {gmt_M_free_options (mode); return (code);}
#define Return(code) {Free_Ctrl (GMT, Ctrl); gmt_end_module (GMT, GMT_cpy); bailout (code);}

GMT_LOCAL void xyz2grd_protect_J (struct GMTAPI_CTRL *API, struct GMT_OPTION *options) {
	if (GMT_Find_Option (API, 'J', options) != NULL) {
		struct GMT_OPTION *opt = GMT_Make_Option (API, 'f', "0f,1f");
		(void)GMT_Append_Option(API, opt, options);
	}
}

EXTERN_MSC int GMT_xyz2grd (void *V_API, int mode, void *args) {
	bool previous_bin_i = false, previous_bin_o = false;
	int error = 0, scol, srow;
	unsigned int zcol, row, col, i, *flag = NULL, n_min = 1, was, n_ave = 0;
	uint64_t n_empty = 0, n_confused = 0;
	uint64_t ij, ij_gmt, n_read, n_filled = 0, n_used = 0, n_req;

	char c, Amode;

	double *in = NULL, wesn[4];

	gmt_grdfloat no_data_f, *data = NULL;

	void * (*save_i) (struct GMT_CTRL *, FILE *, uint64_t *, int *) = NULL;
	int (*save_o) (struct GMT_CTRL *, FILE *, uint64_t, double *, char *);

	struct GMT_GRID *Grid = NULL;
	struct GMT_Z_IO io;
	struct GMT_RECORD *In = NULL;
	struct XYZ2GRD_CTRL *Ctrl = NULL;
	struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL;
	struct GMT_OPTION *options = NULL;
	struct GMTAPI_CTRL *API = gmt_get_api_ptr (V_API);	/* Cast from void to GMTAPI_CTRL pointer */

	/*----------------------- Standard module initialization and parsing ----------------------*/

	if (API == NULL) return (GMT_NOT_A_SESSION);
	if (mode == GMT_MODULE_PURPOSE) return (usage (API, GMT_MODULE_PURPOSE));	/* Return the purpose of program */
	options = GMT_Create_Options (API, mode, args);	if (API->error) return (API->error);	/* Set or get option list */

	if ((error = gmt_report_usage (API, options, 0, usage)) != GMT_NOERROR) bailout (error);	/* Give usage if requested */

	/* Parse the command-line arguments */

	xyz2grd_protect_J(API, options);	/* If -J is used, add a -f0f,1f option to avoid later parsing errors due to -R & -J conflicts */
	if ((GMT = gmt_init_module (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_KEYS, THIS_MODULE_NEEDS, module_kw, &options, &GMT_cpy)) == NULL) bailout (API->error); /* Save current state */
	if (GMT_Parse_Common (API, THIS_MODULE_OPTIONS, options)) Return (API->error);
	Ctrl = New_Ctrl (GMT);	/* Allocate and initialize a new control structure */
	if ((error = parse (GMT, Ctrl, &io, options)) != 0) Return (error);

	/*---------------------------- This is the xyz2grd main code ----------------------------*/

	n_req = (Ctrl->Z.active) ? 1 : ((Ctrl->A.mode == 'n') ? 2 : 3);	/* Required input columns */

	if (Ctrl->S.active) {	/* Just swap data and bail */
		int out_ID;
		unsigned w_mode = GMT_ADD_DEFAULT;

		save_i = GMT->current.io.input;			/* Save previous i/0 parameters */
		save_o = GMT->current.io.output;
		previous_bin_i = GMT->common.b.active[GMT_IN];
		previous_bin_o = GMT->common.b.active[GMT_OUT];
		GMT->current.io.input = gmt_z_input;		/* Override input reader with chosen binary reader for selected type */
		GMT->current.io.output = gmt_z_output;		/* Override output writer with chosen binary writer for selected type */
		GMT->common.b.active[GMT_IN] = io.binary;	/* May have to set input binary as well */
		GMT->common.b.active[GMT_OUT] = io.binary;	/* May have to set output binary as well */
		if ((error = GMT_Set_Columns (API, GMT_IN, 1, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) Return (error);
		/* Initialize the i/o since we are doing record-by-record reading/writing */
		GMT_Report (API, GMT_MSG_INFORMATION, "Swapping data bytes only\n");
		if (Ctrl->S.active) io.swab = true;	/* Need to pass swabbing down to the gut level */

		/* Register the data source */
		if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_IN,  GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {	/* Registers default input sources, unless already set via file */
			Return (API->error);
		}
		if (Ctrl->S.file) {	/* Specified an output file */
			if ((out_ID = GMT_Register_IO (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_NONE, GMT_OUT, NULL, Ctrl->S.file)) == GMT_NOTSET) {
				Return (API->error);
			}
			w_mode = GMT_ADD_EXISTING;
		}
		if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_NONE, GMT_OUT, w_mode, 0, options) != GMT_NOERROR) {	/* Establishes data output to stdout */
			Return (API->error);
		}
		if ((error = GMT_Set_Columns (API, GMT_IN, 1, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {	/* We don't really care or know about columns so must use 1 */
			Return (API->error);
		}
		if ((error = GMT_Set_Columns (API, GMT_OUT, 1, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {	/* We don't really care or know about columns so must use 1 */
			Return (API->error);
		}
		/* Initialize the i/o for doing record-by-record reading/writing */
		if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_OFF) != GMT_NOERROR) {	/* Enables data input and sets access mode */
			Return (API->error);
		}
		if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_OFF) != GMT_NOERROR) {	/* Enables data output and sets access mode */
			Return (API->error);
		}
		if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_POINT) != GMT_NOERROR) {	/* Sets output geometry */
			Return (API->error);
		}
		do {	/* Keep returning records until we reach EOF */
			if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) {	/* Read next record, get NULL if special case */
				if (gmt_M_rec_is_error (GMT)) 		/* Bail if there are any read errors */
					Return (GMT_RUNTIME_ERROR);
				if (gmt_M_rec_is_any_header (GMT)) 	/* Skip all headers */
					continue;
				if (gmt_M_rec_is_eof (GMT)) 		/* Reached end of file */
					break;
				assert (false);						/* Should never get here */
			}

			if (In->data == NULL) {
				gmt_quit_bad_record (API, In);
				Return (API->error);
			}

			/* Data record to process */

			GMT_Put_Record (API, GMT_WRITE_DATA, In);
		} while (true);

		GMT->current.io.input = save_i;			/* Reset input pointer */
		GMT->common.b.active[GMT_IN] = previous_bin_i;	/* Reset input binary */
		GMT->current.io.output = save_o;		/* Reset output pointer */
		GMT->common.b.active[GMT_OUT] = previous_bin_o;	/* Reset output binary */

		if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) {	/* Disables further data input */
			Return (API->error);
		}
		if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) {	/* Disables further data output */
			Return (API->error);
		}
		Return (GMT_NOERROR);	/* We are done here */
	}

	GMT_Report (API, GMT_MSG_INFORMATION, "Processing input table data\n");

	/* PW: This is now done in grdconvert since ESRI Arc Interchange is a recognized format */
	if (Ctrl->E.active && gmt_M_compat_check (GMT, 4)) {	/* Read an ESRI Arc Interchange grid format in ASCII.  This must be a single physical file. */
		uint64_t n_left;
		double value;
		char line[GMT_BUFSIZ];
		FILE *fp = GMT->session.std[GMT_IN];

		if (Ctrl->In.file && (fp = gmt_fopen (GMT, Ctrl->In.file, "r")) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Cannot open file %s\n", Ctrl->In.file);
			Return (GMT_ERROR_ON_FOPEN);
		}

		if ((Grid = gmt_create_grid (GMT)) == NULL) Return (API->error);
		gmt_grd_init (GMT, Grid->header, options, false);
		Grid->header->registration = GMT_GRID_NODE_REG;
		gmt_fgets (GMT, line, GMT_BUFSIZ, fp);
		if (sscanf (line, "%*s %d", &Grid->header->n_columns) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not decode ncols record\n");
			Return (GMT_DATA_READ_ERROR);
		}
		gmt_fgets (GMT, line, GMT_BUFSIZ, fp);
		if (sscanf (line, "%*s %d", &Grid->header->n_rows) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not decode ncols record\n");
			Return (GMT_DATA_READ_ERROR);
		}
		gmt_fgets (GMT, line, GMT_BUFSIZ, fp);
		if (sscanf (line, "%*s %lf", &Grid->header->wesn[XLO]) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not decode xll record\n");
			Return (GMT_DATA_READ_ERROR);
		}
		if (!strncmp (line, "xllcorner", 9U)) Grid->header->registration = GMT_GRID_PIXEL_REG;	/* Pixel grid */
		gmt_fgets (GMT, line, GMT_BUFSIZ, fp);
		if (sscanf (line, "%*s %lf", &Grid->header->wesn[YLO]) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not decode yll record\n");
			Return (GMT_DATA_READ_ERROR);
		}
		if (!strncmp (line, "yllcorner", 9U)) Grid->header->registration = GMT_GRID_PIXEL_REG;	/* Pixel grid */
		gmt_fgets (GMT, line, GMT_BUFSIZ, fp);
		if (sscanf (line, "%*s %lf", &Grid->header->inc[GMT_X]) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Could not decode cellsize record\n");
			Return (GMT_DATA_READ_ERROR);
		}
		Grid->header->inc[GMT_Y] = Grid->header->inc[GMT_X];
		Grid->header->xy_off = 0.5 * Grid->header->registration;
		Grid->header->wesn[XHI] = Grid->header->wesn[XLO] + (Grid->header->n_columns - 1 + Grid->header->registration) * Grid->header->inc[GMT_X];
		Grid->header->wesn[YHI] = Grid->header->wesn[YLO] + (Grid->header->n_rows - 1 + Grid->header->registration) * Grid->header->inc[GMT_Y];
		gmt_set_grddim (GMT, Grid->header);
		if (gmt_M_err_fail (GMT, gmt_grd_RI_verify (GMT, Grid->header, 1), Ctrl->G.file)) Return (GMT_RUNTIME_ERROR);

		GMT_Report (API, GMT_MSG_INFORMATION, "n_columns = %d  n_rows = %d\n", Grid->header->n_columns, Grid->header->n_rows);
		n_left = Grid->header->nm;

		if ((Grid->data = gmt_M_memory_aligned (GMT, NULL, Grid->header->nm, gmt_grdfloat)) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Unable to allocate grid memory\n");
			Return (GMT_MEMORY_ERROR);
		}
		/* ESRI grids are scanline oriented (top to bottom), as are the GMT grids */
		row = col = 0;
		if (fscanf (fp, "%s", line) != 1) {
			GMT_Report (API, GMT_MSG_ERROR, "Unable to read nodata-flag or first data record from ESRI file\n");
			Return (GMT_DATA_READ_ERROR);
		}
		gmt_str_tolower (line);
		if (!strcmp (line, "nodata_value")) {	/* Found the optional nodata word */
			if (fscanf (fp, "%lf", &value) != 1) {
				GMT_Report (API, GMT_MSG_ERROR, "Unable to parse nodata-flag from ESRI file\n");
				Return (GMT_DATA_READ_ERROR);
			}
			if (Ctrl->E.set && !doubleAlmostEqualZero (value, Ctrl->E.nodata)) {
				GMT_Report (API, GMT_MSG_ERROR, "Your -E%g overrides the nodata_value of %g found in the ESRI file\n", Ctrl->E.nodata, value);
			}
			else
				Ctrl->E.nodata = value;
		}
		else {	/* Instead got the very first data value */
			ij = gmt_M_ijp (Grid->header, row, col);
			value = atof (line);
			Grid->data[ij] = (value == Ctrl->E.nodata) ? GMT->session.f_NaN : (gmt_grdfloat) value;
			if (++col == Grid->header->n_columns) col = 0, row++;
			n_left--;
		}
		while (fscanf (fp, "%lf", &value) == 1 && n_left) {
			ij = gmt_M_ijp (Grid->header, row, col);
			Grid->data[ij] = (value == Ctrl->E.nodata) ? GMT->session.f_NaN : (gmt_grdfloat) value;
			if (++col == Grid->header->n_columns) col = 0, row++;
			n_left--;
		}
		gmt_fclose (GMT, fp);
		if (n_left) {
			GMT_Report (API, GMT_MSG_ERROR, "Expected %" PRIu64 " points, found only %" PRIu64 "\n", Grid->header->nm, Grid->header->nm - n_left);
			Return (GMT_RUNTIME_ERROR);
		}
		if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
			Return (API->error);
		}
		Return (GMT_NOERROR);
	}

	/* Here we will read either x,y,z or z data, using -R -I [-r] for sizeing */

	no_data_f = (GMT->common.d.active[GMT_IN]) ? (gmt_grdfloat)GMT->common.d.nan_proxy[GMT_IN] : GMT->session.f_NaN;

	/* Set up and allocate output grid [note: zero padding specified since no BCs required] */
	if ((Grid = GMT_Create_Data (API, GMT_IS_GRID, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, NULL, NULL, \
		GMT_GRID_DEFAULT_REG, 0, NULL)) == NULL) Return (API->error);

	/* See if we have a projection info to add */
	if (GMT->common.J.active)		/* Convert the GMT -J<...> into a proj4 string and save it in the header */
		Grid->header->ProjRefPROJ4 = gmt_export2proj4(GMT);

	Amode = Ctrl->A.active ? Ctrl->A.mode : 'm';

	/* For Amode = 'd' or 'S' we need a second grid, and also require a minimum of 2 points per grid */
	if (Amode == 'd' || Amode == 'S') {
		if ((data = gmt_M_memory_aligned (GMT, NULL, Grid->header->nm, gmt_grdfloat)) == NULL) {
			GMT_Report (API, GMT_MSG_ERROR, "Unable to allocate secondary data memory\n");
			Return (GMT_MEMORY_ERROR);
		}
		n_min = 2;
	}

	if (GMT->common.b.active[GMT_IN] && gmt_M_type (GMT, GMT_IN, GMT_Z) & GMT_IS_RATIME && GMT->current.io.fmt[GMT_IN][GMT_Z].type == GMT_FLOAT) {
		GMT_Report (API, GMT_MSG_WARNING, "Your single precision binary input data are unlikely to hold absolute time coordinates without serious truncation.\n");
		GMT_Report (API, GMT_MSG_WARNING, "You must use double precision when storing absolute time coordinates in binary data tables.\n");
	}

	if (Ctrl->D.active && gmt_decode_grd_h_info (GMT, Ctrl->D.information, Grid->header)) {
		gmt_M_free (GMT, data);
		Return (GMT_PARSE_ERROR);
	}

	GMT_Report (API, GMT_MSG_INFORMATION, "n_columns = %d  n_rows = %d  nm = %" PRIu64 "  size = %" PRIuS "\n", Grid->header->n_columns, Grid->header->n_rows, Grid->header->nm, Grid->header->size);

	if (gmt_M_err_fail (GMT, gmt_set_z_io (GMT, &io, Grid), Ctrl->G.file)) Return (GMT_RUNTIME_ERROR);

	gmt_set_xy_domain (GMT, wesn, Grid->header);	/* May include some padding if gridline-registered */
	if (Ctrl->Z.active && GMT->common.d.active[GMT_IN]) {
		was = GMT->common.d.first_col[GMT_IN];
		if (gmt_M_is_fnan (no_data_f))	/* No point testing since nan_proxy is NaN... */
			GMT->common.d.active[GMT_IN] = false;
		else if (GMT->common.d.first_col[GMT_IN] == GMT_Z)	/* Change the default to the only column available */
			GMT->common.d.first_col[GMT_IN] = GMT_X;
	}

	if (Ctrl->Z.active) {	/* Need to override input method since reading single input column as z (not x,y) */
		zcol = GMT_X;
		save_i = GMT->current.io.input;
		previous_bin_i = GMT->common.b.active[GMT_IN];
		GMT->current.io.input = gmt_z_input;		/* Override and use chosen input mode */
		GMT->common.b.active[GMT_IN] = io.binary;	/* May have to set binary as well */
		GMT->current.io.fmt[GMT_IN][zcol].type = gmt_get_io_type (GMT, Ctrl->Z.type);
	}
	else {
		zcol = GMT_Z;
		if ((flag = gmt_M_memory (GMT, NULL, Grid->header->nm, unsigned int)) == NULL) {	/* No padding needed for flag array */
			GMT_Report (API, GMT_MSG_ERROR, "Unable to allocate flag memory\n");
			gmt_M_free (GMT, data);
			Return (GMT_MEMORY_ERROR);
		}
		gmt_M_memset (Grid->header->pad, 4, unsigned int);	/* Algorithm below expects no padding; we repad at the end */
		GMT->current.setting.io_nan_records = false;	/* Cannot have x,y as NaNs here */
	}

	if ((error = GMT_Set_Columns (API, GMT_IN, (unsigned int)n_req, GMT_COL_FIX_NO_TEXT)) != GMT_NOERROR) {
		gmt_M_free (GMT, data);
		gmt_M_free (GMT, flag);
		Return (error);
	}
	/* Initialize the i/o since we are doing record-by-record reading/writing */
	if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_NOERROR) {
		gmt_M_free (GMT, data);
		gmt_M_free (GMT, flag);
		Return (API->error);	/* Establishes data input */
	}
	if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_NOERROR) {
		gmt_M_free (GMT, data);		gmt_M_free (GMT, flag);
		Return (API->error);	/* Enables data input and sets access mode */
	}

	n_read = ij = 0;
	if (Ctrl->Z.active) {
		size_t nr = 0;
		if (Ctrl->Z.swab) GMT_Report (API, GMT_MSG_INFORMATION, "Binary input data will be byte swapped\n");
		for (i = 0; i < io.skip; i++)
			nr += fread (&c, sizeof (char), 1, API->object[API->current_item[GMT_IN]]->fp);
		GMT_Report (API, GMT_MSG_DEBUG, "Bytes read via -Z: % " PRIuS "\n", nr);
	}

	do {	/* Keep returning records until we reach EOF */
		if ((In = GMT_Get_Record (API, GMT_READ_DATA, NULL)) == NULL) {	/* Read next record, get NULL if special case */
			if (gmt_M_rec_is_error (GMT)) {		/* Bail if there are any read errors */
				gmt_M_free (GMT, data);		gmt_M_free (GMT, flag);
				Return (GMT_RUNTIME_ERROR);
			}
			else if (gmt_M_rec_is_eof (GMT)) 		/* Reached end of file */
				break;
			continue;	/* Go back and read the next record */
		}

		if (In->data == NULL) {
			gmt_quit_bad_record (API, In);
			Return (API->error);
		}

		/* Data record to process */
		n_read++;
		in = In->data;	/* Only need to process numerical part here */

		if (Ctrl->Z.active) {	/* Read separately because of all the possible formats */
			if (ij == io.n_expected) {
				GMT_Report (API, GMT_MSG_ERROR, "More than %" PRIu64 " records, only %" PRIu64 " was expected (aborting)!\n", ij, io.n_expected);
				GMT_Report (API, GMT_MSG_ERROR, "(You are possibly mistaking xyz2grd for an interpolator; see 'greenspline|surface|triangulate' man pages)\n");
				gmt_M_free (GMT, data);		gmt_M_free (GMT, flag);
				Return (GMT_RUNTIME_ERROR);
			}
			ij_gmt = io.get_gmt_ij (&io, Grid, ij);	/* Convert input order to output node (with padding) as per -Z */
			Grid->data[ij_gmt] = (gmt_input_col_is_nan_proxy (GMT, in[zcol], zcol)) ? GMT->session.f_NaN : (gmt_grdfloat)in[zcol];
			ij++;
		}
		else {	/* Get x, y, z */
			if (gmt_M_y_is_outside (GMT, in[GMT_Y],  wesn[YLO], wesn[YHI])) continue;	/* Outside y-range */
			if (gmt_x_is_outside (GMT, &in[GMT_X], wesn[XLO], wesn[XHI])) continue;	/* Outside x-range */

			/* Ok, we are inside the region - process data */

			scol = (int)gmt_M_grd_x_to_col (GMT, in[GMT_X], Grid->header);
			if (scol == -1) scol++, n_confused++;
			col = scol;
			if (col == Grid->header->n_columns) col--, n_confused++;
			srow = (int)gmt_M_grd_y_to_row (GMT, in[GMT_Y], Grid->header);
			if (srow == -1) srow++, n_confused++;
			row = srow;
			if (row == Grid->header->n_rows) row--, n_confused++;
			ij = gmt_M_ij0 (Grid->header, row, col);	/* Because padding is turned off we can use ij for both Grid and flag */
			switch (Amode) {
				case 'f':	/* Want the first value to matter only */
					if (flag[ij] == 0) {	/* Assign first value and that is the end of it */
						Grid->data[ij] = (gmt_grdfloat)in[zcol];
						flag[ij] = (unsigned int)n_read;
					}
					break;
				case 's':	/* Want the last value to matter only */
					Grid->data[ij] = (gmt_grdfloat)in[zcol];	/* Assign last value and that is it */
					flag[ij] = (unsigned int)n_read;
					break;
				case 'l':	/* Keep lowest value */
					if (flag[ij]) {	/* Already assigned the first value */
						if (in[zcol] < (double)Grid->data[ij]) Grid->data[ij] = (gmt_grdfloat)in[zcol];
					}
					else {	/* First time, just assign the current value */
						Grid->data[ij] = (gmt_grdfloat)in[zcol];
					}
					flag[ij]++;
					break;
				case 'u':	/* Keep highest value */
					if (flag[ij]) {	/* Already assigned the first value */
						if (in[zcol] > (double)Grid->data[ij]) Grid->data[ij] = (gmt_grdfloat)in[zcol];
					}
					else {	/* First time, just assign the current value */
						Grid->data[ij] = (gmt_grdfloat)in[zcol];
					}
					flag[ij]++;
					break;
				case 'd':	/* Keep highest and lowest value */
					if (flag[ij]) {	/* Already assigned the first value */
						if (in[zcol] > (double)Grid->data[ij]) Grid->data[ij] = (gmt_grdfloat)in[zcol];
						if (in[zcol] < (double)data[ij]) data[ij] = (gmt_grdfloat)in[zcol];
					}
					else {	/* First time, just assign the current value */
						Grid->data[ij] = data[ij] = (gmt_grdfloat)in[zcol];
					}
					flag[ij]++;
					break;
				case 'S': 	/* Add up squares and means to compute standard deviation */
					data[ij] += (gmt_grdfloat)in[zcol];
					/* Intentionally fall through - This adds up the means; we fall through to next case on purpose to also add up squares */
				case 'r': 	/* Add up squares in case we must rms */
					Grid->data[ij] += (gmt_grdfloat)in[zcol] * (gmt_grdfloat)in[zcol];
					flag[ij]++;
					break;
				default:	/* Add up in case we must sum or mean */
					Grid->data[ij] += (gmt_grdfloat)in[zcol];
					flag[ij]++;
					break;
			}
			n_used++;
		}
	} while (true);

	if (GMT_End_IO (API, GMT_IN, 0) != GMT_NOERROR) {	/* Disables further data input */
		gmt_M_free (GMT, flag);		gmt_M_free (GMT, data);
		Return (API->error);
	}

	if (Ctrl->Z.active) {
		GMT->current.io.input = save_i;	/* Reset pointer */
		GMT->common.b.active[GMT_IN] = previous_bin_i;	/* Reset binary */
		if (GMT->common.d.active[GMT_IN])	/* Reset column */
			GMT->common.d.first_col[GMT_IN] = was;
		if (ij != io.n_expected) {	/* Input amount does not match expectations */
			GMT_Report (API, GMT_MSG_ERROR, "Found %" PRIu64 " records, but %" PRIu64 " was expected (aborting)!\n", ij, io.n_expected);
				GMT_Report (API, GMT_MSG_ERROR, "(You are probably misterpreting xyz2grd with an interpolator; see 'surface' man page)\n");
			gmt_M_free (GMT, flag);		gmt_M_free (GMT, data);
			Return (GMT_RUNTIME_ERROR);
		}
		gmt_check_z_io (GMT, &io, Grid);	/* This fills in missing periodic row or column */
	}
	else {	/* xyz data could have resulted in duplicates */
		if (gmt_M_grd_duplicate_column (GMT, Grid->header, GMT_IN)) {	/* Make sure longitudes got replicated */
			uint64_t ij_west, ij_east;

			for (row = 0; row < Grid->header->n_rows; row++) {	/* For each row, look at west and east bin */
				ij_west = gmt_M_ij0 (Grid->header, row, 0);
				ij_east = gmt_M_ij0 (Grid->header, row, Grid->header->n_columns - 1);

				if (flag[ij_west] && !flag[ij_east]) {		/* Nothing in east bin, just copy from west */
					Grid->data[ij_east] = Grid->data[ij_west];
					flag[ij_east] = flag[ij_west];
					if (n_min == 2) data[ij_east] = data[ij_west];
				}
				else if (flag[ij_east] && !flag[ij_west]) {	/* Nothing in west bin, just copy from east */
					Grid->data[ij_west] = Grid->data[ij_east];
					flag[ij_west] = flag[ij_east];
					if (n_min == 2) data[ij_west] = data[ij_east];
				}
				else {	/* Both have some stuff, consolidate combined value into the west bin, then replicate to the east */
					switch (Amode) {
						case 'f':	/* Keep the first */
							if (flag[ij_east] < flag[ij_west]) Grid->data[ij_west] = Grid->data[ij_east], flag[ij_west] = flag[ij_east];
							break;
						case 's':	/* Keep the last */
							if (flag[ij_east] > flag[ij_west]) Grid->data[ij_west] = Grid->data[ij_east], flag[ij_west] = flag[ij_east];
							break;
						case 'l':	/* Keep the lowest */
							if (Grid->data[ij_east] < Grid->data[ij_west]) Grid->data[ij_west] = Grid->data[ij_east];
							flag[ij_west] += flag[ij_east];
							break;
						case 'd':	/* Keep the lowest in 'data' */
							if (data[ij_east] < data[ij_west]) data[ij_west] = data[ij_east];
							/* Intentionally fall through - since range also needs the highsets */
						case 'u':	/* Keep the highest */
							if (Grid->data[ij_east] > Grid->data[ij_west]) Grid->data[ij_west] = Grid->data[ij_east];
							flag[ij_west] += flag[ij_east];
							break;
						case 'S':	/* Sum up the sums in 'data' */
							data[ij_west] += data[ij_east];
							/* Intentionally fall through */
						default:	/* Add up in case we must sum, rms, mean, or standard deviation */
							Grid->data[ij_west] += Grid->data[ij_east];
							flag[ij_west] += flag[ij_east];
							break;
					}
					/* Replicate: */
					Grid->data[ij_east] = Grid->data[ij_west];
					flag[ij_east] = flag[ij_west];
				}
			}
		}

		for (ij = 0; ij < Grid->header->nm; ij++) {	/* Check if all nodes got one value only */
			if (flag[ij] < n_min) {	/* Cells are not filled enough */
				n_empty++;
				Grid->data[ij] = no_data_f;
			}
			else {	/* Enough values went to this node */
				if (Amode == 'n')
					Grid->data[ij] = (gmt_grdfloat)flag[ij];
				else if (Amode == 'm') {
					Grid->data[ij] /= (gmt_grdfloat)flag[ij];
					if (flag[ij] > 1) n_ave++;
				}
				else if (Amode == 'r')
					Grid->data[ij] = (gmt_grdfloat)sqrt (Grid->data[ij] / (gmt_grdfloat)flag[ij]);
				else if (Amode == 'd')
					Grid->data[ij] -= data[ij];
				else if (Amode == 'S') {
					Grid->data[ij] = (gmt_grdfloat)sqrt ((Grid->data[ij] - data[ij] * data[ij] / (gmt_grdfloat)flag[ij]) / (gmt_grdfloat)(flag[ij] - 1));

				}
				/* Implicit else means return the currently stored value */
				n_filled++;
			}
		}
		if (n_ave) GMT_Report (API, GMT_MSG_INFORMATION, "%" PRIu64 " nodes had multiple entries and were averaged\n", n_ave);
		gmt_M_free (GMT, flag);
		gmt_M_free (GMT, data);

		if (gmt_M_is_verbose (GMT, GMT_MSG_WARNING)) {
			char line[GMT_BUFSIZ], e_value[GMT_LEN32];
			sprintf (line, "%s\n", GMT->current.setting.format_float_out);
			(GMT->common.d.active[GMT_IN]) ? sprintf (e_value, GMT->current.setting.format_float_out, GMT->common.d.nan_proxy[GMT_IN]) : sprintf (e_value, "NaN");
			GMT_Report (API, GMT_MSG_INFORMATION, "Data records read: %" PRIu64 "  used: %" PRIu64 "  nodes filled: %" PRIu64 " nodes empty: %" PRIu64 " [set to %s]\n",
				n_read, n_used, n_filled, n_empty, e_value);
			if (n_confused) GMT_Report (API, GMT_MSG_WARNING, "%" PRIu64 " values gave bad indices: Pixel vs Gridline registration confusion?\n", n_confused);
		}
	}

	gmt_grd_pad_on (GMT, Grid, GMT->current.io.pad);	/* Restore padding */
	if (GMT_Set_Comment (API, GMT_IS_GRID, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, Grid)) {
		gmt_M_free (GMT, flag);		gmt_M_free (GMT, data);
		Return (API->error);
	}
	if (GMT_Write_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->G.file, Grid) != GMT_NOERROR) {
		gmt_M_free (GMT, flag);		gmt_M_free (GMT, data);
		Return (API->error);
	}

	gmt_M_free (GMT, flag);		gmt_M_free (GMT, data);
	Return (GMT_NOERROR);
}
