home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdcontour.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-22  |  29.8 KB  |  833 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdcontour.c    2.46  22 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.  * grdcontour reads a 2-D grd file and contours it. This algoritm handles
  9.  * cases where the old contourm failed, e.g. when a contour line passes
  10.  * exactly through a data point. This is achieved by adding a tiny
  11.  * number to those data points that would have a contour line going
  12.  * through them. This is cheating, but the human eye cannot tell the
  13.  * difference, so who cares?
  14.  * A multitude of options exist and are explained in the usage message.
  15.  * The 2-D grd file format is outlined in the routines read_grd/write_grd.
  16.  *
  17.  * Author:    Paul Wessel
  18.  * Date:    6-JAN-1991-1995
  19.  * Version:    2.0    Based on old v1.x by Paul Wessel & Walter Smith
  20.  *
  21.  */
  22.  
  23. #include "gmt.h"
  24.  
  25. #define MIN_LENGTH 0.01
  26.  
  27. float *grd;
  28. int *edge;
  29. double *x, *y;    /* Arrays holding the contour xy values */
  30.  
  31. struct LABEL {
  32.     double x, y;
  33.     double angle;
  34.     char label[32];
  35.     struct LABEL *next_label, *prev_label;
  36. } *anchor, *old_label;
  37.  
  38. struct SAVE {
  39.     double *x, *y;
  40.     double cval;
  41.     int n, pen;
  42.     BOOLEAN do_it, high;
  43. } *save;
  44.  
  45. main (argc, argv)
  46. int argc;
  47. char **argv; {
  48.     int i, j, n, nm, c, n_edges, label_font_size = 9, br, bg, bb, n_alloc;
  49.     int section, n_contours, smooth_factor = 0, side, id, box = 1, bad;
  50.     int dummy[4], off, nx, ny, half = 5, n_save = 0, got, ncut = 0;
  51.     
  52.     BOOLEAN error = FALSE, first = TRUE, dump = FALSE, anot = FALSE, interior, begin;
  53.     BOOLEAN fix_angle = FALSE, do_ticks = FALSE, multi_segments = FALSE;
  54.     BOOLEAN tick_high = FALSE, tick_low = FALSE, cpt_given = FALSE, tick_label = FALSE;
  55.     
  56.     char *grdfile = 0, dfile[512], line[512], *cont_type, *cont_do_tick, EOL_flag = '>';
  57.     char *cpt_file, *labels;
  58.     
  59.     double c_int = 0.0, c_low, c_high, take_out, aval, a_int = 0.0, tick_gap = 0.2;
  60.     double anot_dist, west = 0.0, east = 0.0, south = 0.0, north = 0.0, tick_length = 0.05;
  61.     double *contour, cval, mult = 1.0, min, max, small, label_angle = 0.0;
  62.     double small_x, small_y, data_west, data_east, data_south, data_north, tmp, *cont_angle;
  63.     
  64.     FILE *fpc = NULL;
  65.     
  66.     struct GRD_HEADER header;
  67.     
  68.     struct PEN pen[2];
  69.  
  70.     argc = gmt_begin (argc, argv);
  71.     
  72.     gmt_init_pen (&pen[0], 1);
  73.     gmt_init_pen (&pen[1], 3);
  74.     c_low = -MAXDOUBLE;
  75.     c_high = MAXDOUBLE;
  76.     dfile[0] = 0;
  77.     
  78.     anot_dist = 4.0; /* Distance in inches between anotations */
  79.     if (!gmtdefs.measure_unit) {    /* change to cm */
  80.         anot_dist = 10.0;
  81.         tick_gap = 0.5;
  82.         tick_length = 0.1;
  83.     }
  84.     br = gmtdefs.page_rgb[0];    /* Default box color is page color */
  85.     bg = gmtdefs.page_rgb[1];
  86.     bb = gmtdefs.page_rgb[2];
  87.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  88.     
  89.     for (i = 1; i < argc; i++) {
  90.         if (argv[i][0] == '-') {
  91.             switch (argv[i][1]) {
  92.         
  93.                 /* Common parameters */
  94.             
  95.                 case 'B':
  96.                 case 'J':
  97.                 case 'K':
  98.                 case 'O':
  99.                 case 'P':
  100.                 case 'R':
  101.                 case 'U':
  102.                 case 'V':
  103.                 case 'X':
  104.                 case 'x':
  105.                 case 'Y':
  106.                 case 'y':
  107.                 case 'c':
  108.                 case '\0':
  109.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  110.                     break;
  111.                 
  112.                 /* Supplemental parameters */
  113.             
  114.                 case 'A':
  115.                     n = sscanf (&argv[i][2], "%lf", &a_int);
  116.                     
  117.                     for (j = 2, bad = 0; argv[i][j] && argv[i][j] != 'f'; j++);
  118.                     if (argv[i][j])    { /* Found font size option */
  119.                         label_font_size = atoi (&argv[i][j+1]);
  120.                         if (label_font_size <= 0) bad++;
  121.                     }
  122.                         
  123.                     for (j = 2; argv[i][j] && argv[i][j] != 'a'; j++);
  124.                     if (argv[i][j])    { /* Found fixed angle option */
  125.                         label_angle = atof (&argv[i][j+1]);
  126.                         fix_angle = TRUE;
  127.                         if (label_angle <= -90.0 || label_angle > 180.0) bad++;
  128.                     }
  129.                         
  130.                     for (j = 2; argv[i][j] && argv[i][j] != '/'; j++);
  131.                     if (argv[i][j] && gmt_getrgb (&argv[i][j+1], &br, &bg, &bb)) bad++;
  132.                     anot = TRUE;
  133.                     if (strchr (argv[i], 'o')) box = 2;
  134.                     if (bad) {
  135.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -A option.  Correct syntax:\n", gmt_program);
  136.                         fprintf (stderr, "\t-A[aint][f<fontsize>][a<fixedangle>][/<r/g/b>][o]\n");
  137.                         error += bad;
  138.                     }
  139.                     break;
  140.                 case 'C':
  141.                     if ((fpc = fopen (&argv[i][2], "r")) == NULL)
  142.                         c_int = atof (&argv[i][2]);
  143.                     else {
  144.                         c_int = 1.0;
  145.                         cpt_given = (int) strstr (&argv[i][2], ".cpt");
  146.                         cpt_file = &argv[i][2];
  147.                     }
  148.                     break;
  149.                 case 'D':
  150.                     dump = TRUE;
  151.                     strcpy (dfile, &argv[i][2]);
  152.                     break;
  153.                 case 'E':
  154.                     sscanf (&argv[i][2], "%lf/%lf", &z_project.view_azimuth, &z_project.view_elevation);
  155.                     break;
  156.                 case 'F':
  157.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  158.                         gmt_pen_syntax ('F');
  159.                         error++;
  160.                     }
  161.                     break;
  162.                 case 'G':
  163.                     sscanf (&argv[i][2], "%lf/%d", &anot_dist, &half);
  164.                     half /= 2;
  165.                     break;
  166.                 case 'L':
  167.                     sscanf (&argv[i][2], "%lf/%lf", &c_low, &c_high);
  168.                     break;
  169.                 case 'M':    /* with -D, create one multiple line segments */
  170.                     multi_segments = 2;
  171.                     if (argv[i][2]) EOL_flag = argv[i][2];
  172.                     break;
  173.                 case 'S':
  174.                     smooth_factor = atoi (&argv[i][2]);
  175.                     break;
  176.                 case 'Q':
  177.                     ncut = atoi (&argv[i][2]);
  178.                     break;
  179.                 case 'T':
  180.                     if (argv[i][2]) {
  181.                         if (argv[i][2] == '+')
  182.                             tick_high = TRUE, j = 1;
  183.                         else if (argv[i][2] == '-')
  184.                             tick_low = TRUE, j = 1;
  185.                         else
  186.                             tick_high = tick_low = TRUE, j = 0;
  187.                         n = sscanf (&argv[i][2+j], "%lf/%lf", &tick_gap, &tick_length);
  188.                         if (n != 2 || tick_gap <= 0.0 || tick_length == 0.0) {
  189.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -T option.  Correct syntax:\n", gmt_program);
  190.                             fprintf (stderr, "\t-T[+|-]<tick_gap>/<tick_length>[:LH], <tick_gap> must be > 0\n");
  191.                             error = TRUE;
  192.                         }
  193.                         do_ticks = TRUE;
  194.                         for (j = 2; argv[i][j] && argv[i][j] != ':'; j++);
  195.                         if (argv[i][j]) {
  196.                             j++;
  197.                             tick_label = TRUE;
  198.                             labels = &argv[i][j];
  199.                         }
  200.                             
  201.                     }
  202.                     break;
  203.                 case 'W':
  204.                     j = (argv[i][2] == 'a' || argv[i][2] == 'c') ? 3 : 2;
  205.                     if (j == 2) {    /* Set both */
  206.                         if (gmt_getpen (&argv[i][j], &pen[0])) {
  207.                             gmt_pen_syntax ('W');
  208.                             error++;
  209.                         }
  210.                         else pen[1] = pen[0];
  211.                     }
  212.                     else {
  213.                         id = (argv[i][2] == 'a') ? 1 : 0;
  214.                         if (gmt_getpen (&argv[i][j], &pen[id])) {
  215.                             gmt_pen_syntax ('W');
  216.                             error++;
  217.                         }
  218.                     }
  219.                     break;
  220.                 case 'Z':
  221.                     mult = atof (&argv[i][2]);
  222.                     break;
  223.                 default:
  224.                     error++;
  225.                     gmt_default_error (argv[i][1]);
  226.                     break;
  227.             }
  228.         }
  229.         else
  230.             grdfile = argv[i];
  231.     }
  232.     
  233.     if (argc == 1 || gmt_quick) {
  234.         fprintf (stderr,"grdcontour %s - Contouring of 2-D gridded data sets\n\n", GMT_VERSION);
  235.         fprintf (stderr, "usage: grdcontour <grdfile> -C<cont_int> -J<params> [-A[<anot_int>][f<fontsize>][/r/g/b][a<angle>][o]]\n");
  236.         fprintf (stderr, "    [-B<tickinfo>] [-D<dumpfile>] [-Eaz/el] [-F<r/g/b>] [-G<gap>/<npoints>] [-K] [-L<Low/high>]\n");
  237.         fprintf (stderr, "    [-M[<flag>]] [-O] [-P] [-Q<cut>] [-R<west/east/south/north>] [-S<smooth>] [-T[+|-][<gap/length>][:LH]] [-U[<label>]] [-V]\n");
  238.         fprintf (stderr, "    [-W<type><pen>] [-X<x_shift>] [-Y<y_shift>] [-Z<fact>] [-c<ncopies>]\n\n");
  239.         
  240.         if (gmt_quick) exit (-1);
  241.         
  242.         fprintf (stderr, "    <grdfile> is 2-D netCDF grdfile to be contoured\n");
  243.         fprintf (stderr, "    -C Contours to be drawn can be specified in one of three ways:\n");
  244.         fprintf (stderr, "       -Fixed contour interval\n");
  245.         fprintf (stderr, "       -Name of file with contour levels in col 1 and C(ont) or A(not) in col 2\n");
  246.         fprintf (stderr, "          [and optionally an individual annotation angle in col 3.]\n");
  247.         fprintf (stderr, "       -Name of cpt-file\n");
  248.         fprintf (stderr, "       If -T is used, only contours with upper case C or A is ticked\n");
  249.         fprintf (stderr, "         [cpt-file contours are set to C unless last column has flags; Use -A to force all to A]\n");
  250.         explain_option ('j');
  251.         fprintf (stderr, "\n\tOPTIONS:\n");
  252.         fprintf (stderr, "    -A Annotation information. [Default is no annotation].\n");
  253.         fprintf (stderr, "       Append f followed by desired font size in points [Default is 9].\n");
  254.         fprintf (stderr, "       Append /r/g/b to change color of text box [Default is %d/%d/%d]\n", br, bg, bb);
  255.         fprintf (stderr, "       Append o to draw outline of text box [Default is no outline]\n");
  256.         fprintf (stderr, "       Append a<angle> to force anotations at this fixed angle [Default follows contour]\n");
  257.         explain_option ('b');
  258.         fprintf (stderr, "    -D to Dump contour lines to individual files (but see -M)\n");
  259.         fprintf (stderr, "    -E set azimuth and elevation of viewpoint for 3-D perspective [180/90]\n");
  260.         fprintf (stderr, "    -F Set color used for Frame and anotation [%d/%d/%d]\n",
  261.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  262.         fprintf (stderr, "    -G Gap between anotations in %s and number of (x,y) values used to find L2 label angle [Default = %lg/10]\n", gmt_unit_names[gmtdefs.measure_unit], anot_dist);
  263.         explain_option ('K');
  264.         fprintf (stderr, "    -L only contour inside this range\n");
  265.         fprintf (stderr, "      -M Used with -D.   Create a single multiple segment file where contours are separated by a record\n");
  266.         fprintf (stderr, "         whose first character is <flag> ['>'].  This header also has the contour level value\n");
  267.         explain_option ('O');
  268.         explain_option ('P');
  269.         fprintf (stderr, "    -Q Do not draw contours with less than <cut> points [Draw all contours]\n");
  270.         explain_option ('R');
  271.         fprintf (stderr, "       [Default is extent of grdfile]\n");
  272.         fprintf (stderr, "    -S will Smooth contours by splining and resampling\n");
  273.         fprintf (stderr, "       at approximately gridsize/<smooth> intervals\n");
  274.         fprintf (stderr, "    -T will embellish innermost, closed contours with ticks pointing in the downward direction\n");
  275.         fprintf (stderr, "       User may specify to tick only highs (-T+) or lows (-T-) [-T means both]\n");
  276.         fprintf (stderr, "       Append spacing/ticklength (in %s) to change defaults [%lg/%lg]\n", gmt_unit_names[gmtdefs.measure_unit], tick_gap, tick_length);
  277.         fprintf (stderr, "       Append :LH to plot the characters L and H in the center of closed contours\n");
  278.         fprintf (stderr, "       for local Lows and Highs (e.g, give :-+ to plot - and + signs)\n");
  279.         explain_option ('U');
  280.         explain_option ('V');
  281.         fprintf (stderr, "    -W sets pen attributes. <type> can be 'a' for anotated contours and\n");
  282.         fprintf (stderr, "       'c' for regular contours.  The default settings are\n");
  283.         fprintf (stderr, "       Contour pen:  width = %d, color = (%d/%d/%d), texture = solid\n", pen[0].width, pen[0].r, pen[0].g, pen[0].b);
  284.         fprintf (stderr, "       Annotate pen: width = %d, color = (%d/%d/%d), texture = solid\n", pen[1].width, pen[1].r, pen[1].g, pen[1].b);
  285.         explain_option ('X');
  286.         fprintf (stderr, "    -Z to multiply data by <fact> before contouring.\n");
  287.         explain_option ('c');
  288.         exit(-1);
  289.     }
  290.     
  291.     if (a_int > 0.0 && (!fpc && c_int == 0.0)) c_int = a_int;
  292.     if (!fpc && c_int <= 0.0) {
  293.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify contour interval, file name with levels, or cpt-file\n", gmt_program);
  294.         error++;
  295.     }
  296.     if (c_low >= c_high) {
  297.         fprintf (stderr, "%s: GMT SYNTAX ERROR -L option: lower limit >= upper!\n", gmt_program);
  298.         error++;
  299.     }
  300.     if (smooth_factor < 0) {
  301.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Smooth_factor must be > 0\n", gmt_program);
  302.         error++;
  303.     }
  304.     if (z_project.view_azimuth > 360.0 || z_project.view_elevation <= 0.0 || z_project.view_elevation > 90.0) {
  305.         fprintf (stderr, "%s: GMT SYNTAX ERROR -E option:  Enter azimuth in 0-360 range, elevation in 0-90 range\n", gmt_program);
  306.         error++;
  307.     }
  308.     if (!grdfile) {
  309.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  310.         error++;
  311.     }
  312.     if (anot_dist <= 0.0 || half <= 0) {
  313.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G option.  Correct syntax:\n", gmt_program);
  314.         fprintf (stderr, "\t-G<anot_dist>/<npoints>, both values must be > 0\n");
  315.         error++;
  316.     }
  317.     if (mult == 0.0) {
  318.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  factor must be nonzero\n", gmt_program);
  319.         error++;
  320.     }
  321.     
  322.     if (error) exit (-1);
  323.  
  324.     if (dump && dfile[0] == 0) strcpy (dfile,"contour");
  325.     
  326.     if (gmtdefs.verbose) fprintf (stderr, "grdcontour: Allocate memory and read data file\n");
  327.     
  328.     if (!strcmp (grdfile,  "=")) {
  329.         fprintf (stderr, "grdcontour: Piping of grdfile not supported!\n");
  330.         exit (-1);
  331.     }
  332.     
  333.     if (read_grd_info (grdfile, &header)) {
  334.         fprintf (stderr, "grdcontour: Error opening file %s\n", grdfile);
  335.         exit (-1);
  336.     }
  337.     off = (header.node_offset) ? 0 : 1;
  338.     
  339.     /* Determine what wesn to pass to map_setup */
  340.  
  341.     if (west == east && south == north) {
  342.         west = header.x_min;
  343.         east = header.x_max;
  344.         south = header.y_min;
  345.         north = header.y_max;
  346.     }
  347.  
  348.     map_setup (west, east, south, north);
  349.  
  350.     /* Determine the wesn to be used to read the grdfile */
  351.  
  352.     grd_setregion (&header, &data_west, &data_east, &data_south, &data_north);
  353.  
  354.     /* Read data */
  355.  
  356.     nx = rint ( (data_east - data_west) / header.x_inc) + off;
  357.     ny = rint ( (data_north - data_south) / header.y_inc) + off;
  358.     nm = nx * ny;
  359.     grd = (float *) memory (CNULL, nm, sizeof (float), "grdcontour");
  360.     if (read_grd (grdfile, &header, grd, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  361.         fprintf (stderr, "grdcontour: Error reading file %s\n", grdfile);
  362.         exit (-1);
  363.     }
  364.  
  365.     anchor = old_label = (struct LABEL *) memory (CNULL, 1, sizeof (struct LABEL), "grdcontour");
  366.     n_edges = header.ny * (int )ceil (header.nx / 16.0);
  367.     edge = (int *) memory (CNULL, n_edges, sizeof (int), "grdcontour");
  368.     
  369.     if (mult != 1.0) {
  370.         if (gmtdefs.verbose) fprintf (stderr, "grdcontour: Multiplying grid by %lg\n", mult);
  371.         for (i = 1; i < nm; i++) grd[i] *= mult;
  372.         header.z_min *= mult;
  373.         header.z_max *= mult;
  374.         if (mult < 0.0) d_swap (header.z_min, header.z_max);
  375.     }
  376.     if (c_low > header.z_min) header.z_min = c_low;
  377.     if (c_high < header.z_max) header.z_max = c_high;
  378.     
  379.     small = c_int * 1.0e-6;
  380.     if (a_int == 0.0) a_int = c_int;
  381.     
  382.     small_x = 0.01 * header.x_inc;    small_y = 0.01 * header.y_inc;
  383.     
  384.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  385.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies, gmtdefs.dpi,
  386.         gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  387.     echo_command (argc, argv);
  388.     if (gmtdefs.unix_time) timestamp (argc, argv);
  389.     if (project_info.three_D) ps_transrotate (-z_project.xmin, -z_project.ymin, 0.0);
  390.     
  391.     if (anot) {    /* Wants anotated contours */
  392.         aval = floor (header.z_min / a_int) * a_int;
  393.         if (aval < header.z_min) aval += a_int;
  394.     }
  395.     else
  396.         aval = header.z_max + 1.0;
  397.     
  398.     n_alloc = GMT_CHUNK;
  399.     contour = (double *) memory (CNULL, n_alloc, sizeof (double), "grdcontour");
  400.     cont_type = memory (CNULL, n_alloc, sizeof (char), "grdcontour");
  401.     cont_do_tick = memory (CNULL, n_alloc, sizeof (char), "grdcontour");
  402.     cont_angle = (double *) memory (CNULL, n_alloc, sizeof (double), "grdcontour");
  403.     
  404.     if (cpt_given) {    /* Presumably got a cpt-file */
  405.         read_cpt (cpt_file);
  406.         n_contours = gmt_n_colors + 1;
  407.         if (n_contours > n_alloc) {
  408.             contour = (double *) memory ((char *)contour, n_contours, sizeof (double), "grdcontour");
  409.             cont_type = memory (cont_type, n_contours, sizeof (char), "grdcontour");
  410.             cont_do_tick = memory (cont_do_tick, n_contours, sizeof (char), "grdcontour");
  411.             cont_angle = (double *) memory ((char *)cont_angle, n_contours, sizeof (double), "grdcontour");
  412.         }
  413.         for (i = 0; i < gmt_n_colors; i++) {
  414.             contour[i] = gmt_lut[i].z_low;
  415.             if (gmt_lut[i].anot)
  416.                 cont_type[i] = 'A';
  417.             else
  418.                 cont_type[i] = (anot) ? 'A' : 'C';
  419.             cont_angle[i] = (fix_angle) ? label_angle : gmt_NaN;
  420.             cont_do_tick[i] = do_ticks;
  421.         }
  422.         contour[gmt_n_colors] = gmt_lut[gmt_n_colors-1].z_high;
  423.         if (gmt_lut[gmt_n_colors-1].anot & 2)
  424.             cont_type[gmt_n_colors] = 'A';
  425.         else
  426.             cont_type[gmt_n_colors] = (anot) ? 'A' : 'C';
  427.         cont_angle[gmt_n_colors] = (fix_angle) ? label_angle : gmt_NaN;
  428.         cont_do_tick[gmt_n_colors] = do_ticks;
  429.     }
  430.     else if (fpc != NULL) {    /* read contour info from file */
  431.         n_contours = 0;
  432.         while (fgets (line, 512, fpc)) {
  433.             if (line[0] == '#') continue;
  434.             got = sscanf (line, "%lf %c %lf", &contour[n_contours], &cont_type[n_contours], &tmp);
  435.             if (cont_type[n_contours] == 0) cont_type[n_contours] = 'C';
  436.             cont_do_tick[n_contours] = (do_ticks && ((cont_type[n_contours] == 'C') || (cont_type[n_contours] == 'A'))) ? 1 : 0;
  437.             cont_angle[n_contours] = (got == 3) ? tmp : gmt_NaN;
  438.             n_contours++;
  439.             if (n_contours == n_alloc) {
  440.                 n_alloc += GMT_CHUNK;
  441.                 contour = (double *) memory ((char *)contour, n_alloc, sizeof (double), "grdcontour");
  442.                 cont_type = memory (cont_type, n_alloc, sizeof (char), "grdcontour");
  443.                 cont_do_tick = memory (cont_do_tick, n_alloc, sizeof (char), "grdcontour");
  444.                 cont_angle = (double *) memory ((char *)cont_angle, n_alloc, sizeof (double), "grdcontour");
  445.             }
  446.         }
  447.         fclose (fpc);
  448.     }
  449.     else {    /* Set up contour intervals automatically from c_int and a_int */
  450.         min = floor (header.z_min / c_int) * c_int; if (min < header.z_min) min += c_int;    
  451.         max = ceil (header.z_max / c_int) * c_int; if (max > header.z_max) max -= c_int;
  452.         for (c = rint(min/c_int), n_contours = 0; c <= rint(max/c_int); c++, n_contours++) {
  453.             if (n_contours == n_alloc) {
  454.                 n_alloc += GMT_CHUNK;
  455.                 contour = (double *) memory ((char *)contour, n_alloc, sizeof (double), "grdcontour");
  456.                 cont_type = memory (cont_type, n_alloc, sizeof (char), "grdcontour");
  457.                 cont_angle = (double *) memory ((char *)cont_angle, n_alloc, sizeof (double), "grdcontour");
  458.             }
  459.             contour[n_contours] = c * c_int;
  460.             if (anot && (contour[n_contours] - aval) > SMALL) aval += a_int;
  461.             cont_type[n_contours] = (fabs (contour[n_contours] - aval) < SMALL) ? 'A' : 'C';
  462.             cont_angle[n_contours] = (fix_angle) ? label_angle : gmt_NaN;
  463.             cont_do_tick[n_contours] = do_ticks;
  464.         }
  465.     }
  466.     contour = (double *) memory ((char *)contour, n_contours, sizeof (double), "grdcontour");
  467.     cont_type = memory (cont_type, n_contours, sizeof (char), "grdcontour");
  468.     cont_do_tick = memory (cont_do_tick, n_contours, sizeof (char), "grdcontour");
  469.     cont_angle = (double *) memory ((char *)cont_angle, n_contours, sizeof (double), "grdcontour");
  470.     
  471.     if (do_ticks) save = (struct SAVE *) memory (CNULL, n_alloc, sizeof (struct SAVE), "grdcontour");
  472.     
  473.     for (c = 0; c < n_contours; c++) {    /* For each contour value cval */
  474.         
  475.         /* Reset markers and set up new zero-contour*/
  476.         
  477.         cval = contour[c];
  478.         if (gmtdefs.verbose) fprintf (stderr, "grdcontour: Tracing the %lg contour\n", cval);
  479.         take_out = (first) ? cval : contour[c] - contour[c-1];
  480.         for (i = 0; i < nm; i++) {
  481.             grd[i] -= take_out;
  482.             if (grd[i] == 0.0) grd[i] += small;
  483.         }
  484.         
  485.         section = 0;
  486.         first = FALSE;
  487.         id = (cont_type[c] == 'A' || cont_type[c] == 'a') ? 1 : 0;
  488.         
  489.         ps_setline (pen[id].width);
  490.         ps_setpaint (pen[id].r, pen[id].g, pen[id].b);
  491.         if (pen[id].texture) ps_setdash (pen[id].texture, pen[id].offset);
  492.         
  493.         side = 0;
  494.         begin = TRUE;
  495.         while (side < 5) {
  496.             while ((n = contours (grd, &header, smooth_factor, gmtdefs.interpolant, &side, edge, begin, &x, &y)) > 0) {
  497.         
  498.                 if (fabs (x[0] - x[n-1]) < small_x && fabs (y[0] - y[n-1]) < small_y) {
  499.                     interior = TRUE;
  500.                     x[n-1] = x[0];    y[n-1] = y[0];
  501.                 }
  502.                 else
  503.                     interior = FALSE;
  504.                 if (dump) dump_contour (x, y, n, cval, section, interior, dfile, &multi_segments, EOL_flag);
  505.                 if (n >= ncut) draw_contour (x, y, n, cval, cont_type[c], cont_angle[c], interior, anot_dist, half);
  506.                 if (cont_do_tick[c] && interior && n >= ncut) {
  507.                     save[n_save].x = (double *) memory (CNULL, n, sizeof (double), "grdcontour");
  508.                     save[n_save].y = (double *) memory (CNULL, n, sizeof (double), "grdcontour");
  509.                     memcpy ((char *)save[n_save].x, (char *)x, n * sizeof (double));
  510.                     memcpy ((char *)save[n_save].y, (char *)y, n * sizeof (double));
  511.                     save[n_save].n = n;
  512.                     save[n_save].pen = pen[id].width;
  513.                     save[n_save].do_it = TRUE;
  514.                     save[n_save].cval = cval;
  515.                     n_save++;
  516.                     if (n_save == n_alloc) {
  517.                         save = (struct SAVE *) memory ((char *)save, n_alloc, sizeof (struct SAVE), "grdcontour");
  518.                         n_alloc += GMT_CHUNK;
  519.                     }
  520.                 }
  521.                 if (n >= ncut) section++;
  522.                 begin = FALSE;
  523.                 free ((char *)x);
  524.                 free ((char *)y);
  525.             }
  526.             begin = FALSE;
  527.         }
  528.     }
  529.     
  530.     if (pen[0].texture || pen[1].texture) ps_setdash (CNULL, 0);
  531.  
  532.     if (do_ticks && n_save) {
  533.         save = (struct SAVE *) memory ((char *)save, n_save, sizeof (struct SAVE), "grdcontour");
  534.         read_grd_info (grdfile, &header);
  535.         read_grd (grdfile, &header, grd, data_west, data_east, data_south, data_north, dummy, FALSE);
  536.         if (mult != 1.0) {
  537.             for (i = 1; i < nm; i++) grd[i] *= mult;
  538.             header.z_min *= mult;
  539.             header.z_max *= mult;
  540.             if (mult < 0.0) d_swap (header.z_min, header.z_max);
  541.         }
  542.         sort_and_plot_ticks (save, n_save, &header, grd, tick_gap, tick_length, tick_low, tick_high, tick_label, labels);
  543.         for (i = 0; i < n_save; i++) {
  544.             free ((char *)save[i].x);
  545.             free ((char *)save[i].y);
  546.         }
  547.         free ((char *)save);
  548.     }
  549.     
  550.     plot_labels (label_font_size, box, br, bg, bb);
  551.     
  552.     if (frame_info.plot) {
  553.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  554.         map_basemap ();
  555.     }
  556.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  557.     
  558.     if (project_info.three_D) ps_rotatetrans (z_project.xmin, z_project.ymin, 0.0);
  559.  
  560.     ps_plotend (gmtdefs.last_page);
  561.     
  562.     free ((char *)grd);
  563.     free ((char *)edge);
  564.     free ((char *)contour);
  565.     free (cont_type);
  566.     free ((char *)cont_angle);
  567.     free (cont_do_tick);
  568.  
  569.     
  570.     if (gmtdefs.verbose) fprintf (stderr, "grdcontour: Done!\n");
  571.     
  572.     gmt_end (argc, argv);
  573. }
  574.  
  575. int dump_contour (xx, yy, nn, cval, id, interior, file, m_segment, flag)
  576. double *xx, *yy, cval;
  577. char *file, flag;
  578. int nn, id;
  579. BOOLEAN interior, *m_segment; {
  580.     int i;
  581.     char fname[512], format[80];
  582.     FILE *fp;
  583.     
  584.     if (nn < 2) return;
  585.     
  586.     sprintf (format, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  587.     if (*m_segment) {
  588.         if (*m_segment == 2) {    /* Must create file the first time around */
  589.             fp = fopen (file, "w");
  590.             *m_segment = TRUE;
  591.         }
  592.         else    /* Later we append to it */
  593.             fp = fopen (file, "a+");
  594.         fprintf (fp, "%c %lg contour\n", flag, cval);
  595.     }
  596.     else {
  597.         if (interior)
  598.             sprintf (fname, "%s_%lg_%d_i.xyz\0", file, cval, id);
  599.         else
  600.             sprintf (fname, "%s_%lg_%d.xyz\0", file, cval, id);
  601.         fp = fopen (fname, "w");
  602.     }
  603.     for (i = 0; i < nn; i++) fprintf (fp, format, xx[i], yy[i], cval);
  604.     fclose (fp);
  605. }
  606.  
  607. int draw_contour (xx, yy, nn, cval, ctype, cangle, closed, gap, half_width)
  608. double *xx, *yy, cval, cangle, gap;
  609. char ctype;
  610. int nn, closed, half_width; {
  611.     int i, j;
  612.     double dist, angle, dx, dy, width, sum_x2, sum_xy, one = 1.0;
  613.     char label[100], format[50];
  614.     struct LABEL *new_label;
  615.     
  616.     if (nn < 2) return;
  617.         
  618.     sprintf (format, "%lg contour\0", cval);
  619.     ps_comment (format);
  620.     
  621.     gmt_n_plot = geo_to_xy_line (xx, yy, nn);
  622.     gmt_plot_line (gmt_x_plot, gmt_y_plot, gmt_pen, gmt_n_plot);
  623.     
  624.     if (ctype == 'A' || ctype == 'a') {    /* Annotated contours */
  625.         if (!gmtdefs.measure_unit) one = 2.5;
  626.         get_format (cval, format);
  627.         sprintf (label, format, cval);
  628.         dist = (closed) ? gap - one : 0.0;    /* Label closed contours longer than 1 inch */
  629.         for (i = 1; i < gmt_n_plot; i++) {
  630.             if (gmt_pen[i] == 3) continue;
  631.  
  632.             dx = gmt_x_plot[i] - gmt_x_plot[i-1];
  633.             if (fabs (dx) > (width = gmt_half_map_width (gmt_y_plot[i-1]))) {
  634.                 width *= 2.0;
  635.                 dx = copysign (width - fabs (dx), -dx);
  636.                 if (gmt_x_plot[i] < width)
  637.                     gmt_x_plot[i-1] -= width;
  638.                 else
  639.                     gmt_x_plot[i-1] += width;
  640.             }
  641.             dy = gmt_y_plot[i] - gmt_y_plot[i-1];
  642.             dist += hypot (dx, dy);
  643.             if (dist > gap) {    /* Time for label */
  644.                 new_label = (struct LABEL *) memory (CNULL, 1, sizeof (struct LABEL), "grdcontour");
  645.                 new_label->x = 0.5 * (gmt_x_plot[i-1] + gmt_x_plot[i]);
  646.                 new_label->y = 0.5 * (gmt_y_plot[i-1] + gmt_y_plot[i]);
  647.                 strcpy (new_label->label, label);
  648.                 sum_x2 = sum_xy = 0.0;
  649.                 for (j = i - half_width; j <= i + half_width; j++) {    /* L2 fit for slope */
  650.                     if (j < 0 || j >= gmt_n_plot) continue;
  651.                     dx = gmt_x_plot[j] - new_label->x;
  652.                     dy = gmt_y_plot[j] - new_label->y;
  653.                     sum_x2 += dx * dx;
  654.                     sum_xy += dx * dy;
  655.                 }
  656.                 if (bad_float (cangle)) {    /* Must calculate label angle */
  657.                     angle = d_atan2 (sum_xy, sum_x2) * R2D;
  658.                     if (angle < 0.0) angle += 360.0;
  659.                     if (angle > 90.0 && angle < 270) angle -= 180.0;
  660.                 }
  661.                 else
  662.                     angle = cangle;
  663.                 new_label->angle = angle;
  664.                 new_label->prev_label = old_label;
  665.                 dist = 0.0;
  666.                 old_label->next_label = new_label;
  667.                 old_label = new_label;
  668.             }
  669.             gmt_x_plot[i-1] = gmt_x_plot[i];    gmt_y_plot[i-1] = gmt_y_plot[i];
  670.         }
  671.     }
  672. }
  673.  
  674. int plot_labels (size, box, r, g, b)
  675. int size, box, r, g, b; {
  676.     /* box = 1: white box only, box = 2: white box + draw outline */
  677.     double dx, dy;
  678.     struct LABEL *this, *old;
  679.     
  680.     dx = 0.5 * gmtdefs.anot_offset;
  681.     dy = 0.05 * gmtdefs.anot_offset;
  682.     box--;
  683.     ps_comment ("Contour annotations:");
  684.     
  685.     ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  686.     
  687.     for (old = anchor; old->next_label; old = old->next_label) {    /* First draw boxes if not 3-D*/
  688.         this = old->next_label;
  689.         gmt_textbox3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, box, dx, dy, r, g, b);
  690.     }
  691.     for (old = anchor; old->next_label; old = old->next_label) { /* Then labels */
  692.         this = old->next_label;
  693.         gmt_text3d (this->x, this->y, project_info.z_level, size, gmtdefs.anot_font, this->label, this->angle, 6, 0);
  694.     }
  695.         
  696.     ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  697.     
  698.     this = anchor;
  699.     while (this) {    /* Free memory */
  700.         old = this;
  701.         this = old->next_label;
  702.         free ((char *)old);
  703.     }
  704. }
  705.  
  706. int sort_and_plot_ticks (save, n, h, grd, tick_gap, tick_length, tick_low, tick_high, tick_label, labels)
  707. struct SAVE save[];
  708. int n;
  709. struct GRD_HEADER *h;
  710. float grd[];
  711. double tick_gap, tick_length;
  712. BOOLEAN tick_low, tick_high, tick_label;
  713. char *labels; {
  714.     int np, i, j, inside, ix, iy, done, n_ticks;
  715.     double x0, y0, add, sign, dx, dy, x_back, y_back, x_front, y_front;
  716.     double xmin, xmax, ymin, ymax, x_mean, y_mean, inc, dist, a, small, *s;
  717.     char txt[2][2];
  718.     
  719.     small = 0.1 * project_info.xmax / h->nx;
  720.     if (tick_label) {
  721.         txt[0][0] = labels[0];    txt[0][1] = 0;
  722.         txt[1][0] = labels[1];    txt[1][1] = 0;
  723.     }
  724.     
  725.     for (i = 0; i < n; i++) {
  726.         np = save[i].n;
  727.         for (j = 0; save[i].do_it && j < n; j++) {
  728.             inside = non_zero_winding (save[j].x[0], save[j].y[0], save[i].x, save[i].y, np);
  729.             if (inside == 2) save[i].do_it = FALSE;
  730.         }
  731.     }
  732.     
  733.     for (i = 0; i < n; i++) {
  734.         if (!save[i].do_it) continue;
  735.         np = save[i].n;
  736.         s = (double *) memory (CNULL, np, sizeof (double), "grdcontour");
  737.         
  738.         xmin = xmax = save[i].x[0];
  739.         ymin = ymax = save[i].y[0];
  740.         for (j = 1; j < np; j++) {
  741.             xmin = MIN (xmin, save[i].x[j]);
  742.             xmax = MAX (xmax, save[i].x[j]);
  743.             ymin = MIN (ymin, save[i].y[j]);
  744.             ymax = MAX (ymax, save[i].y[j]);
  745.         }
  746.         
  747.         ix = rint ((0.5 * (xmin + xmax) - h->x_min) / h->x_inc) + 1;
  748.         iy = rint ((h->y_max - ymax) / h->y_inc) + 1;
  749.         x0 = h->x_min + ix * h->x_inc;
  750.         done = FALSE;
  751.         while (!done && iy < h->ny) {
  752.             y0 = h->y_max - iy * h->y_inc;
  753.             inside = non_zero_winding (x0, y0, save[i].x, save[i].y, np);
  754.             if (inside == 2)
  755.                 done = TRUE;
  756.             else
  757.                 iy++;
  758.         }
  759.         save[i].high = (grd[iy*h->nx+ix] > save[i].cval);
  760.         
  761.         if (save[i].high && !tick_high) continue;    /* Dont tick highs */
  762.         if (!save[i].high && !tick_low) continue;    /* Dont tick lows */
  763.         
  764.         geo_to_xy (save[i].x[0], save[i].y[0], &save[i].x[0], &save[i].y[0]);
  765.         if (tick_label) {
  766.             x_mean = save[i].x[0];
  767.             y_mean = save[i].y[0];
  768.         }
  769.         for (j = 1, s[0] = 0.0; j < np; j++) {
  770.             geo_to_xy (save[i].x[j], save[i].y[j], &save[i].x[j], &save[i].y[j]);
  771.             s[j] = s[j-1] + hypot (save[i].x[j]-save[i].x[j-1], save[i].y[j]-save[i].y[j-1]);
  772.             if (tick_label) {
  773.                 x_mean += save[i].x[j];
  774.                 y_mean += save[i].y[j];
  775.             }
  776.         }
  777.         if (s[np-1] < MIN_LENGTH) {
  778.             free ((char *)s);
  779.             continue;
  780.         }
  781.         dx = save[i].x[1] - save[i].x[0];
  782.         dy = save[i].y[1] - save[i].y[0];
  783.         a = atan2 (dy, dx) + M_PI_2;
  784.         x0 = 0.5 * (save[i].x[0] + save[i].x[1]) + small * cos (a);
  785.         y0 = 0.5 * (save[i].y[0] + save[i].y[1]) + small * sin (a);
  786.         inside = non_zero_winding (x0, y0, save[i].x, save[i].y, np);
  787.         sign = (inside == 2) ? 1.0 : -1.0;
  788.         
  789.         n_ticks = s[np-1] / tick_gap;
  790.         if (n_ticks == 0) {
  791.             free ((char *)s);
  792.             continue;
  793.         }
  794.         inc = s[np-1] / n_ticks;
  795.         if (tick_label) {
  796.             x_mean /= np;
  797.             y_mean /= np;
  798.             gmt_text3d (x_mean, y_mean, project_info.z_level, gmtdefs.anot_font_size, gmtdefs.anot_font, txt[save[i].high], 0.0, 6, 0);
  799.         }
  800.         
  801.         x_back = save[i].x[np-1];
  802.         y_back = save[i].y[np-1];
  803.         dist = 0.0;
  804.         j = 0;
  805.         add = sign * ((save[i].high) ? -1.0 : 1.0);
  806.         ps_setline (save[i].pen);
  807.         while (j < np-1) {
  808.             x_front = save[i].x[j+1];
  809.             y_front = save[i].y[j+1];
  810.             if (s[j] >= dist) {    /* Time for tick */
  811.                 dx = x_front - x_back;
  812.                 dy = y_front - y_back;
  813.                 a = atan2 (dy, dx) + add * M_PI_2;
  814.                 if (project_info.three_D) {
  815.                     xy_do_z_to_xy (save[i].x[j], save[i].y[j], project_info.z_level, &x0, &y0);
  816.                     ps_plot (x0, y0, 3);
  817.                     xy_do_z_to_xy (save[i].x[j] + tick_length * cos (a), save[i].y[j] + tick_length * sin (a), project_info.z_level, &x0, &y0);
  818.                     ps_plot (x0, y0, 2);
  819.                 }
  820.                 else {
  821.                     ps_plot (save[i].x[j], save[i].y[j], 3);
  822.                     ps_plotr (tick_length * cos (a), tick_length * sin (a), 2);
  823.                 }
  824.                 dist += inc;
  825.             }
  826.             x_back = x_front;
  827.             y_back = y_front;
  828.             j++;
  829.         }
  830.         free ((char *)s);
  831.     }
  832. }
  833.