/*--------------------------------------------------------------------
 *	$Id: grdlandmask.c,v 1.3.4.2 2002/02/27 17:41:10 pwessel Exp $
 *
 *	Copyright (c) 1991-2002 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
 *--------------------------------------------------------------------*/
/*
 * grdlandmask defines a grid based on region and xinc/yinc values,
 * reads a shoreline data base, and sets the grid nodes inside, on the
 * boundary, and outside of the polygons to the user-defined values
 * <in>, <on>, and <out>.  These may be any number, including NaN.
 *
 * Author:	P. Wessel
 * Date:	23-Sep-1994
 * Version:	3.0
 * Modified:	24-JUN-1998, for GMT 3.1
 *		18-AUG-1999, for GMT 3.3.2
 *		13-JUL-2000, for GMT 3.3.5
 * Version:	3.4.1
 */
 
#include "gmt.h"

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

struct GMT_SHORE c;

float *data;

main (int argc, char **argv)
{

	int	i, j, k, ij, bin, ind, nm, np, side, i_min, i_max, j_min, j_max, nx1, ny1, np_new;
	int	one_or_zero, base = 3, direction, inside = 1, max_level = MAX_LEVEL, min_level = 0;
	
	BOOLEAN	error = FALSE, pixel = FALSE, dry_wet_only = FALSE, greenwich = FALSE;
	BOOLEAN temp_shift = FALSE, used_polygons;
	
	double	*x, *y, xmin, xmax, ymin, ymax, west_border, east_border, i_dx_inch, i_dy_inch;
	double xinc2, yinc2, min_area = 0.0, i_dx, i_dy, edge = 0.0, del_off, dummy;
	float out_edge_in[5];
	
	char *maskfile = CNULL, res = 'l', line[128], *ptr;
	
	struct GRD_HEADER header;
	struct POL *p;
	
	argc = GMT_begin (argc, argv);
		
	GMT_grd_init (&header, argc, argv, FALSE);
	
	memset ((void *)out_edge_in, 0, (size_t)(5 * sizeof (float)));	/* Default "wet" value = 0 */
	out_edge_in[1] = out_edge_in[3] = 1.0;			/* Default "dry" value = 1 */

	/* Check command line arguments */
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
			
				/* Common parameters */
			
				case 'R':
				case 'V':
				case ':':
				case '\0':
					error += GMT_get_common_args (argv[i], &header.x_min, &header.x_max, &header.y_min, &header.y_max);
					break;
				
				/* Supplemental parameters */
				
				case 'A':
					j = sscanf (&argv[i][2], "%lf/%d/%d", &min_area, &min_level, &max_level);
					if (j == 1) min_level = 0, max_level = MAX_LEVEL;
					break;
				case 'D':
					res = argv[i][2];
					base = GMT_set_resolution (&res, 'D');
					break;
				case 'N':
					strcpy (line, &argv[i][2]);
					if (line[strlen(line)-1] == 'o') { /* Edge is considered outside */
						inside = 2;
						line[strlen(line)-1] = 0;
					}
					ptr = strtok (line, "/");
					j = 0;
					while (j < 5 && ptr) {
						out_edge_in[j] = (ptr[0] == 'N' || ptr[0] == 'n') ? GMT_f_NaN : (float)atof (ptr);
						ptr = strtok (CNULL, "/");
						j++;
					}
					if (!(j == 2 || j == 5)) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Specify 2 or 5 arguments\n", GMT_program);
						exit (EXIT_FAILURE);
					}
					dry_wet_only = (j == 2);
					break;
				case 'F':
					pixel = TRUE;
					break;
				case 'G':
					maskfile = &argv[i][2];
					break;
				case 'I':
					GMT_getinc (&argv[i][2], &header.x_inc, &header.y_inc);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
	}

	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "grdlandmask %s - Create \"wet-dry\" mask grdfile from shoreline data base\n\n", GMT_VERSION);
		fprintf (stderr, "usage: grdlandmask -G<mask_grd_file> -I<xinc[m|c]>[/<yinc>[m|c]] -R<west/east/south/north>\n");
		fprintf (stderr, "\t[-A<min_area>[/<min_level>/<max_level>]] [-D<resolution>] [-F] [-N<maskvalues>[o]] [-V]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\t-G Specify file name for output mask grd file.\n");
		fprintf (stderr, "\t-I sets Increment of the grid; enter xinc, optionally xinc/yinc.\n");
		fprintf (stderr, "\t   Default is yinc = xinc.  Append an m [or c] to indicate minutes [or seconds],\n");
		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] (see pscoast for details)]\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 tasks that need crude continent outlines only\n");
		fprintf (stderr, "\t-F Force pixel registration for output grid [Default is gridline orientation]\n");
		fprintf (stderr, "\t-N gives values to use if a node is outside or inside a feature.\n");
		fprintf (stderr, "\t   Append o to let feature boundary be considered outside [Default is inside].\n");
		fprintf (stderr, "\t   Specify this information using 1 of 2 formats:\n");
		fprintf (stderr, "\t   -N<wet>/<dry>.\n");
		fprintf (stderr, "\t   -N<ocean>/<land>/<lake>/<island>/<pond>.\n");
		fprintf (stderr, "\t   NaN is a valid entry.  Default values are 0/1/0/1/0 (i.e., 0/1)\n");
		GMT_explain_option ('V');
		exit (EXIT_FAILURE);
	}

	if (!project_info.region_supplied) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", GMT_program);
		error++;
	}
	if (header.x_inc <= 0.0 || header.y_inc <= 0.0) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", GMT_program);
		error = TRUE;
	}
	if (!maskfile) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", GMT_program);
		error = TRUE;
	}
	
	if (error) exit (EXIT_FAILURE);

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

	if (header.x_min < 0.0 && header.x_max < 0.0) {	/* Shift longitudes */
		temp_shift = TRUE;
		header.x_min += 360.0;
		header.x_max += 360.0;
	}

	GMT_grd_RI_verify (&header, 1);
	
	if (dry_wet_only) {
		out_edge_in[3] = out_edge_in[1];
		out_edge_in[2] = out_edge_in[4] = out_edge_in[0];
	}

	if (GMT_init_shore (res, &c, header.x_min, header.x_max, header.y_min, header.y_max))  {
		fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", GMT_program, shore_resolution[base]);
		exit (EXIT_FAILURE);
	}
	
	sprintf (line, "%s\n", gmtdefs.d_format);
	if (gmtdefs.verbose && dry_wet_only) {
		fprintf (stderr, "%s: Nodes in water will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[0]);
		fprintf (stderr, "%s: Nodes on land will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[1]);
	}
	else if (gmtdefs.verbose && !dry_wet_only) {
		fprintf (stderr, "%s: Nodes in the oceans will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[0])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[0]);
		fprintf (stderr, "%s: Nodes on land will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[1])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[1]);
		fprintf (stderr, "%s: Nodes in lakes will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[2])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[2]);
		fprintf (stderr, "%s: Nodes in islands will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[3])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[3]);
		fprintf (stderr, "%s: Nodes in ponds will be set to ", GMT_program);
		(GMT_is_fnan (out_edge_in[4])) ? fprintf (stderr, "NaN\n") : fprintf (stderr, line, out_edge_in[4]);
	}
	
	i_dx = 1.0 / header.x_inc;
	i_dy = 1.0 / header.y_inc;
	one_or_zero = (pixel) ? 0 : 1;
	del_off = (pixel) ? 0.5 : 0.0;
	header.nx = irint ((header.x_max - header.x_min) * i_dx) + one_or_zero;
	header.ny = irint ((header.y_max - header.y_min) * i_dy) + one_or_zero;
	header.node_offset = pixel;
	xinc2 = (header.node_offset) ? 0.5 * header.x_inc : 0.0;
	yinc2 = (header.node_offset) ? 0.5 * header.y_inc : 0.0;
	nm = header.nx * header.ny;
	data = (float *) GMT_memory (VNULL, (size_t)nm, sizeof(float), GMT_program);
	/* All data nodes are thus initialized to 0 */
	x = (double *) GMT_memory (VNULL, (size_t)header.nx, sizeof(double), GMT_program);
	y = (double *) GMT_memory (VNULL, (size_t)header.ny, sizeof(double), GMT_program);
	
	nx1 = header.nx - 1;	ny1 = header.ny - 1;

	if (header.x_min < 0.0 && header.x_max > 0.0) {	/* Must shift longitudes */
		greenwich = TRUE;
		edge = header.x_min;
	}
	
	GMT_map_getproject ("x1d");	/* Fake linear projection */
	GMT_map_setup (header.x_min, header.x_max, header.y_min, header.y_max);
	
	/* Fill out gridnode coordinates and apply the implicit linear projection */
	
	for (i = 0; i < header.nx; i++) GMT_geo_to_xy (header.x_min + i * header.x_inc + xinc2, 0.0, &x[i], &dummy);
	for (j = 0; j < header.ny; j++) GMT_geo_to_xy (0.0, header.y_max - j * header.y_inc - yinc2, &dummy, &y[j]);
	i_dx_inch = 1.0 / fabs (x[1] - x[0]);
	i_dy_inch = 1.0 / fabs (y[1] - y[0]);
	
        west_border = floor (project_info.w / c.bsize) * c.bsize;
        east_border = ceil (project_info.e / c.bsize) * c.bsize;	
	for (ind = 0; ind < c.nb; ind++) {	/* Loop over necessary bins only */
		
		bin = c.bins[ind];
		if (gmtdefs.verbose) fprintf (stderr, "%s: Working on block # %5d\r", GMT_program, bin);

		GMT_get_shore_bin (ind, &c, min_area, min_level, max_level);
		
		/* Use polygons, if any.  Go in both directions to cover both land and sea */
		
		used_polygons = FALSE;
		
		for (direction = -1; c.ns > 0 && direction < 2; direction += 2) {
		
			/* Assemble one or more segments into polygons */
		
			np = GMT_assemble_shore (&c, direction, min_level, TRUE, greenwich, west_border, east_border, &p);

			/* Get clipped polygons in x,y inches that can be processed */
			
			np_new = GMT_prep_polygons (&p, np, greenwich, FALSE, 0.0, -1);

			for (k = 0; k < np_new; k++) {
		
				if (p[k].n == 0) continue;

				used_polygons = TRUE;	/* At least som points made it to here */
				
				/* Find min/max of polygon in inches */
			
				xmin = xmax = p[k].lon[0];
				ymin = ymax = p[k].lat[0];
				for (i = 1; i < p[k].n; i++) {
					if (p[k].lon[i] < xmin) xmin = p[k].lon[i];
					if (p[k].lon[i] > xmax) xmax = p[k].lon[i];
					if (p[k].lat[i] < ymin) ymin = p[k].lat[i];
					if (p[k].lat[i] > ymax) ymax = p[k].lat[i];
				}
				i_min = (int)MAX (0, ceil (xmin * i_dx_inch - del_off - GMT_CONV_LIMIT));
				if (i_min > nx1) i_min = 0;
				i_max = (int)MIN (nx1, floor (xmax * i_dx_inch - del_off + GMT_CONV_LIMIT));
				if (i_max <= 0) i_max = nx1;				
				j_min = (int)MAX (0, ceil ((project_info.ymax - ymax) * i_dy_inch - del_off - GMT_CONV_LIMIT));
				j_max = (int)MIN (ny1, floor ((project_info.ymax - ymin) * i_dy_inch - del_off + GMT_CONV_LIMIT));
			
				for (j = j_min; j <= j_max; j++) {
					for (i = i_min; i <= i_max; i++) {
					
						if ((side = GMT_non_zero_winding (x[i], y[j], p[k].lon, p[k].lat, p[k].n)) < inside) continue;	/* Outside */
				
						/* Here, point is inside, we must assign value */
				
						ij = j * header.nx + i;
						if (p[k].level > data[ij]) data[ij] = (float)p[k].level;
					}
				}
			}
		
			GMT_free_polygons (p, np_new);
			GMT_free ((void *)p);
		}
		
		if (!used_polygons) {	/* Loack of polygons or clipping etc resulted in no polygons after all, must deal with background */
		
			k = 10;	/* Initialize to outside range of levels (4 is highest) */
			/* Visit each of the 4 nodes, test if it is inside -R, and if so update lowest level found so far */
			
			if (!GMT_map_outside (c.lon_sw, c.lat_sw)) k = MIN (k, c.node_level[0]);			/* SW */
			if (!GMT_map_outside (c.lon_sw + c.bsize, c.lat_sw)) k = MIN (k, c.node_level[1]);		/* SE */
			if (!GMT_map_outside (c.lon_sw + c.bsize, c.lat_sw - c.bsize)) k = MIN (k, c.node_level[2]);	/* NE */
			if (!GMT_map_outside (c.lon_sw, c.lat_sw - c.bsize)) k = MIN (k, c.node_level[3]);		/* NW */
			
			/* If k is still 10 we must assume this patch should have the min level of the bin */
			
			if (k == 10) k = MIN (MIN (c.node_level[0], c.node_level[1]) , MIN (c.node_level[2], c.node_level[3]));
			
			/* Determine nodes to initialize */
			
			i_min = (int)MAX (0, ceil (fmod (c.lon_sw - header.x_min + 360.0, 360.0) * i_dx - del_off));
			if (i_min > nx1) i_min = 0;
			i_max = (int)MIN (nx1, floor (fmod (c.lon_sw + c.bsize - header.x_min + 360.0, 360.0) * i_dx - del_off));
			if (i_max < 0) i_max = nx1;
			j_min = (int)MAX (0, ceil ((header.y_max - c.lat_sw - c.bsize) * i_dy - del_off));
			j_max = (int)MIN (ny1, floor ((header.y_max - c.lat_sw) * i_dy - del_off));
			for (j = j_min; j <= j_max; j++) for (i = i_min, ij = j * header.nx + i_min; i <= i_max; i++, ij++)
				data[ij] = (float)k;
		}

		GMT_free_shore (&c);
	}
	
	GMT_shore_cleanup (&c);
	if (gmtdefs.verbose) fprintf (stderr, "\n");
	
	for (ij = 0; ij < nm; ij++) {
		k = irint (data[ij]);
		data[ij] = out_edge_in[k];
	}
	
	if (temp_shift) {
		header.x_min -= 360.0;
		header.x_max -= 360.0;
	}

	if (GMT_write_grd (maskfile, &header, data, 0.0, 0.0, 0.0, 0.0, GMT_pad, FALSE)) {
		fprintf (stderr, "%s: Error writing file %s\n", GMT_program, maskfile);
		exit (EXIT_FAILURE);
	}

	if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program);
	
	GMT_free ((void *)data);
	GMT_free ((void *)x);
	GMT_free ((void *)y);
	
	GMT_end (argc, argv);
}