/*--------------------------------------------------------------------
 *    The GMT-system:	@(#)shoredump.c	3.2  05 Sep 1995
 *
 *    shoredump extracts ASCII segments from the shoreline, rivers, and
 *    border databases.
 *
 *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
 *    See README file for copying and redistribution conditions.
 *--------------------------------------------------------------------*/
/*
 * Author:	Paul Wessel
 * Date:	21-JUN-1995
 * Version:	3.0
 *
 */

#include "../../../include/gmt.h"
#include "../../../include/gmt_shore.h"		/* Defines shore structures */

#define NX 10800
#define NY 6336
#define NMASK 8553600   /* NX * NY / 8  for bitmask  */

static unsigned char maskset[8] = {128, 64, 32, 16, 8, 4, 2, 1};
unsigned char mask[NMASK];


struct SHORE c;
struct BR b, r;

char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};

main (argc, argv)
int argc;
char **argv; {
	int i, n, np, ind, bin, base = 3, max_level = MAX_LEVEL, direction = 1, k;
	int blevels[N_BLEVELS], n_blevels = 0, rlevels[N_RLEVELS], n_rlevels = 0;
	
	BOOLEAN	error = FALSE, get_river = FALSE, shift;
	BOOLEAN	greenwich = FALSE, get_shore = FALSE, get_border = FALSE, got_fh[2];
	
	double	west = 0.0, east = 0.0, south = 0.0, north = 0.0, edge = 720.0, left, right;
	double min_area = 0.0;
	
	char res = 'l', format[100];
	
	struct POL *p;

	int	negone = -1, ipixel, jpixel, imgindex;
	double	xpixel, ypixel, radlat, radius;
		
	memset ((char *)rlevels, 0, N_RLEVELS * sizeof (int));
	memset ((char *)blevels, 0, N_BLEVELS * sizeof (int));
	memset (mask, 0, NMASK);	/* Set mask to all off, initially  */

	radius = (NX/2)/M_PI;

	
	argc = gmt_begin (argc, argv);
	
	/* Check and interpret the command line arguments */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch(argv[i][1]) {
		
				/* Common parameters */
				
				case 'R':
				case 'V':
				case '\0':
					error += get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */

				case 'A':
					n = sscanf (&argv[i][2], "%lf/%d", &min_area, &max_level);
					if (n == 1) max_level = MAX_LEVEL;
					break;
				case 'D':
					res = argv[i][2];
					switch (res) {
						case 'f':
							base = 0;
							break;
						case 'h':
							base = 1;
							break;
						case 'i':
							base = 2;
							break;
						case 'l':
							base = 3;
							break;
						case 'c':
							base = 4;
							break;
						default:
							fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Unknown modifier %c\n", gmt_program, argv[i][2]);
							break;
					}
					break;
				case 'N':
					if (!argv[i][2]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR:  -N option takes one argument\n", gmt_program);
						error++;
						continue;
					}
					get_border = TRUE;
					switch (argv[i][2]) {
						case 'a':
							for (k = 0; k < N_BLEVELS; k++) blevels[k] = TRUE;
							break;
						default:
							k = atoi (&argv[i][2]) - 1;
							if (k < 0 || k >= N_BLEVELS) {
								fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Feature not in list!\n", gmt_program);
								error++;
							}
							else
								blevels[k] = TRUE;
							break;
					}
					break;
					
				case 'I':
					if (!argv[i][2]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR:  -I option takes one argument\n", gmt_program);
						error++;
						continue;
					}
					get_river = TRUE;
						
					switch (argv[i][2]) {
						case 'a':
							for (k = 0; k < N_RLEVELS; k++) rlevels[k] = TRUE;
							break;
						case 'r':
							for (k = 0; k < 4; k++) rlevels[k] = TRUE;
							break;
						case 'i':
							for (k = 5; k < 8; k++) rlevels[k] = TRUE;
							break;
						case 'c':
							rlevels[8] = rlevels[9] = TRUE;
							break;
						default:
							k = atoi (&argv[i][2]) - 1;
							if (k < 0 || k >= N_RLEVELS) {
								fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Feature not in list!\n", gmt_program);
								error++;
							}
							else
								rlevels[k] = TRUE;
							break;
					}
					break;
					
				case 'S':
					get_shore = TRUE;
					break;
				default:		/* Options not recognized */
					error = TRUE;
					gmt_default_error (argv[i][1]);
					break;
			}
		}
	}
	
	find_resolutions (got_fh);	/* Check to see if full &| high resolution coastlines are installed */

	if (argc == 1 || gmt_quick) {	/* Display usage */
		fprintf (stderr,"shoredump %s - Extract shorelines, rivers, or borders\n\n", GMT_VERSION);
		fprintf (stderr,"usage: shoredump -R<west>/<east>/<south>/<north> [-A<min_area>[/<max_level>]]\n");
		fprintf (stderr, "	 [-D<resolution>] [-I<feature>] [-N<feature>] [-S] [-V]\n");
		
		if (gmt_quick) exit (-1);
		
		explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "	-A features smaller than <min_area> (in km^2) or of higher level (0-4) than <max_level>\n");
		fprintf (stderr, "	will be skipped [0/4 (4 means lake inside island inside lake)]\n");
		fprintf (stderr, "	-D Choose one of the following resolutions:\n");
		if (got_fh[0]) fprintf (stderr, "	   f - full resolution (may be very slow for large regions)\n");
		if (got_fh[1]) fprintf (stderr, "	   h - high resolution (may be slow for large regions)\n");
		fprintf (stderr, "	   i - intermediate resolution\n");
		fprintf (stderr, "	   l - low resolution [Default]\n");
		fprintf (stderr, "	   c - crude resolution\n");
		fprintf (stderr, "	-I extract rIvers.  Choose features below, repeat -I as many times as needed\n");
		fprintf (stderr, "	       1 = Permanent major rivers\n");
		fprintf (stderr, "	       2 = Additional major rivers\n");
		fprintf (stderr, "	       3 = Additional rivers\n");
		fprintf (stderr, "	       4 = Minor rivers\n");
		fprintf (stderr, "	       5 = Double lined rivers\n");
		fprintf (stderr, "	       6 = Intermittent rivers - major\n");
		fprintf (stderr, "	       7 = Intermittent rivers - additional\n");
		fprintf (stderr, "	       8 = Intermittent rivers - minor\n");
		fprintf (stderr, "	       9 = Major canals\n");
		fprintf (stderr, "	      10 = Minor canals\n");
		fprintf (stderr, "	       a = All rivers and canals (1-10)\n");
		fprintf (stderr, "	       r = All permanent rivers (1-4)\n");
		fprintf (stderr, "	       i = All intermittent rivers (6-8)\n");
		fprintf (stderr, "	       c = All canals (9-10)\n");
		fprintf (stderr, "	-N extract boundaries.  Choose features below, repeat -N as many times as needed\n");
		fprintf (stderr, "	       1 = National boundaries\n");
		fprintf (stderr, "	       2 = State boundaries within the Americas\n");
		fprintf (stderr, "	       3 = Marine boundaries\n");
		fprintf (stderr, "	       a = All boundaries (1-3)\n");
		fprintf (stderr,"	-S extract shorelines [Default].\n");
		explain_option ('V');
		exit (-1);
	}

	/* Check that the options selected are mutually consistant */
	
	if ((get_shore + get_border + get_river) == 0) get_shore = TRUE;	/* Default */

	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
		error++;
	}
	if ((get_shore + get_border + get_river) != 1) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify only one of -I, -N, or -S\n", gmt_program);
		error++;
	}
	
	if (get_river) {
		for (k = 0; k < N_RLEVELS; k++) {
			if (!rlevels[k]) continue;
			rlevels[n_rlevels] = k + 1;
			n_rlevels++;
		}
	}
	
	if (get_border) {
		for (k = 0; k < N_BLEVELS; k++) {
			if (!blevels[k]) continue;
			blevels[n_blevels] = k + 1;
			n_blevels++;
		}
	}

	if (error) exit (-1);

	if (east > 360.0) {
		west -= 360.0;
		east -= 360.0;
	}
	
        map_getproject ("x1d"); /* Fake linear projection */
	map_setup (west, east, south, north);

	if (get_shore && gmt_init_shore (res, &c, project_info.w, project_info.e, project_info.s, project_info.n))  {
		fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", gmt_program, shore_resolution[base]);
		exit (-1);
	}
	
	if (get_border && gmt_init_br ('b', res, &b, project_info.w, project_info.e, project_info.s, project_info.n)) {
		fprintf (stderr, "%s: %s resolution political boundary data base not installed\n", gmt_program, shore_resolution[base]);
		exit (-1);
	}
	
	if (get_river && gmt_init_br ('r', res, &r, project_info.w, project_info.e, project_info.s, project_info.n)) {
		fprintf (stderr, "%s: %s resolution river data base not installed\n", gmt_program, shore_resolution[base]);
		exit (-1);
	}

	if (west < 0.0 && east > 0.0 && project_info.projection > 5) greenwich = TRUE;
	if (gmt_world_map && greenwich)
		edge = project_info.central_meridian;
	else if (!gmt_world_map && greenwich) {
		shift = TRUE;
		edge = west;	if (edge < 0.0) edge += 360.0;
		if (edge > 180.0) edge = 180.0;
	}
	
	if (project_info.w < 0.0 && project_info.e <= 0.0) {	/* Temporarily shift boundaries */
		project_info.w += 360.0;
		project_info.e += 360.0;
	}
	
	sprintf (format, "%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format);

	if (get_shore) {	/* Read shoreline file and extract lines */

		for (ind = 0; ind < c.nb; ind++) {	/* Loop over necessary bins only */
		
			bin = c.bins[ind];
			gmt_get_shore_bin (ind, &c, min_area, max_level);
		
			if (gmtdefs.verbose) fprintf (stderr, "shoredump: Working on shoreline block # %5d\r", bin);

			if (gmt_world_map && greenwich) {
				left = c.bsize * (bin % c.bin_nx);	right = left + c.bsize;
				shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
			}
					
			if (c.ns == 0) continue;
		
			np = gmt_assemble_shore (&c, direction, FALSE, shift, edge, &p);
		
			for (i = 0; i < np; i++) {
				/* Here is where this program differs from shoredump.c  */
				if (p[i].level == 0) {
					/* I think if level is zero it is a sea-level coastline, if not, it isn't.
						So here I am only dumping the zero-level stuff.  But maybe great
						lakes are zero?  Check later with an ascii map.  */
					
					for (k = 0; k < p[i].n; k++) {
						xpixel = p[i].lon[k];
						while (xpixel < 0.0) xpixel += 360.0;
						while (xpixel >= 360.0) xpixel -= 360.0;
						ipixel = floor(30.0 * xpixel);
						if (ipixel == NX) ipixel = 0;
						radlat = M_PI * p[i].lat[k] / 180.0;
						ypixel = NY/2 - radius*log(tan(M_PI_4 + 0.5*radlat));
						if (ypixel < 0.0 || ypixel >= (double)NY) continue;
						jpixel = floor(ypixel);
						imgindex = jpixel * NX + ipixel;
						mask[imgindex/8] |= maskset[imgindex%8];
					}
				}
			}
			
			/* Free up memory */
		
			free_polygons (p, np);
			free ((char *)p);
		
			free_shore (&c);
	
		}
		shore_cleanup (&c);
	}

	else if (get_river) {	/* Read rivers file and plot as lines */
		
		if (gmtdefs.verbose) fprintf (stderr, "shoredump: Extracing Rivers...");

		for (ind = 0; ind < r.nb; ind++) {	/* Loop over necessary bins only */
		
			bin = r.bins[ind];
			gmt_get_br_bin (ind, &r, rlevels, n_rlevels);
			
			if (gmtdefs.verbose) fprintf (stderr, "shoredump: Working on river block # %5d\r", bin);

			if (r.ns == 0) continue;
			
			if (gmt_world_map && greenwich) {
				left = r.bsize * (bin % r.bin_nx);	right = left + r.bsize;
				shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
			}
			
			np = gmt_assemble_br (&r, shift, edge, &p);
			
			for (i = 0; i < np; i++) {
				printf ("> Bin # %d, River segment, level %d\n", bin, p[i].level);
				for (k = 0; k < p[i].n; k++) printf (format, p[i].lon[k], p[i].lat[k]);
			}
			
			/* Free up memory */
		
			free_polygons (p, np);
			free ((char *)p);

			free_br (&r);
		}
		br_cleanup (&r);
	}

	else if (get_border) {	/* Read borders file and plot as lines */
		
		if (gmtdefs.verbose) fprintf (stderr, "shoredump: Extracing Borders...");

		for (ind = 0; ind < b.nb; ind++) {	/* Loop over necessary bins only */
		
			bin = b.bins[ind];
			gmt_get_br_bin (ind, &b, blevels, n_blevels);
			
			if (gmtdefs.verbose) fprintf (stderr, "shoredump: Working on border block # %5d\r", bin);

			if (b.ns == 0) continue;
			
			if (gmt_world_map && greenwich) {
				left = b.bsize * (bin % b.bin_nx);	right = left + b.bsize;
				shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
			}

			np = gmt_assemble_br (&b, shift, edge, &p);
			for (i = 0; i < np; i++) {
				printf ("> Bin # %d, Border segment, level %d\n", bin, p[i].level);
				for (k = 0; k < p[i].n; k++) printf (format, p[i].lon[k], p[i].lat[k]);
			}
			
			/* Free up memory */
		
			free_polygons (p, np);
			free ((char *)p);

			free_br (&b);
		}
		br_cleanup (&b);
	}
	
	if (gmtdefs.verbose) fprintf (stderr, "Done\n");

	for (imgindex = 0; imgindex < NX * NY; imgindex++) {
		if (mask[imgindex/8] & maskset[imgindex%8]) {
			fwrite((char *)&imgindex, 1, 4, stdout);
			fwrite((char *)&negone, 1, 4, stdout);
		}
	}

	gmt_end (argc, argv);
}