/*--------------------------------------------------------------------
 *	$Id: gmtselect.c,v 1.25 2004/06/25 22:39:01 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
 *--------------------------------------------------------------------*/
/*
 * gmtselect is a general-purpose spatial filter.  Data pass or fail based
 * on one or more conditions.  5 conditions may be set:
 *
 *	1. Only data inside a rectangular area may pass
 *	2. Only data within a certain distance from given points may pass
 *	3. Only data within a certain distance from given lines may pass
 *	4. Only data within given polygons may pass
 *	5. Only data within the coastline may pass
 *
 * Distances are calculated in the users units using Euclidian geometry
 * unless a map projection and region (-R -J) are used.  Then, distances
 * are calculated using spherical geometry and converted to km, and any
 * distances given in options or via headers are assumed to be in km.
 *
 * Any one of these conditions may be negated for the opposite result
 * Both binary and ASCII data files are accommodated
 *
 * Author:	Paul Wessel
 * Date:	25-AUG-1998
 * Version:	3.1
 *		3.2, 15-MAR-1999
 *		3.3, 13-APR-1999.	Added quick check for inside a square before checking
 *					for inside a circle (suggested by Joaquim Luis)
 *					Also added option -Cf for faster, approximate distances
 *		3.3.2 8-SEPT-1999.	Sorting -C points in x to speed up search as
 *				 	suggested by Joaquim Luis.
 *		3.3.4 17-FEB-2000.	Polygons containing either S or N pole will work correctly
 *		3.3.5 10-JUL-2000	Added plain -L for periodicity
 * 		3.4
 * Version:	4
 */
 
#include "gmt.h"

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

struct GMT_SHORE c;

struct DATA {
	int      i;
	double   x;
	double   y;
	double   d;
} *data;

int compare_x(const void *point_1, const void *point_2);

main (int argc, char **argv)
{
	int i, j, k, fno, n_read = 0, n_pass = 0, n_files = 0, n_args, n_points, n_lines, n_pols;
	int n_fields, n_expected_fields, out_edge_in[5], ind, bin, last_bin = -1;
	int min_level = 0, max_level = MAX_LEVEL, np[2], base = 3, wd[2], id, this_node;
	int side, is_inside = 1, row, col, n_output = -1;

	
	BOOLEAN error = FALSE, limit_line_distance = FALSE, limit_polygon = FALSE, limit_point_distance = FALSE, limit_point_distance_cart = FALSE;
	BOOLEAN want_inside[6] = {TRUE,TRUE,TRUE,TRUE,TRUE,TRUE}, nofile = TRUE, done = FALSE, first = TRUE, inside, greenwich = FALSE;
	BOOLEAN flat_earth = FALSE, do_project = FALSE, binary = FALSE, limit_coast = FALSE, dry_wet_only, long_verbose = FALSE, limit_z = FALSE;
	
	double west = 0.0, east = 0.0, south = 0.0, north = 0.0, xx, yy, min_area = 0.0;
	double *xp, *yp, *dp, pt_dist = 0.0, ln_dist = 0.0, *in, west_border, east_border;
	double xmin, xmax, ymin, ymax, lon, r, sa, ca, zmin = -DBL_MAX, zmax = DBL_MAX;
	
	char buffer[BUFSIZ], *pt_file, *ln_file, *pol_file, *ptr, res = 'l', za[16], zb[16];
	
	FILE *fp = NULL;
	
	struct GMT_LINES *pol, *line;
	struct POL *p[2];

	PFI near_a_line;

	pt_file = ln_file = pol_file = CNULL;
	memset ((void *)out_edge_in, 0, (size_t)(5 * sizeof (int)));	/* Default for "wet" areas = 0 (outside) */
	out_edge_in[1] = out_edge_in[3] = 1;			/* Default for "dry" areas = 1 (inside) */

	argc = GMT_begin (argc, argv);
	
	for (i = 1; i < argc; i++) {
		if (argv[i][0] == '-') {
			switch (argv[i][1]) {
		
				/* Common parameters */
			
				case 'V':
					if (argv[i][2] == 'l') long_verbose = TRUE;
				case 'H':
				case 'J':
				case 'R':
				case ':':
				case 'b':
				case 'f':
				case '\0':
					error += GMT_get_common_args (argv[i], &west, &east, &south, &north);
					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 'C':
					limit_point_distance = TRUE;
					k = 2;
					if (argv[i][2] == 'f') flat_earth = TRUE, k = 3;
					for (j = k; argv[i][j] && argv[i][j] != '/'; j++);
					if (!argv[i][j]) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -C:  Expects -C[f]dist/file\n", GMT_program);
						error++;
					}
					else {
						pt_file = &argv[i][j+1];
						pt_dist = atof (&argv[i][k]);
					}
					break;
				case 'D':
					res = argv[i][2];
					base = GMT_set_resolution (&res, 'D');
					break;
				case 'L':
					if (argv[i][2]) {	/* Set line options */
						limit_line_distance = TRUE;
						for (j = 2; argv[i][j] && argv[i][j] != '/'; j++);
							if (!argv[i][j]) {
							fprintf (stderr, "%s: GMT SYNTAX ERROR -L:  Expects -Ldist/file\n", GMT_program);
							error++;
						}
						else {
							ln_file = &argv[i][j+1];
							ln_dist = atof (&argv[i][2]);
						}
					}
					else {	/* Obsolete flag to process geographic data (use -f instead)  */
						GMT_io.in_col_type[0] =  GMT_io.out_col_type[0] = GMT_IS_LON;
						GMT_io.in_col_type[1] =  GMT_io.out_col_type[1] = GMT_IS_LAT;
					}
					break;
				case 'F':
					limit_polygon = TRUE;
					pol_file = &argv[i][2];
					break;
				case 'I':
					for (j = 2; argv[i][j]; j++) {
						switch (argv[i][j]) {
							case 'r':
								want_inside[0] = FALSE;
								break;
							case 'c':
								want_inside[1] = FALSE;
								break;
							case 'l':
								want_inside[2] = FALSE;
								break;
							case 'f':
								want_inside[3] = FALSE;
								break;
							case 's':
								want_inside[4] = FALSE;
								break;
							case 'z':
								want_inside[5] = FALSE;
								break;
							default:
								fprintf (stderr, "%s: GMT SYNTAX ERROR -I:  Expects -Icflrsz\n", GMT_program);
								error++;
								break;
						}
					}
					break;
				case 'M':               /* Multiple line segments input */
					GMT_multisegment (&argv[i][2]);
					break;
				case 'N':
					limit_coast = TRUE;
					strcpy (buffer, &argv[i][2]);
					if (buffer[strlen(buffer)-1] == 'o') { /* Edge is considered outside */
						is_inside = 2;
						buffer[strlen(buffer)-1] = 0;
					}
					ptr = strtok (buffer, "/");
					j = 0;
					while (j < 5 && ptr) {
						switch (ptr[0]) {
							case 's':	/* Skip points in this level */
								out_edge_in[j] = 0;
								break;
							case 'k':
								out_edge_in[j] = 1;
								break;
							default:
								fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Bad modifier (use s or k)\n", GMT_program);
								error++;
						}
						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 'Z':
					limit_z = TRUE;
					j = sscanf (&argv[i][2], "%[^/]/%s", za, zb);
					if (j != 2) {
						fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  Specify z_min and z_max\n", GMT_program);
						exit (EXIT_FAILURE);
					}
					if (!(za[0] == '-' && za[1] == '\0')) zmin = atof (za);
					if (!(zb[0] == '-' && zb[1] == '\0')) zmax = atof (zb);
					break;
				default:
					error = TRUE;
					GMT_default_error (argv[i][1]);
					break;
			}
		}
		else
			n_files++;
	}
	
	if (argc == 1 || GMT_quick) {
		fprintf (stderr, "gmtselect %s - Select data subsets based on multiple spatial criteria\n\n", GMT_VERSION);
		fprintf (stderr, "usage: gmtselect <infiles> [-A<min_area>[/<min_level>/<max_level>]]\n\t[-C[f]<dist>/<ptfile>] [-D<resolution>] [-F<polygon>]\n");
		fprintf (stderr, "\t[-J<parameters>] [-L[<dist>/<lfile>]] [-H[<nrec>]] [-I[cflrsz] [-M[<flag>]] [-N<maskvalues>[o]]\n");
		fprintf (stderr, "\t[-R<west/east/south/north>] [-V[l]]\n\t[-:] [-Z<min>/<max>] [-bi[s][<n>]] [-bo[s][<n>]] [-f[i|o]<colinfo>]\n\n");
		
		if (GMT_quick) exit (EXIT_FAILURE);
		
		fprintf (stderr, "\tinfiles (in ASCII or binary) have 2 or more columns with (x,y) or (y,x) in first columns.\n");
		fprintf (stderr, "\t  If no file(s) is given, standard input is read.\n");
		fprintf (stderr, "\n\tOPTIONS:\n");
		fprintf (stderr, "\t-A coastline 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)].  Ignored unless -N is set\n");
		fprintf (stderr, "\t-C pass locations that are within <dist> of any point in ASCII <ptfile>\n");
		fprintf (stderr, "\t   Give 0 distance if 3rd column of ptfile has individual distances.\n");
		fprintf (stderr, "\t   Distances are in users units [or km if -R -J are used].\n");
		fprintf (stderr, "\t   Prepend f to use approximate (flat earth) rather than exact great circle distances.\n");
		fprintf (stderr, "\t-D Choose one of the following resolutions: (Ignored unless -N is set)\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-L Pass locations that are within <dist> of any line in ASCII <linefile>\n");
		fprintf (stderr, "\t   Give 0 distance if 2nd column of segment headers have individual distances.\n");
		fprintf (stderr, "\t   Distances are in users units [or km if -R -J are used].\n");
		fprintf (stderr, "\t-F pass locations that are inside the polygons in the ASCII <polygon> file\n");
		GMT_explain_option ('H');
		fprintf (stderr, "\t-I Used to reverse the tests, i.e. pass locations outside the region\n");
		fprintf (stderr, "\t   Supply a combination of cflrz where each flag means:\n");
		fprintf (stderr, "\t   c will pass locations beyond the minimum distance to the points in -C\n");
		fprintf (stderr, "\t   f will pass locations outside the polygons in -F\n");
		fprintf (stderr, "\t   l will pass locations beyond the minimum distance to the lines in -L\n");
		fprintf (stderr, "\t   r will pass locations outside the region given in -R [and -J]\n");
		fprintf (stderr, "\t   s will pass locations that otherwise would be skipped in -N\n");
		fprintf (stderr, "\t   z will pass locations outside the range given in -Z\n");
		GMT_explain_option ('J');
		GMT_explain_option ('M');
		fprintf (stderr, "\t-N set if a point outside or inside a geographic feature should be s(kipped) or k(ept).\n");
		fprintf (stderr, "\t   Append o to let feature boundary be considered outside [Default is inside].\n");
		fprintf (stderr, "\t   Specify this information with s or k 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   k means keep and s means skip [Default is s/k/s/k/s (i.e., s/k)]\n");
		GMT_explain_option ('R');
		GMT_explain_option ('V');
		fprintf (stderr, "\t   Append l for long verbose, reporting every 1000 points.\n");
		fprintf (stderr, "\t   -Z assumes the 3rd data column contains z-values and we want to keep records with\n");
		fprintf (stderr, "\t      <min> <= z <= <max>.  Use - for <min> or <max> if there is no lower/upper limit.\n");
		GMT_explain_option (':');
		GMT_explain_option ('i');
		GMT_explain_option ('n');
		fprintf (stderr, "\t   Default is 2 input columns\n");
		GMT_explain_option ('o');
		GMT_explain_option ('n');
		GMT_explain_option ('f');
		GMT_explain_option ('.');
		exit (EXIT_FAILURE);
	}
	
	if (limit_z && zmax <= zmin) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR:  -Z must have zmax > zmin!\n", GMT_program);
		error++;
	}
	if (GMT_io.binary[GMT_IN] && gmtdefs.io_header[GMT_IN]) {
		fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", GMT_program);
		error++;
	}
        if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] == 0) GMT_io.ncol[GMT_IN] = 2 + limit_z;
        if (GMT_io.binary[GMT_IN] && GMT_io.ncol[GMT_IN] < (2 + limit_z)) {
                fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data (-bi) must have at least %d columns\n", GMT_program, 2 + limit_z);
		error++;
	}
	if (error) exit (EXIT_FAILURE);

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

	if (!MAPPING && limit_point_distance) limit_point_distance_cart = TRUE;

	if (GMT_io.binary[GMT_IN] && gmtdefs.verbose) {
		char *type[2] = {"double", "single"};
		fprintf (stderr, "%s: Expects %d-column %s-precision binary data\n", GMT_program, GMT_io.ncol[GMT_IN], type[GMT_io.single_precision[GMT_IN]]);
	}

#ifdef SET_IO_MODE
	GMT_setmode (GMT_OUT);
#endif

	if (!project_info.region_supplied && (limit_coast || (GMT_io.in_col_type[0] == GMT_IS_LON && GMT_io.in_col_type[1] == GMT_IS_LAT))) {	/* If we use coastline data or used -fg but didnt give -R we imply -Rg */
		project_info.region_supplied = TRUE;
		west = 0.0;	east = 360.0;	south = -90.0;	north = +90.0;
		GMT_io.in_col_type[0] = GMT_IS_LON;
		GMT_io.in_col_type[1] = GMT_IS_LAT;
	}
	if (project_info.region_supplied) {	/* -R was set directly or indirectly; hence must set -J if not supplied */
		if (project_info.projection < 0) {	/* -J not specifed, set one implicitly */
			if (GMT_io.in_col_type[0] == GMT_IS_LON) {	/* We know we have geographic data (either via -f or -N) */
				project_info.degree[0] = project_info.degree[1] = TRUE;
			}
			/* Supply dummy linear proj */
			project_info.projection = project_info.xyz_projection[0] = project_info.xyz_projection[1] = LINEAR;
			project_info.pars[0] = project_info.pars[1] = 1.0;
		}
		do_project = TRUE;
		if (MAPPING) {
			if (west < 0.0 && east < 0.0) {
				west += 360.0;
				east += 360.0;
			}
			greenwich = (west < 0.0 && east > 0.0);
		}
		GMT_map_setup (west, east, south, north);
	}

	if (limit_coast) {
		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, west, east, south, north))  {
			fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", GMT_program, shore_resolution[base]);
			exit (EXIT_FAILURE);
		}
		west_border = floor (project_info.w / c.bsize) * c.bsize;
		east_border = ceil (project_info.e / c.bsize) * c.bsize;
		wd[0] = 1;	wd[1] = -1;
		np[0] = np[1] = 0;
	}

	binary = (GMT_io.binary[GMT_IN] || GMT_io.binary[GMT_OUT]);

	if (GMT_io.in_col_type[0] & GMT_IS_GEO) {	/* Geographic data */
		GMT_distance_func = (PFD) ((MAPPING) ? ((flat_earth) ? GMT_flatearth_dist : GMT_km_dist) : GMT_km_dist);
		near_a_line = (PFI) GMT_near_a_line_spherical;
	}
	else {
		GMT_distance_func = (PFD) GMT_cartesian_dist;
		near_a_line = (PFI) GMT_near_a_line_cartesian;
	}

	if (limit_point_distance) { 
		n_points = GMT_points_init (pt_file, &xp, &yp, &dp, pt_dist, greenwich);
		
		/* Copy xp into struct data, sort, and copy back */
		
		data = (struct DATA *) GMT_memory (VNULL, (size_t)n_points, sizeof(struct DATA), GMT_program);

		for (k = 0; k < n_points; k++) {
			data[k].i = k;
			data[k].x = xp[k];
			data[k].y = yp[k];
			data[k].d = dp[k];
		}
		/* Sort on x */
		qsort ((void *)data, (size_t)n_points, sizeof (struct DATA), compare_x);
		for (k = 0; k < n_points; k++) {	/* Put back the new order */
			xp[k] = data[k].x;
			yp[k] = data[k].y;
			dp[k] = data[k].d;
		}
		GMT_free ((void *)data);
	}

	if (limit_line_distance)  n_lines  = GMT_lines_init (ln_file, &line, ln_dist, greenwich);
	if (limit_polygon) {
		n_pols = GMT_lines_init (pol_file, &pol, -9999.0, greenwich);
		if (do_project) {	/* Convert all the polygons */
			for (i = 0; i < n_pols; i++) {
				if (pol[i].polar) {	/* Polygon is a polar cap, use simple polar (r-theta) transformation */
					for (j = 0; j < pol[i].np; j++) {
						r = 90.0 - fabs (pol[i].lat[j]);
						sincos (pol[i].lon[j] * D2R, &sa, &ca);
						pol[i].lon[j] = r * ca;
						pol[i].lat[j] = r * sa;
					}
				}
				else {
					for (j = 0; j < pol[i].np; j++) {
						GMT_geo_to_xy (pol[i].lon[j], pol[i].lat[j], &xx, &yy);
						pol[i].lon[j] = xx;
						pol[i].lat[j] = yy;
					}
				}
			}
		}
	}

	/* Now we are ready to take on some input values */
	
	if (n_files > 0)
		nofile = FALSE;
	else
		n_files = 1;
	n_args = (argc > 1) ? argc : 2;
		
	for (fno = 1; !done && fno < n_args; fno++) {	/* Loop over input files, if any */
		if (!nofile && argv[fno][0] == '-') continue;
		
		if (nofile) {	/* Just read standard input */
			fp = GMT_stdin;
			done = TRUE;
			if (gmtdefs.verbose) fprintf (stderr, "%s: Reading from standard input\n", GMT_program);
#ifdef SET_IO_MODE
			GMT_setmode (GMT_IN);
#endif
		}
		else if ((fp = GMT_fopen (argv[fno], GMT_io.r_mode)) == NULL) {
			fprintf (stderr, "%s: Cannot open file %s\n", GMT_program, argv[fno]);
			continue;
		}

		if (!nofile && gmtdefs.verbose) fprintf (stderr, "%s: Working on file %s\n", GMT_program, argv[fno]);
		
		if (gmtdefs.io_header[GMT_IN]) {
			for (i = 0; i < gmtdefs.n_header_recs; i++) {
				GMT_fgets (buffer, BUFSIZ, fp);
				if (first && gmtdefs.io_header[GMT_OUT]) fprintf (GMT_stdout, "%s", buffer);
			}
			first = FALSE;
		}
		n_expected_fields = (GMT_io.ncol[GMT_IN]) ? GMT_io.ncol[GMT_IN] : BUFSIZ;
	
		while ((n_fields = GMT_input (fp, &n_expected_fields, &in)) >= 0 && ! (GMT_io.status & GMT_IO_EOF)) {	/* Not yet EOF */

			while (GMT_io.status & GMT_IO_SEGMENT_HEADER && !(GMT_io.status & GMT_IO_EOF)) {
				GMT_write_segmentheader (GMT_stdout, n_expected_fields);
				n_fields = GMT_input (fp, &n_expected_fields, &in);
			}
			if ((GMT_io.status & GMT_IO_EOF)) continue;	/* At EOF */

			if (GMT_io.status & GMT_IO_MISMATCH) {
				fprintf (stderr, "%s: Mismatch between actual (%d) and expected (%d) fields near line %d\n", GMT_program, n_fields,  n_expected_fields, n_read);
				exit (EXIT_FAILURE);
			}

			n_read++;

			if (limit_z && n_expected_fields < 3) {
				fprintf (stderr, "%s: -Z requires a data file with at least 3 columns; this file only has %d. Exiting.\n", GMT_program, n_fields);
				exit (EXIT_FAILURE);
			}

			if (limit_z) {
				if (GMT_is_dnan (in[2])) continue;	/* cannot keep when no z */
				inside = (in[2] >= zmin && in[2] <= zmax); 
				if (inside != want_inside[5]) continue;
			}

			lon = in[0];
			if (project_info.region_supplied) {
				inside = !GMT_map_outside (lon, in[1]);
				if (inside != want_inside[0]) continue;
			}

			if (limit_point_distance_cart) {
				inside = GMT_near_a_point_cart (lon, in[1], xp, yp, dp, n_points); 
				if (inside != want_inside[1]) continue;
			}
			else if (limit_point_distance) {
				inside = GMT_near_a_point (lon, in[1], xp, yp, dp, n_points);
				if (inside != want_inside[1]) continue;
			}

			if (limit_line_distance) {
				inside = near_a_line (lon, in[1], line, n_lines, FALSE, NULL);
				if (inside != want_inside[2]) continue;
			}

			if (limit_polygon) {
				if (do_project)
					GMT_geo_to_xy (lon, in[1], &xx, &yy);
				else {
					xx = lon;
					yy = in[1];
				}
				i = 0;
				inside = FALSE;
				while (i < n_pols && !inside) {
					if (pol[i].polar == -1 && in[1] <= 0.0) {	/* Point on S hemisphere, check */
						r = 90.0 + in[1];
						sincos (lon * D2R, &sa, &ca);
						inside = (GMT_non_zero_winding (r * ca, r * sa, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					}
					else if (pol[i].polar == +1 && in[1] >= 0.0) {	/* Point on N hemisphere, check */
						r = 90.0 - in[1];
						sincos (lon * D2R, &sa, &ca);
						inside = (GMT_non_zero_winding (r * ca, r * sa, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					}
					else
						inside = (GMT_non_zero_winding (xx, yy, pol[i].lon, pol[i].lat, pol[i].np) >= is_inside);
					i++;
				}
				if (inside != want_inside[3]) continue;
			}

			if (limit_coast) {
				xx = lon;
				while (xx < 0.0) xx += 360.0;
				row = ((int)floor ((90.0 - in[1]) / c.bsize));
				if (row >= c.bin_ny) row = c.bin_ny - 1;	/* Presumably only kicks in for south pole */
				col = (int)floor (xx / c.bsize);
				bin = row * c.bin_nx + col;
				if (bin != last_bin) {	/* Do this upon entering new bin */
					ind = 0;
					while (ind < c.nb && c.bins[ind] != bin) ind++;	/* Set ind to right bin */
					if (ind == c.nb) continue;			/* Bin not among the chosen ones */
					last_bin = bin;
					GMT_get_shore_bin (ind, &c, min_area, min_level, max_level);

					/* Must use polygons.  Go in both directions to cover both land and sea */
					for (id = 0; id < 2; id++) {
						GMT_free_polygons (p[id], np[id]);
						if (np[id]) GMT_free ((void *)p[id]);
						np[id] = GMT_assemble_shore (&c, wd[id],  min_level, TRUE, greenwich, west_border, east_border, &p[id]);
						np[id] = GMT_prep_polygons (&p[id], np[id], greenwich, FALSE, 0.0, -1);
					}
				}

				if (c.ns == 0) {	/* No lines go through, check node level */
					this_node = MIN (MIN (c.node_level[0], c.node_level[1]) , MIN (c.node_level[2], c.node_level[3]));
				}
				else {
					this_node = 0;
					GMT_geo_to_xy (lon, in[1], &xx, &yy);
					for (id = 0; id < 2; id++) {
		
						for (k = 0; k < np[id]; k++) {
		
							if (p[id][k].n == 0) continue;

							/* Find min/max of polygon */
			
							xmin = xmax = p[id][k].lon[0];
							ymin = ymax = p[id][k].lat[0];

							for (i = 1; i < p[id][k].n; i++) {
								if (p[id][k].lon[i] < xmin) xmin = p[id][k].lon[i];
								if (p[id][k].lon[i] > xmax) xmax = p[id][k].lon[i];
								if (p[id][k].lat[i] < ymin) ymin = p[id][k].lat[i];
								if (p[id][k].lat[i] > ymax) ymax = p[id][k].lat[i];
							}
							
							if (yy < ymin || yy > ymax) continue;
							if (xx < xmin || xx > xmax) continue;

							/* Must compare with polygon */

							if ((side = GMT_non_zero_winding (xx, yy, p[id][k].lon, p[id][k].lat, p[id][k].n)) < is_inside) continue;	/* Outside polygon */
				
							/* Here, point is inside, we must assign value */
				
							if (p[id][k].level > this_node) this_node = p[id][k].level;
						}
					}
				}
				inside = out_edge_in[this_node];
				if (inside != want_inside[4]) continue;
			}

			/* Here, we have passed all test and the point is output */

			if (binary) {
				if (n_output < 0) n_output = n_expected_fields;
				GMT_output (GMT_stdout, n_output, in);
			}
			else
				fprintf (GMT_stdout, "%s", GMT_io.current_record);
			n_pass++;
			
			if (long_verbose && n_read%1000 == 0) fprintf (stderr, "%s: Read %d records\r", GMT_program, n_read);
		}
		if (fp != GMT_stdin) GMT_fclose (fp);
	}
	
	if (gmtdefs.verbose) fprintf (stderr, "%s: Read %d records, passed %d records\n", GMT_program, n_read, n_pass);

	if (limit_point_distance) {
		GMT_free ((void *)xp);
		GMT_free ((void *)yp);
		GMT_free ((void *)dp);
	}
	
	GMT_end (argc, argv);
}

int compare_x (const void *point_1, const void *point_2)
{
	struct DATA *p1, *p2;
	
	p1 = (struct DATA *)point_1;
	p2 = (struct DATA *)point_2;

	if (p1->x < p2->x)
		return (-1);
	else if (p1->x > p2->x)
		return (1);
	else
		return (0);
}
