1 /*--------------------------------------------------------------------
2 *
3 * Copyright (c) 1991-2021 by the GMT Team (https://www.generic-mapping-tools.org/team.html)
4 * See LICENSE.TXT file for copying and redistribution conditions.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation; version 3 or any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * Contact info: www.generic-mapping-tools.org
16 *--------------------------------------------------------------------*/
17
18 #include "gmt_dev.h"
19 #include "gmt_internals.h"
20
21 /* Misc functions to find and read DCW polygons.
22 * Some of the countries have state borders too.
23 * The PUBLIC functions are (5):
24 *
25 * gmt_DCW_option : Present the DCW option and modifier usage
26 * gmt_DCW_parse : Parse the DCW option arguments
27 * gmt_DCW_list : List the available polygons and exit
28 * gmt_DCW_operation : Get DCW polygons and operate on them
29 * gmt_DCW_free : Free memory allocated by gmt_DCW_parse
30 *
31 * Author: Paul Wessel
32 * Date: 1-MAY-2013
33 * Version: 5.x
34 *
35 * We expect to find the file dcw-gmt.nc, dcw-countries.txt, and dcw-states.txt
36 * in one of the dirs accessible to GMT or pointed to by the default DIR_DCW.
37 * See separate subversion project DCW for the maintenance of the raw files that
38 * are used to build the netCDF file [svn://gmtserver.soest.hawai.edu/DCW].
39 */
40
41 #define DCW_SITE "ftp://ftp.soest.hawaii.edu/gmt"
42 #define GMT_DCW_N_CONTINENTS 8
43
44 #define DCW_GET_COUNTRY 1 /* Extract countries only */
45 #define DCW_GET_COUNTRY_AND_STATE 2 /* Extract countries and states */
46 #define DCW_DO_OUTLINE 1 /* Draw outline of polygons */
47 #define DCW_DO_FILL 2 /* Fill the polygons */
48
49 struct GMT_DCW_COUNTRY { /* Information per country */
50 char continent[4]; /* 2-char continent code (EU, NA, SA, AF, AU, AN) */
51 char code[4]; /* 2-char country code ISO 3166-1 (e.g., NO, US) */
52 char name[80]; /* Full name of the country */
53 };
54
55 struct GMT_DCW_STATE { /* Information per state */
56 char country[4]; /* 2-char country code ISO 3166-1 (e.g., BR, US) */
57 char code[4]; /* 2/3-char state codes for US, Canada, China, Argentina, Australia, Brazil, Russia (e.g., TX) */
58 char name[80]; /* Full name of the state */
59 };
60
61 struct GMT_DCW_COUNTRY_STATE { /* Information per country with state */
62 char country[4]; /* 2/3-char country code ISO 3166-1 (e.g., BR, US) for countries with states */
63 };
64
65 /* For version 2.0.0 we follow https://en.wikipedia.org/wiki/ISO_3166-2:CN and use 2-char instead of int codes */
66 #define DCW_N_CHINA_PROVINCES 34
67 struct GMT_DCW_CHINA_CODES {
68 unsigned int id;
69 char code[4];
70 };
71
72 /* Compile in read-only structures and arrays with the information */
73
74 static char *GMT_DCW_continents[GMT_DCW_N_CONTINENTS] = {"Africa", "Antarctica", "Asia", "Europe", "Oceania", "North America", "South America", "Miscellaneous"};
75
76 /* Local functions only visible inside this file */
77
gmtdcw_get_path(struct GMT_CTRL * GMT,char * name,char * suffix,char * path)78 GMT_LOCAL bool gmtdcw_get_path (struct GMT_CTRL *GMT, char *name, char *suffix, char *path) {
79 /* This is the order of checking:
80 * 1. Check in GMT->session.DCWDIR, if set
81 * 2. Look via gmt_getsharepath.
82 * 3. Look in userdir/geography/dcw
83 * 4. Try to download from server into (3)
84 * 5. Give up.
85 */
86
87 if (GMT->session.DCWDIR) { /* 1. Check in GMT->session.DCWDIR */
88 sprintf (path, "%s/%s%s", GMT->session.DCWDIR, name, suffix);
89 if (access (path, R_OK) == 0) {
90 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "1. DCW: Read the Digital Chart of the World from %s\n", path);
91 return true;
92 }
93 /* Failed, so remove reference to invalid GMT->session.DCWDIR but don't free
94 * the pointer. this is no leak because the reference still exists
95 * in the previous copy of the current GMT_CTRL struct. */
96 GMT->session.DCWDIR = NULL;
97 }
98 if (gmt_getsharepath (GMT, "dcw", name, suffix, path, R_OK)) {
99 if (access (path, R_OK) == 0) {
100 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "2. DCW: Read the Digital Chart of the World from %s\n", path);
101 return true; /* Found it in share or user somewhere */
102 }
103 }
104 if (GMT->session.USERDIR) { /* Check user dir via remote download */
105 char remote_path[PATH_MAX] = {""};
106 sprintf (path, "%s/geography/dcw/%s%s", GMT->session.USERDIR, name, suffix);
107 if (access (path, R_OK) == 0) { /* Previously downloaded */
108 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "3. DCW: Read the Digital Chart of the World from %s\n", path);
109 return true; /* Found it here */
110 }
111 /* Must download it the first time */
112 if (GMT->current.setting.auto_download == GMT_NO_DOWNLOAD) {
113 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to download the Digital Chart of the World for GMT since GMT_DATA_UPDATE_INTERVAL is off\n");
114 return false;
115 }
116 sprintf (path, "%s/geography/dcw", GMT->session.USERDIR); /* Local directory destination */
117 if (access (path, R_OK) && gmt_mkdir (path)) { /* Must first create the directory */
118 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create GMT directory : %s\n", path);
119 return false;
120 }
121 sprintf (path, "%s/geography/dcw/%s%s", GMT->session.USERDIR, name, suffix); /* Final local path */
122 snprintf (remote_path, PATH_MAX, "%s/geography/dcw/%s%s", gmt_dataserver_url (GMT->parent), name, suffix); /* Unique remote path */
123 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "Downloading %s%s for the first time - be patient\n", name, suffix);
124 if (gmt_download_file (GMT, name, remote_path, path, true)) {
125 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to obtain remote file %s%s\n", name, suffix);
126 return false;
127 }
128 else { /* Successfully downloaded the first time */
129 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "3. DCW: Read the Digital Chart of the World from %s\n", path);
130 return true;
131 }
132 }
133 return (false);
134 }
135
gmtdcw_load_lists(struct GMT_CTRL * GMT,struct GMT_DCW_COUNTRY ** C,struct GMT_DCW_STATE ** S,struct GMT_DCW_COUNTRY_STATE ** CS,unsigned int dim[])136 GMT_LOCAL int gmtdcw_load_lists (struct GMT_CTRL *GMT, struct GMT_DCW_COUNTRY **C, struct GMT_DCW_STATE **S, struct GMT_DCW_COUNTRY_STATE **CS, unsigned int dim[]) {
137 /* Open and read list of countries and states and return via two struct and one char arrays plus dimensions in dim */
138 size_t n_alloc = 300;
139 unsigned int k, n;
140 char path[PATH_MAX] = {""}, line[BUFSIZ] = {""};
141 FILE *fp = NULL;
142 struct GMT_DCW_COUNTRY *Country = NULL;
143 struct GMT_DCW_STATE *State = NULL;
144 struct GMT_DCW_COUNTRY_STATE *Country_State = NULL;
145
146 if (!gmtdcw_get_path (GMT, "dcw-countries", ".txt", path)) return GMT_NOTSET;
147
148 /* Get countries first */
149 if ((fp = fopen (path, "r")) == NULL) {
150 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to open file %s [permission trouble?]\n", path);
151 return GMT_NOTSET;
152 }
153 Country = gmt_M_memory (GMT, NULL, n_alloc, struct GMT_DCW_COUNTRY);
154 k = 0;
155 while ( gmt_fgets (GMT, line, BUFSIZ, fp)) {
156 if (line[0] == '#') continue; /* Skip comments */
157 sscanf (line, "%s %s %[^\n]", Country[k].continent, Country[k].code, Country[k].name);
158 k++;
159 if (k == n_alloc) {
160 n_alloc += 100;
161 Country = gmt_M_memory (GMT, Country, n_alloc, struct GMT_DCW_COUNTRY);
162 }
163 }
164 fclose (fp);
165 dim[0] = k; /* Number of countries */
166 Country = gmt_M_memory (GMT, Country, k, struct GMT_DCW_COUNTRY);
167
168 /* Get states */
169 if (!gmtdcw_get_path (GMT, "dcw-states", ".txt", path)) {
170 gmt_M_free (GMT, Country);
171 return GMT_NOTSET;
172 }
173 if ((fp = fopen (path, "r")) == NULL) {
174 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to open file %s [permission trouble?]\n", path);
175 gmt_M_free (GMT, Country);
176 return GMT_NOTSET;
177 }
178 State = gmt_M_memory (GMT, NULL, n_alloc, struct GMT_DCW_STATE);
179 k = 0; n = 1;
180 while ( gmt_fgets (GMT, line, BUFSIZ, fp)) {
181 if (line[0] == '#') continue; /* Skip comments */
182 sscanf (line, "%s %s %[^\n]", State[k].country, State[k].code, State[k].name);
183 if (k && strcmp (State[k].country, State[k-1].country)) n++; /* New country with states */
184 k++;
185 if (k == n_alloc) {
186 n_alloc += 100;
187 State = gmt_M_memory (GMT, State, n_alloc, struct GMT_DCW_STATE);
188 }
189 }
190 fclose (fp);
191 dim[1] = k; /* Number of states */
192 State = gmt_M_memory (GMT, State, k, struct GMT_DCW_STATE);
193
194 /* Get list of countries with states */
195
196 dim[2] = n; /* Number of countries with states */
197 if (CS) { /* Wants list returned */
198 Country_State = gmt_M_memory (GMT, NULL, n, struct GMT_DCW_COUNTRY_STATE);
199 gmt_M_memcpy (Country_State[0].country, State[0].country, 4, char);
200 for (k = n = 1; k < dim[1]; k++) {
201 if (strcmp (State[k].country, State[k-1].country)) gmt_M_memcpy (Country_State[n++].country, State[k].country, 4, char);
202 }
203 *CS = Country_State;
204 }
205
206 *C = Country;
207 *S = State;
208
209 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "DCW: Found %u countries, %u countries with states, and %u states\n", dim[0], dim[2], dim[1]);
210 return 0;
211 }
212
gmtdcw_comp_countries(const void * p1,const void * p2)213 GMT_LOCAL int gmtdcw_comp_countries (const void *p1, const void *p2) {
214 /* Used to sort countries alphabetically */
215 struct GMT_DCW_COUNTRY *A = (struct GMT_DCW_COUNTRY *)p1;
216 struct GMT_DCW_COUNTRY *B = (struct GMT_DCW_COUNTRY *)p2;
217 return (strcmp (A->code, B->code));
218 }
219
gmtdcw_find_country(char * code,struct GMT_DCW_COUNTRY * list,int n)220 GMT_LOCAL int gmtdcw_find_country (char *code, struct GMT_DCW_COUNTRY *list, int n) {
221 /* Basic binary search for country with given code and an alphabetically sorted list */
222 int low = 0, high = n, mid, last = GMT_NOTSET, way;
223
224 while (low < high) {
225 mid = (low + high) / 2;
226 if (mid == last) return (GMT_NOTSET); /* No such code */
227 way = strcmp (code, list[mid].code);
228 if (way > 0) low = mid;
229 else if (way < 0) high = mid;
230 else return (mid);
231 last = mid;
232 }
233 return (low);
234 }
235
gmtdcw_find_state(struct GMT_CTRL * GMT,char * scode,char * ccode,struct GMT_DCW_STATE * slist,int ns,bool check)236 GMT_LOCAL int gmtdcw_find_state (struct GMT_CTRL *GMT, char *scode, char *ccode, struct GMT_DCW_STATE *slist, int ns, bool check) {
237 /* Return state id given country and state codes using a linear search */
238 int i;
239 static struct GMT_DCW_CHINA_CODES gmtdcw_CN_codes[DCW_N_CHINA_PROVINCES] = {
240 {11, "BJ"}, {12, "TJ"}, {13, "HE"}, {14, "SX"}, {15, "NM"}, {21, "LN"}, {22, "JL"}, {23, "HL"},
241 {31, "SH"}, {32, "JS"}, {33, "ZJ"}, {34, "AH"}, {35, "FJ"}, {36, "JX"}, {37, "SD"}, {41, "HA"},
242 {42, "HB"}, {43, "HN"}, {44, "GD"}, {45, "GX"}, {46, "HI"}, {50, "CQ"}, {51, "SC"}, {52, "GZ"},
243 {53, "YN"}, {54, "XZ"}, {61, "SN"}, {62, "GS"}, {63, "QH"}, {64, "NX"}, {65, "XJ"}, {71, "TW"},
244 {91, "HK"}, {92, "MO"}
245 };
246 if (check && !strncmp (ccode, "CN", 2U) && isdigit (scode[0])) { /* Must switch to 2-character code */
247 unsigned int k = 0, id = atoi (scode);
248 while (k < DCW_N_CHINA_PROVINCES && id > gmtdcw_CN_codes[k].id) k++;
249 if (k == DCW_N_CHINA_PROVINCES) return (GMT_NOTSET); /* No such integer ID found in the list */
250 if (id < gmtdcw_CN_codes[k].id) return (GMT_NOTSET); /* No such integer ID found in the list */
251 GMT_Report (GMT->parent, GMT_MSG_NOTICE, "FYI, Chinese province code %d is deprecated. Use %s instead\n", id, gmtdcw_CN_codes[k].code);
252 scode = gmtdcw_CN_codes[k].code;
253 }
254 for (i = 0; i < ns; i++) if (!strcmp (scode, slist[i].code) && !strcmp (ccode, slist[i].country)) return (i);
255 return (GMT_NOTSET);
256 }
257
gmtdcw_country_has_states(char * code,struct GMT_DCW_COUNTRY_STATE * st_country,unsigned int n)258 GMT_LOCAL bool gmtdcw_country_has_states (char *code, struct GMT_DCW_COUNTRY_STATE *st_country, unsigned int n) {
259 /* Return true if this country has interior state boundaries */
260 unsigned int i;
261 for (i = 0; i < n; i++) if (!strcmp (code, st_country[i].country)) return (true);
262 return (false);
263 }
264
265 /*----------------------------------------------------------|
266 * Public functions that are part of the GMT Devel library |
267 *----------------------------------------------------------|
268 */
269
270 #define GMT_DCW_PLOTTING 1
271 #define GMT_DCW_CLIPPING 2
272
gmt_DCW_operation(struct GMT_CTRL * GMT,struct GMT_DCW_SELECT * F,double wesn[],unsigned int mode)273 struct GMT_DATASET * gmt_DCW_operation (struct GMT_CTRL *GMT, struct GMT_DCW_SELECT *F, double wesn[], unsigned int mode) {
274 /* Given comma-separated names, read the corresponding netCDF variables.
275 * mode = GMT_DCW_REGION : Return the joint w/e/s/n limits
276 * mode = GMT_DCW_PLOT : Plot the polygons [This is actually same as GMT_DCW_EXTRACT internally but plots instead of returning]
277 * mode = GMT_DCW_DUMP : Dump the polygons
278 * mode = GMT_DCW_EXTRACT : Return a dataset structure
279 */
280 int item, ks, retval, ncid, xvarid, yvarid, id;
281 int64_t first, last;
282 size_t np, max_np = 0U, n_alloc;
283 uint64_t k, seg, n_segments;
284 unsigned int n_items = 0, r_item = 0, pos = 0, kk, tbl = 0, j = 0, *order = NULL;
285 unsigned short int *dx = NULL, *dy = NULL;
286 unsigned int GMT_DCW_COUNTRIES = 0, GMT_DCW_STATES = 0, n_bodies[3] = {0, 0, 0}, special = 0;
287 bool done, want_state, outline, fill = false, is_Antarctica = false, hole, new_CN_codes = false;
288 char TAG[GMT_LEN16] = {""}, dim[GMT_LEN16] = {""}, xname[GMT_LEN16] = {""};
289 char yname[GMT_LEN16] = {""}, code[GMT_LEN16] = {""}, state[GMT_LEN16] = {""};
290 char msg[GMT_BUFSIZ] = {""}, path[PATH_MAX] = {""}, list[GMT_BUFSIZ] = {""};
291 char version[GMT_LEN32] = {""}, gmtversion[GMT_LEN32] = {""}, source[GMT_LEN256] = {""}, title[GMT_LEN256] = {""};
292 char label[GMT_LEN256] = {""}, header[GMT_LEN256] = {""};
293 double west, east, south, north, xscl, yscl, out[2], *lon = NULL, *lat = NULL;
294 struct GMT_RANGE *Z = NULL;
295 struct GMT_DATASET *D = NULL;
296 struct GMT_DATASEGMENT *P = NULL, *S = NULL;
297 struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
298 struct GMT_RECORD *Out = NULL;
299 struct GMT_DCW_COUNTRY *GMT_DCW_country = NULL;
300 struct GMT_DCW_STATE *GMT_DCW_state = NULL;
301 struct GMT_FILL *sfill = NULL;
302 struct GMT_PEN *spen = NULL;
303
304 for (j = ks = 0; j < F->n_items; j++) {
305 if (!F->item[j]->codes || F->item[j]->codes[0] == '\0') continue;
306 ks++; /* Gave some codes */
307 }
308 if (ks == 0) return NULL; /* No countries requested */
309 if (mode != GMT_DCW_REGION && F->region && (mode & (GMT_DCW_PLOT+GMT_DCW_DUMP+GMT_DCW_EXTRACT)) == 0) return NULL; /* No plotting/dumping requested, just -R */
310
311 if (mode & GMT_DCW_REGION) { /* Wish to determine region from polygons */
312 if (wesn == NULL) {
313 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Implementation error: Must pass wesn array if mode == %d\n", GMT_DCW_REGION);
314 return NULL;
315 }
316 wesn[XLO] = wesn[XHI] = 0.0; /* Set to zero so it can grow below */
317 wesn[YLO] = +9999.0; wesn[YHI] = -9999.0; /* Initialize so we can shrink it below */
318 }
319
320 if (gmtdcw_load_lists (GMT, &GMT_DCW_country, &GMT_DCW_state, NULL, n_bodies)) /* Something went wrong */
321 return NULL;
322
323 GMT_DCW_COUNTRIES = n_bodies[0];
324 GMT_DCW_STATES = n_bodies[1];
325
326 qsort ((void *)GMT_DCW_country, (size_t)GMT_DCW_COUNTRIES, sizeof (struct GMT_DCW_COUNTRY), gmtdcw_comp_countries); /* Sort on country code */
327
328 n_alloc = n_bodies[0] + n_bodies[1]; /* Presumably max items considered */
329 order = gmt_M_memory (GMT, NULL, n_alloc, unsigned int);
330 for (j = 0; j < F->n_items; j++) {
331 pos = 0;
332 while (gmt_strtok (F->item[j]->codes, ",", &pos, code)) { /* Loop over items */
333 if (code[0] == '=') { /* Must expand a continent into all member countries */
334 for (k = 0; k < GMT_DCW_COUNTRIES; k++) {
335 if (strncmp (GMT_DCW_country[k].continent, &code[1], 2)) continue; /* Not this one */
336 if (n_items) strcat (list, ",");
337 strcat (list, GMT_DCW_country[k].code);
338 order[n_items] = j; /* So we know which color/pen to apply for this item */
339 n_items++;
340 }
341 if (n_items)
342 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Continent code expanded from %s to %s [%d countries]\n", F->item[j]->codes, list, n_items);
343 else
344 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Continent code %s unrecognized\n", code);
345 }
346 else { /* Just append this single one */
347 if (n_items) strcat (list, ",");
348 strcat (list, code);
349 order[n_items] = j; /* So we know which color/pen to apply for this item */
350 n_items++;
351 }
352 }
353 }
354
355 if (n_items)
356 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Requested %d DCW items: %s\n", n_items, list);
357 else {
358 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "No countries selected\n");
359 gmt_M_free (GMT, order);
360 gmt_M_free (GMT, GMT_DCW_country);
361 gmt_M_free (GMT, GMT_DCW_state);
362 return NULL;
363 }
364
365 if (mode & GMT_DCW_PLOT) { /* Plot via psxy instead */
366 /* Because holes in polygons comes last we cannot just plot as we go. Instead, we must assemble
367 * the entire list of polygons for one item, then pass that dataset to psxy for plotting.
368 * So here, that means switch to GMT_DCW_EXTRACT but set a special flag so that we call psxy
369 * and then delete the dataset instead of returning it. */
370 mode -= GMT_DCW_PLOT;
371 mode += GMT_DCW_EXTRACT;
372 special = GMT_DCW_PLOTTING;
373 }
374 else if (mode & (GMT_DCW_CLIP_IN|GMT_DCW_CLIP_OUT)) { /* Lay down clip path via clip instead */
375 /* Because holes in polygons comes last we cannot just set clip path as we go. Instead, we must assemble
376 * the entire list of polygons for one item, then pass that dataset to clip for clipping.
377 * So here, that means switch to GMT_DCW_EXTRACT but set a special flag so that we call clip
378 * and then delete the dataset instead of returning it. */
379 mode -= (mode & GMT_DCW_CLIP_IN) ? GMT_DCW_CLIP_IN : GMT_DCW_CLIP_OUT;
380 mode += GMT_DCW_EXTRACT;
381 special = GMT_DCW_CLIPPING;
382 }
383
384 if (!gmtdcw_get_path (GMT, "dcw-gmt", ".nc", path)) {
385 gmt_M_free (GMT, order);
386 return NULL;
387 }
388
389 if (mode > GMT_DCW_REGION) { /* Wish to get actual polygons */
390 P = GMT_Alloc_Segment (GMT->parent, GMT_NO_STRINGS, 0, 2, NULL, NULL);
391 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Extract polygons from DCW - The Digital Chart of the World\n");
392 }
393
394 if (mode & GMT_DCW_EXTRACT) { /* Plan to return a dataset */
395 uint64_t dim[4] = {n_items, 0, 0, 2}; /* n_items tables whose records (to be allocated) will have 2 columns */
396 if ((D = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_POLY, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL) {
397 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to create empty dataset for DCW polygons\n");
398 gmt_free_segment (GMT, &P);
399 gmt_M_free (GMT, order);
400 return NULL;
401 }
402 }
403
404 if ((retval = gmt_nc_open (GMT, path, NC_NOWRITE, &ncid))) {
405 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot open file %s!\n", path);
406 gmt_free_segment (GMT, &P);
407 gmt_M_free (GMT, order);
408 return NULL;
409 }
410
411 /* Get global attributes */
412 if ((retval = nc_get_att_text (ncid, NC_GLOBAL, "version", version))) {
413 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot obtain DCW attribute version\n");
414 gmt_free_segment (GMT, &P);
415 gmt_M_free (GMT, order);
416 return NULL;
417 }
418 if (version[0] != '1') new_CN_codes = true; /* DCW 2.0.0 or later has new Chinese province codes */
419 if ((retval = nc_get_att_text (ncid, NC_GLOBAL, "title", title))) {
420 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot obtain DCW attribute title\n");
421 gmt_free_segment (GMT, &P);
422 gmt_M_free (GMT, order);
423 return NULL;
424 }
425 if ((retval = nc_get_att_text (ncid, NC_GLOBAL, "source", source))) {
426 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Cannot obtain DCW attribute source\n");
427 gmt_free_segment (GMT, &P);
428 gmt_M_free (GMT, order);
429 return NULL;
430 }
431 if ((retval = nc_get_att_text (ncid, NC_GLOBAL, "gmtversion", gmtversion)) == NC_NOERR)
432 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Found gmtversion string in DCW file: %s\n", gmtversion);
433
434 if (gmt_M_is_verbose (GMT, GMT_MSG_INFORMATION)) {
435 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Using country and state data from dcw-gmt\n");
436 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Title : %s\n", title);
437 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Source : %s\n", source);
438 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Version: %s\n", version);
439 if (gmtversion[0]) GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "DCW version %s requires GMT version %s or later.\n", version, gmtversion);
440 }
441
442 if (gmtversion[0]) { /* The gmtversion attribute was available [starting with DCW 2.0.0] */
443 int maj, min, rel;
444 if (sscanf (gmtversion, "%d.%d.%d", &maj, &min, &rel) != 3) {
445 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Unable to parse minimum GMT version information\n");
446 gmt_free_segment (GMT, &P);
447 gmt_M_free (GMT, order);
448 return NULL;
449 }
450 if (maj > GMT_MAJOR_VERSION || (maj == GMT_MAJOR_VERSION && min > GMT_MINOR_VERSION) || (maj == GMT_MAJOR_VERSION && min == GMT_MINOR_VERSION && rel > GMT_RELEASE_VERSION)) {
451 GMT_Report (GMT->parent, GMT_MSG_ERROR, "This DCW version (%s) requires at least GMT %s; you have %d.%d.%d\n", version, gmtversion, GMT_MAJOR_VERSION, GMT_MINOR_VERSION, GMT_RELEASE_VERSION);
452 gmt_free_segment (GMT, &P);
453 gmt_M_free (GMT, order);
454 return NULL;
455 }
456 }
457
458 if ((mode & GMT_DCW_DUMP) || (mode & GMT_DCW_REGION)) { /* Dump the coordinates to stdout or return -R means setting col types */
459 gmt_set_geographic (GMT, GMT_OUT);
460 }
461
462 if (mode & GMT_DCW_REGION) /* Just update wesn */
463 Z = gmt_M_memory (GMT, NULL, n_items, struct GMT_RANGE);
464
465 pos = item = 0;
466 Out = gmt_new_record (GMT, out, NULL); /* Since we only need to worry about numerics in this module */
467 while (gmt_strtok (list, ",", &pos, code)) { /* Loop over countries */
468 want_state = false;
469 if (code[2] == '.') { /* Requesting a state */
470 gmt_M_memset (state, GMT_LEN16, char);
471 strncpy (state, &code[3], GMT_LEN16-1);
472 code[2] = '\0';
473 want_state = true;
474 }
475 ks = gmtdcw_find_country (code, GMT_DCW_country, GMT_DCW_COUNTRIES);
476 if (ks == GMT_NOTSET) {
477 GMT_Report (GMT->parent, GMT_MSG_WARNING, "No country code matching %s (skipped)\n", code);
478 continue;
479 }
480 k = ks;
481 if (want_state) {
482 item = gmtdcw_find_state (GMT, state, code, GMT_DCW_state, GMT_DCW_STATES, new_CN_codes);
483 if (item == GMT_NOTSET) {
484 GMT_Report (GMT->parent, GMT_MSG_WARNING, "Country %s does not have a state named %s (skipped)\n", code, state);
485 continue;
486 }
487 snprintf (TAG, GMT_LEN16, "%s%s", GMT_DCW_country[k].code, GMT_DCW_state[item].code);
488 if (F->mode & GMT_DCW_ZHEADER)
489 snprintf (msg, GMT_BUFSIZ, "-Z%s %s (%s)\n", TAG, GMT_DCW_state[item].name, GMT_DCW_country[k].name);
490 else
491 snprintf (msg, GMT_BUFSIZ, "%s (%s)\n", GMT_DCW_state[item].name, GMT_DCW_country[k].name);
492 }
493 else {
494 snprintf (TAG, GMT_LEN16, "%s", GMT_DCW_country[k].code);
495 if (F->mode & GMT_DCW_ZHEADER)
496 snprintf (msg, GMT_BUFSIZ, "-Z%s %s\n", TAG, GMT_DCW_country[k].name);
497 else
498 snprintf (msg, GMT_BUFSIZ, "%s\n", GMT_DCW_country[k].name);
499 }
500 if (!strncmp (GMT_DCW_country[k].code, "AQ", 2U)) is_Antarctica = true;
501
502 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, msg);
503 k = strlen (msg) - 1;
504 msg[k] = '\0'; /* Remove the newline for use as segment header */
505
506 /* Open and read the netCDF file */
507
508 snprintf (dim, GMT_LEN16, "%s_length", TAG);
509 if ((retval = nc_inq_dimid (ncid, dim, &id))) {
510 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failure while getting ID for variable %s in %s!\n", dim, path);
511 continue;
512 }
513 if ((retval = nc_inq_dimlen (ncid, id, &np))) {
514 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failure while getting dimension length for variable %s in %s!\n", dim, path);
515 continue;
516 }
517 if (mode > GMT_DCW_REGION && np > max_np) {
518 size_t tmp_size = max_np;
519 gmt_M_malloc2 (GMT, lon, lat, np, &tmp_size, double);
520 gmt_M_malloc2 (GMT, dx, dy, np, &max_np, unsigned short int);
521 if (lon == NULL || lat == NULL || dx == NULL|| dy == NULL) {
522 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Failure while allocating memory!\n");
523 continue;
524 }
525 }
526
527 /* Get the varid of the lon and lat variables, based on their names, and get the data */
528
529 snprintf (xname, GMT_LEN16, "%s_lon", TAG); snprintf (yname, GMT_LEN16, "%s_lat", TAG);
530
531 if ((retval = nc_inq_varid (ncid, xname, &xvarid))) continue;
532 if ((retval = nc_get_att_double (ncid, xvarid, "min", &west))) continue;
533 if ((retval = nc_get_att_double (ncid, xvarid, "max", &east))) continue;
534 if ((retval = nc_get_att_double (ncid, xvarid, "scale", &xscl))) continue;
535 if ((retval = nc_inq_varid (ncid, yname, &yvarid))) continue;
536 if ((retval = nc_get_att_double (ncid, yvarid, "min", &south))) continue;
537 if ((retval = nc_get_att_double (ncid, yvarid, "max", &north))) continue;
538 if ((retval = nc_get_att_double (ncid, yvarid, "scale", &yscl))) continue;
539 if (mode & GMT_DCW_REGION) { /* Just update wesn */
540 Z[r_item].west = west; Z[r_item++].east = east;
541 if (south < wesn[YLO]) wesn[YLO] = south;
542 if (north > wesn[YHI]) wesn[YHI] = north;
543 }
544 if (mode > GMT_DCW_REGION) { /* Need to read the data arrays */
545 if ((retval = nc_get_var_ushort (ncid, xvarid, dx))) continue;
546 if ((retval = nc_get_var_ushort (ncid, yvarid, dy))) continue;
547 }
548 if (mode == GMT_DCW_REGION) continue;
549 xscl = 1.0 / xscl; yscl = 1.0 / yscl;
550 for (k = n_segments = 0; k < np; k++) { /* Unpack */
551 if (dx[k] == 65535U) { /* Start of new segment */
552 n_segments++; /* Count how many segments */
553 lon[k] = GMT->session.d_NaN; /* Flag a segment with lon = NaN */
554 lat[k] = (dy[k] == 1) ? 1.0 : 0.0; /* This is always 0.0 for 1.1.4 and older, which had no holes anyway */
555 }
556 else {
557 lon[k] = dx[k] * xscl + west;
558 lat[k] = dy[k] * yscl + south;
559 }
560 }
561 if (mode & GMT_DCW_EXTRACT) { /* Allocate a table with the right number of segments */
562 gmt_free_table (GMT, D->table[tbl]);
563 D->table[tbl] = gmt_create_table (GMT, n_segments, 0, 2, 0, false);
564 }
565 if (special == GMT_DCW_PLOTTING) { /* Time to consider fill/pen change */
566 outline = (F->item[order[tbl]]->mode & DCW_DO_OUTLINE);
567 fill = (F->item[order[tbl]]->mode & DCW_DO_FILL);
568 spen = (outline) ? &(F->item[order[tbl]]->pen) : NULL;
569 sfill = (fill) ? &(F->item[order[tbl]]->fill) : NULL;
570 }
571
572 /* Extract the pieces into separate segments */
573 k = seg = 0;
574 done = false;
575 while (!done) {
576 first = GMT_NOTSET;
577 while (first == GMT_NOTSET && k < np) { /* Look for next start of segment marker */
578 if (gmt_M_is_dnan (lon[k])) {
579 hole = (lat[k] > 0.0);
580 first = k + 1; /* Start of segment */
581 }
582 k++;
583 }
584 if (first == GMT_NOTSET) { done = true; continue;} /* No more segments */
585 last = GMT_NOTSET;
586 while (last == GMT_NOTSET && k < np) {/* Look for next end of segment marker (or end of line) */
587 if (gmt_M_is_dnan (lon[k])) last = k - 1; /* End of segment */
588 k++;
589 }
590 if (last == GMT_NOTSET) last = np - 1; /* End of last segment */
591 k--; /* Back to last segment marker which will be the next start marker */
592 P->n_rows = last - first + 1; /* Number of points in this segment */
593 P->data[GMT_X] = &lon[first];
594 P->data[GMT_Y] = &lat[first];
595 sprintf (label, " %s Segment %" PRIu64, msg, seg);
596 header[0] = '\0';
597 if (hole)
598 strcat (header, "-Ph");
599
600 if (mode & GMT_DCW_DUMP) { /* Dump the coordinates to stdout */
601 strcat (header, label);
602 strcpy (GMT->current.io.segment_header, header);
603 GMT_Put_Record (GMT->parent, GMT_WRITE_SEGMENT_HEADER, NULL);
604 for (kk = 0; kk < P->n_rows; kk++) {
605 out[GMT_X] = P->data[GMT_X][kk];
606 out[GMT_Y] = P->data[GMT_Y][kk];
607 GMT_Put_Record (GMT->parent, GMT_WRITE_DATA, Out);
608 }
609 seg++;
610 }
611 else if (mode & GMT_DCW_EXTRACT) { /* Attach to dataset */
612 S = D->table[tbl]->segment[seg];
613 SH = gmt_get_DS_hidden (S);
614 if (special == GMT_DCW_PLOTTING) {
615 if (sfill) {
616 strcat (header, " -G"); strcat (header, gmtlib_putfill (GMT, sfill));
617 }
618 else
619 strcat (header, " -G-");
620 if (spen) {
621 strcat (header, " -W"); strcat (header, gmt_putpen (GMT, spen));
622 }
623 else
624 strcat (header, " -W-");
625 }
626 strcat (header, label);
627 S->header = strdup (header);
628 if (hole)
629 SH->pol_mode = GMT_IS_HOLE;
630 S->n_rows = P->n_rows;
631 gmt_M_malloc2 (GMT, S->data[GMT_X], S->data[GMT_Y], S->n_rows, NULL, double);
632 gmt_M_memcpy (S->data[GMT_X], P->data[GMT_X], S->n_rows, double);
633 gmt_M_memcpy (S->data[GMT_Y], P->data[GMT_Y], S->n_rows, double);
634 SH->alloc_mode[GMT_X] = SH->alloc_mode[GMT_Y] = GMT_ALLOC_INTERNALLY; /* Allocated in GMT */
635 seg++;
636 }
637 }
638 tbl++;
639 }
640
641 gmt_nc_close (GMT, ncid);
642 gmt_M_free (GMT, GMT_DCW_country);
643 gmt_M_free (GMT, GMT_DCW_state);
644 gmt_M_free (GMT, Out);
645
646 if (mode & GMT_DCW_REGION) {
647 gmt_find_range (GMT, Z, n_items, &wesn[XLO], &wesn[XHI]);
648 gmt_M_free (GMT, Z);
649 GMT->current.io.geo.range = GMT_IGNORE_RANGE; /* Override this setting explicitly */
650 gmt_extend_region (GMT, wesn, F->adjust, F->inc);
651 if (is_Antarctica) { /* Must override to include pole and full longitude range */
652 wesn[YLO] = -90.0; /* Since it is a South polar cap */
653 wesn[XLO] = 0.0;
654 wesn[XHI] = 360.0;
655 }
656 /* Do basic sanity checks */
657 if (wesn[YLO] < -90.0) wesn[YLO] = -90.0;
658 if (wesn[YHI] > +90.0) wesn[YHI] = +90.0;
659 if (fabs (wesn[XHI] - wesn[XLO]) > 360.0) {
660 wesn[XLO] = 0.0;
661 wesn[XHI] = 360.0;
662 }
663 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Region implied by DCW polygons is %g/%g/%g/%g\n", wesn[XLO], wesn[XHI], wesn[YLO], wesn[YHI]);
664 }
665 gmt_M_free (GMT, order);
666 if (mode > GMT_DCW_REGION) {
667 gmt_M_free (GMT, dx);
668 gmt_M_free (GMT, dy);
669 gmt_M_free (GMT, lon);
670 gmt_M_free (GMT, lat);
671 P->data[GMT_X] = P->data[GMT_Y] = NULL;
672 gmt_free_segment (GMT, &P);
673 }
674
675 if (D) gmt_set_dataset_minmax (GMT, D); /* Update stats */
676
677 if (special) { /* Plot via psxy or clip via psclip, then free dataset */
678 char cmd[GMT_BUFSIZ] = {""}, in_string[GMT_VF_LEN] = {""};
679 static char *module[2] = {"psxy", "psclip"};
680 /* Get a virtual file for the current DCW dataset */
681 if (GMT_Open_VirtualFile (GMT->parent, GMT_IS_DATASET, GMT_IS_POLY, GMT_IN|GMT_IS_REFERENCE, D, in_string) == GMT_NOTSET) {
682 return (NULL);
683 }
684 /* All pen and fill settings are passed via segment headers, so this part is common to both psxy and psclip: */
685 snprintf (cmd, GMT_BUFSIZ, "-R -J -O -K %s", in_string);
686 if (F->mode & GMT_DCW_CLIP_OUT) /* Set the outside clip flag */
687 strcat (cmd, " -N"); /* Select the outside clip flag */
688 strcat (cmd, " --GMT_HISTORY=readonly"); /* Ignore history on exit */
689 GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "Calling %s with args %s\n", module[special-1], cmd);
690 if (GMT_Call_Module (GMT->parent, module[special-1], GMT_MODULE_CMD, cmd) != GMT_OK) {
691 return (NULL);
692 }
693 /* Close the virtual file and destroy the resource */
694 GMT_Close_VirtualFile (GMT->parent, in_string);
695 GMT_Destroy_Data (GMT->parent, &D);
696 }
697
698 return (D);
699 }
700
gmt_DCW_list(struct GMT_CTRL * GMT,struct GMT_DCW_SELECT * F)701 unsigned int gmt_DCW_list (struct GMT_CTRL *GMT, struct GMT_DCW_SELECT *F) {
702 /* Write to stdout the available countries [and optionally states], then make calling program exit */
703 unsigned int list_mode, i, j, k, kk, GMT_DCW_COUNTRIES = 0, GMT_DCW_STATES = 0, GMT_DCW_N_COUNTRIES_WITH_STATES = 0, n_bodies[3] = {0, 0, 0};
704 bool search = false;
705 char string[GMT_LEN128] = {""};
706 struct GMT_DCW_COUNTRY *GMT_DCW_country = NULL;
707 struct GMT_DCW_STATE *GMT_DCW_state = NULL;
708 struct GMT_DCW_COUNTRY_STATE *GMT_DCW_country_with_state = NULL;
709 struct GMT_RECORD *Out = NULL;
710 struct GMTAPI_CTRL *API = GMT->parent;
711
712 list_mode = F->mode;
713 if ((list_mode & GMT_DCW_LIST) == 0) return 0;
714 if (gmtdcw_load_lists (GMT, &GMT_DCW_country, &GMT_DCW_state, &GMT_DCW_country_with_state, n_bodies)) return 0; /* Something went wrong */
715 GMT_DCW_COUNTRIES = n_bodies[0];
716 GMT_DCW_STATES = n_bodies[1];
717 GMT_DCW_N_COUNTRIES_WITH_STATES = n_bodies[2];
718 GMT_Report (API, GMT_MSG_INFORMATION, "List of ISO 3166-1 alpha-2 codes for DCW supported countries:\n\n");
719 for (k = 0; k < F->n_items; k++) {
720 if (!F->item[k]->codes || F->item[k]->codes[0] == '\0') continue;
721 search = true; /* Gave some codes */
722 }
723
724 /* Initialize rec-by-rec output */
725 if (GMT_Set_Columns (API, GMT_OUT, 0, GMT_COL_FIX) != GMT_NOERROR) {
726 return (API->error);
727 }
728 if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_TEXT, GMT_OUT, GMT_ADD_DEFAULT, 0, F->options) != GMT_NOERROR) { /* Establishes data output */
729 return (API->error);
730 }
731 if (GMT_Begin_IO (API, GMT_IS_DATASET, GMT_OUT, GMT_HEADER_ON) != GMT_NOERROR) { /* Enables data output and sets access mode */
732 return (API->error);
733 }
734 if (GMT_Set_Geometry (API, GMT_OUT, GMT_IS_TEXT) != GMT_NOERROR) { /* Sets output geometry */
735 return (API->error);
736 }
737
738 Out = gmt_new_record (GMT, NULL, string); /* Since we only need to worry about text in this module */
739
740 for (i = k = 0; i < GMT_DCW_COUNTRIES; i++) {
741 if (search) { /* Listed continent(s) */
742 bool found = false;
743 for (kk = 0; kk < F->n_items; kk++) {
744 if (F->item[kk]->codes[0] == '=') {
745 if (strstr (F->item[kk]->codes, GMT_DCW_country[i].continent))
746 found = true;
747 }
748 else if (strncmp (F->item[kk]->codes, GMT_DCW_country[i].code, 2U) == 0)
749 found = true;
750 }
751 if (!found) continue;
752 }
753 if (F->n_items == 0 && (i == 0 || strcmp (GMT_DCW_country[i].continent, GMT_DCW_country[i-1].continent)) ) {
754 sprintf (string, "%s [%s]", GMT_DCW_continents[k++], GMT_DCW_country[i].continent);
755 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
756 }
757 if ((list_mode & 2) == 0) {
758 sprintf (string, "%s\t%s", GMT_DCW_country[i].code, GMT_DCW_country[i].name);
759 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
760 }
761 if ((list_mode & 2) && gmtdcw_country_has_states (GMT_DCW_country[i].code, GMT_DCW_country_with_state, GMT_DCW_N_COUNTRIES_WITH_STATES)) {
762 for (j = 0; j < GMT_DCW_STATES; j++) {
763 if (!strcmp (GMT_DCW_country[i].code, GMT_DCW_state[j].country)) {
764 sprintf (string, "%s.%s\t%s", GMT_DCW_country[i].code, GMT_DCW_state[j].code, GMT_DCW_state[j].name);
765 GMT_Put_Record (API, GMT_WRITE_DATA, Out);
766 }
767 }
768 }
769 }
770 if (GMT_End_IO (API, GMT_OUT, 0) != GMT_NOERROR) { /* Disables further data output */
771 return (API->error);
772 }
773 gmt_M_free (GMT, Out);
774 gmt_M_free (GMT, GMT_DCW_country);
775 gmt_M_free (GMT, GMT_DCW_state);
776 gmt_M_free (GMT, GMT_DCW_country_with_state);
777 return ((list_mode & GMT_DCW_LIST));
778 }
779
gmt_DCW_option(struct GMTAPI_CTRL * API,char option,unsigned int plot)780 void gmt_DCW_option (struct GMTAPI_CTRL *API, char option, unsigned int plot) {
781 /* Show the usage */
782 char *action[2] = {"extract", "plot"};
783 char *action2[2] = {"extracting", "plotting"};
784 char *usage[2] = {"Extract clipping polygons", "Apply different fill or outlines"};
785 GMT_Usage (API, 1, "\n-%c%s", option, DCW_OPT);
786 GMT_Usage (API, -2, "%s for specified list of countries. "
787 "Based on closed polygons from the Digital Chart of the World (DCW). "
788 "Append comma-separated list of ISO 3166 codes for countries to %s, i.e., "
789 "<code1>,<code2>,... etc., using the 2-character country codes. "
790 "To select a state of a country (if available), append .state, e.g, US.TX for Texas. "
791 "To select a whole continent, use =AF|AN|AS|EU|OC|NA|SA as <code>. Some modifiers:", usage[plot], action[plot]);
792 if (plot == 1) {
793 GMT_Usage (API, 3, "+c Set clip paths for the inside area [none].");
794 GMT_Usage (API, 3, "+C Set clip paths for the outside area [none].");
795 GMT_Usage (API, 3, "+g Fill polygons using given <fill> to fill [none].");
796 }
797 GMT_Usage (API, 3, "+l Just list the countries and their codes [no %s takes place].", action2[plot]);
798 GMT_Usage (API, 3, "+L List states/territories for Argentina, Australia, Brazil, Canada, China, India, Russia and the US. "
799 "Select =<continent>+l|L to only list countries from that continent or <code>+L for that country(repeatable).");
800 if (plot == 1)
801 GMT_Usage (API, 3, "+p Draw outline using given <pen> [none].");
802 GMT_Usage (API, 3, "+z Add -Z<countrycode> to multisegment headers if extracting polygons.");
803 if (plot == 1) {
804 GMT_Usage (API, -2, "Note: One of +c|C|p|g must be specified to plot; if -M is in effect we just get the data. "
805 "Repeat -%c to give different groups of items their own pen/fill settings.", option);
806 }
807 }
808
gmt_DCW_parse(struct GMT_CTRL * GMT,char option,char * args,struct GMT_DCW_SELECT * F)809 unsigned int gmt_DCW_parse (struct GMT_CTRL *GMT, char option, char *args, struct GMT_DCW_SELECT *F) {
810 /* Parse the F option in pscoast */
811 unsigned int n_errors = 0, pos = 0;
812 char p[GMT_BUFSIZ] = {""}, *c = NULL, *a = NULL, *q = NULL;
813 struct GMT_DCW_ITEM *this_item = NULL;
814
815 if ((a = strchr (args, '+'))) a[0] = '\0'; /* Temporarily chop off modifiers */
816 this_item = gmt_M_memory (GMT, NULL, 1, struct GMT_DCW_ITEM); /* New item to fill */
817 this_item->codes = strdup (args);
818 if (a) a[0] = '+'; /* Reset modifiers */
819
820 /* If +g is used with patterns and +r<dpi> is appended then there is conflict with +r for the deprecated region modification.
821 * We avoid this by checking for this case and replacing +r with @r to avoid the strtok splitting off that modifier. */
822
823 if (a && (c = strchr (a, '+'))) { /* Handle modifiers */
824 if ((q = strstr (c, "+g")) && strchr ("Pp", q[2]) && strstr (&q[3], "+r")) { /* There is a +r<dpi> that follows a +g pattern modifier */
825 char *t = &q[3]; /* First character of pattern name or number */
826 while (t[0] != '+') t++; /* Wind to next modifier or run out of chars */
827 if (t[0] == '+' && t[1] == 'r') { /* Found a +r<value> */
828 char *r = t++; /* Now t is at the 'r' */
829 t++; /* Now t is at first char afterwards */
830 while (t[0] && isdigit (t[0])) t++; /* Wind pass all integers */
831 if (t[0] == '\0' || t[0] == '+') { /* The modifier could be +r<dpi> or +r<inc>, assume dpi */
832 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Option -%c: Ambiguous modifier +r<val>; could be dpi of the pattern or (a deprecated) region increment - choosing dpi.\n", option);
833 GMT_Report (GMT->parent, GMT_MSG_DEBUG, "If you meant the region modifier then place it before the +g pattern specification.\n", option);
834 r[0] = GMT_ASCII_US; /* Change +r<dpi> to ASCII31<dpi> to pass strtok splitting */
835 }
836 /* Else it is taken to be a deprecated region increment */
837 }
838 }
839 while ((gmt_strtok (c, "+", &pos, p))) {
840 switch (p[0]) {
841 /* Listings*/
842 case 'R': case 'e': case 'r': /* Get region from polygon(s) BB */
843 F->region = true;
844 n_errors += gmt_parse_region_extender (GMT, option, p, &(F->adjust), F->inc); /* Possibly extend the final region before reporting */
845 break;
846 case 'l': /* Country list */
847 F->mode = DCW_GET_COUNTRY;
848 F->mode |= GMT_DCW_LIST;
849 break;
850 case 'L': /* Country and state list */
851 F->mode = DCW_GET_COUNTRY_AND_STATE;
852 F->mode |= GMT_DCW_LIST;
853 break;
854 case 'c': /* Set up a clip path around the selections instead */
855 F->mode |= GMT_DCW_CLIP_IN;
856 break;
857 case 'C': /* Set up a clip path outside the selections instead */
858 F->mode |= GMT_DCW_CLIP_OUT;
859 break;
860 case 'p':
861 if (gmt_getpen (GMT, &p[1], &(this_item->pen))) { /* Error decoding pen */
862 gmt_pen_syntax (GMT, option, NULL, " ", NULL, 0);
863 n_errors++;
864 }
865 this_item->mode |= DCW_DO_OUTLINE;
866 F->mode |= GMT_DCW_PLOT;
867 break;
868 case 'g':
869 if ((q = strchr (p, GMT_ASCII_US))) q[0] = '+'; /* Restore +r<dpi> */
870 if (gmt_getfill (GMT, &p[1], &(this_item->fill))) {
871 gmt_fill_syntax (GMT, option, NULL, " ");
872 n_errors++;
873 }
874 this_item->mode |= DCW_DO_FILL;
875 F->mode |= GMT_DCW_PLOT;
876 break;
877 case 'z': /* Add country code to -Z<code> in segment header */
878 F->mode |= GMT_DCW_ZHEADER;
879 break;
880 default:
881 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Unrecognized modifier +%s.\n", option, p);
882 n_errors++;
883 break;
884 }
885 }
886 }
887 if ((F->mode & (GMT_DCW_PLOT|DCW_DO_FILL)) && (F->mode & (GMT_DCW_CLIP_IN|GMT_DCW_CLIP_OUT))) {
888 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Cannot mix clipping and plotting!\n", option);
889 n_errors++;
890 }
891 if ((F->mode & (GMT_DCW_CLIP_IN|GMT_DCW_CLIP_OUT)) && (F->mode & GMT_DCW_ZHEADER)) {
892 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: Cannot mix clipping and setting header codes!\n", option);
893 n_errors++;
894 }
895 if (this_item->codes[0] == '\0' && !(F->mode & (DCW_GET_COUNTRY+DCW_GET_COUNTRY_AND_STATE))) { /* Did not give +l or +L, and no codes */
896 GMT_Report (GMT->parent, GMT_MSG_ERROR, "Option -%c: No country codes given\n", option);
897 n_errors++;
898 }
899 F->item = gmt_M_memory (GMT, F->item, F->n_items+1, struct GMT_DCW_ITEM *); /* Add one more pointer space to the structure (NULL first time) */
900 F->item[F->n_items] = this_item;
901 F->n_items++;
902 return (n_errors);
903 }
904
gmt_DCW_free(struct GMT_CTRL * GMT,struct GMT_DCW_SELECT * F)905 void gmt_DCW_free (struct GMT_CTRL *GMT, struct GMT_DCW_SELECT *F) {
906 /* Free what we might have allocated during parsing */
907 unsigned int k;
908 if (F->n_items == 0) return; /* Nothing to free */
909 for (k = 0; k < F->n_items; k++) {
910 gmt_M_str_free (F->item[k]->codes);
911 gmt_M_free (GMT, F->item[k]);
912 }
913 gmt_M_free (GMT, F->item);
914 }
915