/*--------------------------------------------------------------------
 *	$Id: pscoast.c,v 1.27 2004/05/23 05:52:07 pwessel Exp $
 *
 *	Copyright (c) 1991-2004 by P. Wessel and W. H. F. Smith
 *	See COPYING 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 General Public License as published by
 *	the Free Software Foundation; version 2 of the License.
 *
 *	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 General Public License for more details.
 *
 *	Contact info: gmt.soest.hawaii.edu
 *--------------------------------------------------------------------*/
/*
 * pscoast can draw shorelines, paint the landward and/or seaward side of
 * the shoreline, paint lakes, draw political borders, and rivers.  Three
 * data sets are considered: shores, borders, and rivers.  Each of these
 * data sets come in 5 different resolutions.  The lower resolution files
 * are derived from the full resolution data using the Douglas-Peucker (DP)
 * line reduction algorithm.  By giving a tolerance in km, the algorithm
 * will remove points that do not depart more than the tolerance from the
 * straight line that would connect the points if the point in question was
 * removed.  The resolutions are:
 *
 * full		- The complete World Vector Shoreline + CIA data base
 * high		- DP reduced with tolerance = 0.2 km
 * intermediate	- DP reduced with tolerance = 1 km
 * low		- DP reduced with tolerance = 5 km
 * crude	- DP reduced with tolerance = 25 km
 *
 * If selected, pscoast will open the desired binned shoreline file and fill
 * the landmasses and/or oceans/lakes as shaded/colored/textured polygons. The
 * shoreline may also be drawn.  Political boundaries and rivers may be plotted,
 * too.  If only land is filled, then the 'wet' areas remain transparent.
 * If only oceans are filled, then the 'dry' areas remain transparent.  This
 * allows the user to overlay land or ocean without wiping out the plot.
 * For more information about the binned polygon file, see the GMT Technical
 * Reference and Cookbook.
 * Optionally, the user may choose to issue clip paths rather than paint the
 * polygons.  That way one may clip subsequent images to be visible only
 * inside or outside the coastline.
 * 
 *
 * Author:	Paul Wessel
 * Date:	24-JUN-2000
 * Version:	4
 *
 */

#include "gmt.h"

#ifdef DEBUG
int debug_bin = -1;
#endif

void recursive_painting (int k0, int np, struct POL *p, int level, int start_level, struct GMT_FILL fill[]);
void recursive_clipping (int k0, int np, struct POL *p, int level, int start_level, int *n_pieces);

main (int argc, char **argv)
{
	int i, n, np, ind, bin, base = 3, anti_bin = -1, level, flag;
	int level_to_be_painted, max_level = MAX_LEVEL, direction, np_new, k, last_k;
	int blevels[N_BLEVELS], n_blevels = 0, rlevels[N_RLEVELS], n_rlevels = 0, bin_trouble;
	int start_direction, stop_direction, lowest_level = 1, min_level = 0, n_pieces = 0;
	
	BOOLEAN	error = FALSE, set_lake_fill = FALSE, draw_river = FALSE, shift = FALSE, need_coast_base, recursive;
	BOOLEAN	greenwich = FALSE, draw_coast = FALSE, fill_ocean = FALSE, fill_land = FALSE, possibly_donut_hell = FALSE;
	BOOLEAN clobber_background, draw_border = FALSE, paint_polygons = FALSE, donut, fill_in_use = FALSE;
	BOOLEAN donut_hell = FALSE, world_map_save, did_clip = FALSE;
	BOOLEAN clip_ocean = FALSE, clip_land = FALSE, end_of_clip = FALSE, clipping = FALSE, dumping = FALSE;
	
	double	west = 0.0, east = 0.0, south = 0.0, north = 0.0, bin_x[5], bin_y[5];
	double west_border, east_border, anti_lon, anti_lat, edge = 720.0, left, right;
	double *xtmp, *ytmp, min_area = 0.0, step = 0.01;
	double anti_x, anti_y, x_0, y_0, x_c, y_c, dist, out[2];
	
	char res = 'l', key[5], *string, comment[128];
	char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};

	struct GMT_PEN coast, pen, bpen[N_BLEVELS], rpen[N_RLEVELS];
	struct GMT_FILL fill[5];	/* Colors for (0) water, (1) land, (2) lakes, (3) islands in lakes, (4) lakes in islands in lakes */
	struct GMT_SHORE c;
	struct GMT_BR b, r;
	struct POL *p;
	struct MAP_SCALE ms;
	struct MAP_ROSE mr;
		
	GMT_init_fill (&fill[0], 255, 255, 255);	/* Default Ocean color */
	GMT_init_fill (&fill[1], 0, 0, 0);		/* Default Land color */
	GMT_init_fill (&fill[2], 255, 255, 255);	/* Default Lake color */
	fill[3] = fill[1];
	fill[4] = fill[2];
	GMT_init_pen (&coast, GMT_PENWIDTH);
	GMT_init_pen (&pen, GMT_PENWIDTH);
	for (k = 0; k < N_RLEVELS; k++) GMT_init_pen (&rpen[k], GMT_PENWIDTH);
	for (k = 0; k < N_BLEVELS; k++) GMT_init_pen (&bpen[k], GMT_PENWIDTH);
	
	memset ((void *)rlevels, 0, (size_t)(N_RLEVELS * sizeof (int)));
	memset ((void *)blevels, 0, (size_t)(N_BLEVELS * sizeof (int)));
	ms.plot = mr.plot = FALSE;
	
	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 'B':	/* Get tickmark information */
				case 'J':
				case 'K':
				case 'O':
				case 'P':
				case 'R':
				case 'U':
				case 'V':
				case 'X':
				case 'x':
				case 'Y':
				case 'y':
				case 'c':
				case 'b':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					break;
				
				/* Supplemental parameters */

				case 'A':
					n = sscanf (&argv[i][2], "%lf/%d/%d", &min_area, &min_level, &max_level);
					if (n == 1) min_level = 0, max_level = MAX_LEVEL;
					break;
				case 'C':	/* Lake color */
					if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[2])) {
						GMT_fill_syntax ('C');
						error++;
					}
					set_lake_fill = TRUE;
					fill[4] = fill[2];
					break;
				case 'D':
					res = argv[i][2];
					base = GMT_set_resolution (&res, 'D');
					break;
				case 'E':
					sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
					break;
				case 'G':		/* Set Gray shade, pattern, or clipping */
					if (argv[i][2] == 'c')
						clip_land = TRUE;
					else {
						fill_land = TRUE;
						if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[1])) {
							GMT_fill_syntax ('G');
							error++;
						}
						fill[3] = fill[1];
					}
					break;
				case 'L':
					error += GMT_getscale (&argv[i][2], &ms);
					break;
				case 'M':
					GMT_multisegment (&argv[i][2]);
					dumping = TRUE;
					break;

				case 'N':
					if (!argv[i][2]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR:  -N option takes at least one argument\n", GMT_program);
						error++;
						continue;
					}
					draw_border = TRUE;
					strncpy(key, "", 5);
					for (k = 2; argv[i][k] && argv[i][k] != '/'; k++);
					if (argv[i][k] && argv[i][k] == '/') {	/* Specified pen args */
						strncpy (key, &argv[i][2], (size_t)k-2);
						string = &argv[i][k+1];
					}
					else {	/* No pen specified, use default */
						strcpy (key, &argv[i][2]);
						string = CNULL;
					}
					GMT_init_pen (&pen, GMT_PENWIDTH);
					if (string && GMT_getpen (string, &pen)) {
						GMT_pen_syntax ('N');
						error++;
					}
						
					switch (key[0]) {
						case 'a':
							for (k = 0; k < N_BLEVELS; k++) blevels[k] = TRUE;
							for (k = 0; k < N_BLEVELS; k++) bpen[k] = pen;
							break;
						default:
							k = atoi (key) - 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;
								bpen[k] = pen;
							}
							break;
					}
					break;
					
				case 'I':
					if (!argv[i][2]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR:  -I option takes at least one argument\n", GMT_program);
						error++;
						continue;
					}
					draw_river = TRUE;
					strncpy(key, "", 5);
					for (k = 2; argv[i][k] && argv[i][k] != '/'; k++);
					if (argv[i][k] && argv[i][k] == '/') {	/* Specified pen args */
						strncpy (key, &argv[i][2], (size_t)k-2);
						string = &argv[i][k+1];
					}
					else {	/* No pen specified, use default */
						strcpy (key, &argv[i][2]);
						string = CNULL;
					}
					GMT_init_pen (&pen, GMT_PENWIDTH);
					if (string && GMT_getpen (string, &pen)) {
						GMT_pen_syntax ('I');
						error++;
					}
						
					switch (key[0]) {
						case 'a':
							for (k = 0; k < N_RLEVELS; k++) rlevels[k] = TRUE;
							for (k = 0; k < N_RLEVELS; k++) rpen[k] = pen;
							break;
						case 'r':
							for (k = 0; k < 4; k++) rlevels[k] = TRUE;
							for (k = 0; k < 4; k++) rpen[k] = pen;
							break;
						case 'i':
							for (k = 4; k < 7; k++) rlevels[k] = TRUE;
							for (k = 4; k < 7; k++) rpen[k] = pen;
							break;
						case 'c':
							for (k = 7; k < 10; k++) rlevels[k] = TRUE;
							for (k = 7; k < 10; k++) rpen[k] = pen;
							break;
						default:
							k = atoi (key) - 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;
								rpen[k] = pen;
							}
							break;
					}
					break;
					
				case 'S':		/* Set ocean color if needed */
					if (argv[i][2] == 'c')
						clip_ocean = TRUE;
					else {
						if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[0])) {
							GMT_fill_syntax ('S');
							error++;
						}
						fill_ocean = TRUE;
					}
					break;
				case 'T':
					error += GMT_getrose (&argv[i][2], &mr);
					break;
				case 'W':
					if (argv[i][2] && GMT_getpen (&argv[i][2], &coast)) {
						GMT_pen_syntax ('W');
						error++;
					}
					draw_coast = TRUE;
					break;
				case 'Q':
					end_of_clip = TRUE;
					break;

				case '-':	/* Undocumented option to increase path-fix resolution */
					step = atof (&argv[i][2]);
					break;
#ifdef DEBUG
				case '+':
					debug_bin = atoi (&argv[i][2]);
					break;
#endif
				default:		/* Options not recognized */
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
	}
	
	if (argc == 1 || GMT_quick) {	/* Display usage */
		fprintf (stderr,"pscoast %s - Plot continents, shorelines, rivers, and borders on maps\n\n", GMT_VERSION);
		fprintf (stderr,"usage: pscoast -J<params> -R<west>/<east>/<south>/<north>\n\t[-A<min_area>[/<min_level>/<max_level>]]\n");
		fprintf (stderr, "\t[-B<tickinfo>] [-C[<fill>]] [-D<resolution>] [-Eaz/el] [-G[<fill>]] [-I<feature>[/<pen>]] [-K]\n");
		fprintf (stderr, "\t[-L[f][x]<lon0>/<lat0>[/<slon>]/<slat>/<length>[m|n|k][:label:<just>][+p<pen>][+f<fill>]]\n");
		fprintf (stderr, "\t[-M[<flag>]] [-N<feature>[/<pen>]] [-O] [-P] [-Q] [-S<fill>]\n");
		fprintf (stderr, "\t[-T[f|m][x]<lon0>/<lat0>/<size>[/<info>][:w,e,s,n:][+<gints>[/<mints>]]] [-U[dx/dy/][label]]\n");
		fprintf (stderr, "\t[-V] [-W[<pen>]] [-X<x_shift>] [-Y<y_shift>] [-bo[s][<n>]] [-c<ncopies>]\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		GMT_explain_option ('j');
		GMT_explain_option ('R');
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A features smaller than <min_area> (in km^2) or of levels (0-4) outside the min-max levels\n");
		fprintf (stderr, "\t   will be skipped [0/4 (4 means lake inside island inside lake)]\n");
		GMT_explain_option ('b');
		fprintf (stderr, "\t-C sets separate color for lakes. [Default is same as ocean (white)]. Set\n");
		fprintf (stderr, "\t   1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
		fprintf (stderr, "\t   2) p[or P]<dpi>/<pattern> for predefined patterns (1-90).\n");
		fprintf (stderr, "\t-D Choose one of the following resolutions:\n");
		fprintf (stderr, "\t   f - full resolution (may be very slow for large regions)\n");
		fprintf (stderr, "\t   h - high resolution (may be slow for large regions)\n");
		fprintf (stderr, "\t   i - intermediate resolution\n");
		fprintf (stderr, "\t   l - low resolution [Default]\n");
		fprintf (stderr, "\t   c - crude resolution, for busy plots that need crude continent outlines only\n");
		fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
		fprintf (stderr, "\t-G Paint or clip \"dry\" areas.  Append desired effect as:\n");
		fprintf (stderr, "\t   1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
		fprintf (stderr, "\t   2) p[or P]<dpi>/<pattern> for predefined patterns (1-90).\n");
		fprintf (stderr, "\t   3) c to issue clip paths for land areas.\n");
		fprintf (stderr, "\t-I draws rIvers.  Specify feature and optionally append pen [Default is solid line of unit thickness].\n");
		fprintf (stderr, "\t   Choose from the features below.  Repeat the -I option as many times as needed\n");
		fprintf (stderr, "\t      1 = Permanent major rivers\n");
		fprintf (stderr, "\t      2 = Additional major rivers\n");
		fprintf (stderr, "\t      3 = Additional rivers\n");
		fprintf (stderr, "\t      4 = Minor rivers\n");
		fprintf (stderr, "\t      5 = Intermittent rivers - major\n");
		fprintf (stderr, "\t      6 = Intermittent rivers - additional\n");
		fprintf (stderr, "\t      7 = Intermittent rivers - minor\n");
		fprintf (stderr, "\t      8 = Major canals\n");
		fprintf (stderr, "\t      9 = Minor canals\n");
		fprintf (stderr, "\t     10 = Irrigation canals\n");
		fprintf (stderr, "\t      a = All rivers and canals (1-10)\n");
		fprintf (stderr, "\t      r = All permanent rivers (1-4)\n");
		fprintf (stderr, "\t      i = All intermittent rivers (5-7)\n");
		fprintf (stderr, "\t      c = All canals (8-10)\n");
		GMT_explain_option ('K');
		fprintf (stderr, "\t-L draws a simple map scaLe centered on <lon0>/<lat0>.  Use -Lx to specify Cartesian\n");
		fprintf (stderr, "\t   coordinates instead.  Scale is calculated at latitude <slat>; optionally give longitude\n");
		fprintf (stderr, "\t   <slon> [Default is central longitude].  <length> is in km, or [nautical] miles if [n] m\n");
		fprintf (stderr, "\t   is appended.  -Lf draws a \"fancy\" scale [Default is plain]. By default, the label is set\n");
		fprintf (stderr, "\t   to the distance unit and placed on top [t].  Use the :label:<just> mechanism to specify another\n");
		fprintf (stderr, "\t   label (or - to keep the default) and placement (t,b,l,r, and u - to use the label as a unit.\n");
		fprintf (stderr, "\t   Append +p<pen> and/or +f<fill> to draw/paint a rectangle beneath the scale [no rectangle]\n");
		fprintf (stderr, "\t-M Dump a multisegment ascii (or binary, see -bo) file to standard output.  No plotting occurs\n");
		fprintf (stderr, "\t   Specify any combination of -W, -I, -N.  Optionally, append 1-char\n");
		fprintf (stderr, "\t   segment header <flag> [%c]\n", GMT_io.EOF_flag);
		fprintf (stderr, "\t-N draws boundaries.  Specify feature and optionally append pen [Default is solid line of unit thickness]\n");
		fprintf (stderr, "\t   Choose from the features below.  Repeat the -N option as many times as needed\n");
		fprintf (stderr, "\t     1 = National boundaries\n");
		fprintf (stderr, "\t     2 = State boundaries within the Americas\n");
		fprintf (stderr, "\t     3 = Marine boundaries\n");
		fprintf (stderr, "\t     a = All boundaries (1-3)\n");
		GMT_explain_option ('O');
		GMT_explain_option ('P');
		fprintf (stderr,"\t-Q means terminate previously set clip-paths\n");
		fprintf (stderr, "\t-S Paint or clip \"wet\" areas.  Append desired effect as:\n");
		fprintf (stderr, "\t   1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
		fprintf (stderr, "\t   2) p[or P]<dpi>/<pattern> for predefined patterns (1-90).\n");
		fprintf (stderr, "\t   3) c to issue clip paths for water areas.\n");
		fprintf (stderr, "\t-T draws a north-poinTing rose centered on <lon0>/<lat0>.  Use -Tx to specify Cartesian\n");
		fprintf (stderr, "\t   coordinates instead.  -Tf draws a \"fancy\" rose [Default is plain].  Give rose diameter\n");
		fprintf (stderr, "\t   <size> and optionally the west, east, south, north labels desired [W,E,S,N].\n");
		fprintf (stderr, "\t   For fancy rose, specify as <info> the kind you want: 1 draws E-W, N-S directions [Default],\n");
		fprintf (stderr, "\t   2 adds NW-SE and NE-SW, while 3 adds WNW-ESE, NNW-SSE, NNE-SSW, and ENE-WSW.\n");
		fprintf (stderr, "\t   For Magnetic compass rose, specify -Tm.  Use the optional <info> = <dec>/<dlabel> (where <dec> is\n");
		fprintf (stderr, "\t   the magnetic declination and <dlabel> is a label for the magnetic compass needle) to plot\n");
		fprintf (stderr, "\t   directions to both magnetic and geographic north [Default is just geographic].\n");
		fprintf (stderr, "\t   If the North label = \'*\' then a north star is plotted instead of the label.\n");
		fprintf (stderr, "\t   Append +<gints>/<mints> to override default annotation/tick interval(s) [10/5/1/30/5/1].\n");
		GMT_explain_option ('U');
		GMT_explain_option ('V');
		fprintf (stderr,"\t-W draw shorelines.  Append pen [Default is solid line of unit thickness].\n");
		GMT_explain_option ('X');
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('c');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}

	if (set_lake_fill && !(fill_land || fill_ocean || draw_coast)) {	/* Just lakes, fix -A */
		if (min_level < 2) min_level = 2;
	}
		
	/* Check that the options selected are mutually consistent */
	
	clipping = (clip_land || clip_ocean);
	if (!(end_of_clip || project_info.region_supplied)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (coast.width < 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  Pen thickness cannot be negative\n", GMT_program);
		error++;
	}
	if (!(fill_land || fill_ocean || set_lake_fill || draw_coast || draw_border || draw_river || clipping || end_of_clip)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -C, -G, -S, -I, -N, and -W\n", GMT_program);
		error++;
	}
	if ((fill_land + fill_ocean + set_lake_fill) > 0 && (clip_land + clip_ocean) > 0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must choose between clipping and painting\n", GMT_program);
		error++;
	}
	if (clip_land && clip_ocean) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must choose between clipping land OR water\n", GMT_program);
		error++;
	}
	if (dumping && (clipping + fill_land + fill_ocean + set_lake_fill)) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must choose between dumping and clipping/plotting\n", GMT_program);
		error++;
	}
	if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program);
		error++;
	}
	if (clipping && project_info.projection == AZ_EQDIST && fabs (project_info.w - project_info.e) == 360.0 && (project_info.n - project_info.s) == 180.0) {
		fprintf (stderr, "%s: -JE not implemented for global clipping - I quit\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (clipping && (draw_border || draw_river || draw_coast)) {
		fprintf (stderr, "%s: Cannot do clipping AND draw coastlines, rivers, or borders\n", GMT_program);
		exit (EXIT_FAILURE);
	}
	if (dumping && !(draw_border || draw_river || draw_coast)) {
		fprintf (stderr, "%s: Must specify at least one of -I, -N, and -W\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (draw_river) {
		for (k = 0; k < N_RLEVELS; k++) {
			if (!rlevels[k]) continue;
			if (rpen[k].width < 0.0) {
				fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Pen thickness cannot be negative\n", GMT_program);
				error++;
			}
			rlevels[n_rlevels] = k + 1;
			rpen[n_rlevels] = rpen[k];
			n_rlevels++;
		}
	}
	
	if (draw_border) {
		for (k = 0; k < N_BLEVELS; k++) {
			if (!blevels[k]) continue;
			if (bpen[k].width < 0.0) {
				fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Pen thickness cannot be negative\n", GMT_program);
				error++;
			}
			blevels[n_blevels] = k + 1;
			bpen[n_blevels] = bpen[k];
			n_blevels++;
		}
	}

	if (error) exit (EXIT_FAILURE);

	GMT_put_history (argc, argv);	/* Update .gmtcommands4 */

	need_coast_base = (fill_land || fill_ocean || set_lake_fill || draw_coast || clip_land || clip_ocean);
	clobber_background = (fill_land && fill_ocean);
	recursive = (((fill_land + (fill_ocean || set_lake_fill)) == 1) || clipping);
	paint_polygons = (fill_land || fill_ocean || set_lake_fill) || clipping;

	if (east > 360.0) {
		west -= 360.0;
		east -= 360.0;
	}
	
	if (end_of_clip && project_info.projection == -1) {	/* Fake area and linear projection */
		GMT_map_getproject ("x1d");
		east = north = 1.0;
	}
	else if (dumping) {
		GMT_map_getproject ("x1d");
	}

	GMT_map_setup (west, east, south, north);

	world_map_save = GMT_world_map;
	
	if (need_coast_base && 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]);
		need_coast_base = FALSE;
	}
	
	if (draw_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]);
		draw_border = FALSE;
	}
	
	if (draw_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]);
		draw_river = FALSE;
	}
	
	if (! (need_coast_base || draw_border || draw_river || end_of_clip)) {
		fprintf (stderr, "%s: No databases available - aborts\n", GMT_program);
		exit (EXIT_FAILURE);
	}

	if (dumping) {
		if (!GMT_io.binary[1]) fprintf (GMT_stdout, "# Data from the %s resolution GMT shoreline, borders, and rivers database\n", shore_resolution[base]);
	}
	else {
		if (end_of_clip)
			gmtdefs.page_orientation |= GMT_CLIP_OFF;	/* Signal that this program terminates clipping that initiated prior to this process */
		else if (clipping)
			gmtdefs.page_orientation |= GMT_CLIP_ON;	/* Signal that this program initiates clipping that wil outlive this process */
			
		ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
		gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
		gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, gmtdefs.encoding.name, GMT_epsinfo(argv[0]));
		GMT_echo_command (argc, argv);
		if (gmtdefs.unix_time) GMT_timestamp (argc, argv);
		if (!end_of_clip && base > 0) ps_command ("2 setlinejoin");	/* Avoid excessive mitering for cruder lines */
		 
        	if (end_of_clip) {  /* Just undo previous clip-path */
              		ps_clipoff ();

			GMT_map_basemap (); /* Basemap needed */

			ps_plotend (gmtdefs.last_page);

			if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program);

			GMT_end (argc, argv);
		}

		if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
	}

	if (fill_ocean && !set_lake_fill) fill[2] = fill[4] = fill[0];	/* Use ocean fill by default for lakes */
	
	for (i = 0; i < 5; i++) if (fill[i].use_pattern) {	/* Must initialize patterns */
		ps_imagefill_init (fill[i].pattern_no, fill[i].pattern, fill[i].inverse, fill[i].dpi, fill[i].colorize, fill[i].f_rgb, fill[i].b_rgb);
		fill_in_use = TRUE;
	}
	
	if (fill_in_use && !clobber_background) {	/* Force clobber when fill is used since our routine cannot deal with clipped fills */
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Pattern fill requires oceans to be painted first\n", GMT_program);
		clobber_background = TRUE;
		recursive = FALSE;
	}
		
	if (clipping && frame_info.plot) {	/* Want basemap drawn before clipping is activated */
		/* project_info.w = west_border;
		project_info.e = east_border; */
		GMT_map_basemap ();
	}

	if (project_info.projection == AZ_EQDIST && fabs (project_info.w - project_info.e) == 360.0 && (project_info.n - project_info.s) == 180.0) {
		possibly_donut_hell = TRUE;
		anti_lon = project_info.central_meridian + 180.0;
		if (anti_lon >= 360.0) anti_lon -= 360.0;
		anti_lat = -project_info.pole;
		anti_bin = (int)floor ((90.0 - anti_lat) / c.bsize) * c.bin_nx + (int)floor (anti_lon / c.bsize);
		GMT_geo_to_xy (anti_lon, anti_lat, &anti_x, &anti_y);
		GMT_geo_to_xy (project_info.central_meridian, project_info.pole, &x_0, &y_0);
		if (fill_land && gmtdefs.verbose) {
			fprintf (stderr, "%s: Warning: Fill continent option (-G) may not work for this projection.\n", GMT_program);
			fprintf (stderr, "If the antipole (%g/%g) is in the ocean then chances are good\n", anti_lon, anti_lat);
			fprintf (stderr, "Else: avoid projection center coordinates that are exact multiples of %g degrees\n", c.bsize);
		}
	}

	if (possibly_donut_hell && paint_polygons && !clobber_background) {	/* Force clobber when donuts may be called for now */
		if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: -JE requires oceans to be painted first\n", GMT_program);
		clobber_background = TRUE;
		recursive = FALSE;
	}
		
	if (clobber_background) {	/* Paint entire map as ocean first, then lay land on top */
		n = GMT_map_clip_path (&xtmp, &ytmp, &donut);
		if (donut) {
			ps_line (xtmp, ytmp, n, 1, TRUE, FALSE);
			ps_line (&xtmp[n], &ytmp[n], n, -1, TRUE, FALSE);
			GMT_fill (xtmp, ytmp, n, &fill[0], -1);
		}
		else 
			GMT_fill (xtmp, ytmp, n, &fill[0], FALSE);
		GMT_free ((void *)xtmp);
		GMT_free ((void *)ytmp);
		fill_ocean = FALSE; 
	}
	if (clipping) {
		start_direction = (clip_land) ? 1 : -1;
		stop_direction  = (clip_land) ? 1 : -1; 
		level_to_be_painted = (clip_land) ? 1 : 0;
	}
	else if (recursive) {
		start_direction = (fill_land) ? 1 : -1;
		stop_direction  = (fill_land) ? 1 : -1; 
		level_to_be_painted = (fill_land) ? 1 : 0;
	}
	else {
		start_direction = -1;
		stop_direction = 1;
		lowest_level = (fill_ocean) ? 0 : 1;
	}

	if (west < 0.0 && east > 0.0 && project_info.projection > 5) greenwich = TRUE;
	if ((360.0 - fabs (project_info.e - project_info.w) ) < c.bsize) {
		GMT_world_map = TRUE;
		if (project_info.projection == GNOMONIC) GMT_world_map = FALSE;
	}
	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;
	}
	
	/* west_border = project_info.w;	east_border = project_info.e; */
	if (project_info.w < 0.0 && project_info.e <= 0.0) {	/* Temporarily shift boundaries */
		project_info.w += 360.0;
		project_info.e += 360.0;
		if (project_info.central_meridian < 0.0) project_info.central_meridian += 360.0;
	}
	/* west_border = project_info.w;	east_border = project_info.e; */
	west_border = floor (project_info.w / c.bsize) * c.bsize;
	east_border = ceil (project_info.e / c.bsize) * c.bsize;
		
	if (!dumping && draw_coast) GMT_setpen (&coast);

	for (ind = 0; need_coast_base && ind < c.nb; ind++) {	/* Loop over necessary bins only */
		
		bin = c.bins[ind];
#ifdef DEBUG
		if (debug_bin >= 0 && bin != debug_bin) continue;
#endif
		GMT_get_shore_bin (ind, &c, min_area, min_level, max_level);
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Working on block # %5d\r", GMT_program, bin);
		if (!dumping) {
			sprintf (comment, "Bin # %d", bin);
			ps_comment (comment);
		}

		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);
		}

		bin_trouble = (anti_bin == bin) ? anti_bin : -1;
			
		if (possibly_donut_hell) {
			bin_x[0] = bin_x[3] = bin_x[4] = c.lon_sw;
			bin_x[1] = bin_x[2] = c.lon_sw + c.bsize;
			bin_y[0] = bin_y[1] = bin_y[4] = c.lat_sw;
			bin_y[2] = bin_y[3] = c.lat_sw + c.bsize;
			GMT_geo_to_xy (c.lon_sw + 0.5 * c.bsize, c.lat_sw + 0.5 * c.bsize, &x_c, &y_c);
			dist = hypot (x_c - x_0, y_c - y_0);
			donut_hell = GMT_non_zero_winding (anti_lon, anti_lat, bin_x, bin_y, 5);
			if (!donut_hell) donut_hell = (dist > 0.8 * project_info.r);
		}
		
		for (direction = start_direction; paint_polygons && direction <= stop_direction; direction += 2) {
		
			/* Assemble one or more segments into polygons */
		
	/*		if ((np = GMT_assemble_shore (&c, direction, TRUE, shift, edge, &p)) == 0) continue; */
			if ((np = GMT_assemble_shore (&c, direction, min_level, TRUE, shift, west_border, east_border, &p)) == 0) continue;

			/* Get clipped polygons in x,y inches that can be plotted */
			
			np_new = GMT_prep_polygons (&p, np, greenwich, donut_hell, step, bin_trouble);
		
			if (clipping) {
				for (k = level_to_be_painted; k < MAX_LEVEL - 1; k++) {
					recursive_clipping (-1, np_new, p, k, k, &n_pieces);
					did_clip = TRUE;
				}
				
				for (k = 0; k < np_new; k++) {	/* Do any remaining interior polygons */
					if (p[k].n == 0) continue;
						if (p[k].level%2 == level_to_be_painted || p[k].level > 2) {
							flag = (n_pieces == 0) ? 1 : 0;
							n_pieces++;
							ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag);
						}
				}
			}
			else if (!recursive) {	/* Simply points all polygons as is */
				for (k = 0; k < np_new; k++) {
					if (p[k].n == 0 || p[k].level < lowest_level) continue;
					level = p[k].level;
					if (donut_hell && GMT_non_zero_winding (anti_x, anti_y, p[k].lon, p[k].lat, p[k].n)) {	/* Antipode inside polygon, must do donut */
						n = GMT_map_clip_path (&xtmp, &ytmp, &donut);
						ps_line (xtmp, ytmp, n, 1, TRUE, FALSE);
						ps_line (p[k].lon, p[k].lat, p[k].n, -1, TRUE, FALSE);
						GMT_fill (xtmp, ytmp, n, &fill[level], -1);
						GMT_free ((void *)xtmp);
						GMT_free ((void *)ytmp);
					}
					else
						GMT_fill (p[k].lon, p[k].lat, p[k].n, &fill[level], FALSE);
				}
			}
			else {	/* Must avoid pointing anything but the polygons inside */
			
				for (k = level_to_be_painted; k < MAX_LEVEL - 1; k++) {
					recursive_painting (-1, np_new, p, k, k, fill);
				}
				
				for (k = 0; k < np_new; k++) {	/* Do any remaining interior polygons */
					if (p[k].n == 0) continue;
					if (p[k].level%2 == level_to_be_painted || p[k].level > 2)
						GMT_fill (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], FALSE);
				}
				
			}
			
			GMT_free_polygons (p, np_new);
			GMT_free ((void *)p);
		}
		
		if (draw_coast && c.ns) {	/* Draw or dump shorelines, no need to assemble polygons */
			/* if ((np = GMT_assemble_shore (&c, direction, FALSE, shift, edge, &p)) == 0) continue; */
			if ((np = GMT_assemble_shore (&c, direction, min_level, FALSE, shift, west_border, east_border, &p)) == 0) continue;
		
			for (i = 0; i < np; i++) {
				if (dumping) {
					sprintf (GMT_io.segment_header, "%c Shore Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level);
					GMT_write_segmentheader (GMT_stdout, 2);
					for (k = 0; k < p[i].n; k++) {
						out[0] = p[i].lon[k];
						out[1] = p[i].lat[k];
						GMT_output (GMT_stdout, 2, out);
					}
				}
				else {
					if (donut_hell) p[i].n = GMT_fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step);
					GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
					if (GMT_n_plot) GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot);
				}
			}
			
			GMT_free_polygons (p, np);
			GMT_free ((void *)p);
		}
		
		GMT_free_shore (&c);
	
	}
	if (need_coast_base) GMT_shore_cleanup (&c);

	if (clipping && did_clip) ps_clipon ((double *)NULL, (double *)NULL, 0, GMT_no_rgb, 2);	/* Dummy to end path */

	if (gmtdefs.verbose) fprintf (stderr, "\n");
	
	if (draw_river) {	/* Read rivers file and plot as lines */
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Adding Rivers...", GMT_program);
		if (!dumping) ps_comment ("Start of River segments");
		last_k = -1;

		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 (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);
			}
			
			if ((np = GMT_assemble_br (&r, shift, edge, &p)) == 0) continue;
			
			for (i = 0; i < np; i++) {
				if (dumping) {
					sprintf (GMT_io.segment_header, "%c River Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level);
					GMT_write_segmentheader (GMT_stdout, 2);
					for (k = 0; k < p[i].n; k++) {
						out[0] = p[i].lon[k];
						out[1] = p[i].lat[k];
						GMT_output (GMT_stdout, 2, out);
					}
				}
				else {
					GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
					if (!GMT_n_plot) continue;

					k = p[i].level - 1;
					if (k != last_k) {
						GMT_setpen (&rpen[k]);
						last_k = k;
					}
					GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot);
				}
			}
			
			/* Free up memory */
		
			GMT_free_br (&r);
			for (k = 0; k < np; k++) {
				GMT_free ((void *)p[k].lon);
				GMT_free ((void *)p[k].lat);
			}
			GMT_free ((void *)p);
		}
		GMT_br_cleanup (&r);
		
		if (gmtdefs.verbose) fprintf (stderr, "\n");
	}

	if (draw_border) {	/* Read borders file and plot as lines */
		
		if (gmtdefs.verbose) fprintf (stderr, "%s: Adding Borders...", GMT_program);
		if (!dumping) ps_comment ("Start of Border segments");
		
		/* Must resample borders because some points may be too far apart and look like 'jumps' */
		
		step = MAX (fabs(project_info.w - project_info.e), fabs (project_info.n - project_info.s)) / 4.0;

		last_k = -1;

		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 (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);
			}

			if ((np = GMT_assemble_br (&b, shift, edge, &p)) == 0) continue;

			for (i = 0; i < np; i++) {
				if (dumping) {
					sprintf (GMT_io.segment_header, "%c Border Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level);
					GMT_write_segmentheader (GMT_stdout, 2);
					for (k = 0; k < p[i].n; k++) {
						out[0] = p[i].lon[k];
						out[1] = p[i].lat[k];
						GMT_output (GMT_stdout, 2, out);
					}
				}
				else {
					p[i].n = GMT_fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step);
					GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
					if (!GMT_n_plot) continue;

					k = p[i].level - 1;
					if (k != last_k) {
						GMT_setpen (&bpen[k]);
						last_k = k;
					}
					GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot);
				}
			}
			
			/* Free up memory */
		
			GMT_free_br (&b);
			for (k = 0; k < np; k++) {
				GMT_free ((void *)p[k].lon);
				GMT_free ((void *)p[k].lat);
			}
			GMT_free ((void *)p);
		}
		GMT_br_cleanup (&b);
		
		if (gmtdefs.verbose) fprintf (stderr, "\n");
	}
	
	if (!dumping) {
		if (base > 0) ps_command ("0 setlinejoin");	/* Reset mitering */

		if (frame_info.plot) {
			/* project_info.w = west_border;
			project_info.e = east_border; */
			GMT_world_map = world_map_save;
			GMT_map_basemap ();
		}
	
		if (ms.plot) GMT_draw_map_scale (&ms);
		if (mr.plot) GMT_draw_map_rose (&mr);

		if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
		ps_plotend (gmtdefs.last_page);
	}

	GMT_end (argc, argv);
}

void recursive_painting (int k0, int np, struct POL p[], int level, int start_level, struct GMT_FILL fill[]) {
	int k, flag;
	BOOLEAN go;
	char text[64];
	if (level > MAX_LEVEL) return;
	flag = (level == start_level) ? 1 : 0;
	for (k = k0 + 1; k < np; k++) {
		if (p[k].n == 0 || p[k].level < level) continue;
		go = (p[k].level == level) ? ((k0 == -1) ? TRUE : GMT_non_zero_winding (p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) : FALSE;
		if (go) {
			sprintf (text, "Polygon %d", k);
			ps_comment (text);
			ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag);
			recursive_painting (k, np, p, level+1, start_level, fill);
			p[k].n = 0;
		}
		
		if (p[k].level == start_level) {	/* Done nesting, time to paint */
			ps_clipon (p[k].lon, p[k].lat, 0, GMT_no_rgb, 4);	/* Dummy to end path */
			GMT_fill (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], -1);
		}
	}
}

void recursive_clipping (int k0, int np, struct POL *p, int level, int start_level, int *n_pieces)
{
	int k, flag;
	BOOLEAN go;
	char text[64];
	if (level > MAX_LEVEL) return;
	for (k = k0 + 1; k < np; k++) {
		if (p[k].n == 0 || p[k].level < level) continue;
		go = (p[k].level == level) ? ((k0 == -1) ? TRUE : GMT_non_zero_winding (p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) : FALSE;
		if (go) {
			sprintf (text, "Polygon %d", k);
			ps_comment (text);
			flag = (*n_pieces == 0) ? 1 : 0;
			(*n_pieces)++;
			ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag);
			recursive_clipping (k, np, p, level+1, start_level, n_pieces);
			p[k].n = 0;
		}
	}
}
