/*-------------------------------------------------------------------- * $Id: pscoast.c 12849 2014-02-04 00:18:26Z pwessel $ * * 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 *--------------------------------------------------------------------*/ /* * pscoast can draw shorelines, paint the landward and/or seaward side of * the shoreline, paint lakes, draw political borders, and rivers. Three * data sets are considered: shores, borders, and rivers. Each of these * data sets come in 5 different resolutions. The lower resolution files * are derived from the full resolution data using the Douglas-Peucker (DP) * line reduction algorithm. By giving a tolerance in km, the algorithm * will remove points that do not depart more than the tolerance from the * straight line that would connect the points if the point in question was * removed. The resolutions are: * * full - The complete World Vector Shoreline + CIA data base * high - DP reduced with tolerance = 0.2 km * intermediate - DP reduced with tolerance = 1 km * low - DP reduced with tolerance = 5 km * crude - DP reduced with tolerance = 25 km * * If selected, pscoast will open the desired binned shoreline file and fill * the landmasses and/or oceans/lakes as shaded/colored/textured polygons. The * shoreline may also be drawn. Political boundaries and rivers may be plotted, * too. If only land is filled, then the 'wet' areas remain transparent. * If only oceans are filled, then the 'dry' areas remain transparent. This * allows the user to overlay land or ocean without wiping out the plot. * For more information about the binned polygon file, see the GMT Technical * Reference and Cookbook. * Optionally, the user may choose to issue clip paths rather than paint the * polygons. That way one may clip subsequent images to be visible only * inside or outside the coastline. * * * Author: Paul Wessel * Date: 1-JAN-2010 * Version: 5 * */ #define THIS_MODULE_NAME "pscoast" #define THIS_MODULE_LIB "core" #define THIS_MODULE_PURPOSE "Plot continents, countries, shorelines, rivers, and borders on maps" #include "gmt_dev.h" #define GMT_PROG_OPTIONS "->BJKOPRUVXYbcptxy" GMT_OPT("EZ") #define LAKE 0 #define RIVER 1 /* Control structure for pscoast */ struct PSCOAST_CTRL { struct A { /* -A[//][+ag|i|s][+r|l][+p] */ bool active; struct GMT_SHORE_SELECT info; } A; struct C { /* -C */ bool active; struct GMT_FILL fill[2]; /* lake and riverlake fill */ } 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 F { /* -F */ bool active; struct GMT_DCW_SELECT info; } F; struct G { /* -G */ bool active; bool clip; struct GMT_FILL fill; } G; struct I { /* -I[/] */ bool active; unsigned int use[GSHHS_N_RLEVELS]; unsigned int list[GSHHS_N_RLEVELS]; unsigned int n_rlevels; struct GMT_PEN pen[GSHHS_N_RLEVELS]; } I; struct L { /* -L */ bool active; struct GMT_MAP_SCALE item; } L; struct M { /* -M */ bool active; } M; struct N { /* -N[/] */ bool active; unsigned int use[GSHHS_N_BLEVELS]; unsigned int list[GSHHS_N_BLEVELS]; unsigned int n_blevels; struct GMT_PEN pen[GSHHS_N_BLEVELS]; } N; struct Q { /* -Q */ bool active; } Q; struct S { /* -S */ bool active; bool clip; struct GMT_FILL fill; } S; struct T { /* -L */ bool active; struct GMT_MAP_ROSE item; } T; struct W { /* -W[/] */ bool active; bool use[GSHHS_MAX_LEVEL]; struct GMT_PEN pen[GSHHS_MAX_LEVEL]; } W; #ifdef DEBUG struct DBG { /* -+ */ bool active; int bin; } debug; #endif }; void *New_pscoast_Ctrl (struct GMT_CTRL *GMT) { /* Allocate and initialize a new control structure */ unsigned int k; struct PSCOAST_CTRL *C = GMT_memory (GMT, NULL, 1, struct PSCOAST_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 */ if (GMT->current.map.frame.paint) /* Default Ocean color = Frame background color */ C->S.fill = GMT->current.map.frame.fill; else GMT_init_fill (GMT, &C->S.fill, 1.0, 1.0, 1.0); /* Default Ocean color = white */ C->C.fill[LAKE] = C->C.fill[RIVER] = C->S.fill; /* Default Lake/Riverlake color = Ocean color */ GMT_init_fill (GMT, &C->G.fill, 0.0, 0.0, 0.0); /* Default Land color = black */ for (k = 0; k < GSHHS_N_RLEVELS; k++) C->I.pen[k] = GMT->current.setting.map_default_pen; /* Default river pens */ for (k = 0; k < GSHHS_N_BLEVELS; k++) C->N.pen[k] = GMT->current.setting.map_default_pen; /* Default border pens */ for (k = 0; k < GSHHS_MAX_LEVEL; k++) C->W.pen[k] = GMT->current.setting.map_default_pen; /* Default coastline pens */ GMT_memset (&C->L.item, 1, struct GMT_MAP_SCALE); GMT_memset (&C->T.item, 1, struct GMT_MAP_ROSE); return (C); } void Free_pscoast_Ctrl (struct GMT_CTRL *GMT, struct PSCOAST_CTRL *C) { /* Deallocate control structure */ if (!C) return; if (C->F.active && C->F.info.codes) free ((void *)C->F.info.codes); GMT_free (GMT, C); } int GMT_pscoast_usage (struct GMTAPI_CTRL *API, int level) { /* This displays the pscoast synopsis and optionally full usage information */ 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: pscoast %s %s [%s]\n", GMT_B_OPT, GMT_J_OPT, GMT_A_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [-C[/]]\n\t[-D][+] [%s] -G[]]\n", GMT_Rgeoz_OPT, DCW_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[-I[/]] [%s] [-K]\n\t[-L%s]\n", GMT_Jz_OPT, GMT_SCALE); GMT_Message (API, GMT_TIME_NONE, "\t[-M] [-N[/]] [-O] [-P] [-Q] [-S]\n"); GMT_Message (API, GMT_TIME_NONE, "\t[-T%s]\n\t[%s] [%s] [-W[/][]]\n", GMT_TROSE, GMT_U_OPT, GMT_V_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s] [%s]\n", GMT_X_OPT, GMT_Y_OPT, GMT_bo_OPT); GMT_Message (API, GMT_TIME_NONE, "\t[%s] [%s]\n\t[%s] [%s]\n", GMT_c_OPT, GMT_p_OPT, GMT_t_OPT, GMT_colon_OPT); #ifdef DEBUG GMT_Message (API, GMT_TIME_NONE, "\t[-+]\n"); #endif if (level == GMT_SYNOPSIS) return (EXIT_FAILURE); GMT_Option (API, "J-Z,G"); GMT_Message (API, GMT_TIME_NONE, "\n\tOPTIONS:\n"); GMT_Option (API, "A,B-"); GMT_fill_syntax (API->GMT, 'C', "Set separate color for lakes and riverlakes [Default is same as ocean]."); GMT_Message (API, GMT_TIME_NONE, "\t Alternatively, set custom fills below. Repeat the -C option as needed.\n"); GMT_Message (API, GMT_TIME_NONE, "\t l = Lakes.\n"); GMT_Message (API, GMT_TIME_NONE, "\t r = River-lakes.\n"); GMT_Message (API, GMT_TIME_NONE, "\t-D Choose one of the following resolutions:\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 busy plots 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_DCW_option (API, 'F', 1U); GMT_fill_syntax (API->GMT, 'G', "Paint or clip \"dry\" areas."); GMT_Message (API, GMT_TIME_NONE, "\t 6) c to issue clip paths for land areas.\n"); GMT_pen_syntax (API->GMT, 'I', "Draw rivers. Specify feature and optionally append pen [Default for all levels: %s]."); GMT_Message (API, GMT_TIME_NONE, "\t Choose from the features below. Repeat the -I option as many times as needed.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 0 = Double-lined rivers (river-lakes).\n"); GMT_Message (API, GMT_TIME_NONE, "\t 1 = Permanent major rivers.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 2 = Additional major rivers.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 3 = Additional rivers.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 4 = Minor rivers.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 5 = Intermittent rivers - major.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 6 = Intermittent rivers - additional.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 7 = Intermittent rivers - minor.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 8 = Major canals.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 9 = Minor canals.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 10 = Irrigation canals.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Alternatively, specify preselected river groups:\n"); GMT_Message (API, GMT_TIME_NONE, "\t a = All rivers and canals (0-10).\n"); GMT_Message (API, GMT_TIME_NONE, "\t A = All rivers and canals except river-lakes (1-10).\n"); GMT_Message (API, GMT_TIME_NONE, "\t r = All permanent rivers (0-4).\n"); GMT_Message (API, GMT_TIME_NONE, "\t R = All permanent rivers except river-lakes (1-4).\n"); GMT_Message (API, GMT_TIME_NONE, "\t i = All intermittent rivers (5-7).\n"); GMT_Message (API, GMT_TIME_NONE, "\t c = All canals (8-10).\n"); GMT_Option (API, "K"); GMT_mapscale_syntax (API->GMT, 'L', "Draws a simple map scale centered on /"); GMT_Message (API, GMT_TIME_NONE, "\t-M Dump a multisegment ASCII (or binary, see -bo) file to standard output. No plotting occurs.\n"); GMT_Message (API, GMT_TIME_NONE, "\t Specify one of -F, -I, -N, or -W.\n"); GMT_pen_syntax (API->GMT, 'N', "Draw boundaries. Specify feature and optionally append pen [Default for all levels: %s]."); GMT_Message (API, GMT_TIME_NONE, "\t Choose from the features below. Repeat the -N option as many times as needed.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 1 = National boundaries.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 2 = State boundaries within the Americas.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 3 = Marine boundaries.\n"); GMT_Message (API, GMT_TIME_NONE, "\t a = All boundaries (1-3).\n"); GMT_Option (API, "O,P"); GMT_Message (API, GMT_TIME_NONE, "\t-Q Terminate previously set clip-paths.\n"); GMT_fill_syntax (API->GMT, 'S', "Paint of clip \"wet\" areas."); GMT_Message (API, GMT_TIME_NONE, "\t 6) c to issue clip paths for water areas.\n"); GMT_maprose_syntax (API->GMT, 'T', "Draw a north-pointing map rose centered on /."); GMT_Option (API, "U,V"); GMT_pen_syntax (API->GMT, 'W', "Draw shorelines. Append pen [Default for all levels: %s]."); GMT_Message (API, GMT_TIME_NONE, "\t Alternatively, set custom pens below. Repeat the -W option as many times as needed.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 1 = Coastline.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 2 = Lake shores.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 3 = Island in lakes shores.\n"); GMT_Message (API, GMT_TIME_NONE, "\t 4 = Lake in island in lake shores.\n"); GMT_Message (API, GMT_TIME_NONE, "\t When feature-specific pens are used, those not set are deactivated.\n"); GMT_Option (API, "X,bo,c,p,t"); #ifdef DEBUG GMT_Message (API, GMT_TIME_NONE, "\t-+ Print only a single bin (debug option).\n"); #endif GMT_Option (API, "."); return (EXIT_FAILURE); } int GMT_pscoast_parse (struct GMT_CTRL *GMT, struct PSCOAST_CTRL *Ctrl, struct GMT_OPTION *options) { /* This parses the options provided to pscoast and sets parameters in Ctrl. * Note Ctrl has already been initialized and non-zero default values set. * 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, n_files = 0, k; int ks; bool clipping; struct GMT_OPTION *opt = NULL; struct GMTAPI_CTRL *API = GMT->parent; struct GMT_PEN pen; char *string = NULL; for (opt = options; opt; opt = opt->next) { /* Process all the options given */ switch (opt->option) { case '<': /* Input files */ n_files++; break; /* Processes program-specific parameters */ case 'A': Ctrl->A.active = true; GMT_set_levels (GMT, opt->arg, &Ctrl->A.info); break; case 'C': /* Lake colors */ Ctrl->C.active = true; if ((opt->arg[0] == 'l' || opt->arg[0] == 'r') && opt->arg[1] == '/') { /* Specific lake or river-lake fill */ k = (opt->arg[0] == 'l') ? LAKE : RIVER; if (opt->arg[2] && GMT_getfill (GMT, &opt->arg[2], &Ctrl->C.fill[k])) { GMT_fill_syntax (GMT, 'C', " "); n_errors++; } } else if (opt->arg[0]) { if (GMT_getfill (GMT, opt->arg, &Ctrl->C.fill[LAKE])) { GMT_fill_syntax (GMT, 'C', " "); n_errors++; } Ctrl->C.fill[RIVER] = Ctrl->C.fill[LAKE]; } break; case 'D': Ctrl->D.active = true; Ctrl->D.set = (opt->arg[0]) ? opt->arg[0] : 'l'; Ctrl->D.force = (opt->arg[1] == '+'); break; case 'F': Ctrl->F.active = true; n_errors += GMT_DCW_parse (GMT, opt->option, opt->arg, &Ctrl->F.info); break; case 'G': /* Set Gray shade, pattern, or clipping */ Ctrl->G.active = true; if (opt->arg[0] == 'c' && !opt->arg[1]) Ctrl->G.clip = true; else if (GMT_getfill (GMT, opt->arg, &Ctrl->G.fill)) { GMT_fill_syntax (GMT, 'G', " "); n_errors++; } break; case 'I': Ctrl->I.active = true; if (!opt->arg[0]) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -I option takes at least one argument\n"); n_errors++; continue; } pen = GMT->current.setting.map_default_pen; /* Set default pen */ if ((string = strchr (opt->arg, '/'))) { /* Get specified pen */ if (GMT_getpen (GMT, ++string, &pen)) { /* Error decoding pen */ GMT_pen_syntax (GMT, 'I', " "); n_errors++; } } switch (opt->arg[0]) { case 'a': for (k = 0; k < GSHHS_N_RLEVELS; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; case 'A': for (k = 1; k < GSHHS_N_RLEVELS; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; case 'r': for (k = 0; k < GSHHS_RIVER_INTERMITTENT; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; case 'R': for (k = 1; k < GSHHS_RIVER_INTERMITTENT; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; case 'i': for (k = GSHHS_RIVER_INTERMITTENT; k < GSHHS_RIVER_CANALS; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; case 'c': for (k = GSHHS_RIVER_CANALS; k < GSHHS_N_RLEVELS; k++) Ctrl->I.use[k] = 1, Ctrl->I.pen[k] = pen; break; default: ks = atoi (opt->arg); if (ks < 0 || ks >= GSHHS_N_RLEVELS) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -I option: Feature not in list!\n"); n_errors++; } else Ctrl->I.use[ks] = 1, Ctrl->I.pen[ks] = pen; break; } break; case 'L': Ctrl->L.active = true; n_errors += GMT_getscale (GMT, 'L', opt->arg, &Ctrl->L.item); break; case 'm': if (GMT_compat_check (GMT, 4)) /* Warn and fall through */ GMT_Report (API, GMT_MSG_COMPAT, "Warning: -m option is deprecated and reverted back to -M.\n"); else { n_errors += GMT_default_error (GMT, opt->option); break; } case 'M': Ctrl->M.active = true; break; case 'N': Ctrl->N.active = true; if (!opt->arg[0]) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error: -N option takes at least one argument\n"); n_errors++; continue; } pen = GMT->current.setting.map_default_pen; /* Set default pen */ if ((string = strchr (opt->arg, '/'))) { /* Get specified pen */ if (GMT_getpen (GMT, ++string, &pen)) { /* Error decoding pen */ GMT_pen_syntax (GMT, 'N', " "); n_errors++; } } switch (opt->arg[0]) { case 'a': for (k = 0; k < GSHHS_N_BLEVELS; k++) Ctrl->N.use[k] = 1, Ctrl->N.pen[k] = pen; break; default: ks = opt->arg[0] - '1'; if (ks < 0 || ks >= GSHHS_N_BLEVELS) { GMT_Report (API, GMT_MSG_NORMAL, "Syntax error -N option: Feature not in list!\n"); n_errors++; } else Ctrl->N.use[ks] = 1, Ctrl->N.pen[ks] = pen; break; } break; case 'Q': Ctrl->Q.active = true; break; case 'S': /* Set ocean color if needed */ Ctrl->S.active = true; if (opt->arg[0] == 'c' && opt->arg[1] == '\0') Ctrl->S.clip = true; else if (GMT_getfill (GMT, opt->arg, &Ctrl->S.fill)) { GMT_fill_syntax (GMT, 'S', " "); n_errors++; } break; case 'T': Ctrl->T.active = true; n_errors += GMT_getrose (GMT, 'T', opt->arg, &Ctrl->T.item); break; case 'W': Ctrl->W.active = true; /* Want to draw shorelines */ if ((opt->arg[0] >= '1' && opt->arg[0] <= '4') && opt->arg[1] == '/') { /* Specific pen for this feature */ k = (int)(opt->arg[0] - '1'); if (opt->arg[2] && GMT_getpen (GMT, &opt->arg[2], &Ctrl->W.pen[k])) { GMT_pen_syntax (GMT, 'W', " "); n_errors++; } Ctrl->W.use[k] = true; } else if (opt->arg[0]) { /* Same pen for all features */ Ctrl->W.use[0] = true; if (GMT_getpen (GMT, opt->arg, &Ctrl->W.pen[0])) { GMT_pen_syntax (GMT, 'W', " "); n_errors++; } else { for (k = 1; k < GSHHS_MAX_LEVEL; k++) Ctrl->W.pen[k] = Ctrl->W.pen[0], Ctrl->W.use[k] = true; } } else { /* Accept default pen for all features */ for (k = 0; k < GSHHS_MAX_LEVEL; k++) Ctrl->W.use[k] = true; } break; #ifdef DEBUG case '+': Ctrl->debug.active = true; Ctrl->debug.bin = atoi (opt->arg); break; #endif default: /* Report bad options */ n_errors += GMT_default_error (GMT, opt->option); break; } } if (GMT_DCW_list (GMT, Ctrl->F.info.mode)) return 1; if (Ctrl->C.active && !(Ctrl->G.active || Ctrl->S.active || Ctrl->W.active)) { /* Just lakes, fix -A */ if (Ctrl->A.info.low < 2) Ctrl->A.info.low = 2; } if (!GMT->common.R.active && Ctrl->F.active && Ctrl->M.active && !Ctrl->F.info.region) Ctrl->F.info.region = true; /* For -M and -F with no plotting, get -R from pols */ if (!GMT->common.R.active && Ctrl->F.info.region) { /* Must pick up region from chosen polygons */ (void) GMT_DCW_operation (GMT, &Ctrl->F.info, GMT->common.R.wesn, GMT_DCW_REGION); GMT->common.R.active = true; if (!GMT->common.J.active && !Ctrl->M.active) { /* No plotting or no dumping means just return the -R string */ char record[GMT_BUFSIZ] = {"-R"}, text[GMT_LEN64] = {""}; size_t i, j; GMT_ascii_format_col (GMT, text, GMT->common.R.wesn[XLO], GMT_OUT, GMT_X); strcat (record, text); strcat (record, "/"); GMT_ascii_format_col (GMT, text, GMT->common.R.wesn[XHI], GMT_OUT, GMT_X); strcat (record, text); strcat (record, "/"); GMT_ascii_format_col (GMT, text, GMT->common.R.wesn[YLO], GMT_OUT, GMT_Y); strcat (record, text); strcat (record, "/"); GMT_ascii_format_col (GMT, text, GMT->common.R.wesn[YHI], GMT_OUT, GMT_Y); strcat (record, text); /* Remove any white space due to selected formatting */ for (i = j = 2; i < strlen (record); i++) { if (record[i] == ' ') continue; /* Skip spaces */ record[j++] = record[i]; } record[j] = '\0'; if (GMT_Init_IO (API, GMT_IS_TEXTSET, GMT_IS_NONE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_OK) { /* Establishes data output */ return (API->error); } if (GMT_Begin_IO (API, GMT_IS_TEXTSET, GMT_OUT, GMT_HEADER_OFF) != GMT_OK) { return (API->error); } GMT_Put_Record (API, GMT_WRITE_TEXT, record); /* Write text record to output destination */ if (GMT_End_IO (API, GMT_OUT, 0) != GMT_OK) { /* Disables further data output */ return (API->error); } return 1; /* To exit, basically */ } } /* Check that the options selected are mutually consistent */ clipping = (Ctrl->G.clip || Ctrl->S.clip); if (Ctrl->M.active) { /* Need -R only */ n_errors += GMT_check_condition (GMT, !GMT->common.R.active, "Syntax error: Must specify -R option\n"); } else if (!Ctrl->Q.active) { /* Need -R -J */ n_errors += GMT_check_condition (GMT, !GMT->common.R.active, "Syntax error: Must specify -R option\n"); n_errors += GMT_check_condition (GMT, !GMT->common.J.active, "Syntax error: Must specify a map projection with the -J option\n"); } for (k = 0; k < GSHHS_MAX_LEVEL; k++) { n_errors += GMT_check_condition (GMT, Ctrl->W.pen[k].width < 0.0, "Syntax error -W option: Pen thickness for feature %d cannot be negative\n", k); } n_errors += GMT_check_condition (GMT, n_files > 0, "Syntax error: No input files allowed\n"); n_errors += GMT_check_condition (GMT, (Ctrl->G.active + Ctrl->S.active + Ctrl->C.active) > 1 && clipping, "Syntax error: Cannot combine -C, -G, -S while clipping\n"); n_errors += GMT_check_condition (GMT, Ctrl->G.clip && Ctrl->S.clip, "Syntax error: Must choose between clipping land OR water\n"); n_errors += GMT_check_condition (GMT, Ctrl->M.active && (Ctrl->G.active || Ctrl->S.active || Ctrl->C.active), "Syntax error: Must choose between dumping and clipping/plotting\n"); n_errors += GMT_check_condition (GMT, clipping && GMT->current.proj.projection == GMT_AZ_EQDIST && fabs (GMT->common.R.wesn[XLO] - GMT->common.R.wesn[XHI]) == 360.0 && (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO]) == 180.0, "-JE not implemented for global clipping - I quit\n"); n_errors += GMT_check_condition (GMT, clipping && (Ctrl->N.active || Ctrl->I.active || Ctrl->W.active), "Cannot do clipping AND draw coastlines, rivers, or borders\n"); n_errors += GMT_check_condition (GMT, !(Ctrl->G.active || Ctrl->S.active || Ctrl->C.active || Ctrl->F.active || Ctrl->W.active || Ctrl->N.active || Ctrl->I.active || Ctrl->Q.active), "Syntax error: Must specify at least one of -C, -G, -S, -I, -N, -Q and -W\n"); n_errors += GMT_check_condition (GMT, Ctrl->M.active && (Ctrl->F.active + Ctrl->N.active + Ctrl->I.active + Ctrl->W.active) != 1, "Syntax error -M: Must specify one of -F, -I, -N, and -W\n"); if (Ctrl->I.active) { /* Generate list of desired river features in sorted order */ for (k = Ctrl->I.n_rlevels = 0; k < GSHHS_N_RLEVELS; k++) { if (!Ctrl->I.use[k]) continue; n_errors += GMT_check_condition (GMT, Ctrl->I.pen[k].width < 0.0, "Syntax error -I option: Pen thickness cannot be negative\n"); Ctrl->I.list[Ctrl->I.n_rlevels++] = k; /* Since goes from 0-10 */ } } if (Ctrl->N.active) { /* Generate list of desired border features in sorted order */ for (k = Ctrl->N.n_blevels = 0; k < GSHHS_N_BLEVELS; k++) { if (!Ctrl->N.use[k]) continue; n_errors += GMT_check_condition (GMT, Ctrl->N.pen[k].width < 0.0, "Syntax error -N option: Pen thickness cannot be negative\n"); Ctrl->N.list[Ctrl->N.n_blevels++] = k + 1; /* Add one so we get range 1-3 */ } } return (n_errors); } bool add_this_polygon_to_path (struct GMT_CTRL *GMT, int k0, struct GMT_GSHHS_POL *p, int level, int k) { /* Determines if we should add the current polygon pol[k] to the growing path we are constructing */ int j; if (p[k].level != level) return (false); /* Sorry, this polygon does not have the correct level */ if (k0 == -1) return (true); /* The start level is always used */ /* Must make sure the polygon is inside the mother polygon. Because numerical noise can trick us we try up to two points */ if (GMT_non_zero_winding (GMT, p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) return (true); /* OK, we're in! */ /* First point was not inside. Test another point to see if the first failure was just an accident */ j = p[k].n / 2; /* We pick the second point from the other half of the polygon */ if (GMT_non_zero_winding (GMT, p[k].lon[j], p[k].lat[j], p[k0].lon, p[k0].lat, p[k0].n)) return (true); /* One of those rare occasions when we almost missed a polygon */ return (false); /* No, we do not want to use this polygon */ } void recursive_path (struct GMT_CTRL *GMT, struct PSL_CTRL *PSL, int k0, int np, struct GMT_GSHHS_POL p[], int level, struct GMT_FILL fill[]) { /* To avoid painting where no paint should go we assemble all the swiss cheese polygons by starting * with the lowest level (0) and looking for polygons inside polygons of the same level. We do this * using a recursive algorithm */ int k; if (level > GSHHS_MAX_LEVEL) return; for (k = k0 + 1; k < np; k++) { if (p[k].n == 0 || p[k].level < level) continue; if (add_this_polygon_to_path (GMT, k0, p, level, k)) { /* Add this to the current path */ PSL_comment (PSL, "Polygon %d\n", k); PSL_plotline (PSL, p[k].lon, p[k].lat, p[k].n, PSL_MOVE); #ifdef DEBUGX GMT_Message (API, GMT_TIME_NONE, "#\n"); for (i = 0; i < p[k].n; i++) GMT_Message (API, GMT_TIME_NONE, "%g\t%g\n", p[k].lon[i], p[k].lat[i]); #endif recursive_path (GMT, PSL, k, np, p, level + 1, fill); p[k].n = 0; /* Mark as used */ if (k0 == -1 && fill) { /* At start level: done nesting, time to paint the assembled swiss cheese polygon */ GMT_setfill (GMT, &fill[p[k].fid], false); PSL_command (PSL, "FO\n"); } } } } #define bailout(code) {GMT_Free_Options (mode); return (code);} #define Return(code) {Free_pscoast_Ctrl (GMT, Ctrl); GMT_end_module (GMT, GMT_cpy); bailout (code);} int GMT_pscoast (void *V_API, int mode, void *args) { /* High-level function that implements the pscoast task */ int i, np, ind, bin = 0, base, anti_bin = -1, np_new, k, last_k, err, bin_trouble, error, n; int level_to_be_painted = 0, direction, start_direction, stop_direction, last_pen_level; bool shift = false, need_coast_base, recursive; bool greenwich = false, possibly_donut_hell = false, fill_in_use = false; bool clobber_background, paint_polygons = false, donut; bool donut_hell = false, world_map_save, clipping; double bin_x[5], bin_y[5], out[2], *xtmp = NULL, *ytmp = NULL; double west_border, east_border, anti_lon = 0.0, anti_lat = -90.0, edge = 720.0; double left, right, anti_x = 0.0, anti_y = 0.0, x_0 = 0.0, y_0 = 0.0, x_c, y_c, dist; char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"}; struct GMT_FILL fill[6]; /* Colors for (0) water, (1) land, (2) lakes, (3) islands in lakes, (4) lakes in islands in lakes, (5) riverlakes */ struct GMT_SHORE c; struct GMT_BR b, r; struct GMT_GSHHS_POL *p = NULL; struct PSCOAST_CTRL *Ctrl = NULL; /* Control structure specific to program */ struct GMT_CTRL *GMT = NULL, *GMT_cpy = NULL; /* General GMT interal parameters */ struct GMT_OPTION *options = NULL; struct PSL_CTRL *PSL = NULL; /* General PSL interal parameters */ 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_pscoast_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_pscoast_usage (API, GMT_USAGE)); /* Return the usage message */ if (options->option == GMT_OPT_SYNOPSIS) bailout (GMT_pscoast_usage (API, GMT_SYNOPSIS)); /* Return the synopsis */ /* Parse the command-line arguments; return if errors are encountered */ 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_pscoast_Ctrl (GMT); /* Allocate and initialize defaults in a new control structure */ if ((error = GMT_pscoast_parse (GMT, Ctrl, options))) Return (error); /*---------------------------- This is the pscoast main code ----------------------------*/ /* Check and interpret the command line arguments */ clipping = (Ctrl->G.clip || Ctrl->S.clip); if (Ctrl->D.force) Ctrl->D.set = GMT_shore_adjust_res (GMT, Ctrl->D.set); base = GMT_set_resolution (GMT, &Ctrl->D.set, 'D'); fill[0] = Ctrl->S.fill; fill[1] = fill[3] = Ctrl->G.fill; fill[2] = fill[4] = (Ctrl->C.active) ? Ctrl->C.fill[LAKE] : Ctrl->S.fill; fill[5] = (Ctrl->C.active) ? Ctrl->C.fill[RIVER] : fill[2]; need_coast_base = (Ctrl->G.active || Ctrl->S.active || Ctrl->C.active || Ctrl->W.active); if (Ctrl->Q.active) need_coast_base = false; /* Since we just end clipping */ clobber_background = (Ctrl->G.active && Ctrl->S.active); recursive = (Ctrl->G.active != (Ctrl->S.active || Ctrl->C.active) || clipping); paint_polygons = (Ctrl->G.active || Ctrl->S.active || Ctrl->C.active); if (GMT->common.R.wesn[XLO] > 360.0) { GMT->common.R.wesn[XLO] -= 360.0; GMT->common.R.wesn[XHI] -= 360.0; } if (Ctrl->Q.active && !GMT->common.J.active) { /* Fake area and linear projection */ GMT_parse_common_options (GMT, "J", 'J', "x1d"); /* Fake linear projection */ GMT->common.R.wesn[XLO] = GMT->common.R.wesn[YLO] = 0.0; GMT->common.R.wesn[XHI] = GMT->common.R.wesn[YHI] = 1.0; } else if (Ctrl->M.active && !GMT->common.J.active) GMT_parse_common_options (GMT, "J", 'J', "x1d"); /* Fake linear projection */ if (GMT_err_pass (GMT, GMT_map_setup (GMT, GMT->common.R.wesn), "")) Return (GMT_RUNTIME_ERROR); if (!GMT_is_geographic (GMT, GMT_IN)) { GMT_Report (API, GMT_MSG_NORMAL, "Error: You must use a map projection or -Jx|X...d[/...d] for geographic data\n"); Return (EXIT_FAILURE); } world_map_save = GMT->current.map.is_world; if (need_coast_base && 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]); need_coast_base = false; } if (Ctrl->N.active && GMT_init_br (GMT, 'b', Ctrl->D.set, &b, GMT->common.R.wesn)) { GMT_Report (API, GMT_MSG_NORMAL, "%s resolution political boundary data base not installed\n", shore_resolution[base]); Ctrl->N.active = false; } if (need_coast_base) GMT_Report (API, GMT_MSG_VERBOSE, "GSHHG version %s\n%s\n%s\n", c.version, c.title, c.source); if (Ctrl->I.active && GMT_init_br (GMT, 'r', Ctrl->D.set, &r, GMT->common.R.wesn)) { GMT_Report (API, GMT_MSG_NORMAL, "%s resolution river data base not installed\n", shore_resolution[base]); Ctrl->I.active = false; } if (!(need_coast_base || Ctrl->F.active || Ctrl->N.active || Ctrl->I.active || Ctrl->Q.active)) { GMT_Report (API, GMT_MSG_NORMAL, "No databases available - aborts\n"); Return (EXIT_FAILURE); } if (Ctrl->M.active) { /* Dump linesegments to stdout; no plotting takes place */ int id = 0; char header[GMT_BUFSIZ] = {""}, *kind[3] = {"Coastlines", "Political boundaries", "Rivers"}, *version = NULL, *title = NULL, *source = NULL; GMT_set_geographic (GMT, GMT_OUT); /* Output lon/lat */ if (Ctrl->N.active) id = 1; if (Ctrl->I.active) id = 2; GMT_set_segmentheader (GMT, GMT_OUT, true); /* Turn on segment headers on output */ if ((error = GMT_set_cols (GMT, GMT_OUT, 2)) != GMT_OK) { Return (error); } if (GMT_Init_IO (API, GMT_IS_DATASET, GMT_IS_LINE, GMT_OUT, GMT_ADD_DEFAULT, 0, options) != GMT_OK) { /* Establishes data output */ 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); } if (Ctrl->W.active) { version = c.version; title = c.title; source = c.source; } else if (Ctrl->N.active) { version = b.version; title = b.title; source = b.source; } else { version = r.version; title = r.title; source = r.source; } if (Ctrl->F.active) { sprintf (header, "# Country polygons extracted from the DCW database\n"); GMT_Put_Record (API, GMT_WRITE_TEXT, header); } else { sprintf (header, "# %s extracted from the %s resolution GSHHG version %s database\n", kind[id], shore_resolution[base], version); GMT_Put_Record (API, GMT_WRITE_TEXT, header); sprintf (header, "# %s\n# %s\n", title, source); GMT_Put_Record (API, GMT_WRITE_TEXT, header); } } else { if (Ctrl->Q.active) GMT->current.ps.nclip = -1; /* Signal that this program terminates polygon clipping that initiated prior to this process */ else if (clipping) GMT->current.ps.nclip = +1; /* Signal that this program initiates new clipping that will outlive this process */ PSL = GMT_plotinit (GMT, options); if (Ctrl->Q.active) { /* Just undo previous clip-path */ PSL_endclipping (PSL, 1); GMT_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); GMT_plotcanvas (GMT); /* Fill canvas if requested */ GMT_map_basemap (GMT); /* Basemap needed */ GMT_plane_perspective (GMT, -1, 0.0); GMT_plotend (GMT); GMT_Report (API, GMT_MSG_VERBOSE, "Done!\n"); Return (GMT_OK); } GMT_plane_perspective (GMT, GMT->current.proj.z_project.view_plane, GMT->current.proj.z_level); GMT_plotcanvas (GMT); /* Fill canvas if requested */ } for (i = 0; i < 5; i++) if (fill[i].use_pattern) fill_in_use = true; if (fill_in_use && !clobber_background) { /* Force clobber when fill is used since our routine cannot deal with clipped fills */ GMT_Report (API, GMT_MSG_VERBOSE, "Warning: Pattern fill requires oceans to be painted first\n"); clobber_background = true; recursive = false; } if (clipping) GMT_map_basemap (GMT); if (GMT->current.proj.projection == GMT_AZ_EQDIST && GMT_360_RANGE (GMT->common.R.wesn[XLO], GMT->common.R.wesn[XHI]) && GMT_180_RANGE (GMT->common.R.wesn[YHI], GMT->common.R.wesn[YLO])) { possibly_donut_hell = true; anti_lon = GMT->current.proj.central_meridian + 180.0; if (anti_lon >= 360.0) anti_lon -= 360.0; anti_lat = -GMT->current.proj.pole; anti_bin = irint (floor ((90.0 - anti_lat) / c.bsize)) * c.bin_nx + irint (floor (anti_lon / c.bsize)); GMT_geo_to_xy (GMT, anti_lon, anti_lat, &anti_x, &anti_y); GMT_geo_to_xy (GMT, GMT->current.proj.central_meridian, GMT->current.proj.pole, &x_0, &y_0); if (Ctrl->G.active) GMT_Report (API, GMT_MSG_VERBOSE, "Warning: Fill/clip continent option (-G) may not work for this projection.\nIf the antipole (%g/%g) is in the ocean then chances are good\nElse: avoid projection center coordinates that are exact multiples of %g degrees\n", anti_lon, anti_lat, c.bsize); } if (possibly_donut_hell && paint_polygons && !clobber_background) { /* Force clobber when donuts may be called for now */ GMT_Report (API, GMT_MSG_VERBOSE, "Warning: -JE requires oceans to be painted first\n"); clobber_background = true; recursive = false; } if (clobber_background) { /* Paint entire map as ocean first, then lay land on top */ n = (int)GMT_map_clip_path (GMT, &xtmp, &ytmp, &donut); GMT_setfill (GMT, &Ctrl->S.fill, false); if (donut) { /* If donut, then the path consists of two path of np points */ PSL_plotline (PSL, xtmp, ytmp, n, PSL_MOVE + PSL_CLOSE); PSL_plotpolygon (PSL, &xtmp[n], &ytmp[n], n); } else PSL_plotpolygon (PSL, xtmp, ytmp, n); GMT_free (GMT, xtmp); GMT_free (GMT, ytmp); level_to_be_painted = 1; } if (recursive) { start_direction = stop_direction = (Ctrl->G.active) ? 1 : -1; level_to_be_painted = (Ctrl->G.active) ? 1 : 0; } else { start_direction = -1; stop_direction = 1; level_to_be_painted = (Ctrl->S.active) ? 0 : 1; } if (GMT->common.R.wesn[XLO] < 0.0 && GMT->common.R.wesn[XHI] > 0.0 && !GMT_IS_LINEAR (GMT)) greenwich = true; if (need_coast_base && (360.0 - fabs (GMT->common.R.wesn[XHI] - GMT->common.R.wesn[XLO]) ) < c.bsize) { GMT->current.map.is_world = true; if (GMT->current.proj.projection == GMT_GNOMONIC || GMT->current.proj.projection == GMT_GENPER) GMT->current.map.is_world = false; if (GMT_IS_AZIMUTHAL (GMT)) GMT->current.map.is_world = false; } if (GMT->current.map.is_world && greenwich) edge = GMT->current.proj.central_meridian; else if (!GMT->current.map.is_world && greenwich) { shift = true; edge = GMT->common.R.wesn[XLO]; if (edge < 0.0) edge += 360.0; if (edge > 180.0) edge = 180.0; } if (GMT->common.R.wesn[XLO] < 0.0 && GMT->common.R.wesn[XHI] <= 0.0) { /* Temporarily shift boundaries */ GMT->common.R.wesn[XLO] += 360.0; GMT->common.R.wesn[XHI] += 360.0; if (GMT->current.proj.central_meridian < 0.0) GMT->current.proj.central_meridian += 360.0; } if (need_coast_base) { west_border = floor (GMT->common.R.wesn[XLO] / c.bsize) * c.bsize; east_border = ceil (GMT->common.R.wesn[XHI] / c.bsize) * c.bsize; } if (!Ctrl->M.active && Ctrl->W.active) GMT_setpen (GMT, &Ctrl->W.pen[0]); last_pen_level = -1; if (clipping) PSL_beginclipping (PSL, xtmp, ytmp, 0, GMT->session.no_rgb, 1); /* Start clippath */ for (ind = 0; need_coast_base && ind < c.nb; ind++) { /* Loop over necessary bins only */ bin = c.bins[ind]; #ifdef DEBUG if (Ctrl->debug.active && bin != Ctrl->debug.bin) continue; #endif 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); } GMT_Report (API, GMT_MSG_VERBOSE, "Working on bin # %5d\r", bin); if (!Ctrl->M.active) PSL_comment (PSL, "Bin # %d\n", bin); if (GMT->current.map.is_world && greenwich) { left = c.bsize * (bin % c.bin_nx); right = left + c.bsize; shift = ((left - edge) <= 180.0 && (right - edge) > 180.0); } bin_trouble = (anti_bin == bin) ? anti_bin : -1; if (possibly_donut_hell) { /* Fill a 5-point polygon with bin corners */ bin_x[0] = bin_x[3] = bin_x[4] = c.lon_sw; bin_x[1] = bin_x[2] = c.lon_sw + c.bsize; bin_y[0] = bin_y[1] = bin_y[4] = c.lat_sw; bin_y[2] = bin_y[3] = c.lat_sw + c.bsize; /* Determine (x,y) coordinate of projected bin center */ GMT_geo_to_xy (GMT, c.lon_sw + 0.5 * c.bsize, c.lat_sw + 0.5 * c.bsize, &x_c, &y_c); dist = hypot (x_c - x_0, y_c - y_0); /* Distance from projection center to mid-point of current bin */ /* State donut_hell if we are close (within 20%) to the horizon or if the bin contains the antipole. * This is simply done with a Cartesian inside/outside test, which is adequate. * the 0.8 factor is arbitrary of course [PW] */ donut_hell = (dist > 0.8 * GMT->current.proj.r || GMT_non_zero_winding (GMT, anti_lon, anti_lat, bin_x, bin_y, 5)); } for (direction = start_direction; paint_polygons && direction <= stop_direction; direction += 2) { /* Assemble one or more segments into polygons */ if ((np = GMT_assemble_shore (GMT, &c, direction, true, west_border, east_border, &p)) == 0) continue; /* Get clipped polygons in x,y inches that can be plotted */ np_new = GMT_prep_shore_polygons (GMT, &p, np, donut_hell, 0.0, bin_trouble); if (clipping) { for (k = level_to_be_painted; k < GSHHS_MAX_LEVEL - 1; k++) recursive_path (GMT, PSL, -1, np_new, p, k, NULL); for (k = 0; k < np_new; k++) { /* Do any remaining interior polygons */ if (p[k].n == 0) continue; if (p[k].level % 2 == level_to_be_painted || p[k].level > 2) { PSL_plotline (PSL, p[k].lon, p[k].lat, p[k].n, PSL_MOVE); } } } else if (recursive) { /* Must avoid pointing anything but the polygons inside */ for (k = level_to_be_painted; k < GSHHS_MAX_LEVEL - 1; k++) recursive_path (GMT, PSL, -1, np_new, p, k, fill); for (k = 0; k < np_new; k++) { /* Do any remaining interior polygons */ if (p[k].n == 0) continue; if (p[k].level % 2 == level_to_be_painted || p[k].level > 2) { GMT_setfill (GMT, &fill[p[k].fid], false); PSL_plotpolygon (PSL, p[k].lon, p[k].lat, p[k].n); } } } else { /* Simply paints all polygons as is */ for (k = 0; k < np_new; k++) { if (p[k].n == 0 || p[k].level < level_to_be_painted) continue; if (donut_hell && GMT_non_zero_winding (GMT, anti_x, anti_y, p[k].lon, p[k].lat, p[k].n)) { /* Antipode inside polygon, must do donut */ n = (int)GMT_map_clip_path (GMT, &xtmp, &ytmp, &donut); GMT_setfill (GMT, &fill[p[k].fid], false); PSL_plotline (PSL, xtmp, ytmp, n, PSL_MOVE + PSL_CLOSE); PSL_plotpolygon (PSL, p[k].lon, p[k].lat, p[k].n); GMT_free (GMT, xtmp); GMT_free (GMT, ytmp); } else { GMT_setfill (GMT, &fill[p[k].fid], false); PSL_plotpolygon (PSL, p[k].lon, p[k].lat, p[k].n); } } } GMT_free_shore_polygons (GMT, p, np_new); GMT_free (GMT, p); } if (Ctrl->W.active && c.ns) { /* Draw or dump shorelines, no need to assemble polygons */ if ((np = GMT_assemble_shore (GMT, &c, 1, false, west_border, east_border, &p)) == 0) continue; for (i = 0; i < np; i++) { if (Ctrl->M.active) { sprintf (GMT->current.io.segment_header, "Shore Bin # %d, Level %d", bin, p[i].level); GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); for (k = 0; k < p[i].n; k++) { out[GMT_X] = p[i].lon[k]; out[GMT_Y] = p[i].lat[k]; GMT_Put_Record (API, GMT_WRITE_DOUBLE, out); } } else if (Ctrl->W.use[p[i].level-1]) { if (donut_hell) p[i].n = (int)GMT_fix_up_path (GMT, &p[i].lon, &p[i].lat, p[i].n, 0.0, 0); GMT->current.plot.n = GMT_geo_to_xy_line (GMT, p[i].lon, p[i].lat, p[i].n); if (!GMT->current.plot.n) continue; k = p[i].level - 1; if (k != last_pen_level) { GMT_setpen (GMT, &Ctrl->W.pen[k]); last_pen_level = k; } GMT_plot_line (GMT, GMT->current.plot.x, GMT->current.plot.y, GMT->current.plot.pen, GMT->current.plot.n); } } GMT_free_shore_polygons (GMT, p, np); GMT_free (GMT, p); } GMT_free_shore (GMT, &c); } if (need_coast_base) { GMT_Report (API, GMT_MSG_VERBOSE, "Working on bin # %5ld\n", bin); GMT_shore_cleanup (GMT, &c); } (void)GMT_DCW_operation (GMT, &Ctrl->F.info, NULL, Ctrl->M.active ? GMT_DCW_DUMP : GMT_DCW_PLOT); if (clipping) PSL_beginclipping (PSL, xtmp, ytmp, 0, GMT->session.no_rgb, 2); /* End clippath */ if (Ctrl->I.active) { /* Read rivers file and plot as lines */ GMT_Report (API, GMT_MSG_VERBOSE, "Adding Rivers..."); if (!Ctrl->M.active) PSL_comment (PSL, "Start of River segments\n"); last_k = -1; for (ind = 0; ind < r.nb; ind++) { /* Loop over necessary bins only */ bin = r.bins[ind]; GMT_get_br_bin (GMT, ind, &r, Ctrl->I.list, Ctrl->I.n_rlevels); if (r.ns == 0) continue; if (GMT->current.map.is_world && greenwich) { left = r.bsize * (bin % r.bin_nx); right = left + r.bsize; shift = ((left - edge) <= 180.0 && (right - edge) > 180.0); } if ((np = GMT_assemble_br (GMT, &r, shift, edge, &p)) == 0) { GMT_free_br (GMT, &r); continue; } for (i = 0; i < np; i++) { if (Ctrl->M.active) { sprintf (GMT->current.io.segment_header, "River Bin # %d, Level %d", bin, p[i].level); GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); for (k = 0; k < p[i].n; k++) { out[GMT_X] = p[i].lon[k]; out[GMT_Y] = p[i].lat[k]; GMT_Put_Record (API, GMT_WRITE_DOUBLE, out); } } else { GMT->current.plot.n = GMT_geo_to_xy_line (GMT, p[i].lon, p[i].lat, p[i].n); if (!GMT->current.plot.n) continue; k = p[i].level; if (k != last_k) { GMT_setpen (GMT, &Ctrl->I.pen[k]); last_k = k; } GMT_plot_line (GMT, GMT->current.plot.x, GMT->current.plot.y, GMT->current.plot.pen, GMT->current.plot.n); } } /* Free up memory */ GMT_free_br (GMT, &r); GMT_free_shore_polygons (GMT, p, np); GMT_free (GMT, p); } GMT_br_cleanup (GMT, &r); } if (Ctrl->N.active) { /* Read borders file and plot as lines */ double step; GMT_Report (API, GMT_MSG_VERBOSE, "Adding Borders..."); if (!Ctrl->M.active) PSL_comment (PSL, "Start of Border segments\n"); /* Must resample borders because some points may be too far apart and look like 'jumps' */ step = MAX (fabs(GMT->common.R.wesn[XLO] - GMT->common.R.wesn[XHI]), fabs (GMT->common.R.wesn[YHI] - GMT->common.R.wesn[YLO])) / 4.0; last_k = -1; for (ind = 0; ind < b.nb; ind++) { /* Loop over necessary bins only */ bin = b.bins[ind]; GMT_get_br_bin (GMT, ind, &b, Ctrl->N.list, Ctrl->N.n_blevels); if (b.ns == 0) continue; if (GMT->current.map.is_world && greenwich) { left = b.bsize * (bin % b.bin_nx); right = left + b.bsize; shift = ((left - edge) <= 180.0 && (right - edge) > 180.0); } if ((np = GMT_assemble_br (GMT, &b, shift, edge, &p)) == 0) { GMT_free_br (GMT, &b); continue; } for (i = 0; i < np; i++) { if (Ctrl->M.active) { sprintf (GMT->current.io.segment_header, "Border Bin # %d, Level %d", bin, p[i].level); GMT_Put_Record (API, GMT_WRITE_SEGMENT_HEADER, NULL); for (k = 0; k < p[i].n; k++) { out[GMT_X] = p[i].lon[k]; out[GMT_Y] = p[i].lat[k]; GMT_Put_Record (API, GMT_WRITE_DOUBLE, out); } } else { p[i].n = (int)GMT_fix_up_path (GMT, &p[i].lon, &p[i].lat, p[i].n, step, 0); GMT->current.plot.n = GMT_geo_to_xy_line (GMT, p[i].lon, p[i].lat, p[i].n); if (!GMT->current.plot.n) continue; k = p[i].level - 1; if (k != last_k) { GMT_setpen (GMT, &Ctrl->N.pen[k]); last_k = k; } GMT_plot_line (GMT, GMT->current.plot.x, GMT->current.plot.y, GMT->current.plot.pen, GMT->current.plot.n); } } /* Free up memory */ GMT_free_br (GMT, &b); GMT_free_shore_polygons (GMT, p, np); GMT_free (GMT, p); } GMT_br_cleanup (GMT, &b); } if (!Ctrl->M.active) { if (GMT->current.map.frame.init) { GMT->current.map.is_world = world_map_save; GMT_map_basemap (GMT); } if (Ctrl->L.active) GMT_draw_map_scale (GMT, &Ctrl->L.item); if (Ctrl->T.active) GMT_draw_map_rose (GMT, &Ctrl->T.item); GMT_plane_perspective (GMT, -1, 0.0); GMT_plotend (GMT); } else if (GMT_End_IO (API, GMT_OUT, 0) != GMT_OK) { Return (API->error); /* Disables further data output */ } GMT_Report (API, GMT_MSG_VERBOSE, "Done\n"); Return (GMT_OK); }