/*-------------------------------------------------------------------- * $Id: pscoast.c,v 1.3.4.3 2002/02/27 17:58:55 pwessel Exp $ * * Copyright (c) 1991-2002 by P. Wessel and W. H. F. Smith * See COPYING 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 General Public License as published by * the Free Software Foundation; version 2 of the License. * * 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 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: 24-JUN-2000 * Version: 3.4.1 * */ #include "gmt.h" #ifdef DEBUG int debug_bin = -1; #endif void recursive_painting (int k0, int np, struct POL *p, int level, int start_level, struct GMT_FILL fill[]); void recursive_clipping (int k0, int np, struct POL *p, int level, int start_level, int *n_pieces); main (int argc, char **argv) { int i, n, np, ind, bin, base = 3, anti_bin = -1, level, flag; int level_to_be_painted, max_level = MAX_LEVEL, direction, np_new, k, last_k; int blevels[N_BLEVELS], n_blevels = 0, rlevels[N_RLEVELS], n_rlevels = 0, bin_trouble; int start_direction, stop_direction, lowest_level = 1, min_level = 0, n_pieces = 0; BOOLEAN error = FALSE, set_lake_fill = FALSE, draw_river = FALSE, shift = FALSE, need_coast_base, recursive; BOOLEAN greenwich = FALSE, draw_coast = FALSE, fill_ocean = FALSE, fill_land = FALSE, possibly_donut_hell = FALSE; BOOLEAN clobber_background, draw_border = FALSE, paint_polygons = FALSE, donut, fill_in_use = FALSE; BOOLEAN draw_scale = FALSE, gave_xy = FALSE, fancy = FALSE, donut_hell = FALSE, world_map_save, did_clip = FALSE; BOOLEAN clip_ocean = FALSE, clip_land = FALSE, end_of_clip = FALSE, clipping = FALSE, dumping = FALSE; double west = 0.0, east = 0.0, south = 0.0, north = 0.0, bin_x[5], bin_y[5]; double west_border, east_border, anti_lon, anti_lat, edge = 720.0, left, right; double *xtmp, *ytmp, min_area = 0.0, x0, y0, scale_lat, length, step = 0.01; double anti_x, anti_y, x_0, y_0, x_c, y_c, dist, out[2]; char res = 'l', key[5], *string, measure = '\0', comment[128]; char *shore_resolution[5] = {"full", "high", "intermediate", "low", "crude"}; struct GMT_PEN coast, pen, bpen[N_BLEVELS], rpen[N_RLEVELS]; struct GMT_FILL fill[5]; /* Colors for (0) water, (1) land, (2) lakes, (3) islands in lakes, (4) lakes in islands in lakes */ struct GMT_SHORE c; struct GMT_BR b, r; struct POL *p; GMT_init_fill (&fill[0], 255, 255, 255); /* Default Ocean color */ GMT_init_fill (&fill[1], 0, 0, 0); /* Default Land color */ GMT_init_fill (&fill[2], 255, 255, 255); /* Default Lake color */ fill[3] = fill[1]; fill[4] = fill[2]; GMT_init_pen (&coast, GMT_PENWIDTH); GMT_init_pen (&pen, GMT_PENWIDTH); for (k = 0; k < N_RLEVELS; k++) GMT_init_pen (&rpen[k], GMT_PENWIDTH); for (k = 0; k < N_BLEVELS; k++) GMT_init_pen (&bpen[k], GMT_PENWIDTH); memset ((void *)rlevels, 0, (size_t)(N_RLEVELS * sizeof (int))); memset ((void *)blevels, 0, (size_t)(N_BLEVELS * sizeof (int))); argc = GMT_begin (argc, argv); /* Check and interpret the command line arguments */ for (i = 1; i < argc; i++) { if (argv[i][0] == '-') { switch(argv[i][1]) { /* Common parameters */ case 'B': /* Get tickmark information */ case 'J': case 'K': case 'O': case 'P': case 'R': case 'U': case 'V': case 'X': case 'x': case 'Y': case 'y': case 'c': case '\0': error += GMT_get_common_args (argv[i], &west, &east, &south, &north); break; /* Supplemental parameters */ case 'b': error += GMT_io_selection (&argv[i][2]); break; case 'A': n = sscanf (&argv[i][2], "%lf/%d/%d", &min_area, &min_level, &max_level); if (n == 1) min_level = 0, max_level = MAX_LEVEL; break; case 'C': /* Lake color */ if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[2])) { GMT_fill_syntax ('C'); error++; } set_lake_fill = TRUE; fill[4] = fill[2]; break; case 'D': res = argv[i][2]; base = GMT_set_resolution (&res, 'D'); break; case 'E': sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation); break; case 'G': /* Set Gray shade, pattern, or clipping */ if (argv[i][2] == 'c') clip_land = TRUE; else { fill_land = TRUE; if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[1])) { GMT_fill_syntax ('G'); error++; } fill[3] = fill[1]; } break; case 'L': error += GMT_getscale (&argv[i][2], &x0, &y0, &scale_lat, &length, &measure, &fancy, &gave_xy); draw_scale = TRUE; break; case 'M': GMT_multisegment (&argv[i][2]); dumping = TRUE; break; case 'N': if (!argv[i][2]) { fprintf (stderr, "%s: GMT SYNTAX ERROR: -N option takes two arguments\n", GMT_program); error++; continue; } draw_border = TRUE; memset ((void *)key, 0, (size_t)5); for (k = 2; argv[i][k] && argv[i][k] != '/'; k++); if (argv[i][k] && argv[i][k] == '/') { /* Specified pen args */ strncpy (key, &argv[i][2], (size_t)k-2); string = &argv[i][k+1]; } else { /* No pen specified, use default */ strcpy (key, &argv[i][2]); string = CNULL; } GMT_init_pen (&pen, GMT_PENWIDTH); if (string && GMT_getpen (string, &pen)) { GMT_pen_syntax ('N'); error++; } switch (key[0]) { case 'a': for (k = 0; k < N_BLEVELS; k++) blevels[k] = TRUE; for (k = 0; k < N_BLEVELS; k++) bpen[k] = pen; break; default: k = atoi (key) - 1; if (k < 0 || k >= N_BLEVELS) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Feature not in list!\n", GMT_program); error++; } else { blevels[k] = TRUE; bpen[k] = pen; } break; } break; case 'I': if (!argv[i][2]) { fprintf (stderr, "%s: GMT SYNTAX ERROR: -I option takes two arguments\n", GMT_program); error++; continue; } draw_river = TRUE; memset ((void *)key, 0, (size_t)5); for (k = 2; argv[i][k] && argv[i][k] != '/'; k++); if (argv[i][k] && argv[i][k] == '/') { /* Specified pen args */ strncpy (key, &argv[i][2], (size_t)k-2); string = &argv[i][k+1]; } else { /* No pen specified, use default */ strcpy (key, &argv[i][2]); string = CNULL; } GMT_init_pen (&pen, GMT_PENWIDTH); if (string && GMT_getpen (string, &pen)) { GMT_pen_syntax ('I'); error++; } switch (key[0]) { case 'a': for (k = 0; k < N_RLEVELS; k++) rlevels[k] = TRUE; for (k = 0; k < N_RLEVELS; k++) rpen[k] = pen; break; case 'r': for (k = 0; k < 4; k++) rlevels[k] = TRUE; for (k = 0; k < 4; k++) rpen[k] = pen; break; case 'i': for (k = 4; k < 7; k++) rlevels[k] = TRUE; for (k = 4; k < 7; k++) rpen[k] = pen; break; case 'c': for (k = 7; k < 10; k++) rlevels[k] = TRUE; for (k = 7; k < 10; k++) rpen[k] = pen; break; default: k = atoi (key) - 1; if (k < 0 || k >= N_RLEVELS) { fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Feature not in list!\n", GMT_program); error++; } else { rlevels[k] = TRUE; rpen[k] = pen; } break; } break; case 'S': /* Set ocean color if needed */ if (argv[i][2] == 'c') clip_ocean = TRUE; else { if (argv[i][2] && GMT_getfill (&argv[i][2], &fill[0])) { GMT_fill_syntax ('S'); error++; } fill_ocean = TRUE; } break; case 'W': if (argv[i][2] && GMT_getpen (&argv[i][2], &coast)) { GMT_pen_syntax ('W'); error++; } draw_coast = TRUE; break; case 'Q': end_of_clip = TRUE; break; case '-': /* Undocumented option to increase path-fix resolution */ step = atof (&argv[i][2]); break; #ifdef DEBUG case '+': debug_bin = atoi (&argv[i][2]); break; #endif default: /* Options not recognized */ error = TRUE; GMT_default_error (argv[i][1]); break; } } } if (argc == 1 || GMT_quick) { /* Display usage */ fprintf (stderr,"pscoast %s - Plot continents, shorelines, rivers, and borders on maps\n\n", GMT_VERSION); fprintf (stderr,"usage: pscoast -J -R/// [-A[//]] [-B]\n"); fprintf (stderr, "\t[-C[]] [-D] [-Eaz/el] [-G[]] [-I[/]]\n"); fprintf (stderr, "\t[-K] [-L[f][x]///[m|n|k]] [-M[]] [-N[/]] [-O] [-P] [-Q]\n"); fprintf (stderr, "\t[-S] [-U[dx/dy/][label]] [-V] [-W[]] [-X] [-Y] [-bo[s][]] [-c]\n"); if (GMT_quick) exit (EXIT_FAILURE); GMT_explain_option ('j'); GMT_explain_option ('R'); fprintf (stderr, "\n\tOPTIONS:\n"); fprintf (stderr, "\t-A features smaller than (in km^2) or of levels (0-4) outside the min-max levels\n"); fprintf (stderr, "\t will be skipped [0/4 (4 means lake inside island inside lake)]\n"); GMT_explain_option ('b'); fprintf (stderr, "\t-C sets separate color for lakes. [Default is same as ocean (white)]. Set\n"); fprintf (stderr, "\t 1) (each 0-255) for color or (0-255) for gray-shade [0].\n"); fprintf (stderr, "\t 2) p[or P]/ for predefined patterns (1-90).\n"); fprintf (stderr, "\t-D Choose one of the following resolutions:\n"); fprintf (stderr, "\t f - full resolution (may be very slow for large regions)\n"); fprintf (stderr, "\t h - high resolution (may be slow for large regions)\n"); fprintf (stderr, "\t i - intermediate resolution\n"); fprintf (stderr, "\t l - low resolution [Default]\n"); fprintf (stderr, "\t c - crude resolution, for busy plots that need crude continent outlines only\n"); fprintf (stderr, "\t-E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n"); fprintf (stderr, "\t-G Paint or clip \"dry\" areas. Append desired effect as:\n"); fprintf (stderr, "\t 1) (each 0-255) for color or (0-255) for gray-shade [0].\n"); fprintf (stderr, "\t 2) p[or P]/ for predefined patterns (1-90).\n"); fprintf (stderr, "\t 3) c to issue clip paths for land areas.\n"); fprintf (stderr, "\t-I draws rIvers. Specify feature and optionally append pen [Default is solid line of unit thickness].\n"); fprintf (stderr, "\t Choose from the features below. Repeat the -I option as many times as needed\n"); fprintf (stderr, "\t 1 = Permanent major rivers\n"); fprintf (stderr, "\t 2 = Additional major rivers\n"); fprintf (stderr, "\t 3 = Additional rivers\n"); fprintf (stderr, "\t 4 = Minor rivers\n"); fprintf (stderr, "\t 5 = Intermittent rivers - major\n"); fprintf (stderr, "\t 6 = Intermittent rivers - additional\n"); fprintf (stderr, "\t 7 = Intermittent rivers - minor\n"); fprintf (stderr, "\t 8 = Major canals\n"); fprintf (stderr, "\t 9 = Minor canals\n"); fprintf (stderr, "\t 10 = Irrigation canals\n"); fprintf (stderr, "\t a = All rivers and canals (1-10)\n"); fprintf (stderr, "\t r = All permanent rivers (1-4)\n"); fprintf (stderr, "\t i = All intermittent rivers (5-7)\n"); fprintf (stderr, "\t c = All canals (8-10)\n"); GMT_explain_option ('K'); fprintf (stderr, "\t-L draws a simple map scaLe centered on /. Use -Lx to specify cartesian coordinates instead.\n"); fprintf (stderr, "\t Scale is calculated at latitude . is in km, or [nautical] miles if [n] m is appended\n"); fprintf (stderr, "\t -Lf draws a \"fancy\" scale [Default is plain]\n"); fprintf (stderr, "\t-M Dump a multisegment ascii (or binary, see -bo) file to standard output. No plotting occurs\n"); fprintf (stderr, "\t Specify any combination of -W, -I, -N. Optionally, append 1-char\n"); fprintf (stderr, "\t segment header [%c]\n", GMT_io.EOF_flag); fprintf (stderr, "\t-N draws boundaries. Specify feature and optionally append pen [Default is solid line of unit thickness]\n"); fprintf (stderr, "\t Choose from the features below. Repeat the -N option as many times as needed\n"); fprintf (stderr, "\t 1 = National boundaries\n"); fprintf (stderr, "\t 2 = State boundaries within the Americas\n"); fprintf (stderr, "\t 3 = Marine boundaries\n"); fprintf (stderr, "\t a = All boundaries (1-3)\n"); GMT_explain_option ('O'); GMT_explain_option ('P'); fprintf (stderr,"\t-Q means terminate previously set clip-paths\n"); fprintf (stderr, "\t-S Paint or clip \"wet\" areas. Append desired effect as:\n"); fprintf (stderr, "\t 1) (each 0-255) for color or (0-255) for gray-shade [0].\n"); fprintf (stderr, "\t 2) p[or P]/ for predefined patterns (1-90).\n"); fprintf (stderr, "\t 3) c to issue clip paths for water areas.\n"); GMT_explain_option ('U'); GMT_explain_option ('V'); fprintf (stderr,"\t-W draw shorelines. Append pen [Default is solid line of unit thickness].\n"); GMT_explain_option ('X'); GMT_explain_option ('c'); GMT_explain_option ('o'); GMT_explain_option ('.'); exit (EXIT_FAILURE); } if (set_lake_fill && !(fill_land || fill_ocean || draw_coast)) { /* Just lakes, fix -A */ if (min_level < 2) min_level = 2; } /* Check that the options selected are mutually consistant */ clipping = (clip_land || clip_ocean); if (!(end_of_clip || project_info.region_supplied)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify -R option\n", GMT_program); error++; } if (coast.width < 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -W option: Pen thickness cannot be negative\n", GMT_program); error++; } if (!(fill_land || fill_ocean || set_lake_fill || draw_coast || draw_border || draw_river || clipping || end_of_clip)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must specify at least one of -C, -G, -S, -I, -N, and -W\n", GMT_program); error++; } if ((fill_land + fill_ocean + set_lake_fill) > 0 && (clip_land + clip_ocean) > 0) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must choose between clipping and painting\n", GMT_program); error++; } if (clip_land && clip_ocean) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must choose between clipping land OR water\n", GMT_program); error++; } if (dumping && (clipping + fill_land + fill_ocean + set_lake_fill)) { fprintf (stderr, "%s: GMT SYNTAX ERROR: Must choose between dumping and clipping/plotting\n", GMT_program); error++; } if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -E option: Enter azimuth in 0-360 range, elevation in 0-90 range\n", GMT_program); error++; } if (clipping && project_info.projection == AZ_EQDIST && fabs (project_info.w - project_info.e) == 360.0 && (project_info.n - project_info.s) == 180.0) { fprintf (stderr, "%s: -JE not implemented for global clipping - I quit\n", GMT_program); exit (EXIT_FAILURE); } if (clipping && (draw_border || draw_river || draw_coast)) { fprintf (stderr, "%s: Cannot do clipping AND draw coastlines, rivers, or borders\n", GMT_program); exit (EXIT_FAILURE); } if (dumping && !(draw_border || draw_river || draw_coast)) { fprintf (stderr, "%s: Must specify at least one of -I, -N, and -W\n", GMT_program); exit (EXIT_FAILURE); } if (draw_river) { for (k = 0; k < N_RLEVELS; k++) { if (!rlevels[k]) continue; if (rpen[k].width < 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -I option: Pen thickness cannot be negative\n", GMT_program); error++; } rlevels[n_rlevels] = k + 1; rpen[n_rlevels] = rpen[k]; n_rlevels++; } } if (draw_border) { for (k = 0; k < N_BLEVELS; k++) { if (!blevels[k]) continue; if (bpen[k].width < 0.0) { fprintf (stderr, "%s: GMT SYNTAX ERROR -N option: Pen thickness cannot be negative\n", GMT_program); error++; } blevels[n_blevels] = k + 1; bpen[n_blevels] = bpen[k]; n_blevels++; } } if (error) exit (EXIT_FAILURE); GMT_put_history (argc, argv); /* Update .gmtcommands */ need_coast_base = (fill_land || fill_ocean || set_lake_fill || draw_coast || clip_land || clip_ocean); clobber_background = (fill_land && fill_ocean); recursive = (((fill_land + (fill_ocean || set_lake_fill)) == 1) || clipping); paint_polygons = (fill_land || fill_ocean || set_lake_fill) || clipping; if (east > 360.0) { west -= 360.0; east -= 360.0; } if (end_of_clip && project_info.projection == -1) { /* Fake area and linear projection */ GMT_map_getproject ("x1d"); east = north = 1.0; } else if (dumping) { GMT_map_getproject ("x1d"); } GMT_map_setup (west, east, south, north); world_map_save = GMT_world_map; if (need_coast_base && GMT_init_shore (res, &c, project_info.w, project_info.e, project_info.s, project_info.n)) { fprintf (stderr, "%s: %s resolution shoreline data base not installed\n", GMT_program, shore_resolution[base]); need_coast_base = FALSE; } if (draw_border && GMT_init_br ('b', res, &b, project_info.w, project_info.e, project_info.s, project_info.n)) { fprintf (stderr, "%s: %s resolution political boundary data base not installed\n", GMT_program, shore_resolution[base]); draw_border = FALSE; } if (draw_river && GMT_init_br ('r', res, &r, project_info.w, project_info.e, project_info.s, project_info.n)) { fprintf (stderr, "%s: %s resolution river data base not installed\n", GMT_program, shore_resolution[base]); draw_river = FALSE; } if (! (need_coast_base || draw_border || draw_river || end_of_clip)) { fprintf (stderr, "%s: No databases available - aborts\n", GMT_program); exit (EXIT_FAILURE); } if (dumping) { if (!GMT_io.binary[1]) fprintf (GMT_stdout, "# Data from the %s resolution GMT shoreline, borders, and rivers database\n", shore_resolution[base]); } else { ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin, gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi, GMT_INCH, gmtdefs.paper_width, gmtdefs.page_rgb, GMT_epsinfo(argv[0])); GMT_echo_command (argc, argv); if (gmtdefs.unix_time) GMT_timestamp (argc, argv); if (!end_of_clip && base > 0) ps_command ("2 setlinejoin"); /* Avoid excessive mitering for cruder lines */ if (end_of_clip) { /* Just undo previous clip-path */ ps_clipoff (); if (frame_info.plot) GMT_map_basemap (); /* Basemap needed */ ps_plotend (gmtdefs.last_page); if (gmtdefs.verbose) fprintf (stderr, "%s: Done!\n", GMT_program); GMT_end (argc, argv); } if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0); } if (fill_ocean && !set_lake_fill) fill[2] = fill[4] = fill[0]; /* Use ocean fill by default for lakes */ for (i = 0; i < 5; i++) if (fill[i].use_pattern) { /* Must initialize patterns */ ps_imagefill_init (fill[i].pattern_no, fill[i].pattern, fill[i].inverse, fill[i].dpi, fill[i].colorize, fill[i].f_rgb, fill[i].b_rgb); fill_in_use = TRUE; } if (fill_in_use && !clobber_background) { /* Force clobber when fill is used since our routine cannot deal with clipped fills */ if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: Pattern fill requires oceans to be painted first\n", GMT_program); clobber_background = TRUE; recursive = FALSE; } if (clipping && frame_info.plot) { /* Want basemap drawn before clipping is activated */ /* project_info.w = west_border; project_info.e = east_border; */ GMT_map_basemap (); } if (project_info.projection == AZ_EQDIST && fabs (project_info.w - project_info.e) == 360.0 && (project_info.n - project_info.s) == 180.0) { possibly_donut_hell = TRUE; anti_lon = project_info.central_meridian + 180.0; if (anti_lon >= 360.0) anti_lon -= 360.0; anti_lat = -project_info.pole; anti_bin = (int)floor ((90.0 - anti_lat) / c.bsize) * c.bin_nx + (int)floor (anti_lon / c.bsize); GMT_geo_to_xy (anti_lon, anti_lat, &anti_x, &anti_y); GMT_geo_to_xy (project_info.central_meridian, project_info.pole, &x_0, &y_0); if (fill_land && gmtdefs.verbose) { fprintf (stderr, "%s: Warning: Fill continent option (-G) may not work for this projection.\n", GMT_program); fprintf (stderr, "If the antipole (%lg/%lg) is in the ocean then chances are good\n", anti_lon, anti_lat); fprintf (stderr, "Else: avoid projection center coordinates that are exact multiples of %lg degrees\n", c.bsize); } } if (possibly_donut_hell && !clobber_background) { /* Force clobber when donuts may be called for now */ if (gmtdefs.verbose) fprintf (stderr, "%s: Warning: -JE requires oceans to be painted first\n", GMT_program); clobber_background = TRUE; recursive = FALSE; } if (clobber_background) { /* Paint entire map as ocean first, then lay land on top */ n = GMT_map_clip_path (&xtmp, &ytmp, &donut); if (donut) { ps_line (xtmp, ytmp, n, 1, TRUE, FALSE); ps_line (&xtmp[n], &ytmp[n], n, -1, TRUE, FALSE); GMT_fill (xtmp, ytmp, n, &fill[0], -1); } else GMT_fill (xtmp, ytmp, n, &fill[0], FALSE); GMT_free ((void *)xtmp); GMT_free ((void *)ytmp); fill_ocean = FALSE; } if (clipping) { start_direction = (clip_land) ? 1 : -1; stop_direction = (clip_land) ? 1 : -1; level_to_be_painted = (clip_land) ? 1 : 0; } else if (recursive) { start_direction = (fill_land) ? 1 : -1; stop_direction = (fill_land) ? 1 : -1; level_to_be_painted = (fill_land) ? 1 : 0; } else { start_direction = -1; stop_direction = 1; lowest_level = (fill_ocean) ? 0 : 1; } if (west < 0.0 && east > 0.0 && project_info.projection > 5) greenwich = TRUE; if ((360.0 - fabs (project_info.e - project_info.w) ) < c.bsize) { GMT_world_map = TRUE; if (project_info.projection == GNOMONIC) GMT_world_map = FALSE; } if (GMT_world_map && greenwich) edge = project_info.central_meridian; else if (!GMT_world_map && greenwich) { shift = TRUE; edge = west; if (edge < 0.0) edge += 360.0; if (edge > 180.0) edge = 180.0; } /* west_border = project_info.w; east_border = project_info.e; */ if (project_info.w < 0.0 && project_info.e <= 0.0) { /* Temporarily shift boundaries */ project_info.w += 360.0; project_info.e += 360.0; if (project_info.central_meridian < 0.0) project_info.central_meridian += 360.0; } /* west_border = project_info.w; east_border = project_info.e; */ west_border = floor (project_info.w / c.bsize) * c.bsize; east_border = ceil (project_info.e / c.bsize) * c.bsize; if (!dumping && draw_coast) GMT_setpen (&coast); for (ind = 0; need_coast_base && ind < c.nb; ind++) { /* Loop over necessary bins only */ bin = c.bins[ind]; #ifdef DEBUG if (debug_bin >= 0 && bin != debug_bin) continue; #endif GMT_get_shore_bin (ind, &c, min_area, min_level, max_level); if (gmtdefs.verbose) fprintf (stderr, "%s: Working on block # %5d\r", GMT_program, bin); if (!dumping) { sprintf (comment, "Bin # %d", bin); ps_comment (comment); } if (GMT_world_map && 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) { 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; GMT_geo_to_xy (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); donut_hell = GMT_non_zero_winding (anti_lon, anti_lat, bin_x, bin_y, 5); if (!donut_hell) donut_hell = (dist > 0.8 * project_info.r); } for (direction = start_direction; paint_polygons && direction <= stop_direction; direction += 2) { /* Assemble one or more segments into polygons */ /* if ((np = GMT_assemble_shore (&c, direction, TRUE, shift, edge, &p)) == 0) continue; */ if ((np = GMT_assemble_shore (&c, direction, min_level, TRUE, shift, west_border, east_border, &p)) == 0) continue; /* Get clipped polygons in x,y inches that can be plotted */ np_new = GMT_prep_polygons (&p, np, greenwich, donut_hell, step, bin_trouble); if (clipping) { for (k = level_to_be_painted; k < MAX_LEVEL - 1; k++) { recursive_clipping (-1, np_new, p, k, k, &n_pieces); did_clip = TRUE; } 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) { flag = (n_pieces == 0) ? 1 : 0; n_pieces++; ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag); } } } else if (!recursive) { /* Simply points all polygons as is */ for (k = 0; k < np_new; k++) { if (p[k].n == 0 || p[k].level < lowest_level) continue; level = p[k].level; if (donut_hell && GMT_non_zero_winding (anti_x, anti_y, p[k].lon, p[k].lat, p[k].n)) { /* Antipode inside polygon, must do donut */ n = GMT_map_clip_path (&xtmp, &ytmp, &donut); ps_line (xtmp, ytmp, n, 1, TRUE, FALSE); ps_line (p[k].lon, p[k].lat, p[k].n, -1, TRUE, FALSE); GMT_fill (xtmp, ytmp, n, &fill[level], -1); GMT_free ((void *)xtmp); GMT_free ((void *)ytmp); } else GMT_fill (p[k].lon, p[k].lat, p[k].n, &fill[level], FALSE); } } else { /* Must avoid pointing anything but the polygons inside */ for (k = level_to_be_painted; k < MAX_LEVEL - 1; k++) { recursive_painting (-1, np_new, p, k, 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_fill (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], FALSE); } } GMT_free_polygons (p, np_new); GMT_free ((void *)p); } if (draw_coast && c.ns) { /* Draw or dump shorelines, no need to assemble polygons */ /* if ((np = GMT_assemble_shore (&c, direction, FALSE, shift, edge, &p)) == 0) continue; */ if ((np = GMT_assemble_shore (&c, direction, min_level, FALSE, shift, west_border, east_border, &p)) == 0) continue; for (i = 0; i < np; i++) { if (dumping) { sprintf (GMT_io.segment_header, "%c Shore Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level); GMT_write_segmentheader (GMT_stdout, 2); for (k = 0; k < p[i].n; k++) { out[0] = p[i].lon[k]; out[1] = p[i].lat[k]; GMT_output (GMT_stdout, 2, out); } } else { if (donut_hell) p[i].n = GMT_fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step); GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n); if (GMT_n_plot) GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot); } } GMT_free_polygons (p, np); GMT_free ((void *)p); } GMT_free_shore (&c); } if (need_coast_base) GMT_shore_cleanup (&c); if (clipping && did_clip) ps_clipon (p[k].lon, p[k].lat, 0, GMT_no_rgb, 2); /* Dummy to end path */ if (gmtdefs.verbose) fprintf (stderr, "\n"); if (draw_river) { /* Read rivers file and plot as lines */ if (gmtdefs.verbose) fprintf (stderr, "%s: Adding Rivers...", GMT_program); if (!dumping) ps_comment ("Start of River segments"); last_k = -1; for (ind = 0; ind < r.nb; ind++) { /* Loop over necessary bins only */ bin = r.bins[ind]; GMT_get_br_bin (ind, &r, rlevels, n_rlevels); if (r.ns == 0) continue; if (GMT_world_map && 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 (&r, shift, edge, &p)) == 0) continue; for (i = 0; i < np; i++) { if (dumping) { sprintf (GMT_io.segment_header, "%c River Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level); GMT_write_segmentheader (GMT_stdout, 2); for (k = 0; k < p[i].n; k++) { out[0] = p[i].lon[k]; out[1] = p[i].lat[k]; GMT_output (GMT_stdout, 2, out); } } else { GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n); if (!GMT_n_plot) continue; k = p[i].level - 1; if (k != last_k) { GMT_setpen (&rpen[k]); last_k = k; } GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot); } } /* Free up memory */ GMT_free_br (&r); for (k = 0; k < np; k++) { GMT_free ((void *)p[k].lon); GMT_free ((void *)p[k].lat); } GMT_free ((void *)p); } GMT_br_cleanup (&r); if (gmtdefs.verbose) fprintf (stderr, "\n"); } if (draw_border) { /* Read borders file and plot as lines */ if (gmtdefs.verbose) fprintf (stderr, "%s: Adding Borders...", GMT_program); if (!dumping) ps_comment ("Start of Border segments"); /* Must resample borders because some points may be too far apart and look like 'jumps' */ step = MAX (fabs(project_info.w - project_info.e), fabs (project_info.n - project_info.s)) / 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 (ind, &b, blevels, n_blevels); if (b.ns == 0) continue; if (GMT_world_map && 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 (&b, shift, edge, &p)) == 0) continue; for (i = 0; i < np; i++) { if (dumping) { sprintf (GMT_io.segment_header, "%c Border Bin # %d, Level %d\n", GMT_io.EOF_flag, bin, p[i].level); GMT_write_segmentheader (GMT_stdout, 2); for (k = 0; k < p[i].n; k++) { out[0] = p[i].lon[k]; out[1] = p[i].lat[k]; GMT_output (GMT_stdout, 2, out); } } else { p[i].n = GMT_fix_up_path (&p[i].lon, &p[i].lat, p[i].n, greenwich, step); GMT_n_plot = GMT_geo_to_xy_line (p[i].lon, p[i].lat, p[i].n); if (!GMT_n_plot) continue; k = p[i].level - 1; if (k != last_k) { GMT_setpen (&bpen[k]); last_k = k; } GMT_plot_line (GMT_x_plot, GMT_y_plot, GMT_pen, GMT_n_plot); } } /* Free up memory */ GMT_free_br (&b); for (k = 0; k < np; k++) { GMT_free ((void *)p[k].lon); GMT_free ((void *)p[k].lat); } GMT_free ((void *)p); } GMT_br_cleanup (&b); if (gmtdefs.verbose) fprintf (stderr, "\n"); } if (!dumping) { if (base > 0) ps_command ("0 setlinejoin"); /* Reset mitering */ if (frame_info.plot) { /* project_info.w = west_border; project_info.e = east_border; */ GMT_world_map = world_map_save; GMT_map_basemap (); } if (draw_scale) GMT_draw_map_scale (x0, y0, scale_lat, length, measure, gave_xy, fancy); if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0); ps_plotend (gmtdefs.last_page); } GMT_end (argc, argv); } void recursive_painting (int k0, int np, struct POL p[], int level, int start_level, struct GMT_FILL fill[]) { int k, flag; BOOLEAN go; char text[64]; if (level > MAX_LEVEL) return; flag = (level == start_level) ? 1 : 0; for (k = k0 + 1; k < np; k++) { if (p[k].n == 0 || p[k].level < level) continue; go = (p[k].level == level) ? ((k0 == -1) ? TRUE : GMT_non_zero_winding (p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) : FALSE; if (go) { sprintf (text, "Polygon %d", k); ps_comment (text); ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag); recursive_painting (k, np, p, level+1, start_level, fill); p[k].n = 0; } if (p[k].level == start_level) { /* Done nesting, time to paint */ ps_clipon (p[k].lon, p[k].lat, 0, GMT_no_rgb, 4); /* Dummy to end path */ GMT_fill (p[k].lon, p[k].lat, p[k].n, &fill[p[k].level], -1); } } } void recursive_clipping (int k0, int np, struct POL *p, int level, int start_level, int *n_pieces) { int k, flag; BOOLEAN go; char text[64]; if (level > MAX_LEVEL) return; for (k = k0 + 1; k < np; k++) { if (p[k].n == 0 || p[k].level < level) continue; go = (p[k].level == level) ? ((k0 == -1) ? TRUE : GMT_non_zero_winding (p[k].lon[0], p[k].lat[0], p[k0].lon, p[k0].lat, p[k0].n)) : FALSE; if (go) { sprintf (text, "Polygon %d", k); ps_comment (text); flag = (*n_pieces == 0) ? 1 : 0; (*n_pieces)++; ps_clipon (p[k].lon, p[k].lat, p[k].n, GMT_no_rgb, flag); recursive_clipping (k, np, p, level+1, start_level, n_pieces); p[k].n = 0; } } }