/*----------------------------------------------------------------------*\
 | Program to create a subset of etopo5, and optionally add a constant	|
 | to every cell value.													|
 |																		|
 | Usage: subset_etopo5 [lat1] [lon1] [lat2] [lon2]						|
 |						[-add offset]									|
 |						[-verbose]										|
 |						[-byteorder natural | LSBfirst | MSBfirst]		|
 |						[-output output_file]							|
 |						input_file										|
 |																		|
 | where	lat1 and lat2 are decimal degrees (possibly including the	|
 |			decimal point) followed immediately by 'N' or 'S' with no	|
 |			intervening spaces.  Only positive values are permitted.	|
 |																		|
 |			lon1 and lon2 are decimal degrees (possibly including the	|
 |			decimal point) followed immediately by 'E' or 'W' with no	|
 |			intervening spaces.  Only positive values are permitted.	|
 |																		|
 |			offset is an integer										|
 |																		|
 |			the argument to -byteorder is one of the text strings		|
 |			"natural", "LSBfirst", or "MSBfirst"						|
 |																		|
 | lat1, lon1, lat2, and lon2 define the (rectangular) subset that will	|
 | be output.  The subset includes the row containing min(lat1,lat2)	|
 | and the column containing min(lon1,lon2); it does not include the	|
 | row containing max(lat1,lat2), nor does it include the column		|
 | containing max(lon1,lon2).											|
 |																		|
 | If only one latitude is specified, the other is assumed to be 90S.	|
 | If only one longitude is specified, the other is assumed to be 180E.	|
 | If no latitudes or longitudes are specified, the entire image will	|
 | be output.  This allows you to add a constant to etopo5 or to alter	|
 | its byte ordering.													|
 |																		|
 | -add may be abbreviated as "-a"										|
 | -byteorder may be abbreviated as "-byte"								|
 | -output may be abbreviated as "-o"									|
 |																		|
 | If -byteorder is specified, it must be followed with a valid string,	|
 | and the strings have the following meanings:							|
 |		LSBfirst: The output file will be written little-endian.		|
 |		MSBfirst: The output file will be written big-endian.			|
 |		natural: The output file will be written using the byte order	|
 |				 that is natural for the current machine.				|
 |																		|
 | Peter N. Schweitzer (U.S. Geological Survey, Reston VA 22092)		|
\*----------------------------------------------------------------------*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>

#define WIDTH	360*12
#define HEIGHT	180*12
#define MAXVAL	65535

#define nint(x) ((x - floor(x) < 0.5) ? (int)(x) : (int)(x + 0.5))

enum byte_order_t {
	LITTLE_ENDIAN,
	BIG_ENDIAN
	} byte_order;

static enum byte_order_t get_byte_order (void) {
	short w = 255;
	unsigned char *s = (unsigned char *) &w;
	if (*s) return (LITTLE_ENDIAN);
	else return (BIG_ENDIAN);
	}

main (int argc, char *argv[]) {
	short i,j;
	int x,y;
	FILE *in, *out;
	char *input_file, *output_file;
	int x1,y1,x2,y2;
	double v;
	char *s;
	int add = 0;
	int verbose = 0;

	long k,n;
	short *depth;
	short *p;
	char scratch [WIDTH * sizeof(short)];
	enum byte_order_t byte_order = get_byte_order();
	int swap_bytes = 0;

	if (argc <= 1) {
		fprintf (stderr,"Usage: %s [lat1] [lon1] [lat2] [lon2] [-add offset] [-byteorder natural | LSBfirst | MSBfirst] [-output output_file] [-verbose] input_file\n",argv[0]);
		exit (0);
		}

	input_file  = NULL;
	output_file = NULL;
	x1 = y1 = x2 = y2 = -999;

	for (i=0; i < argc; i++)
		if (isdigit (*argv[i])) {
			v = strtod (argv[i],&s);
			switch (toupper(*s)) {
				case 'N':
					v = 12.0 * (90.0 - v);
					y = nint (v);
					if (y1 == -999) y1 = y;
					else
						if (y2 == -999) y2 = y;
						else {
							fprintf (stderr,"Error: you specified more than two latitude values.\n");
							fprintf (stderr,"       The third (argument %d) was \"%s\"\n",i,argv[i]);
							exit (1);
							}
					break;
				case 'S':
					v = 12 * (90 + v);
					y = nint (v);
					if (y1 == -999) y1 = y;
					else
						if (y2 == -999) y2 = y;
						else {
							fprintf (stderr,"Error: you specified more than two latitude values.\n");
							fprintf (stderr,"       The third (argument %d) was \"%s\"\n",i,argv[i]);
							exit (1);
							}
					break;
				case 'E':
					v = 12.0 * (180.0 + v);
					x = nint (v);
					if (x > WIDTH) x -= WIDTH;
					if (x1 == -999) x1 = x;
					else
						if (x2 == -999) x2 = x;
						else {
							fprintf (stderr,"Error: you specified more than two longitude values.\n");
							fprintf (stderr,"       The third (argument %d) was \"%s\"\n",i,argv[i]);
							exit (1);
							}
					break;
				case 'W':
					v = 12.0 * (180.0 - v);
					x = nint (v);
					if (x > WIDTH) x -= WIDTH;
					if (x1 == -999) x1 = x;
					else
						if (x2 == -999) x2 = x;
						else {
							fprintf (stderr,"Error: you specified more than two longitude values.\n");
							fprintf (stderr,"       The third (argument %d) was \"%s\"\n",i,argv[i]);
							exit (1);
							}
					break;
				default:
					fprintf (stderr,"Error: geographic coordinates must include N,S,E, or W\n");
					fprintf (stderr,"       argument %d (\"%s\") does not.\n",i,argv[i]);
					exit (1);
					break;
				}
			}
		else
			if (memcmp (argv[i],"-o",2) == 0) {
				i++;
				output_file = argv[i];
				}
			else
				if (memcmp (argv[i],"-add",4) == 0) {
					i++;
					add = atoi (argv[i]);
					}
				else
					if (memcmp (argv[i],"-byte",5) == 0) {
						i++;
						if (strcmp (argv[i],"natural") == 0) swap_bytes = 0;
						else
							if (strcmp (argv[i],"LSBfirst") == 0)
								swap_bytes = ((byte_order == LITTLE_ENDIAN) ? 0 : 1);
							else
								if (strcmp (argv[i],"MSBfirst") == 0)
									swap_bytes = ((byte_order == LITTLE_ENDIAN) ? 1 : 0);
								else {
									fprintf (stderr,"Error: -byteorder must be natural, LSBfirst, or MSBfirst; not \"%s\"\n",argv[i]);
									exit (1);
									}
						}
				    else
				    	if (memcmp (argv[i],"-v",2) == 0)
				    		verbose = 1;
				    	else
							input_file = argv[i];

	if (x1 > x2) { x = x1; x1 = x2; x2 = x; }
	if (y1 > y2) { y = y1; y1 = y2; y2 = y; }

	if (x1 == -999) x1 = 0;
	if (x2 == -999) x2 = WIDTH;
	if (y1 == -999) y1 = 0;
	if (y2 == -999) y2 = HEIGHT;

	if (!input_file) input_file = "etopo5.img";
	if (!output_file) output_file = "subset.img";

	if (verbose) {
		fprintf (stderr,"x1 = %d, y1 = %d, x2 = %d, y2 = %d\n",x1,y1,x2,y2);
		fprintf (stderr,"add %d to all values\n",add);
		fprintf (stderr,"input_file  = \"%s\"\n",input_file);
		fprintf (stderr,"output_file = \"%s\"\n",output_file);
		fprintf (stderr,"              bytes will ");
		if (!swap_bytes) fprintf (stderr,"not ");
		fprintf (stderr,"be swapped in the output.\n");
		fprintf (stderr,"Resulting image will have width %d and height %d\n",x2-x1,y2-y1);
		}

	if (in = fopen (input_file,"rb")) {
		if (depth = (short *) malloc (WIDTH * sizeof (short))) {
			if (out = fopen (output_file,"wb")) {
				for (y=0; y < y2; y++)
					if ((n = fread (scratch,sizeof(short),WIDTH,in)) == WIDTH) {
						if (y >= y1) {
							if (byte_order != LITTLE_ENDIAN)
								swab (scratch,(char *)depth,n*sizeof(short));
							else
								memcpy (depth,scratch,n*sizeof(short));
							if (add) for (x=0, p=depth; x < WIDTH; x++) *p++ += add;
							if (swap_bytes) {
								swab ((char *)depth,scratch,n*sizeof(short));
								memcpy (depth,scratch,n*sizeof(short));
								}
							p = depth + x1;
							fwrite (p,sizeof(short),x2-x1,out);
							}
						}
					else
						fprintf (stderr,"Error: only %ld of %ld words read from input file\n",n,WIDTH);
				fclose (out);
				}
				else
					fprintf (stderr,"Error: could not create output file %s\n",output_file);
			free (depth);
			}
		else
			fprintf (stderr,"Error: could not allocate array for input data\n");
		fclose (in);
		}
	else
		fprintf (stderr,"Error: could not open input file %s\n",input_file);

	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/