/* srtm30_imgbm.c <srtm30.DTM> > srtm.imgbm
 *
 * read a srtm30 file and write average elevation in imgbm format to stdout
 *
 * 14 may 97, whfs
 * 24 may 07, dts - modified for 1 minute processing
 *
 * 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 srtm30 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 
 *
 * Changed to omit zero elevation
 *
 */

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

#define NX 21600	/* Img dimensions  */
#define NY 17280
#define MINUTE 1
#define	NX_IN 4800	/* SRTM30 dimensions, except Antarctica  */
#define NY_IN 6000
#define	NX_IN_A 7200	/* SRTM30 dimensions, Antarctica only  */
#define NY_IN_A 3600
#define NX_OUT_MAX 3600	/* Max possible pixels of output (Antarctica)  */
#define SRTM30_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	bufin[NX_IN_A],buffer[NX_IN_A];	/* input buffer for a row of input  */

main(argc, argv)
int	argc;
char	**argv;
{
	int	north, west, parse_srtm30_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_srtm30_filename(argv[1], &north, &west))) {
		fprintf(stderr,"usage:  srtm30_imgbm <srtm3.DTM filename>  >  srtm30.imgbm\n");
		exit(-1);
	}

/* radous of the earh in pixels */
	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/(2*MINUTE);

/*  compute the column indices for the output file */
	k = (west < 0) ? NX + west * 60/MINUTE : west * 60/MINUTE;
	for (i_out = 0; i_out < nx_out; i_out++) {
		io[i_out] = (k + i_out)%NX;
	}

/* loop over input rorw */
	for (j_in = 0; j_in < ny_in; j_in++) {

/* zero the rows of the putput */
		for (i_out = 0; i_out < nx_out; i_out++) {
			n[i_out] = 0;
			a[i_out] = 0.0;
		}

/* compute the row indices of the output file and read a row of input */
		dlat = north - (j_in + 0.5) / 120.0;
		if ( (fread((char *)bufin, 2, nx_in, fp)) != nx_in) {
			fprintf(stderr,"srtm30_imgbm  ERROR reading input row %ld\n", j_in);
			exit(-1);
		}
                swap16(bufin,buffer,nx_in);
		if (fabs(dlat) > 80.738000) 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);

/* average the input line to create an output line */
		for (i_in = 0; i_in < nx_in; i_in++) {
			if (buffer[i_in] == SRTM30_NAN || buffer[i_in] == 0) continue;
			i_out = i_in/(2*MINUTE);
			a[i_out] += buffer[i_in];
			n[i_out]++;
		}

/* write the output line and zero the memory */
		for (i_out = 0; i_out < nx_out; i_out++) {
                        if (n[i_out] <= MINUTE) continue;
			a[i_out] /= n[i_out];
			imgpair.index = jpixel*NX + io[i_out];
			imgpair.depth = (int)floor(a[i_out]+.5);
			fwrite((char *)&imgpair, sizeof(struct IMGPAIR), 1, stdout);
		}
	}
	fclose(fp);
	exit(0);
}

int	parse_srtm30_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' || fname[0] == 'W') {
		lon = -i;
	}
	else if (fname[0] == 'e' || fname[0] == 'E') {
		lon = i;
	}
	else {
		return(-1);
	}

	strncpy(s1, &fname[5], 2);
	s1[2] = '\0';
	i = atoi(s1);
	if (fname[4] == 's' || fname[4] == 'S') {
		lat = -i;
	}
	else if (fname[4] == 'n' || 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);
}
/************************************************************************
* Creator: David T. Sandwell    Scripps Institution of Oceanography    *
* Date   : 09/12/93             Copyright, David T. Sandwell           *
************************************************************************/

swap16( in, out, n )    /* Swaps 2 bytes within each 16-bit word of
                                array in. */
char    *in;            /* Input array */
char    *out;           /* Output array */
int     n;              /* # of short integers to swap */

{
register char *ip, *op; /* Local register variables */

        if( n > 0 )     /* Make sure n is positive */
        {
                ip = in+2;      /* Load the pointers into temporary registers */
                op = out;
                while( n-- )
                {
                        *op++ = *--ip;  /* Do the swap */
                        *op++ = *--ip;
                        ip += 4;
                }
        }
}