home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / pscoast.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-11-03  |  31.1 KB  |  853 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)pscoast.c    3.14  30 Jun 1995
  3.  *
  4.  *    Copyright (c) 1991-1995 by P. Wessel and W. H. F. Smith
  5.  *    See README file for copying and redistribution conditions.
  6.  *--------------------------------------------------------------------*/
  7. /*
  8.  * pscoast can draw shorelines, paint the landward and/or seaward side of
  9.  * the shoreline, paint lakes, draw political borders, and rivers.  Three
  10.  * data sets are considered: shores, borders, and rivers.  Each of these
  11.  * data sets come in 5 different resolutions.  The lower resolution files
  12.  * are derived from the full resolution data using the Douglas-Peucker (DP)
  13.  * line reduction algorithm.  By giving a tolerance in km, the algorithm
  14.  * will remove points that do not depart more than the tolerance from the
  15.  * straight line that would connect the points if the point in question was
  16.  * removed.  The resolutions are:
  17.  *
  18.  * full        - The complete World Vector Shoreline + CIA data base
  19.  * high        - DP reduced with tolerance = 0.2 km
  20.  * intermediate    - DP reduced with tolerance = 1 km
  21.  * low        - DP reduced with tolerance = 5 km
  22.  * crude    - DP reduced with tolerance = 25 km
  23.  *
  24.  * If selected, pscoast will open the desired binned shoreline file and fill
  25.  * the landmasses and/or oceans/lakes as shaded/colored/textured polygons. The
  26.  * shoreline may also be drawn.  Political boundaries and rivers may be plotted,
  27.  * too.  If only land is filled, then the 'wet' areas remain transparent.
  28.  * If only oceans are filled, then the 'dry' areas remain transparent.  This
  29.  * allows the user to overlay land or ocean without wiping out the plot.
  30.  * For more information about the binned polygon file, see the GMT Technical
  31.  * Reference and Cookbook.
  32.  * 
  33.  *
  34.  * Author:    Paul Wessel
  35.  * Date:    26-AUG-1994
  36.  * Version:    3.0    Rewritten from 2.1.4 to handle high-resolution
  37.  *            WVS + CIA reformatted and binned database
  38.  *
  39.  */
  40.  
  41. #include "gmt.h"
  42. #include "gmt_shore.h"        /* Defines shore structures */
  43.  
  44. struct SHORE c;
  45. struct BR b, r;
  46.  
  47. char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"};
  48.  
  49. main (argc, argv)
  50. int argc;
  51. char **argv; {
  52.     int i, j, n, np, ind, bin, base = 3, n_use, start, anti_bin = -1, level;
  53.     int level_to_be_painted, max_level = MAX_LEVEL, direction, np_new, k, last_k;
  54.     int blevels[N_BLEVELS], n_blevels = 0, rlevels[N_RLEVELS], n_rlevels = 0;
  55.     
  56.     BOOLEAN    error = FALSE, set_lake_fill = FALSE, draw_river = FALSE, shift, need_coast_base;
  57.     BOOLEAN    greenwich = FALSE, draw_coast = FALSE, fill_ocean = FALSE, fill_land = FALSE, possibly_donut_hell = FALSE;
  58.     BOOLEAN clobber_background, draw_border = FALSE, paint_polygons = FALSE, donut, fill_in_use = FALSE;
  59.     BOOLEAN miles = FALSE, draw_scale = FALSE, gave_xy = FALSE, fancy = FALSE, donut_hell = FALSE, got_fh[2];
  60.     
  61.     double    west = 0.0, east = 0.0, south = 0.0, north = 0.0, bin_x[5], bin_y[5];
  62.     double west_border, east_border, anti_lon, anti_lat, edge = 720.0, left, right;
  63.     double *xtmp, *ytmp, min_area = 0.0, x0, y0, scale_lat, length, step = 0.01;
  64.     double anti_x, anti_y, x_0, y_0, x_c, y_c, dist;
  65.     
  66.     char res = 'l', key[5], *string, txt_a[32], txt_b[32], txt_c[32];
  67.     
  68.     struct PEN coast, pen, bpen[N_BLEVELS], rpen[N_RLEVELS];
  69.     struct FILL fill[5];    /* Colors for (0) water, (1) land, (2) lakes, (3) islands in lakes, (4) lakes in islands in lakes */
  70.     
  71.     struct POL *p;
  72.         
  73.     gmt_init_fill (&fill[0], 255, 255, 255);    /* Default Ocean color */
  74.     gmt_init_fill (&fill[1], 0, 0, 0);        /* Default Land color */
  75.     gmt_init_fill (&fill[2], 255, 255, 255);    /* Default Lake color */
  76.     fill[3] = fill[1];
  77.     fill[4] = fill[2];
  78.     gmt_init_pen (&coast, 1);
  79.     gmt_init_pen (&pen, 1);
  80.     for (k = 0; k < N_RLEVELS; k++) gmt_init_pen (&rpen[k], 1);
  81.     for (k = 0; k < N_BLEVELS; k++) gmt_init_pen (&bpen[k], 1);
  82.     
  83.     memset ((char *)rlevels, 0, N_RLEVELS * sizeof (int));
  84.     memset ((char *)blevels, 0, N_BLEVELS * sizeof (int));
  85.     
  86.     argc = gmt_begin (argc, argv);
  87.     
  88.     /* Check and interpret the command line arguments */
  89.     
  90.     for (i = 1; i < argc; i++) {
  91.         if (argv[i][0] == '-') {
  92.             switch(argv[i][1]) {
  93.         
  94.                 /* Common parameters */
  95.                 
  96.                 case 'B':    /* Get tickmark information */
  97.                 case 'J':
  98.                 case 'K':
  99.                 case 'O':
  100.                 case 'P':
  101.                 case 'R':
  102.                 case 'U':
  103.                 case 'V':
  104.                 case 'X':
  105.                 case 'x':
  106.                 case 'Y':
  107.                 case 'y':
  108.                 case 'c':
  109.                 case '\0':
  110.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  111.                     break;
  112.                 
  113.                 /* Supplemental parameters */
  114.  
  115.                 case 'A':
  116.                     n = sscanf (&argv[i][2], "%lf/%d", &min_area, &max_level);
  117.                     if (n == 1) max_level = MAX_LEVEL;
  118.                     break;
  119.                 case 'C':    /* Lake color */
  120.                     if (argv[i][2] && gmt_getfill (&argv[i][2], &fill[2])) {
  121.                         gmt_fill_syntax ('C');
  122.                         error++;
  123.                     }
  124.                     set_lake_fill = TRUE;
  125.                     fill[4] = fill[2];
  126.                     break;
  127.                 case 'D':
  128.                     res = argv[i][2];
  129.                     switch (res) {
  130.                         case 'f':
  131.                             base = 0;
  132.                             break;
  133.                         case 'h':
  134.                             base = 1;
  135.                             break;
  136.                         case 'i':
  137.                             base = 2;
  138.                             break;
  139.                         case 'l':
  140.                             base = 3;
  141.                             break;
  142.                         case 'c':
  143.                             base = 4;
  144.                             break;
  145.                         default:
  146.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -D option:  Unknown modifier %c\n", gmt_program, argv[i][2]);
  147.                             break;
  148.                     }
  149.                     break;
  150.                 case 'E':
  151.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  152.                     break;
  153.                 case 'F':
  154.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  155.                         gmt_pen_syntax ('F');
  156.                         error++;
  157.                     }
  158.                     break;
  159.                 case 'G':        /* Set Gray shade or pattern */
  160.                     fill_land = TRUE;
  161.                     if (argv[i][2] && gmt_getfill (&argv[i][2], &fill[1])) {
  162.                         gmt_fill_syntax ('G');
  163.                         error++;
  164.                     }
  165.                     fill[3] = fill[1];
  166.                     break;
  167.                 case 'L':
  168.                     j = 2;
  169.                     if (argv[i][j] == 'f') fancy = TRUE, j++;
  170.                     if (argv[i][j] == 'x') gave_xy = TRUE, j++;
  171.                     k = sscanf (&argv[i][j], "%[^/]/%[^/]/%[^/]/%lf", txt_a, txt_b, txt_c, &length);
  172.                     x0 = ddmmss_to_degree (txt_a);
  173.                     y0 = ddmmss_to_degree (txt_b);
  174.                     scale_lat = ddmmss_to_degree (txt_c);
  175.                     if (k != 4) {
  176.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -L option:  Correct syntax\n", gmt_program);
  177.                         fprintf (stderr, "    -L[f][x]<x0>/<y0>/<lat>/<length>[m], append m for miles [Default is km]\n");
  178.                         error++;
  179.                     }
  180.                     miles = (strchr (&argv[i][2], 'm') != CNULL);
  181.                     draw_scale = TRUE;
  182.                     break;
  183.                 case 'N':
  184.                     if (!argv[i][2]) {
  185.                         fprintf (stderr, "%s: GMT SYNTAX ERROR:  -N option takes two arguments\n", gmt_program);
  186.                         error++;
  187.                         continue;
  188.                     }
  189.                     draw_border = TRUE;
  190.                     for (k = 2; argv[i][k] && argv[i][k] != '/'; k++);
  191.                     if (argv[i][k] && argv[i][k] == '/') {    /* Specified pen args */
  192.                         strncpy (key, &argv[i][2], k-2);
  193.                         string = &argv[i][k+1];
  194.                     }
  195.                     else {    /* No pen specified, use default */
  196.                         strcpy (key, &argv[i][2]);
  197.                         string = CNULL;
  198.                     }
  199.                     gmt_init_pen (&pen, 1);
  200.                     if (string && gmt_getpen (string, &pen)) {
  201.                         gmt_pen_syntax ('N');
  202.                         error++;
  203.                     }
  204.                         
  205.                     switch (key[0]) {
  206.                         case 'a':
  207.                             for (k = 0; k < N_BLEVELS; k++) blevels[k] = TRUE;
  208.                             for (k = 0; k < N_BLEVELS; k++) bpen[k] = pen;
  209.                             break;
  210.                         default:
  211.                             k = atoi (key) - 1;
  212.                             if (k < 0 || k >= N_BLEVELS) {
  213.                                 fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Feature not in list!\n", gmt_program);
  214.                                 error++;
  215.                             }
  216.                             else {
  217.                                 blevels[k] = TRUE;
  218.                                 bpen[k] = pen;
  219.                             }
  220.                             break;
  221.                     }
  222.                     break;
  223.                     
  224.                 case 'I':
  225.                     if (!argv[i][2]) {
  226.                         fprintf (stderr, "%s: GMT SYNTAX ERROR:  -I option takes two arguments\n", gmt_program);
  227.                         error++;
  228.                         continue;
  229.                     }
  230.                     draw_river = TRUE;
  231.                     for (k = 2; argv[i][k] && argv[i][k] != '/'; k++);
  232.                     if (argv[i][k] && argv[i][k] == '/') {    /* Specified pen args */
  233.                         strncpy (key, &argv[i][2], k-2);
  234.                         string = &argv[i][k+1];
  235.                     }
  236.                     else {    /* No pen specified, use default */
  237.                         strcpy (key, &argv[i][2]);
  238.                         string = CNULL;
  239.                     }
  240.                     gmt_init_pen (&pen, 1);
  241.                     if (string && gmt_getpen (string, &pen)) {
  242.                         gmt_pen_syntax ('I');
  243.                         error++;
  244.                     }
  245.                         
  246.                     switch (key[0]) {
  247.                         case 'a':
  248.                             for (k = 0; k < N_RLEVELS; k++) rlevels[k] = TRUE;
  249.                             for (k = 0; k < N_RLEVELS; k++) rpen[k] = pen;
  250.                             break;
  251.                         case 'r':
  252.                             for (k = 0; k < 4; k++) rlevels[k] = TRUE;
  253.                             for (k = 0; k < 4; k++) rpen[k] = pen;
  254.                             break;
  255.                         case 'i':
  256.                             for (k = 5; k < 8; k++) rlevels[k] = TRUE;
  257.                             for (k = 5; k < 8; k++) rpen[k] = pen;
  258.                             break;
  259.                         case 'c':
  260.                             rlevels[8] = rlevels[9] = TRUE;
  261.                             rpen[8] = rpen[9] = pen;
  262.                             break;
  263.                         default:
  264.                             k = atoi (key) - 1;
  265.                             if (k < 0 || k >= N_RLEVELS) {
  266.                                 fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Feature not in list!\n", gmt_program);
  267.                                 error++;
  268.                             }
  269.                             else {
  270.                                 rlevels[k] = TRUE;
  271.                                 rpen[k] = pen;
  272.                             }
  273.                             break;
  274.                     }
  275.                     break;
  276.                     
  277.                 case 'S':        /* Set ocean color if needed */
  278.                     if (argv[i][2] && gmt_getfill (&argv[i][2], &fill[0])) {
  279.                         gmt_fill_syntax ('S');
  280.                         error++;
  281.                     }
  282.                     fill_ocean = TRUE;
  283.                     break;
  284.                 case 'W':
  285.                     if (argv[i][2] && gmt_getpen (&argv[i][2], &coast)) {
  286.                         gmt_pen_syntax ('W');
  287.                         error++;
  288.                     }
  289.                     draw_coast = TRUE;
  290.                     break;
  291.                 case '-':    /* Undocumented option to increase path-fix resolution */
  292.                     step = atof (&argv[i][2]);
  293.                     break;
  294.  
  295.                 default:        /* Options not recognized */
  296.                     error = TRUE;
  297.                     gmt_default_error (argv[i][1]);
  298.                     break;
  299.             }
  300.         }
  301.     }
  302.     
  303.     find_resolutions (got_fh);    /* Check to see if full &| high resolution coastlines are installed */
  304.  
  305.     if (argc == 1 || gmt_quick) {    /* Display usage */
  306.         fprintf (stderr,"pscoast %s - Plot continents, shorelines, rivers, and borders on maps\n\n", GMT_VERSION);
  307.         fprintf (stderr,"usage: pscoast -J<params> -R<west>/<east>/<south>/<north> [-A<min_area>[/<max_level>]] [-B<tickinfo>]\n");
  308.         fprintf (stderr, "     [-C[<fill>]] [-D<resolution>] [-Eaz/el] [-F<r/g/b] [-G[<fill>]] [-I<feature>[/<pen>]]\n");
  309.         fprintf (stderr, "     [-K] [-L[f][x]<lon0>/<lat0>/<slat>/<length>[m]] [-N<feature>[/<pen>]] [-O] [-P]\n");
  310.         fprintf (stderr, "     [-S<fill>] [-U[dx/dy/][label]] [-V] [-W[<pen>]] [-X<x_shift>] [-Y<y_shift>] [-c<ncopies>]\n");
  311.         
  312.         if (gmt_quick) exit (-1);
  313.         
  314.         explain_option ('j');
  315.         explain_option ('R');
  316.         fprintf (stderr, "\n\tOPTIONS:\n");
  317.         fprintf (stderr, "    -A features smaller than <min_area> (in km^2) or of higher level (0-4) than <max_level>\n");
  318.         fprintf (stderr, "    will be skipped [0/4 (4 means lake inside island inside lake)]\n");
  319.         explain_option ('b');
  320.         fprintf (stderr, "    -C sets separate color for lakes. [Default is same as ocean (white)]. Set\n");
  321.         fprintf (stderr, "       1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
  322.         fprintf (stderr, "       2) p[or P]<iconsize>/<pattern> for predefined patterns (0-31).\n");
  323.         fprintf (stderr, "    -D Choose one of the following resolutions:\n");
  324.         if (got_fh[0]) fprintf (stderr, "       f - full resolution (may be very slow for large regions)\n");
  325.         if (got_fh[1]) fprintf (stderr, "       h - high resolution (may be slow for large regions)\n");
  326.         fprintf (stderr, "       i - intermediate resolution\n");
  327.         fprintf (stderr, "       l - low resolution [Default]\n");
  328.         fprintf (stderr, "       c - crude resolution, for busy plots that need crude continent outlines only\n");
  329.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  330.         fprintf (stderr, "      -F Set color used for Frame and anotation [%d/%d/%d]\n",
  331.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  332.         fprintf (stderr, "    -G Paint \"dry\" areas.  Append desired paint [Default is black] as:\n");
  333.         fprintf (stderr, "       1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
  334.         fprintf (stderr, "       2) p[or P]<iconsize>/<pattern> for predefined patterns (0-31).\n");
  335.         fprintf (stderr, "    -I draws rIvers.  Specify feature and optionally append pen [Default is solid line of unit thickness].\n");
  336.         fprintf (stderr, "       Choose from the features below.  Repeat the -I option as many times as needed\n");
  337.         fprintf (stderr, "           1 = Permanent major rivers\n");
  338.         fprintf (stderr, "           2 = Additional major rivers\n");
  339.         fprintf (stderr, "           3 = Additional rivers\n");
  340.         fprintf (stderr, "           4 = Minor rivers\n");
  341.         fprintf (stderr, "           5 = Double lined rivers\n");
  342.         fprintf (stderr, "           6 = Intermittent rivers - major\n");
  343.         fprintf (stderr, "           7 = Intermittent rivers - additional\n");
  344.         fprintf (stderr, "           8 = Intermittent rivers - minor\n");
  345.         fprintf (stderr, "           9 = Major canals\n");
  346.         fprintf (stderr, "          10 = Minor canals\n");
  347.         fprintf (stderr, "           a = All rivers and canals (1-10)\n");
  348.         fprintf (stderr, "           r = All permanent rivers (1-4)\n");
  349.         fprintf (stderr, "           i = All intermittent rivers (6-8)\n");
  350.         fprintf (stderr, "           c = All canals (9-10)\n");
  351.         explain_option ('K');
  352.         fprintf (stderr, "      -L draws a simple map scaLe centered on <lon0>/<lat0>.  Use -Lx to specify cartesian coordinates instead.\n");
  353.         fprintf (stderr, "         Scale is calculated at latitude <slat>. <length> is in km, or miles is m is appended\n");
  354.         fprintf (stderr, "         -Lf draws a \"fancy\" scale [Default is plain]\n");
  355.         fprintf (stderr, "    -N draws boundaries.  Specify feature and optionally append pen [Default is solid line of unit thickness]\n");
  356.         fprintf (stderr, "       Choose from the features below.  Repeat the -N option as many times as needed\n");
  357.         fprintf (stderr, "           1 = National boundaries\n");
  358.         fprintf (stderr, "           2 = State boundaries within the Americas\n");
  359.         fprintf (stderr, "           3 = Marine boundaries\n");
  360.         fprintf (stderr, "           a = All boundaries (1-3)\n");
  361.         explain_option ('O');
  362.         explain_option ('P');
  363.         fprintf (stderr, "    -S Paint \"wet\" areas.  Append desired paint [Default is white] as:\n");
  364.         fprintf (stderr, "       1) <r/g/b> (each 0-255) for color or <gray> (0-255) for gray-shade [0].\n");
  365.         fprintf (stderr, "       2) p[or P]<iconsize>/<pattern> for predefined patterns (0-31).\n");
  366.         explain_option ('U');
  367.         explain_option ('V');
  368.         fprintf (stderr,"    -W draw shorelines.  Append pen [Default is solid line of unit thickness].\n");
  369.         explain_option ('X');
  370.         explain_option ('c');
  371.         explain_option ('.');
  372.         exit (-1);
  373.     }
  374.  
  375.     /* Check that the options selected are mutually consistant */
  376.     
  377.     if (!project_info.region_supplied) {
  378.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  379.         error++;
  380.     }
  381.     if (coast.width < 0) {
  382.         fprintf (stderr, "%s: GMT SYNTAX ERROR -W option:  Pen thickness cannot be negative\n", gmt_program);
  383.         error++;
  384.     }
  385.     if (!(fill_land || fill_ocean || draw_coast || draw_border || draw_river)) {
  386.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -G, -S, -I, -N, and -W\n", gmt_program);
  387.         error++;
  388.     }
  389.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  390.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  391.         error++;
  392.     }
  393.     
  394.     if (draw_river) {
  395.         for (k = 0; k < N_RLEVELS; k++) {
  396.             if (!rlevels[k]) continue;
  397.             if (rpen[k].width < 0) {
  398.                 fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Pen thickness cannot be negative\n", gmt_program);
  399.                 error++;
  400.             }
  401.             rlevels[n_rlevels] = k + 1;
  402.             rpen[n_rlevels] = rpen[k];
  403.             n_rlevels++;
  404.         }
  405.     }
  406.     
  407.     if (draw_border) {
  408.         for (k = 0; k < N_BLEVELS; k++) {
  409.             if (!blevels[k]) continue;
  410.             if (bpen[k].width < 0) {
  411.                 fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Pen thickness cannot be negative\n", gmt_program);
  412.                 error++;
  413.             }
  414.             blevels[n_blevels] = k + 1;
  415.             bpen[n_blevels] = bpen[k];
  416.             n_blevels++;
  417.         }
  418.     }
  419.  
  420.     if (error) exit (-1);
  421.  
  422.     need_coast_base = (fill_land || fill_ocean || draw_coast);
  423.     clobber_background = (fill_land && fill_ocean);
  424.     paint_polygons = (fill_land || fill_ocean);
  425.  
  426.     if (east > 360.0) {
  427.         west -= 360.0;
  428.         east -= 360.0;
  429.     }
  430.     
  431.     map_setup (west, east, south, north);
  432.  
  433.     if (need_coast_base && gmt_init_shore (res, &c, project_info.w, project_info.e, project_info.s, project_info.n))  {
  434.         fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", gmt_program, shore_resolution[base]);
  435.         need_coast_base = FALSE;
  436.     }
  437.     
  438.     if (draw_border && gmt_init_br ('b', res, &b, project_info.w, project_info.e, project_info.s, project_info.n)) {
  439.         fprintf (stderr, "%s: %s resolution political boundary data base not installed\n", gmt_program, shore_resolution[base]);
  440.         draw_border = FALSE;
  441.     }
  442.     
  443.     if (draw_river && gmt_init_br ('r', res, &r, project_info.w, project_info.e, project_info.s, project_info.n)) {
  444.         fprintf (stderr, "%s: %s resolution river data base not installed\n", gmt_program, shore_resolution[base]);
  445.         draw_river = FALSE;
  446.     }
  447.     
  448.     if (! (need_coast_base || draw_border || draw_river)) {
  449.         fprintf (stderr, "%s: No databases available - aborts\n", gmt_program);
  450.         exit (-1);
  451.     }
  452.     
  453.     
  454.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  455.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  456.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo(argv[0]));
  457.     echo_command (argc, argv);
  458.     if (gmtdefs.unix_time) timestamp (argc, argv);
  459.     if (base > 0) ps_command ("2 setlinejoin");    /* Avoid excessive mitering for cruder lines */
  460.          
  461.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  462.     
  463.     if (fill_ocean && !set_lake_fill) fill[2] = fill[4] = fill[0];    /* Use ocean fill by default for lakes */
  464.     
  465.     for (i = 0; i < 5; i++) if (fill[i].use_pattern) {    /* Must initialize patterns */
  466.         ps_imagefill_init (fill[i].pattern_no, fill[i].pattern, fill[i].inverse, fill[i].icon_size);
  467.         fill_in_use = TRUE;
  468.     }
  469.     
  470.     if (fill_in_use && !clobber_background) {    /* Force clobber when fill is used for now */
  471.         fprintf (stderr, "%s: Warning: Pattern fill requires oceans to be painted first\n", gmt_program);
  472.         clobber_background = TRUE;
  473.     }
  474.         
  475.     if (project_info.projection == AZ_EQDIST && fabs (project_info.w - project_info.e) == 360.0 && (project_info.n - project_info.s) == 180.0) {
  476.         possibly_donut_hell = TRUE;
  477.         anti_lon = project_info.central_meridian + 180.0;
  478.         if (anti_lon >= 360.0) anti_lon -= 360.0;
  479.         anti_lat = -project_info.pole;
  480.         anti_bin = floor ((90.0 - anti_lat) / c.bsize) * c.bin_nx + floor (anti_lon / c.bsize);
  481.         geo_to_xy (anti_lon, anti_lat, &anti_x, &anti_y);
  482.         geo_to_xy (project_info.central_meridian, project_info.pole, &x_0, &y_0);
  483.         if (fill_land) {
  484.             fprintf (stderr, "%s: Warning: Fill continent option (-G) may not work for this projection.\n", gmt_program);
  485.             fprintf (stderr, "If the antipole (%lg/%lg) is in the ocean then chances are good\n", anti_lon, anti_lat);
  486.             fprintf (stderr, "Else: avoid projection center coordinates that are exact multiples of %lg degrees\n", c.bsize);
  487.         }
  488.     }
  489.  
  490.     if (possibly_donut_hell && !clobber_background) {    /* Force clobber when donuts may be called for now */
  491.         fprintf (stderr, "%s: Warning: -JE requires oceans to be painted first\n", gmt_program);
  492.         clobber_background = TRUE;
  493.     }
  494.         
  495.     if (clobber_background) {    /* Paint entire map as ocean first, then lay land on top */
  496.         n = map_clip_path (&xtmp, &ytmp, &donut);
  497.         if (donut) {
  498.             ps_line (xtmp, ytmp, n, 1, TRUE, FALSE);
  499.             ps_line (&xtmp[n], ytmp[n], n, -1, TRUE, FALSE);
  500.             dump_paint (xtmp, ytmp, n, &fill[0], -1);
  501.         }
  502.         else 
  503.             dump_paint (xtmp, ytmp, n, &fill[0], FALSE);
  504.         free ((char *)xtmp);
  505.         free ((char *)ytmp);
  506.         fill_ocean = FALSE;
  507.         direction = 1;
  508.     }
  509.     else if (fill_land) {
  510.         direction = 1;
  511.         level_to_be_painted = 1;
  512.     }
  513.     else {
  514.         direction = -1;
  515.         level_to_be_painted = 0;
  516.     }
  517.  
  518.     if (west < 0.0 && east > 0.0 && project_info.projection > 5) greenwich = TRUE;
  519.     if (gmt_world_map && greenwich)
  520.         edge = project_info.central_meridian;
  521.     else if (!gmt_world_map && greenwich) {
  522.         shift = TRUE;
  523.         edge = west;    if (edge < 0.0) edge += 360.0;
  524.         if (edge > 180.0) edge = 180.0;
  525.     }
  526.     
  527.     west_border = project_info.w;    east_border = project_info.e;
  528.     if (project_info.w < 0.0 && project_info.e <= 0.0) {    /* Temporarily shift boundaries */
  529.         project_info.w += 360.0;
  530.         project_info.e += 360.0;
  531.     }
  532.         
  533.     if (draw_coast) {
  534.         ps_setpaint (coast.r, coast.g, coast.b);
  535.         ps_setline (coast.width);
  536.         if (coast.texture) ps_setdash (coast.texture, coast.offset);
  537.     }
  538.  
  539.     for (ind = 0; need_coast_base && ind < c.nb; ind++) {    /* Loop over necessary bins only */
  540.         
  541.         bin = c.bins[ind];
  542.         gmt_get_shore_bin (ind, &c, min_area, max_level);
  543.         
  544.         if (gmtdefs.verbose) fprintf (stderr, "pscoast: Working on block # %5d\r", bin);
  545.         printf ("%% Bin # %d\n", bin);
  546.  
  547.         if (gmt_world_map && greenwich) {
  548.             left = c.bsize * (bin % c.bin_nx);    right = left + c.bsize;
  549.             shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
  550.         }
  551.         
  552.         if (anti_bin == bin)
  553.             anti_bin = bin;
  554.             
  555.         if (possibly_donut_hell) {
  556.             bin_x[0] = bin_x[3] = bin_x[4] = c.lon_sw;
  557.             bin_x[1] = bin_x[2] = c.lon_sw + c.bsize;
  558.             bin_y[0] = bin_y[1] = bin_y[4] = c.lat_sw;
  559.             bin_y[2] = bin_y[3] = c.lat_sw + c.bsize;
  560.             geo_to_xy (c.lon_sw + 0.5 * c.bsize, c.lat_sw + 0.5 * c.bsize, &x_c, &y_c);
  561.             dist = hypot (x_c - x_0, y_c - y_0);
  562.             donut_hell = non_zero_winding (anti_lon, anti_lat, bin_x, bin_y, 5);
  563.             if (!donut_hell) donut_hell = (dist > 0.8 * project_info.r);
  564.         }
  565.         
  566.         if (paint_polygons && (np = gmt_assemble_shore (&c, direction, TRUE, shift, edge, &p))) {
  567.         
  568.             /* Have assembled segments into polygons */
  569.         
  570.             np_new = np;
  571.         
  572.             for (k = 0; k < np; k++) {
  573.             
  574.                 if (donut_hell) p[k].n = fix_up_path (&p[k].lon, &p[k].lat, p[k].n, greenwich, step);
  575.                 
  576.                 /* Clip polygon against map boundary if necessary and return plot x,y in inches */
  577.                 
  578.                 if ((n = clip_to_map (p[k].lon, p[k].lat, p[k].n, &xtmp, &ytmp)) == 0) {    /* Completely outside */
  579.                     p[k].n = 0;
  580.                     continue;
  581.                 }
  582.             
  583.                 /* Must check if polygon must be split and partially plotted at both edges of map */
  584.                 
  585.                 if (will_it_wrap (xtmp, ytmp, n, &start)) {    /* Polygon does indeed wrap */
  586.                 
  587.                     /* First truncate agains left border */
  588.                         
  589.                     gmt_n_plot = truncate_left (xtmp, ytmp, n, start);
  590.                     n_use = compact_line (gmt_x_plot, gmt_y_plot, gmt_n_plot, FALSE, 0);
  591.                     if (project_info.three_D) two_d_to_three_d (gmt_x_plot, gmt_y_plot, gmt_n_plot);
  592.                     p[k].lon = (double *) memory ((char *)p[k].lon, n_use, sizeof (double), "pscoast");
  593.                     p[k].lat = (double *) memory ((char *)p[k].lat, n_use, sizeof (double), "pscoast");
  594.                     memcpy ((char *)p[k].lon, (char *)gmt_x_plot, n_use * sizeof (double));
  595.                     memcpy ((char *)p[k].lat, (char *)gmt_y_plot, n_use * sizeof (double));
  596.                     p[k].n = n_use;
  597.                                 
  598.                     /* Then truncate agains right border */
  599.                         
  600.                     gmt_n_plot = truncate_right (xtmp, ytmp, n, start);
  601.                     n_use = compact_line (gmt_x_plot, gmt_y_plot, gmt_n_plot, FALSE, 0);
  602.                     if (project_info.three_D) two_d_to_three_d (gmt_x_plot, gmt_y_plot, gmt_n_plot);
  603.                     p = (struct POL *) memory ((char *)p, np_new + 1, sizeof (struct POL), "pscoast");
  604.                     p[np_new].lon = (double *) memory (CNULL, n_use, sizeof (double), "pscoast");
  605.                     p[np_new].lat = (double *) memory (CNULL, n_use, sizeof (double), "pscoast");
  606.                     memcpy ((char *)p[np_new].lon, (char *)gmt_x_plot, n_use * sizeof (double));
  607.                     memcpy ((char *)p[np_new].lat, (char *)gmt_y_plot, n_use * sizeof (double));
  608.                     p[np_new].n = n_use;
  609.                     p[np_new].interior = p[k].interior;
  610.                     p[np_new].level = p[k].level;
  611.                     np_new++;
  612.                 }
  613.                 else {
  614.                     n_use = compact_line (xtmp, ytmp, n, FALSE, 0);
  615.                     if (project_info.three_D) two_d_to_three_d (xtmp, ytmp, n_use);
  616.                     if (anti_bin == bin && step == 0.0) {    /* Must warn for donut effect */
  617.                         fprintf (stderr, "%s: GMT Warning: Antipodal bin # %d not filled!\n", gmt_program, anti_bin);
  618.                         free ((char *)xtmp);
  619.                         free ((char *)ytmp);
  620.                         continue;
  621.                     }
  622.                     
  623.                     else {
  624.                         p[k].lon = (double *) memory ((char *)p[k].lon, n_use, sizeof (double), "pscoast");
  625.                         p[k].lat = (double *) memory ((char *)p[k].lat, n_use, sizeof (double), "pscoast");
  626.                         memcpy ((char *)p[k].lon, (char *)xtmp, n_use * sizeof (double));
  627.                         memcpy ((char *)p[k].lat, (char *)ytmp, n_use * sizeof (double));
  628.                         p[k].n = n_use;
  629.                     }
  630.                 }
  631.                 
  632.                 free ((char *)xtmp);
  633.                 free ((char *)ytmp);
  634.             }
  635.         
  636.             /* Now we have clipped polygons in x,y inches that can be plotted */
  637.         
  638.             if (clobber_background) {    /* Simply points all polygons as is */
  639.                 for (k = 0; k < np_new; k++) {
  640.                     level = MIN (p[k].level, 2);
  641.                     if (donut_hell && non_zero_winding (anti_x, anti_y, p[k].lon, p[k].lat, p[k].n)) {    /* Antipode inside polygon, must do donut */
  642.                         n = map_clip_path (&xtmp, &ytmp, &donut);
  643.                         ps_line (xtmp, ytmp, n, 1, TRUE, FALSE);
  644.                         ps_line (p[k].lon, p[k].lat, p[k].n, -1, TRUE, FALSE);
  645.                         dump_paint (xtmp, ytmp, n, &fill[level], -1);
  646.                         free ((char *)xtmp);
  647.                         free ((char *)ytmp);
  648.                     }
  649.                     else
  650.                         dump_paint (p[k].lon, p[k].lat, p[k].n, &fill[level], FALSE);
  651.                 }
  652.             }
  653.             else {    /* Must avoid pointing anything but the polygons inside */
  654.             
  655.                 for (k = level_to_be_painted; k < MAX_LEVEL - 1; k++) {
  656.                     recursive_painting (-1, np_new, p, k, k, fill);
  657.                 }
  658.                 
  659.                 for (k = 0; k < np_new; k++) {    /* Do any remaining interior polygons */
  660.                     if (p[k].n == 0) continue;
  661.                     if (p[k].level%2 == level_to_be_painted)
  662.                         dump_paint (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], FALSE);
  663.                 }
  664.                 
  665.             }
  666.             
  667.             free_polygons (p, np_new);
  668.             free ((char *)p);
  669.         }
  670.         
  671.         if (draw_coast && c.ns) {    /* Draw shorelines, no need to assemble polygons */
  672.         
  673.             np = gmt_assemble_shore (&c, direction, FALSE, shift, edge, &p);
  674.         
  675.             for (i = 0; i < np; i++) {
  676.                 if (donut_hell) p[i].n = fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step);
  677.                 gmt_n_plot = geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
  678.                 if (gmt_n_plot) gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  679.             }
  680.             
  681.             free_polygons (p, np);
  682.             free ((char *)p);
  683.         }
  684.         
  685.         free_shore (&c);
  686.     
  687.     }
  688.     if (need_coast_base) shore_cleanup (&c);
  689.  
  690.     if (gmtdefs.verbose) fprintf (stderr, "\n");
  691.     
  692.     if (draw_river) {    /* Read rivers file and plot as lines */
  693.         
  694.         if (gmtdefs.verbose) fprintf (stderr, "pscoast: Adding Rivers...");
  695.         ps_comment ("Start of River segments");
  696.         last_k = -1;
  697.  
  698.         for (ind = 0; ind < r.nb; ind++) {    /* Loop over necessary bins only */
  699.         
  700.             bin = r.bins[ind];
  701.             gmt_get_br_bin (ind, &r, rlevels, n_rlevels);
  702.             
  703.             if (r.ns == 0) continue;
  704.             
  705.             if (gmt_world_map && greenwich) {
  706.                 left = r.bsize * (bin % r.bin_nx);    right = left + r.bsize;
  707.                 shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
  708.             }
  709.             
  710.             np = gmt_assemble_br (&r, shift, edge, &p);
  711.             
  712.             for (i = 0; i < np; i++) {
  713.                 gmt_n_plot = geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
  714.                 if (!gmt_n_plot) continue;
  715.  
  716.                 k = p[i].level - 1;
  717.                 if (k != last_k) {
  718.                     ps_setline (rpen[k].width);
  719.                     ps_setpaint (rpen[k].r, rpen[k].g, rpen[k].b);
  720.                     if (rpen[k].texture) ps_setdash (rpen[k].texture, rpen[k].offset);
  721.                     last_k = k;
  722.                 }
  723.                 gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  724.             }
  725.             
  726.             /* Free up memory */
  727.         
  728.             free_br (&r);
  729.             for (k = 0; k < np; k++) {
  730.                 free ((char *)p[k].lon);
  731.                 free ((char *)p[k].lat);
  732.             }
  733.             free ((char *)p);
  734.         }
  735.         br_cleanup (&r);
  736.         ps_setdash (CNULL, 0);
  737.         
  738.         if (gmtdefs.verbose) fprintf (stderr, "\n");
  739.     }
  740.  
  741.     if (draw_border) {    /* Read borders file and plot as lines */
  742.         
  743.         if (gmtdefs.verbose) fprintf (stderr, "pscoast: Adding Borders...");
  744.         ps_comment ("Start of Border segments");
  745.         
  746.         /* Must resample borders because some points may be too far apart and look like 'jumps' */
  747.         
  748.         step = MAX (fabs(project_info.w - project_info.e), fabs (project_info.n - project_info.s)) / 4.0;
  749.  
  750.         last_k = -1;
  751.  
  752.         for (ind = 0; ind < b.nb; ind++) {    /* Loop over necessary bins only */
  753.         
  754.             bin = b.bins[ind];
  755.             gmt_get_br_bin (ind, &b, blevels, n_blevels);
  756.             
  757.             if (b.ns == 0) continue;
  758.             
  759.             if (gmt_world_map && greenwich) {
  760.                 left = b.bsize * (bin % b.bin_nx);    right = left + b.bsize;
  761.                 shift = ((left - edge) <= 180.0 && (right - edge) > 180.0);
  762.             }
  763.  
  764.             np = gmt_assemble_br (&b, shift, edge, &p);
  765.             for (i = 0; i < np; i++) {
  766.                 p[i].n = fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step);
  767.                 gmt_n_plot = geo_to_xy_line (p[i].lon, p[i].lat, p[i].n);
  768.                 if (!gmt_n_plot) continue;
  769.  
  770.                 k = p[i].level - 1;
  771.                 if (k != last_k) {
  772.                     ps_setline (bpen[k].width);
  773.                     ps_setpaint (bpen[k].r, bpen[k].g, bpen[k].b);
  774.                     if (bpen[k].texture) ps_setdash (bpen[k].texture, bpen[k].offset);
  775.                     last_k = k;
  776.                 }
  777.                 gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  778.             }
  779.             
  780.             /* Free up memory */
  781.         
  782.             free_br (&b);
  783.             for (k = 0; k < np; k++) {
  784.                 free ((char *)p[k].lon);
  785.                 free ((char *)p[k].lat);
  786.             }
  787.             free ((char *)p);
  788.         }
  789.         br_cleanup (&b);
  790.         ps_setdash (CNULL, 0);
  791.         
  792.         if (gmtdefs.verbose) fprintf (stderr, "\n");
  793.     }
  794.     
  795.     if (base > 0) ps_command ("0 setlinejoin");    /* Reset mitering */
  796.  
  797.     if (frame_info.plot) {
  798.         project_info.w = west_border;
  799.         project_info.e = east_border;
  800.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  801.         map_basemap ();
  802.     }
  803.     
  804.     if (draw_scale) draw_map_scale (x0, y0, scale_lat, length, miles, gave_xy, fancy);
  805.  
  806.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  807.  
  808.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  809.     ps_plotend (gmtdefs.last_page);
  810.     
  811.     gmt_end (argc, argv);
  812. }
  813.  
  814. int dump_paint (x, y, n, f, status)
  815. double x[], y[];
  816. int n;
  817. struct FILL *f;
  818. BOOLEAN status; {
  819.     if (f->use_pattern)
  820.         ps_imagefill (x, y, n, f->pattern_no, f->pattern, f->inverse, f->icon_size, status);
  821.     else
  822.         ps_polygon (x, y, n, f->r, f->g, f->b, status);
  823. }
  824.  
  825. int recursive_painting (k0, np, p, level, start_level, fill)
  826. int k0, np, level, start_level;
  827. struct POL p[];
  828. struct FILL fill[]; {
  829.     int k, flag;
  830.     BOOLEAN go;
  831.     char text[50];
  832.     if (level > MAX_LEVEL) return;
  833.     flag = (level == start_level) ? 1 : 0;
  834.     for (k = k0 + 1; k < np; k++) {
  835.         if (p[k].n == 0 || p[k].level < level) continue;
  836.         go = (p[k].level == level) ? ((k0 == -1) ? TRUE : non_zero_winding (p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) : FALSE;
  837.         if (go) {
  838.             sprintf (text, "Polygon %d\0", k);
  839.             ps_comment (text);
  840.             ps_clipon (p[k].lon, p[k].lat, p[k].n, -1, -1, -1, flag);
  841.             recursive_painting (k, np, p, level+1, start_level, fill);
  842.             p[k].n = 0;
  843.         }
  844.         
  845.         if (p[k].level == start_level) {    /* Done nesting, time to paint */
  846.             ps_clipon (p[k].lon, p[k].lat, 0, -1, -1, -1, 4);    /* Dummy to end path */
  847.             dump_paint (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], -1);
  848.         }
  849.     }
  850.     
  851.     return;
  852. }
  853.