/* gtopo30_imgbm.c <gtopo30.DTM> > gtopo.imgbm
 *
 * read a gtopo30 file and write average elevation in imgbm format to stdout
 *
 * 14 may 97, whfs
 *
 * Assumes for now, that -9999 is always the NaN value, and that the file
 * names and structures are as given in the web page documentation.
 *
 * Also, no provision is made for the possible case that an img pixel
 * straddles data from two different gtopo30 files.  This could happen
 * (straddle in lat, never in lon), so after all outputs from this program
 * are assembled, they should be run through imgbm_combine 
 *
 * No implementation of masks yet.  Could do that later.
 *
 */

#include <stdio.h>
#include <math.h>

#define NX 10800	/* Img dimensions  */
#define NY 6336
#define	NX_IN 4800	/* GTOPO30 dimensions, except Antarctica  */
#define NY_IN 6000
#define	NX_IN_A 7200	/* GTOPO30 dimensions, Antarctica only  */
#define NY_IN_A 3600
#define NX_OUT_MAX 1800	/* Max possible pixels of output (Antarctica)  */
#define GTOPO30_NAN -9999	/* hardwired for now; could read HDR file to get it  */

struct IMGPAIR {
        int     index;
        int     depth;
} imgpair;

double	a[NX_OUT_MAX];	/* Average elevation.  */
int	n[NX_OUT_MAX];	/* N points averaged  */
int	io[NX_OUT_MAX];	/* i index of imgfile for output  */
short	buffer[NX_IN_A];	/* input buffer for a row of input  */


main(argc, argv)
int	argc;
char	**argv;
{
	
	int	north, west, parse_gtopo30_filename();
	int	nx_in, nx_out, ny_in, i_in, j_in, i_out, k;
	int	jpixel, last_jpixel;
	double	dlat, rlat, ypixel, radius;
	FILE	*fp;




	if (argc != 2 || (fp = fopen(argv[1], "r")) == NULL || (parse_gtopo30_filename(argv[1], &north, &west))) {
		fprintf(stderr,"usage:  gtopo30_imgbm <gtopo3.DTM filename>  >  gtopo30.imgbm\n");
		exit(-1);
	}

	radius = (NX/2)/M_PI;

	nx_in = (north == -60) ? NX_IN_A : NX_IN;
	ny_in = (north == -60) ? NY_IN_A : NY_IN;
	nx_out = nx_in / 4;

	k = (west < 0) ? NX + 30 * west : 30 * west;
	for (i_out = 0; i_out < nx_out; i_out++) {
		io[i_out] = (k + i_out)%NX;
	}

	last_jpixel = -1;
	for (j_in = 0; j_in < ny_in; j_in++) {
		dlat = north - (j_in + 0.5) / 120.0;
		if ( (fread((char *)buffer, 2, nx_in, fp)) != nx_in) {
			fprintf(stderr,"gtopo30_imgbm  ERROR reading input row %ld\n", j_in);
			exit(-1);
		}
		if (fabs(dlat) > 72.005977) continue;
		rlat = M_PI * dlat / 180.0;
		ypixel = NY/2 - radius*log(tan(M_PI_4 + 0.5*rlat));
		if (ypixel < 0.0 || ypixel >= (double)NY) continue;
		jpixel = floor(ypixel);
		if (jpixel != last_jpixel) {
			/* Starting a new row.  */
			if (last_jpixel != -1) {
				/* We have a row done already.  Write out and zero memory. */
				k = jpixel * NX;
				for (i_out = 0; i_out < nx_out; i_out++) {
					if (n[i_out] > 0) {
						if (n[i_out] > 1) a[i_out] /= n[i_out];
						imgpair.index = k + io[i_out];
						imgpair.depth = nint(a[i_out]);
						fwrite((char *)&imgpair, sizeof(struct IMGPAIR), 1, stdout);
					}
					n[i_out] = 0;
					a[i_out] = 0.0;
				}
			}
			else {
				/* Initialize memory for first time  */
				for (i_out = 0; i_out < nx_out; i_out++) {
					n[i_out] = 0;
					a[i_out] = 0.0;
				}
			}
			last_jpixel = jpixel;
		}
		/* Load input into averages  */
		for (i_in = 0; i_in < nx_in; i_in++) {
			if (buffer[i_in] == GTOPO30_NAN) continue;
			i_out = i_in/4;
			a[i_out] += buffer[i_in];
			n[i_out]++;
		}
	}
	fclose(fp);

	/* Write out final row.  */
	for (i_out = 0; i_out < nx_out; i_out++) {
		if (n[i_out] > 0) {
			if (n[i_out] > 1) a[i_out] /= n[i_out];
			imgpair.index = k + io[i_out];
			imgpair.depth = nint(a[i_out]);
			fwrite((char *)&imgpair, sizeof(struct IMGPAIR), 1, stdout);
		}
	}
	exit(0);

}

int	parse_gtopo30_filename(fname, n, w)
char	fname[];
int	*n, *w;
{
	char	s1[4];
	int	i, lon, lat;

	if (strlen(fname) < 11) return(-1);

	strncpy(s1, &fname[1], 3);
	s1[3] = '\0';
	i = atoi(s1);
	if (fname[0] == 'W') {
		lon = -i;
	}
	else if (fname[0] == 'E') {
		lon = i;
	}
	else {
		return(-1);
	}

	strncpy(s1, &fname[5], 2);
	s1[2] = '\0';
	i = atoi(s1);
	if (fname[4] == 'S') {
		lat = -i;
	}
	else if (fname[4] == 'N') {
		lat = i;
	}
	else {
		return(-1);
	}

	if (lat != 90 && lat != 40 && lat != -10 && lat != -60) return(-1);

	i = abs(lon);
	if (lat == -60) {
		if (i != 0 && i != 60 && i != 120 && lon != -180) return(-1);
	}
	else {
		if (i != 20 && i != 60 && i != 100 && i != 140 && lon != -180) return(-1);
	}
	*n = lat;
	*w = lon;
	return(0);
}