Skip to content

Commit 9d13399

Browse files
authored
Let gmt select -C accept a single point as an alternative to a file (#4009)
* Let gmt seelct -C accept a single point instead of file If we only have a single point then it is a hassle to make a file in order to pass it to gmt select. Closes #4008. * Update gmtselect.rst * Update gmtselect.c
1 parent 9ebe118 commit 9d13399

File tree

3 files changed

+43
-21
lines changed

3 files changed

+43
-21
lines changed

doc/rst/source/gmtselect.rst

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ Synopsis
1414

1515
**gmt select** [ *table* ]
1616
[ |SYN_OPT-Area| ]
17-
[ |-C|\ *pointfile*\ **+d**\ *dist* ]
17+
[ |-C|\ *pointfile*\ \|\ *lon*/*lat*\ **+d**\ *dist* ]
1818
[ |-D|\ *resolution*\ [**+f**] ]
1919
[ |-E|\ [**fn**] ]
2020
[ |-F|\ *polygonfile* ]
@@ -71,12 +71,14 @@ Optional Arguments
7171

7272
.. _-C:
7373

74-
**-C**\ *pointfile*\ **+d**\ *dist*
74+
**-C**\ *pointfile*\ \|\ *lon*/*lat*\ **+d**\ *dist*
7575
Pass all records whose location is within *dist* of any of the
7676
points in the ASCII file *pointfile*. If *dist* is zero then the 3rd
7777
column of *pointfile* must have each point's individual radius of
78-
influence. Distances are Cartesian and in user units; specify
79-
**-fg** to indicate spherical distances and append a distance unit
78+
influence. . If you only have a single point then you can specify
79+
*lon*/*lat* instead of *pointfile*. Distances are Cartesian and in
80+
user units; specify **-fg** to indicate spherical distances and
81+
append a distance unit, even if the distance specified is 0.
8082
(see `Units`_). Alternatively, if **-R** and **-J** are used then
8183
geographic coordinates are projected to map coordinates (in cm,
8284
inch, or points, as determined by :term:`PROJ_LENGTH_UNIT`) before

src/gmt_io.c

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5182,8 +5182,6 @@ char *gmt_get_filename (struct GMTAPI_CTRL *API, const char* filename, const cha
51825182
if ((c = gmt_first_modifier (API->GMT, file, mods)))
51835183
c[0] = '\0'; /* Begone with you */
51845184
}
5185-
if (file[0] == ' ')
5186-
c++;
51875185
clean_file = strdup (file);
51885186

51895187
GMT_Report (API, GMT_MSG_DEBUG, "gmt_get_filename: In: %s Out: %s\n", filename, clean_file);

src/gmtselect.c

Lines changed: 37 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@
2222
* Version: 6 API
2323
*
2424
* Brief synopsis: gmtselect is a general-purpose spatial filter. Data pass
25-
* or fail basedon one or more conditions. Seven conditions may be set:
25+
* or fail based on one or more conditions. Seven conditions may be set:
2626
*
2727
* 1. Only data inside a rectangular area may pass
2828
* 2. Only data within a certain distance from given points may pass
@@ -87,8 +87,10 @@ struct GMTSELECT_CTRL { /* All control options for this program (except common a
8787
} A;
8888
struct GMTSELECT_C { /* [-C[-|=|+]<dist>/<ptfile>] */
8989
bool active;
90+
bool point; /* True if we got a single point */
9091
int mode; /* Form of distance calculation (can be negative) */
9192
double dist; /* Radius of influence for each point */
93+
double lon, lat; /* Coordinates give for a single point */
9294
char unit; /* Unit name */
9395
char *file; /* Name of file with points */
9496
} C;
@@ -179,7 +181,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
179181
const char *name = gmt_show_name_and_purpose (API, THIS_MODULE_LIB, THIS_MODULE_CLASSIC_NAME, THIS_MODULE_PURPOSE);
180182
if (level == GMT_MODULE_PURPOSE) return (GMT_NOERROR);
181183
GMT_Message (API, GMT_TIME_NONE, "usage: %s [<table>] [%s]\n", name, GMT_A_OPT);
182-
GMT_Message (API, GMT_TIME_NONE, "\t[-C<ptfile>+d%s] [-D<resolution>][+f] [-E[f][n]] [-F<polygon>] [-G<gridmask>] [%s]\n",
184+
GMT_Message (API, GMT_TIME_NONE, "\t[-C<ptfile|lon/lat>+d%s] [-D<resolution>][+f] [-E[f][n]] [-F<polygon>] [-G<gridmask>] [%s]\n",
183185
GMT_DIST_OPT, GMT_J_OPT);
184186
GMT_Message (API, GMT_TIME_NONE, "\t[-I[cfglrsz] [-L<lfile>+d%s[+p]] [-N<info>] [%s]\n\t[%s] [-Z<min>[/<max>][+c<col>][+a][+i]] [%s] "
185187
"[%s]\n\t[%s] [%s] [%s] [%s]\n\t[%s] [%s]\n\t[%s] [%s] [%s] [%s] [%s]\n\n",
@@ -194,6 +196,7 @@ static int usage (struct GMTAPI_CTRL *API, int level) {
194196
GMT_Message (API, GMT_TIME_NONE, "\t (ignored unless -N is set).\n");
195197
gmt_dist_syntax (API->GMT, 'C', "Pass locations that are within <dist> of any point in the ASCII <ptfile>.");
196198
GMT_Message (API, GMT_TIME_NONE, "\t Give distance as 0 if 3rd column of <ptfile> has individual distances.\n");
199+
GMT_Message (API, GMT_TIME_NONE, "\t For a single point you can instead specify <lon>/<lat>+d[unit].\n");
197200
GMT_Message (API, GMT_TIME_NONE, "\t Use -R -J to compute mapped Cartesian distances in cm, inch, m, or points [%s].\n",
198201
API->GMT->session.unit_name[API->GMT->current.setting.proj_length_unit]);
199202
GMT_Message (API, GMT_TIME_NONE, "\t-D Choose one of the following resolutions: (Ignored unless -N is set).\n");
@@ -336,7 +339,7 @@ static int parse (struct GMT_CTRL *GMT, struct GMTSELECT_CTRL *Ctrl, struct GMT_
336339
Ctrl->A.active = true;
337340
n_errors += gmt_set_levels (GMT, opt->arg, &Ctrl->A.info);
338341
break;
339-
case 'C': /* Near a point test Syntax -C<pfile>+d<distance> */
342+
case 'C': /* Near a point test Syntax -C<pfile>+d<distance> or -C<lon/lat>+d<distance> */
340343
Ctrl->C.active = true;
341344
if ((c = strstr (opt->arg, "+d")) == NULL) { /* Must be old syntax or error */
342345
n_errors += gmtselect_old_C_parser (API, opt->arg, Ctrl);
@@ -347,8 +350,16 @@ static int parse (struct GMT_CTRL *GMT, struct GMTSELECT_CTRL *Ctrl, struct GMT_
347350
GMT_Report (API, GMT_MSG_ERROR, "Option -C: No file given\n");
348351
n_errors++;
349352
}
350-
else
353+
else {
351354
Ctrl->C.file = gmt_get_filename (API, opt->arg, "d");
355+
if (gmt_count_char (GMT, Ctrl->C.file, '/') == 1 && gmt_access (GMT, Ctrl->C.file, R_OK)) { /* Check if we got a lon/lat point */
356+
if ((j = sscanf (Ctrl->C.file, "%[^/]/%s", za, zb)) != 2) continue; /* No, strange... */
357+
n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_X), gmt_scanf_arg (GMT, za, gmt_M_type (GMT, GMT_IN, GMT_X), false, &Ctrl->C.lon), za);
358+
n_errors += gmt_verify_expectations (GMT, gmt_M_type (GMT, GMT_IN, GMT_Y), gmt_scanf_arg (GMT, zb, gmt_M_type (GMT, GMT_IN, GMT_Y), false, &Ctrl->C.lat), zb);
359+
if (n_errors == 0)
360+
Ctrl->C.point = true;
361+
}
362+
}
352363
Ctrl->C.mode = gmt_get_distance (GMT, &c[2], &(Ctrl->C.dist), &(Ctrl->C.unit));
353364
break;
354365
case 'D': /* Set GSHHS resolution */
@@ -507,7 +518,7 @@ static int parse (struct GMT_CTRL *GMT, struct GMTSELECT_CTRL *Ctrl, struct GMT_
507518
n_errors += gmt_M_check_condition (GMT, Ctrl->C.mode == -1, "Option -C: Unrecognized distance unit\n");
508519
n_errors += gmt_M_check_condition (GMT, Ctrl->C.mode == -2, "Option -C: Unable to decode distance\n");
509520
n_errors += gmt_M_check_condition (GMT, Ctrl->C.mode == -3, "Option -C: Distance is negative\n");
510-
n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && gmt_access (GMT, Ctrl->C.file, R_OK),
521+
n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->C.point && gmt_access (GMT, Ctrl->C.file, R_OK),
511522
"Option -C: Cannot read file %s!\n", Ctrl->C.file);
512523
n_errors += gmt_M_check_condition (GMT, Ctrl->F.active && gmt_access (GMT, Ctrl->F.file, R_OK),
513524
"Option -F: Cannot read file %s!\n", Ctrl->F.file);
@@ -640,17 +651,28 @@ EXTERN_MSC int GMT_gmtselect (void *V_API, int mode, void *args) {
640651

641652
gmt_disable_bghi_opts (GMT); /* Do not want any -b -g -h -i to affect the reading from -C,-F,-L files */
642653

643-
if (Ctrl->C.active) { /* Initialize point structure used in test for proximity to points [use Ctrl->C.dist ]*/
644-
if ((Cin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_IO_ASCII, NULL, Ctrl->C.file, NULL)) == NULL) {
645-
Return (API->error);
646-
}
647-
if (Cin->n_columns < 2) { /* Trouble */
648-
GMT_Report (API, GMT_MSG_ERROR, "Option -C: %s does not have at least 2 columns with coordinates\n", Ctrl->C.file);
649-
Return (GMT_RUNTIME_ERROR);
654+
if (Ctrl->C.active) { /* Initialize point structure used in test for proximity to points [use Ctrl->C.dist] */
655+
if (Ctrl->C.point) {
656+
uint64_t dim[4] = {1, 1, 1, 2};
657+
if ((Cin = GMT_Create_Data (API, GMT_IS_DATASET, GMT_IS_POINT, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
658+
Return (API->error);
659+
}
660+
/* Place the given point in the array */
661+
Cin->table[0]->segment[0]->data[GMT_X][0] = Ctrl->C.lon;
662+
Cin->table[0]->segment[0]->data[GMT_Y][0] = Ctrl->C.lat;
650663
}
651-
if (Ctrl->C.dist == 0.0 && Cin->n_columns <= 2) { /* Trouble */
652-
GMT_Report (API, GMT_MSG_ERROR, "Option -C: %s does not have a 3rd column with distances, yet -C0/<file> was given\n", Ctrl->C.file);
653-
Return (GMT_RUNTIME_ERROR);
664+
else { /* Read from file */
665+
if ((Cin = GMT_Read_Data (API, GMT_IS_DATASET, GMT_IS_FILE, GMT_IS_POINT, GMT_IO_ASCII, NULL, Ctrl->C.file, NULL)) == NULL) {
666+
Return (API->error);
667+
}
668+
if (Cin->n_columns < 2) { /* Trouble */
669+
GMT_Report (API, GMT_MSG_ERROR, "Option -C: %s does not have at least 2 columns with coordinates\n", Ctrl->C.file);
670+
Return (GMT_RUNTIME_ERROR);
671+
}
672+
if (Ctrl->C.dist == 0.0 && Cin->n_columns <= 2) { /* Trouble */
673+
GMT_Report (API, GMT_MSG_ERROR, "Option -C: %s does not have a 3rd column with distances, yet -C0/<file> was given\n", Ctrl->C.file);
674+
Return (GMT_RUNTIME_ERROR);
675+
}
654676
}
655677
point = Cin->table[0]; /* Can only be one table since we read a single file */
656678
for (seg = 0; seg < point->n_segments; seg++) {

0 commit comments

Comments
 (0)