/*-------------------------------------------------------------------- * $Id: gmtselect.c 12822 2014-01-31 23:39:56Z remko $ * * Copyright (c) 1991-2014 by P. Wessel, W. H. F. Smith, R. Scharroo, J. Luis and F. Wobbe * See LICENSE.TXT 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 Lesser General Public License as published by * the Free Software Foundation; version 3 or any later version. * * 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 Lesser General Public License for more details. * * Contact info: gmt.soest.hawaii.edu *--------------------------------------------------------------------*/ /* * API functions to support the gmtselect application. * * Author: Paul Wessel * Date: 1-JAN-2010 * Version: 5 API * * Brief synopsis: gmtselect is a general-purpose spatial filter. Data pass * or fail basedon one or more conditions. Six 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 * 6. Only data with z-values within specified range 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 */ #define THIS_MODULE_NAME "gmtselect" #define THIS_MODULE_LIB "core" #define THIS_MODULE_PURPOSE "Select data table subsets based on multiple spatial criteria" #include "gmt_dev.h" #define GMT_PROG_OPTIONS "-:>JRVabfghios" GMT_OPT("HMm") #define GMTSELECT_N_TESTS 6 /* Number of specific tests available */ #define GMTSELECT_N_CLASSES (GSHHS_MAX_LEVEL + 1) /* Number of bands separated by the levels */ #define F_ITEM 0 #define N_ITEM 1 struct GMTSELECT_DATA { /* Used for temporary storage when sorting data on x coordinate */ double x, y, d; }; struct GMTSELECT_ZLIMIT { /* Used to hold info for each -Z option given */ unsigned int col; /* Column to test */ bool equal; /* Just check if z == min withing 5 ULps */ double min; /* Smallest z-value to pass through, for this column */ double max; /* Largest z-value to pass through, for this column */ }; struct GMTSELECT_CTRL { /* All control options for this program (except common args) */ /* active is true if the option has been activated */ struct A { /* -A[//] */ bool active; struct GMT_SHORE_SELECT info; } A; struct C { /* [-C[-|=|+][unit]/] */ bool active; int mode; /* Form of distance calculation (can be negative) */ double dist; /* Radius of influence for each point */ char unit; /* Unit name */ char *file; /* Name of file with points */ } C; struct D { /* -D */ bool active; bool force; /* if true, select next highest level if current set is not avaialble */ char set; /* One of f, h, i, l, c */ } D; struct E { /* -E , = combination or f,n */ bool active; unsigned int inside[2]; /* if 2, then a point exactly on a polygon boundary is considered OUTSIDE, else 1 */ } E; struct L { /* -L[p][-|=|+][unit]/ */ bool active; unsigned int end_mode; /* Controls what happens beyond segment endpoints */ int mode; /* Form of distance calculation (can be negative) */ double dist; /* Distance of influence for each line */ char unit; /* Unit name */ char *file; /* Name of file with lines */ } L; struct F { /* -F */ bool active; char *file; /* Name of file with polygons */ } F; struct I { /* -Icflrsz */ bool active; bool pass[GMTSELECT_N_TESTS]; /* One flag for each setting */ } I; struct N { /* -N */ bool active; unsigned int mode; /* 1 if dry/wet only, 0 if 5 mask levels */ bool mask[GMTSELECT_N_CLASSES]; /* Mask for each level */ } N; struct Z { /* -Z/[+c] */ bool active; unsigned int n_tests; /* How many of these tests did we get */ unsigned int max_col; /* The largest column we encountered */ struct GMTSELECT_ZLIMIT *limit; } Z; struct dbg { /* -+step */ bool active; double step; } dbg; }; void *New_gmtselect_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ int i; struct GMTSELECT_CTRL *C; C = GMT_memory (GMT, NULL, 1, struct GMTSELECT_CTRL); /* Initialize values whose defaults are not 0/false/NULL */ C->A.info.high = GSHHS_MAX_LEVEL; /* Include all GSHHS levels */ C->D.set = 'l'; /* Low-resolution coastline data */ C->E.inside[F_ITEM] = C->E.inside[N_ITEM] = GMT_ONEDGE; /* Default is that points on a boundary are inside */ for (i = 0; i < GMTSELECT_N_TESTS; i++) C->I.pass[i] = true; /* Default is to pass if we are inside */ GMT_memset (C->N.mask, GMTSELECT_N_CLASSES, bool); /* Default for "wet" areas = false (outside) */ C->N.mask[1] = C->N.mask[3] = true; /* Default for "dry" areas = true (inside) */ C->Z.max_col = 1; /* Minimum number of columns to expect is 1 unless "x,y" data are implied [2] */ return (C); } void Free_gmtselect_Ctrl (struct GMT_CTRL *GMT, struct GMTSELECT_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->C.file) free (C->C.file); if (C->F.file) free (C->F.file); if (C->L.file) free (C->L.file); if (C->Z.n_tests) GMT_free (GMT, C->Z.limit); GMT_free (GMT, C); } int compare_x (const void *point_1, const void *point_2) { const struct GMTSELECT_DATA *p1 = point_1, *p2 = point_2; if (p1->x < p2->x) return (-1); if (p1->x > p2->x) return (1); return (0); } int GMT_gmtselect_usage (struct GMTAPI_CTRL *API, int level) { GMT_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_NAME, THIS_MODULE_PURPOSE); if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR); GMT_Message (API, GMT_TIME_NONE, "usage: gmtselect [] [%s]\n", GMT_A_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[-C%s/] [-D][+] [-E[f][n]] [-F] [%s]\n", GMT_DIST_OPT, GMT_J_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[-I[cflrsz] [-L[p]%s/] [-N]\n\t[%s] [%s] [%s]\n\t[-Z[/][+c]] [%s] [%s]\n\t[%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s]\n\n", GMT_DIST_OPT, GMT_Rgeo_OPT, GMT_V_OPT, GMT_a_OPT, GMT_b_OPT, GMT_f_OPT, GMT_g_OPT, GMT_h_OPT, GMT_i_OPT, GMT_o_OPT, GMT_s_OPT, GMT_colon_OPT); if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Option (API, "<,A"); GMT_Message (API, GMT_TIME_NONE, "\t (ignored unless -N is set).\n"); GMT_dist_syntax (API->GMT, 'C', "Pass locations that are within of any point in the ASCII ."); GMT_Message (API, GMT_TIME_NONE, "\t Give distance as 0 if 3rd column of has individual distances.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use -R -J to compute mapped Cartesian distances in cm, inch, m, or points [%s].\n", API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]); GMT_Message (API, GMT_TIME_NONE, "\t-D Choose one of the following resolutions: (Ignored unless -N is set).\n"); GMT_Message (API, GMT_TIME_NONE, "\t f - full resolution (may be very slow for large regions).\n"); GMT_Message (API, GMT_TIME_NONE, "\t h - high resolution (may be slow for large regions).\n"); GMT_Message (API, GMT_TIME_NONE, "\t i - intermediate resolution.\n"); GMT_Message (API, GMT_TIME_NONE, "\t l - low resolution [Default].\n"); GMT_Message (API, GMT_TIME_NONE, "\t c - crude resolution, for tasks that need crude continent outlines only.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append + to use a lower resolution should the chosen one not be available [abort].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-E Indicate if points exactly on a polygon boundary are inside or outside.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append f and/or n to modify the -F option or -N option, respectively,\n"); GMT_Message (API, GMT_TIME_NONE, "\t to consider such points to be outside the feature [inside].\n"); GMT_Message (API, GMT_TIME_NONE, "\t-F Pass locations that are inside the polygons in the ASCII file.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-I Reverse the tests, i.e., pass locations outside the region.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Supply a combination of cflrz where each flag means:\n"); GMT_Message (API, GMT_TIME_NONE, "\t c will pass locations beyond the minimum distance to the points in -C.\n"); GMT_Message (API, GMT_TIME_NONE, "\t f will pass locations outside the polygons in -F.\n"); GMT_Message (API, GMT_TIME_NONE, "\t l will pass locations beyond the minimum distance to the lines in -L.\n"); GMT_Message (API, GMT_TIME_NONE, "\t r will pass locations outside the region given in -R [and -J].\n"); GMT_Message (API, GMT_TIME_NONE, "\t s will pass locations that otherwise would be skipped in -N.\n"); GMT_Message (API, GMT_TIME_NONE, "\t z will pass locations outside the range given in -Z.\n"); GMT_Option (API, "J"); GMT_dist_syntax (API->GMT, 'L', "Pass locations that are within of any line in ASCII ."); GMT_Message (API, GMT_TIME_NONE, "\t Give distance as 0 if 2nd column of segment headers have individual distances.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Use -R -J to compute mapped Cartesian distances in cm, inch, or points [%s].\n", API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]); GMT_Message (API, GMT_TIME_NONE, "\t Optionally, use -Lp to exclude points projecting beyond a line's endpoints.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-N Set if a point outside or inside a geographic feature should be s(kipped) or k(ept).\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append o to let feature boundary be considered outside [Default is inside].\n"); GMT_Message (API, GMT_TIME_NONE, "\t Specify this information with s or k using 1 of 2 formats:\n"); GMT_Message (API, GMT_TIME_NONE, "\t -N/.\n"); GMT_Message (API, GMT_TIME_NONE, "\t -N////.\n"); GMT_Message (API, GMT_TIME_NONE, "\t k means keep and s means skip [Default is s/k/s/k/s (i.e., s/k)].\n"); GMT_Option (API, "R,V"); GMT_Message (API, GMT_TIME_NONE, "\t-Z Assume the 3rd data column contains z-values and we want to keep records with\n"); GMT_Message (API, GMT_TIME_NONE, "\t <= z <= . Use - for or if there is no lower/upper limit.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Append +cto select another column than the third [2].\n"); GMT_Message (API, GMT_TIME_NONE, "\t If is not given we pass records whose z equal .\n"); GMT_Message (API, GMT_TIME_NONE, "\t The -Z option is repeatable.\n"); GMT_Option (API, "a,bi0"); GMT_Message (API, GMT_TIME_NONE, "\t Default is 2 input columns (3 if -Z is used).\n"); GMT_Option (API, "bo,f,g,h,i,o,s,:,."); return (EXIT_FAILURE); } int GMT_gmtselect_parse (struct GMT_CTRL *GMT, struct GMTSELECT_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to gmtselect and sets parameters in CTRL. * Any GMT common options will override values set previously by other commands. * It also replaces any file names specified as input or output with the data ID * returned when registering these sources/destinations with the API. */ unsigned int n_errors = 0, pos, j, k, col, n_z_alloc = 0; char ptr[GMT_BUFSIZ] = {""}, buffer[GMT_BUFSIZ] = {""}, za[GMT_LEN64] = {""}, zb[GMT_LEN64] = {""}, *c = NULL; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; bool fix = false; for (opt = options; opt; opt = opt->next) { switch (opt->option) { case '<': /* Skip input files */ if (!GMT_check_filearg (GMT, '<', opt->arg, GMT_IN)) n_errors++; break; /* Processes program-specific parameters */ case 'A': /* Limit GSHHS features */ Ctrl->A.active = true; GMT_set_levels (GMT, opt->arg, &Ctrl->A.info); break; case 'C': /* Near a point test */ Ctrl->C.active = true; if (opt->arg[0] == 'f' && GMT_compat_check (GMT, 4)) { GMT_Report (API, GMT_MSG_COMPAT, "Warning: Option -Cf is deprecated; use -C- instead\n"); opt->arg[0] = '-'; fix = true; } for (j = 0; opt->arg[j] && opt->arg[j] != '/'; j++); if (opt->arg[j]) { Ctrl->C.file = strdup (&opt->arg[j+1]); opt->arg[j] = '\0'; /* Chop off the /filename part */ Ctrl->C.mode = GMT_get_distance (GMT, opt->arg, &(Ctrl->C.dist), &(Ctrl->C.unit)); opt->arg[j] = '/'; /* Restore the /filename part */ } else { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -C option: Expects -C%s/\n", GMT_DIST_OPT); n_errors++; } if (fix) opt->arg[0] = 'f'; break; case 'D': /* Set GSHHS resolution */ Ctrl->D.active = true; Ctrl->D.set = opt->arg[0]; if (opt->arg[1] == '+') Ctrl->D.force = true; break; case 'E': /* On-boundary selection */ Ctrl->E.active = true; for (j = 0; opt->arg[j]; j++) { switch (opt->arg[j]) { case 'f': Ctrl->E.inside[F_ITEM] = GMT_INSIDE; break; case 'n': Ctrl->E.inside[N_ITEM] = GMT_INSIDE; break; default: GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -E option: Expects -Ef, -En, or -Efn\n"); n_errors++; break; } } break; case 'F': /* Inside/outside polygon test */ if ((Ctrl->F.active = GMT_check_filearg (GMT, 'F', opt->arg, GMT_IN))) Ctrl->F.file = strdup (opt->arg); else n_errors++; break; case 'I': /* Invert these tests */ Ctrl->I.active = true; for (j = 0; opt->arg[j]; j++) { switch (opt->arg[j]) { case 'r': Ctrl->I.pass[0] = false; break; case 'c': Ctrl->I.pass[1] = false; break; case 'l': Ctrl->I.pass[2] = false; break; case 'f': Ctrl->I.pass[3] = false; break; case 's': Ctrl->I.pass[4] = false; break; case 'z': Ctrl->I.pass[5] = false; break; default: GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -I option: Expects -Icflrsz\n"); n_errors++; break; } } break; case 'L': /* Near a line test */ Ctrl->L.active = true; k = 0; if (opt->arg[k] == 'p') { /* Disallow points beyond endpoints */ Ctrl->L.end_mode = 10; k++; } for (j = k; opt->arg[j] && opt->arg[j] != '/'; j++); /* Find the first slash */ if (!opt->arg[j]) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -L option: Expects -L[p]%s/\n", GMT_DIST_OPT); n_errors++; } else { if (GMT_check_filearg (GMT, 'L', &opt->arg[j+1], GMT_IN)) Ctrl->L.file = strdup (&opt->arg[j+1]); else n_errors++; opt->arg[j] = '\0'; /* Chop off the /filename part */ Ctrl->L.mode = GMT_get_distance (GMT, &opt->arg[k], &(Ctrl->L.dist), &(Ctrl->L.unit)); opt->arg[j] = '/'; /* Restore the /filename part */ } break; case 'N': /* Inside/outside GSHHS land */ Ctrl->N.active = true; strncpy (buffer, opt->arg, GMT_BUFSIZ); if (buffer[strlen(buffer)-1] == 'o' && GMT_compat_check (GMT, 4)) { /* Edge is considered outside */ GMT_Report (API, GMT_MSG_COMPAT, "Warning: Option -N...o is deprecated; use -E instead\n"); Ctrl->E.active = true; Ctrl->E.inside[N_ITEM] = GMT_INSIDE; buffer[strlen(buffer)-1] = 0; } j = pos = 0; while (j < GMTSELECT_N_CLASSES && (GMT_strtok (buffer, "/", &pos, ptr))) { switch (ptr[0]) { case 's': /* Skip points in this level */ Ctrl->N.mask[j] = false; break; case 'k': Ctrl->N.mask[j] = true; break; default: GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -N option: Bad modifier (use s or k)\n"); n_errors++; } j++; } if (!(j == 2 || j == GMTSELECT_N_CLASSES)) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -N option: Specify 2 or 5 arguments\n"); n_errors++; } Ctrl->N.mode = (j == 2); break; case 'Z': /* Test column-ranges */ Ctrl->Z.active = true; col = ((c = strstr (opt->arg, "+c"))) ? atoi (&c[2]) : GMT_Z; if (c) c[0] = '\0'; j = sscanf (opt->arg, "%[^/]/%s", za, zb); if (c) c[0] = '+'; if (j < 1) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -Z option: Specify z_min [and z_max]\n"); n_errors++; } if (Ctrl->Z.n_tests == n_z_alloc) Ctrl->Z.limit = GMT_memory (GMT, Ctrl->Z.limit, n_z_alloc += 8, struct GMTSELECT_ZLIMIT); Ctrl->Z.limit[Ctrl->Z.n_tests].min = -DBL_MAX; Ctrl->Z.limit[Ctrl->Z.n_tests].max = +DBL_MAX; if (!(za[0] == '-' && za[1] == '\0')) n_errors += GMT_verify_expectations (GMT, GMT->current.io.col_type[GMT_IN][GMT_Z], GMT_scanf_arg (GMT, za, GMT->current.io.col_type[GMT_IN][GMT_Z], &Ctrl->Z.limit[Ctrl->Z.n_tests].min), za); if (j == 1) Ctrl->Z.limit[Ctrl->Z.n_tests].equal = true; else { if (!(zb[0] == '-' && zb[1] == '\0')) n_errors += GMT_verify_expectations (GMT, GMT->current.io.col_type[GMT_IN][GMT_Z], GMT_scanf_arg (GMT, zb, GMT->current.io.col_type[GMT_IN][GMT_Z], &Ctrl->Z.limit[Ctrl->Z.n_tests].max), zb); } n_errors += GMT_check_condition (GMT, Ctrl->Z.limit[Ctrl->Z.n_tests].max <= Ctrl->Z.limit[Ctrl->Z.n_tests].min, "Syntax error: -Z must have zmax > zmin!\n"); Ctrl->Z.limit[Ctrl->Z.n_tests].col = col; if (col > Ctrl->Z.max_col) Ctrl->Z.max_col = col; Ctrl->Z.n_tests++; break; #ifdef DEBUG case '+': /* Undocumented option to increase path-fix resolution */ Ctrl->dbg.active = true; Ctrl->dbg.step = atof (opt->arg); break; #endif default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } if (Ctrl->Z.max_col == 1 && (Ctrl->C.active || Ctrl->E.active || Ctrl->F.active || Ctrl->L.active || Ctrl->N.active || GMT->common.R.active)) Ctrl->Z.max_col = 2; if (Ctrl->Z.n_tests) Ctrl->Z.limit = GMT_memory (GMT, Ctrl->Z.limit, Ctrl->Z.n_tests, struct GMTSELECT_ZLIMIT); n_errors += GMT_check_condition (GMT, Ctrl->C.mode == -1, "Syntax error -C: Unrecognized distance unit\n"); n_errors += GMT_check_condition (GMT, Ctrl->C.mode == -2, "Syntax error -C: Unable to decode distance\n"); n_errors += GMT_check_condition (GMT, Ctrl->C.mode == -3, "Syntax error -C: Distance is negative\n"); n_errors += GMT_check_condition (GMT, Ctrl->C.active && GMT_access (GMT, Ctrl->C.file, R_OK), "Syntax error -C: Cannot read file %s!\n", Ctrl->C.file); n_errors += GMT_check_condition (GMT, Ctrl->F.active && GMT_access (GMT, Ctrl->F.file, R_OK), "Syntax error -F: Cannot read file %s!\n", Ctrl->F.file); n_errors += GMT_check_condition (GMT, Ctrl->L.active && GMT_access (GMT, Ctrl->L.file, R_OK), "Syntax error -L: Cannot read file %s!\n", Ctrl->L.file); n_errors += GMT_check_condition (GMT, Ctrl->L.mode == -1, "Syntax error -L: Unrecognized distance unit\n"); n_errors += GMT_check_condition (GMT, Ctrl->L.mode == -2, "Syntax error -L: Unable to decode distance\n"); n_errors += GMT_check_condition (GMT, Ctrl->L.mode == -3, "Syntax error -L: Distance is negative\n"); n_errors += GMT_check_condition (GMT, !Ctrl->N.active && (Ctrl->A.active || Ctrl->D.active), "Syntax error: -A and -D requires -N!\n"); n_errors += GMT_check_condition (GMT, Ctrl->L.active && Ctrl->C.active && !(Ctrl->C.mode == Ctrl->L.mode && Ctrl->C.unit == Ctrl->L.unit), "Syntax error: If both -C and -L are used they must use the same distance unit and calculation mode\n"); n_errors += GMT_check_binary_io (GMT, Ctrl->Z.max_col); n_errors += GMT_check_condition (GMT, Ctrl->Z.n_tests > 1 && Ctrl->I.active && !Ctrl->I.pass[5], "Syntax error: -Iz can only be used with one -Z range\n"); return (n_errors ? GMT_PARSE_ERROR : GMT_OK); } /* Must free allocated memory before returning */ #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_gmtselect_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_gmtselect (void *V_API, int mode, void *args) { int err; /* Required by GMT_err_fail */ unsigned int base = 3, np[2] = {0, 0}, r_mode; unsigned int side, col, id; int n_fields, ind, wd[2] = {0, 0}, n_minimum = 2, bin, last_bin = INT_MAX, error = 0; bool inside, need_header = false, shuffle, just_copy_record = false, pt_cartesian = false; bool output_header = false, do_project = false, no_resample = false, keep; uint64_t k, row, seg, n_read = 0, n_pass = 0; double xx, yy, *in = NULL; double west_border = 0.0, east_border = 0.0, xmin, xmax, ymin, ymax, lon; char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"}; struct GMT_DATATABLE *pol = NULL, *line = NULL, *point = NULL; struct GMT_GSHHS_POL *p[2] = {NULL, NULL}; struct GMT_SHORE c; struct GMT_DATASET *Cin = NULL, *Lin = NULL, *Fin = NULL; struct GMTSELECT_CTRL *Ctrl = NULL; struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; struct GMT_OPTION *options = NULL; struct GMTAPI_CTRL *API = GMT_get_API_ptr (V_API); /* Cast from void to GMTAPI_CTRL pointer */ /*----------------------- Standard module initialization and parsing ----------------------*/ if (API == NULL) return (GMT_NOT_A_SESSION); if (mode == GMT_MODULE_PURPOSE) return (GMT_gmtselect_usage (API, GMT_MODULE_PURPOSE)); /* Return the purpose of program */ options = GMT_Create_Options (API, mode, args); if (API->error) return (API->error); /* Set or get option list */ if (!options || options->option == GMT_OPT_USAGE) bailout (GMT_gmtselect_usage (API, GMT_USAGE));/* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_gmtselect_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments */ GMT = GMT_begin_module (API, THIS_MODULE_LIB, THIS_MODULE_NAME, &GMT_cpy); /* Save current state */ if (GMT_Parse_Common (API, GMT_PROG_OPTIONS, options)) Return (API->error); Ctrl = New_gmtselect_Ctrl (GMT); /* Allocate and initialize a new control structure */ if ((error = GMT_gmtselect_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the gmtselect main code ----------------------------*/ GMT_Report (API, GMT_MSG_VERBOSE, "Processing input table data\n"); if (Ctrl->C.active && !GMT_is_geographic (GMT, GMT_IN)) pt_cartesian = true; shuffle = (GMT->current.setting.io_lonlat_toggle[GMT_IN] != GMT->current.setting.io_lonlat_toggle[GMT_OUT]); /* Must rewrite output record */ n_minimum = Ctrl->Z.max_col; /* Minimum number of columns in ASCII input */ if (!GMT->common.R.active && Ctrl->N.active) { /* If we use coastline data or used -fg but didnt give -R we implicitly set -Rg */ GMT->common.R.active = true; GMT->common.R.wesn[XLO] = 0.0; GMT->common.R.wesn[XHI] = 360.0; GMT->common.R.wesn[YLO] = -90.0; GMT->common.R.wesn[YHI] = +90.0; GMT_set_geographic (GMT, GMT_IN); } if (GMT->common.R.active) { /* -R was set directly or indirectly; hence must set -J if not supplied */ if (!GMT->common.J.active) { /* -J not specified, set one implicitly */ /* Supply dummy linear proj */ GMT->current.proj.projection = GMT->current.proj.xyz_projection[GMT_X] = GMT->current.proj.xyz_projection[GMT_Y] = GMT_LINEAR; GMT->current.proj.pars[0] = GMT->current.proj.pars[1] = 1.0; GMT->common.J.active = no_resample = true; } else do_project = true; /* Only true when the USER selected -J, not when we supply a dummy -Jx1d */ Ctrl->dbg.step = 0.01; if (GMT_is_geographic (GMT, GMT_IN)) { while (GMT->common.R.wesn[XLO] < 0.0 && GMT->common.R.wesn[XHI] < 0.0) { /* Make all-negative longitude range positive instead */ GMT->common.R.wesn[XLO] += 360.0; GMT->common.R.wesn[XHI] += 360.0; } } GMT_err_fail (GMT, GMT_map_setup (GMT, GMT->common.R.wesn), ""); if (no_resample) GMT->current.map.parallel_straight = GMT->current.map.meridian_straight = 2; /* No resampling along bin boundaries */ } if (do_project) GMT_Report (API, GMT_MSG_VERBOSE, "Warning: -J means all data will be projected before tests are applied\n"); if (Ctrl->N.active) { /* Set up GSHHS */ if (Ctrl->D.force) Ctrl->D.set = GMT_shore_adjust_res (GMT, Ctrl->D.set); if (Ctrl->D.active) base = GMT_set_resolution (GMT, &Ctrl->D.set, 'D'); if (Ctrl->N.mode) { /* Post-process -N choice */ Ctrl->N.mask[3] = Ctrl->N.mask[1]; Ctrl->N.mask[2] = Ctrl->N.mask[4] = Ctrl->N.mask[0]; } if (GMT_init_shore (GMT, Ctrl->D.set, &c, GMT->common.R.wesn, &Ctrl->A.info)) { GMT_Report (API, GMT_MSG_NORMAL, "%s resolution shoreline data base not installed\n", shore_resolution[base]); Return (EXIT_FAILURE); } GMT_Report (API, GMT_MSG_LONG_VERBOSE, "GSHHG version %s\n%s\n%s\n", c.version, c.title, c.source); west_border = floor (GMT->common.R.wesn[XLO] / c.bsize) * c.bsize; east_border = ceil (GMT->common.R.wesn[XHI] / c.bsize) * c.bsize; wd[0] = 1; wd[1] = -1; np[0] = np[1] = 0; } just_copy_record = (GMT_is_ascii_record (GMT) && !shuffle); /* Initiate pointer to distance calculation function */ if (GMT_is_geographic (GMT, GMT_IN) && !do_project) { /* Geographic data and no -R -J conversion */ if (Ctrl->C.active) GMT_init_distaz (GMT, Ctrl->C.unit, Ctrl->C.mode, GMT_MAP_DIST); else if (Ctrl->L.active) GMT_init_distaz (GMT, Ctrl->L.unit, Ctrl->L.mode, GMT_MAP_DIST); } else if (do_project) /* Lon/lat projected via -R -J */ GMT_init_distaz (GMT, 'Z', 0, GMT_MAP_DIST); /* Compute r-squared instead of r after projection to avoid hypot */ else /* Cartesian data */ GMT_init_distaz (GMT, 'R', 0, GMT_MAP_DIST); /* Compute r-squared instead of r to avoid hypot */ if (Ctrl->C.active) { /* Initialize point structure used in test for proximity to points [use Ctrl->C.dist ]*/ if ((Cin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_IO_ASCII, NULL, Ctrl->C.file, NULL)) == NULL) { Return (API->error); } if (Cin->n_columns < 2) { /* Trouble */ GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -C option: %s does not have at least 2 columns with coordinates\n", Ctrl->C.file); Return (EXIT_FAILURE); } if (Ctrl->C.dist == 0.0 && Cin->n_columns <= 2) { /* Trouble */ GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -C option: %s does not have a 3rd column with distances, yet -C0/ was given\n", Ctrl->C.file); Return (EXIT_FAILURE); } point = Cin->table[0]; /* Can only be one table since we read a single file */ for (seg = 0; seg < point->n_segments; seg++) { if (Cin->n_columns == 2) point->segment[seg]->dist = Ctrl->C.dist; if (do_project) { /* Convert all the points using the map projection */ for (row = 0; row < point->segment[seg]->n_rows; row++) { GMT_geo_to_xy (GMT, point->segment[seg]->coord[GMT_X][row], point->segment[seg]->coord[GMT_Y][row], &xx, &yy); point->segment[seg]->coord[GMT_X][row] = xx; point->segment[seg]->coord[GMT_Y][row] = yy; } pt_cartesian = true; /* Well, now it is */ } } if (pt_cartesian) { /* Speed up testing by sorting points on the x-coordinate first */ struct GMTSELECT_DATA *data = NULL; /* Used for temporary storage when sorting data on x coordinate */ /* Copy xp into struct data, sort, and copy back */ data = GMT_memory (GMT, NULL, point->n_records, struct GMTSELECT_DATA); for (seg = k = 0; seg < point->n_segments; seg++) { for (row = 0; row < point->segment[seg]->n_rows; row++, k++) { data[k].x = point->segment[seg]->coord[GMT_X][row]; data[k].y = point->segment[seg]->coord[GMT_Y][row]; data[k].d = (Ctrl->C.dist == 0.0) ? point->segment[seg]->coord[GMT_Z][row] : Ctrl->C.dist; } } /* Sort on x to speed up inside testing */ qsort (data, point->n_records, sizeof (struct GMTSELECT_DATA), compare_x); for (seg = k = 0; seg < point->n_segments; seg++) { /* Put back the new order */ for (row = 0; row < point->segment[seg]->n_rows; row++, k++) { point->segment[seg]->coord[GMT_X][row] = data[k].x; point->segment[seg]->coord[GMT_Y][row] = data[k].y; if (Ctrl->C.dist == 0.0) point->segment[seg]->coord[GMT_Z][row] = data[k].d ; } } GMT_free (GMT, data); } } if (Ctrl->L.active) { /* Initialize lines structure used in test for proximity to lines [use Ctrl->L.dist, ] */ if ((Lin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_LINE, GMT_IO_ASCII, NULL, Ctrl->L.file, NULL)) == NULL) { Return (API->error); } if (Lin->n_columns < 2) { /* Trouble */ GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -L option: %s does not have at least 2 columns with coordinates\n", Ctrl->L.file); Return (EXIT_FAILURE); } line = Lin->table[0]; /* Can only be one table since we read a single file */ for (seg = 0; seg < line->n_segments; seg++) { if (Ctrl->L.dist > 0.0) line->segment[seg]->dist = Ctrl->L.dist; /* Only override when nonzero */ if (do_project) { /* Convert all the line points using the map projection */ for (row = 0; row < line->segment[seg]->n_rows; row++) { GMT_geo_to_xy (GMT, line->segment[seg]->coord[GMT_X][row], line->segment[seg]->coord[GMT_Y][row], &xx, &yy); line->segment[seg]->coord[GMT_X][row] = xx; line->segment[seg]->coord[GMT_Y][row] = yy; } } } } if (Ctrl->F.active) { /* Initialize polygon structure used in test for polygon in/out test */ GMT_skip_xy_duplicates (GMT, true); /* Avoid repeating x/y points in polygons */ if ((Fin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POLY, GMT_IO_ASCII, NULL, Ctrl->F.file, NULL)) == NULL) { Return (API->error); } GMT_skip_xy_duplicates (GMT, false); /* Reset */ if (Fin->n_columns < 2) { /* Trouble */ GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -F option: %s does not have at least 2 columns with coordinates\n", Ctrl->F.file); Return (EXIT_FAILURE); } pol = Fin->table[0]; /* Can only be one table since we read a single file */ if (do_project) { /* Convert all the polygons points using the map projection */ for (seg = 0; seg < pol->n_segments; seg++) { for (row = 0; row < pol->segment[seg]->n_rows; row++) { GMT_geo_to_xy (GMT, pol->segment[seg]->coord[GMT_X][row], pol->segment[seg]->coord[GMT_Y][row], &xx, &yy); pol->segment[seg]->coord[GMT_X][row] = xx; pol->segment[seg]->coord[GMT_Y][row] = yy; } } } } /* Gather input/output file names (or stdin/out) and enable i/o */ if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_IN, GMT_ADD_DEFAULT, 0, options) != GMT_OK) { /* Establishes data input */ Return (API->error); } if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_POINT, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_OK) { /* Establishes data output */ Return (API->error); } if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_IN, GMT_HEADER_ON) != GMT_OK) { /* Enables data input and sets access mode */ Return (API->error); } if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_OK) { /* Enables data output and sets access mode */ Return (API->error); } /* Now we are ready to take on some input values */ GMT->common.b.ncol[GMT_OUT] = UINT_MAX; /* Flag to have it reset to GMT->common.b.ncol[GMT_IN] when writing */ r_mode = (just_copy_record) ? GMT_READ_MIXED : GMT_READ_DOUBLE; GMT_set_segmentheader (GMT, GMT_OUT, false); /* Since processing of -C|L|F files might have turned it on [should be determined below] */ do { /* Keep returning records until we reach EOF */ if ((in = GMT_Get_Record (API, r_mode, &n_fields)) == NULL) { /* Read next record, get NULL if special case */ if (GMT_REC_IS_ERROR (GMT)) /* Bail if there are any read errors */ Return (GMT_RUNTIME_ERROR); if (GMT_REC_IS_TABLE_HEADER (GMT)) { /* Echo table headers */ GMT_Put_Record (API, GMT_WRITE_TABLE_HEADER, NULL); continue; } if (GMT_REC_IS_EOF (GMT)) /* Reached end of file */ break; else if (GMT_REC_IS_SEGMENT_HEADER (GMT)) { output_header = true; need_header = GMT->current.io.multi_segments[GMT_OUT]; /* Only need to break up segments */ continue; } } /* Data record to process */ n_read++; if (n_read%1000 == 0) GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Read %" PRIu64 " records, passed %" PRIu64 "records\r", n_read, n_pass); if (n_fields < n_minimum) { /* Bad number of columns */ if (Ctrl->Z.active) GMT_Report (API, GMT_MSG_NORMAL, "-Z requires a data file with at least %u columns; this file only has %d near line %" PRIu64 ". Exiting.\n", n_minimum, n_fields, n_read); else GMT_Report (API, GMT_MSG_NORMAL, "Data file must have at least 2 columns; this file only has %d near line %" PRIu64 ". Exiting.\n", n_fields, n_read); Return (EXIT_FAILURE); } if (Ctrl->Z.active) { /* Apply z-range test(s) */ for (k = 0, keep = true; keep && k < Ctrl->Z.n_tests; k++) { col = Ctrl->Z.limit[k].col; /* Shorthand notation */ if (GMT_is_dnan (in[col])) keep = false; /* Cannot keep when no z-value */ if (Ctrl->Z.limit[k].equal) inside = doubleAlmostEqualZero (in[col], Ctrl->Z.limit[k].min); else inside = (in[col] >= Ctrl->Z.limit[k].min && in[col] <= Ctrl->Z.limit[k].max); if (inside != Ctrl->I.pass[5]) keep = false; } if (!keep) { output_header = need_header; continue;} } lon = in[GMT_X]; /* Use copy since we may have to wrap 360 */ if (GMT->common.R.active) { /* Apply region test */ inside = !GMT_map_outside (GMT, lon, in[GMT_Y]); if (inside != Ctrl->I.pass[0]) { output_header = need_header; continue;} } if (do_project) /* First project the input point */ GMT_geo_to_xy (GMT, lon, in[GMT_Y], &xx, &yy); else { xx = lon; yy = in[GMT_Y]; } if (Ctrl->C.active) { /* Check for distance to points */ inside = GMT_near_a_point (GMT, xx, yy, point, Ctrl->C.dist); if (inside != Ctrl->I.pass[1]) { output_header = need_header; continue;} } if (Ctrl->L.active) { /* Check for distance to lines */ inside = GMT_near_lines (GMT, xx, yy, line, Ctrl->L.end_mode, NULL, NULL, NULL); if (inside != Ctrl->I.pass[2]) { output_header = need_header; continue;} } if (Ctrl->F.active) { /* Check if we are in/out-side polygons */ if (do_project) /* Projected lon/lat; temporary reset input type for GMT_inonout to do Cartesian mode */ GMT_set_cartesian (GMT, GMT_IN); inside = 0; for (seg = 0; seg < pol->n_segments && !inside; seg++) { /* Check each polygon until we find that our point is inside */ if (GMT_polygon_is_hole (pol->segment[seg])) continue; /* Holes are handled within GMT_inonout */ inside = (GMT_inonout (GMT, xx, yy, pol->segment[seg]) >= Ctrl->E.inside[F_ITEM]); } if (do_project) /* Reset input type for GMT_inonout to do Cartesian mode */ GMT_set_geographic (GMT, GMT_IN); if (inside != Ctrl->I.pass[3]) { output_header = need_header; continue;} } if (Ctrl->N.active) { /* Check if on land or not */ int brow, i, this_node; xx = lon; while (xx < 0.0) xx += 360.0; brow = irint (floor ((90.0 - in[GMT_Y]) / c.bsize)); if (brow >= c.bin_ny) brow = c.bin_ny - 1; /* Presumably only kicks in for south pole */ col = urint (floor (xx / c.bsize)); bin = brow * 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_free_shore (GMT, &c); /* Free previously allocated arrays */ if ((err = GMT_get_shore_bin (GMT, ind, &c))) { GMT_Report (API, GMT_MSG_NORMAL, "%s [%s resolution shoreline]\n", GMT_strerror(err), shore_resolution[base]); Return (EXIT_FAILURE); } /* Must use polygons. Go in both directions to cover both land and sea */ for (id = 0; id < 2; id++) { GMT_free_shore_polygons (GMT, p[id], np[id]); if (np[id]) GMT_free (GMT, p[id]); np[id] = GMT_assemble_shore (GMT, &c, wd[id], true, west_border, east_border, &p[id]); np[id] = GMT_prep_shore_polygons (GMT, &p[id], np[id], !no_resample, Ctrl->dbg.step, -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 (GMT, lon, in[GMT_Y], &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; holes are handled explicitly via the levels */ if ((side = GMT_non_zero_winding (GMT, xx, yy, p[id][k].lon, p[id][k].lat, p[id][k].n)) < Ctrl->E.inside[N_ITEM]) 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 = Ctrl->N.mask[this_node]; if (inside != Ctrl->I.pass[4]) { output_header = need_header; continue;} } /* Here, we have passed all active test and the point is to be output */ if (output_header) { /* First output point for this segment - write the header */ GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); output_header = false; } if (just_copy_record) GMT_Put_Record (API, GMT_WRITE_TEXT, NULL); else GMT_Put_Record (API, GMT_WRITE_DOUBLE, in); n_pass++; } while (true); if (GMT_End_IO (API, GMT_IN, 0) != GMT_OK) { /* Disables further data input */ Return (API->error); } if (GMT_End_IO (API, GMT_OUT, 0) != GMT_OK) { /* Disables further data output */ Return (API->error); } GMT_Report (API, GMT_MSG_VERBOSE, "Read %" PRIu64 " records, passed %" PRIu64" records\n", n_read, n_pass); if (Ctrl->N.active) { GMT_free_shore (GMT, &c); GMT_shore_cleanup (GMT, &c); for (id = 0; id < 2; id++) { GMT_free_shore_polygons (GMT, p[id], np[id]); if (np[id]) GMT_free (GMT, p[id]); } } Return (GMT_OK); }