home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdvector.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-23  |  12.2 KB  |  386 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdvector.c    3.15  23 Jul 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.    grdvector reads 2 grdfiles that contains the 2 components of a vector
  9.    field (cartesian or polar) and plots vectors at the grid positions.
  10.    This is basically a short-hand for using grd2xyz | psxy -SV and is
  11.    more convenient for such plots on a grid.
  12.  
  13.  
  14.  */
  15.  
  16.  
  17. #include "gmt.h"
  18.  
  19. float *r, *theta;
  20.  
  21. main (argc, argv)
  22. int argc;
  23. char **argv; {
  24.  
  25.     int    i, j, n = 0, nm, nx, ny, ij, i0, j0, di, dj, dummy[4];
  26.     
  27.     BOOLEAN convert_angles = FALSE, get_rgb = FALSE, cartesian = TRUE, shrink = FALSE, set_fill = FALSE;
  28.     BOOLEAN error = FALSE, center = FALSE, outline = FALSE, azimuth = FALSE, stick_plot = TRUE, inc_set = FALSE;
  29.     BOOLEAN clip = TRUE;
  30.     
  31.     char *file[2], *cpt;
  32.     
  33.     double dx2, dy2, v_width = 0.03, h_length = 0.12, h_width = 0.1;
  34.     double v_w, h_l, h_w, v_shrink, v_norm = 0.0, tmp, x, y, plot_x, plot_y, x_off, y_off;
  35.     double west, east, south, north, off, dir, length, x2, y2, scale = 1.0;
  36.     double data_west, data_east, data_south, data_north, value;
  37.     
  38.     struct GRD_HEADER h[2];
  39.     struct FILL fill;
  40.     struct PEN pen;
  41.  
  42.     gmt_init_pen (&pen, 1);
  43.     gmt_init_fill (&fill, -1, -1, -1);
  44.     west = east = south = north = 0.0;
  45.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  46.     di = dj = 1;
  47.     i0 = j0 = 0;
  48.     dx2 = dy2 = 0.0;
  49.     
  50.     argc = gmt_begin (argc, argv);
  51.  
  52.     for (i = 1; i < argc; i++) {
  53.         if (argv[i][0] == '-') {
  54.             switch (argv[i][1]) {
  55.                 /* Common parameters */
  56.             
  57.                 case 'B':
  58.                 case 'J':
  59.                 case 'K':
  60.                 case 'O':
  61.                 case 'P':
  62.                 case 'R':
  63.                 case 'U':
  64.                 case 'V':
  65.                 case 'X':
  66.                 case 'x':
  67.                 case 'Y':
  68.                 case 'y':
  69.                 case 'c':
  70.                 case '\0':
  71.                     error += get_common_args (argv[i], &west, &east, &south, &north);
  72.                     break;
  73.                 
  74.                 /* Supplemental parameters */
  75.             
  76.                 case 'A':
  77.                     cartesian = FALSE;
  78.                     break;
  79.                 case 'C':    /* Vary symbol color with z */
  80.                     cpt = &argv[i][2];
  81.                     get_rgb = TRUE;
  82.                     break;
  83.                 case 'E':
  84.                     center = TRUE;
  85.                     break;
  86.                 case 'F':
  87.                     if (gmt_getrgb (&argv[i][2], &gmtdefs.basemap_frame_rgb[0], &gmtdefs.basemap_frame_rgb[1], &gmtdefs.basemap_frame_rgb[2])) {
  88.                         gmt_pen_syntax ('F');
  89.                         error++;
  90.                     }
  91.                     break;
  92.                 case 'G':        /* Set Gray shade for polygon */
  93.                     if (gmt_getfill (&argv[i][2], &fill)) {
  94.                         gmt_fill_syntax ('G');
  95.                         error++;
  96.                     }
  97.                     set_fill = TRUE;
  98.                     break;
  99.                 case 'I':    /* Only use gridnodes dx2,dy2 apart */
  100.                     gmt_getinc (&argv[i][2], &dx2, &dy2);
  101.                     inc_set = TRUE;
  102.                     break;
  103.                 case 'N':    /* Do not clip at border */
  104.                     clip = FALSE;
  105.                     break;
  106.                 case 'Q':
  107.                     if (argv[i][2] && argv[i][3] != 'n') {
  108.                         if (sscanf (&argv[i][2], "%lf/%lf/%lf", &v_width, &h_length, &h_width) != 3) {
  109.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -Q option:  Could not decode arrowwidth/headlength/headwidth\n", gmt_program);
  110.                             error++;
  111.                         }
  112.                     }
  113.                     for (j = 2; argv[i][j] && argv[i][j] != 'n'; j++);
  114.                     if (argv[i][j]) {    /* Normalize option used */
  115.                         v_norm = atof (&argv[i][j+1]);
  116.                         if (v_norm > 0.0) {
  117.                             v_shrink = 1.0 / v_norm;
  118.                             shrink = TRUE;
  119.                         }
  120.                         else {
  121.                             fprintf (stderr, "%s: GMT SYNTAX ERROR -Qn option:  No reference length given\n", gmt_program);
  122.                             error++;
  123.                         }
  124.                     }
  125.                     stick_plot = FALSE;
  126.                     break;
  127.                 case 'S':
  128.                     scale = atof (&argv[i][2]);
  129.                     break;
  130.                 case 'T':
  131.                     convert_angles = TRUE;
  132.                     break;
  133.                 case 'W':        /* Set line attributes */
  134.                     if (argv[i][2] && gmt_getpen (&argv[i][2], &pen)) {
  135.                         gmt_pen_syntax ('W');
  136.                         error++;
  137.                     }
  138.                     outline = TRUE;
  139.                     break;
  140.                 case 'Z':
  141.                     azimuth = TRUE;
  142.                     break;
  143.                 default:
  144.                     error = TRUE;
  145.                     gmt_default_error (argv[i][1]);
  146.                     break;
  147.             }
  148.         }
  149.         else if (n < 2)
  150.             file[n++] = argv[i];
  151.         else
  152.             n++;
  153.     }
  154.     
  155.     if (argc == 1 || gmt_quick) {
  156.         fprintf (stderr, "grdvector %s - Plot vector fields from grdfiles\n\n", GMT_VERSION);
  157.         fprintf (stderr, "usage: grdvector compx.grd compy.grd -J<params> -R<west/east/south/north> [-A]\n");
  158.         fprintf (stderr, "\t[-B<tickinfo>] [-C<cpt>] [-E] [-F<r/g/b>] [-G<fill>] [-K] [-O] [-P] [-Q<params>] [-N] [-S<scale>] [-T]\n");
  159.         fprintf (stderr, "\t[-U[<label>]] [-V] [-W<pen>] [-X<x_shift>] [-Y<y_shift>] [-Z] [-c<ncopies>] [-:]\n\n");
  160.         
  161.         if (gmt_quick) exit (-1);
  162.         
  163.         fprintf (stderr, "\tcompx & compy are grdfiles with the 2 vector components.\n");
  164.         explain_option ('j');
  165.         fprintf (stderr, "\n\tOPTIONS:\n");
  166.         fprintf (stderr, "\t-A means grdfiles have polar (r, theta) components [Default is Cartesian]\n");
  167.         explain_option ('b');
  168.         fprintf (stderr, "    -C Use cpt-file to assign colors based on vector length\n");
  169.         fprintf (stderr, "    -E cEnter vectors on grid nodes [Default draws from grid node]\n");
  170.         fprintf (stderr, "      -F Set color used for Frame and anotation [%d/%d/%d]\n",
  171.             gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  172.         fprintf (stderr, "    -G Specify color. fill can be either <r/g/b> (each 0-255) for color or (0-255) for gray [0].\n");
  173.         fprintf (stderr, "       Default is no fill (vector outlines only)\n");
  174.         explain_option ('K');
  175.         fprintf (stderr, "    -N Do Not clip vectors that exceed the map boundaries [Default will clip]\n");
  176.         explain_option ('O');
  177.         explain_option ('P');
  178.         fprintf (stderr, "    -Q Select vector plot [Default is stick-plot].\n");
  179.         fprintf (stderr, "       Optionally, specify vector parameters\n");
  180.         fprintf (stderr, "       <params> are arrowwidth/headlength/headwidth [Default is 0.03/0.12/0.09]\n");
  181.         fprintf (stderr, "       Append n<size> which will cause vectors shorter than <size> to be\n");
  182.         fprintf (stderr, "         scaled down\n");
  183.         explain_option ('R');
  184.         fprintf (stderr, "    -S sets scale for vector length in units per inch [1]\n");
  185.         fprintf (stderr, "    -T means aximuth should be convered to angles based on map projection\n");
  186.         explain_option ('U');
  187.         explain_option ('V');
  188.         fprintf (stderr, "    -W sets pen attributes [width = %d, color = (%d/%d/%d), texture = solid line].\n", 
  189.             pen.width, pen.r, pen.g, pen.b);
  190.         explain_option ('X');
  191.         fprintf (stderr, "    -Z means the angles provided are azimuths rahter than direction\n");
  192.         explain_option ('c');
  193.         explain_option ('.');
  194.         exit (-1);
  195.     }
  196.  
  197.     if (inc_set && (dx2 <= 0.0 || dy2 <= 0.0)) {
  198.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option:  Must specify positive increments\n", gmt_program);
  199.         error++;
  200.     }
  201.     if (scale == 0.0) {
  202.         fprintf (stderr, "%s: GMT SYNTAX ERROR -S option:  Scale must be nonzero\n", gmt_program);
  203.         error++;
  204.     }
  205.     if (azimuth && cartesian) {
  206.         fprintf (stderr, "%s: GMT SYNTAX ERROR -Z option:  Azimuths not valid input for Cartesian data\n", gmt_program);
  207.         error++;
  208.     }
  209.     if (get_rgb && !cpt) {
  210.         fprintf (stderr, "%s: GMT SYNTAX ERROR -C option:  Must specify a color palette table\n", gmt_program);
  211.         error++;
  212.     }
  213.     if (!(set_fill || outline)) {
  214.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify at least one of -G, -W\n", gmt_program);
  215.         error++;
  216.     }
  217.     if (n != 2) {
  218.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify two input grdfiles\n", gmt_program);
  219.         error++;
  220.     }
  221.     
  222.     if (error) exit (-1);
  223.  
  224.     if (get_rgb) read_cpt (cpt);
  225.  
  226.     if (!(strcmp (file[0], "=") || strcmp (file[1], "="))) {
  227.         fprintf (stderr, "grdvector: Piping of grdfiles not supported!\n");
  228.         exit (-1);
  229.     }
  230.  
  231.     if (read_grd_info (file[0], &h[0])) {
  232.         fprintf (stderr, "grdvector: Error opening file %s\n", file[0]);
  233.         exit (-1);
  234.     }
  235.     
  236.     if (read_grd_info (file[1], &h[1])) {
  237.         fprintf (stderr, "grdvector: Error opening file %s\n", file[1]);
  238.         exit (-1);
  239.     }
  240.     
  241.     if (!(h[0].nx == h[1].nx && h[0].ny == h[1].ny && h[0].x_min == h[1].x_min && h[0].y_min == h[1].y_min 
  242.         && h[0].x_inc == h[1].x_inc && h[0].y_inc == h[1].y_inc)) {
  243.         fprintf (stderr, "grdvector: files %s and %s does not match!\n", file[0], file[1]);
  244.         exit (-1);
  245.     }
  246.     off = (h[0].node_offset) ? 0 : 1;
  247.         
  248.     /* Determine what wesn to pass to map_setup */
  249.  
  250.     if (west == east && south == north) {
  251.         west = h[0].x_min;
  252.         east = h[0].x_max;
  253.         south = h[0].y_min;
  254.         north = h[0].y_max;
  255.     }
  256.  
  257.     map_setup (west, east, south, north);
  258.  
  259.     /* Determine the wesn to be used to read the grdfile */
  260.  
  261.     grd_setregion (&h[0], &data_west, &data_east, &data_south, &data_north);
  262.  
  263.     /* Read data */
  264.  
  265.     nx = rint ( (data_east - data_west) / h[0].x_inc) + off;
  266.     ny = rint ( (data_north - data_south) / h[0].y_inc) + off;
  267.     nm = nx * ny;
  268.     r = (float *) memory (CNULL, nm, sizeof (float), "grdvector");
  269.     theta = (float *) memory (CNULL, nm, sizeof (float), "grdvector");
  270.     if (read_grd (file[0], &h[0], r, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  271.         fprintf (stderr, "grdvector: Error reading file %s\n", file[0]);
  272.         exit (-1);
  273.     }
  274.     if (read_grd (file[1], &h[1], theta, data_west, data_east, data_south, data_north, dummy, FALSE)) {
  275.         fprintf (stderr, "grdvector: Error reading file %s\n", file[1]);
  276.         exit (-1);
  277.     }
  278.     scale = 1.0 / scale;
  279.     
  280.     ps_plotinit (CNULL, gmtdefs.overlay, gmtdefs.page_orientation, gmtdefs.x_origin, gmtdefs.y_origin,
  281.         gmtdefs.global_x_scale, gmtdefs.global_y_scale, gmtdefs.n_copies,
  282.         gmtdefs.dpi, gmtdefs.measure_unit, gmtdefs.paper_width, gmtdefs.page_rgb, gmt_epsinfo (argv[0]));
  283.     echo_command (argc, argv);
  284.     if (gmtdefs.unix_time) timestamp (argc, argv);
  285.  
  286.     ps_setline (pen.width);
  287.     if (pen.texture) ps_setdash (pen.texture, pen.offset);
  288.     ps_setpaint (pen.r, pen.g, pen.b);
  289.  
  290.         if (clip) map_clip_on (-1, -1, -1, 3);
  291.  
  292.     if (dx2 != 0.0 && dy2 != 0.0) {
  293.         dj = rint (dy2 / h[0].y_inc);
  294.         di = rint (dx2 / h[0].x_inc);
  295.         tmp = ceil (h[0].y_max / dy2) * dy2;
  296.         if (tmp > h[0].y_max) tmp -= dy2;
  297.         j0 = rint ((h[0].y_max - tmp) / h[0].y_inc);
  298.         tmp = floor (h[0].x_min / dx2) * dx2;
  299.         if (tmp < h[0].x_min) tmp += dx2;
  300.         i0 = rint ((tmp - h[0].x_min) / h[0].x_inc);
  301.     }
  302.     
  303.     for (j = j0; j < h[1].ny; j += dj) {
  304.         y = h[0].y_max - j *h[0].y_inc;
  305.         for (i = i0; i < h[1].nx; i += di) {
  306.         
  307.             ij = j * h[0].nx + i;
  308.             value = r[ij];
  309.             
  310.             if (cartesian) {
  311.                 length = hypot (theta[ij], r[ij]);
  312.                 if (length == 0.0) continue;
  313.                 dir = R2D * atan2 (theta[ij], r[ij]);
  314.                 if (r[ij] < 0) dir += M_PI;
  315.                 value = copysign (length, r[ij]);
  316.                 r[ij] = length;
  317.                 theta[ij] = dir;
  318.             }
  319.             else if (r[ij] < 0.0) {
  320.                 r[ij] = -r[ij];
  321.                 theta[ij] += M_PI;
  322.             }
  323.             else if (r[ij] == 0.0) continue;
  324.             
  325.             if (get_rgb) get_rgb24 (value, &fill.r, &fill.g, &fill.b);
  326.             
  327.             x = h[0].x_min + i * h[0].x_inc;
  328.             geo_to_xy (x, y, &plot_x, &plot_y);
  329.             
  330.             if (convert_angles) {
  331.                 if (!azimuth) theta[ij] = 90.0 - theta[ij];
  332.                 azim_2_angle (x, y, 0.1, (double)theta[ij], &tmp);
  333.                 theta[ij] = tmp * D2R;
  334.             }
  335.             else
  336.                 theta[ij] *= D2R;
  337.             
  338.             r[ij] *= scale;
  339.             
  340.             x2 = plot_x + r[ij] * cos (theta[ij]);
  341.             y2 = plot_y + r[ij] * sin (theta[ij]);
  342.             
  343.             if (center) {
  344.                 x_off = 0.5 * (x2 - plot_x);
  345.                 y_off = 0.5 * (y2 - plot_y);
  346.                 plot_x -= x_off;
  347.                 plot_y -= y_off;
  348.                 x2 -= x_off;
  349.                 y2 -= y_off;
  350.             }
  351.             
  352.             if (stick_plot) {
  353.                 if (get_rgb) ps_setpaint (fill.r, fill.g, fill.b);
  354.                 ps_plot (plot_x, plot_y, 3);
  355.                 ps_plot (x2, y2, 2);
  356.                 continue;
  357.             }
  358.             
  359.             if (shrink && r[ij] < v_norm) {    /* Scale arrow attributes down with length */
  360.                 v_w = v_width * r[ij] * v_shrink;
  361.                 h_l = h_length * r[ij] * v_shrink;
  362.                 h_w = h_width * r[ij] * v_shrink;
  363.                 ps_vector (plot_x, plot_y, x2, y2, v_w, h_l, h_w, gmtdefs.vector_shape, fill.r, fill.g, fill.b, outline);
  364.             }
  365.             else    /* Leave as specified */
  366.                 ps_vector (plot_x, plot_y, x2, y2, v_width, h_length, h_width, gmtdefs.vector_shape, fill.r, fill.g, fill.b, outline);
  367.         }
  368.     }
  369.     
  370.     if (pen.texture) ps_setdash (CNULL, 0);
  371.  
  372.         if (clip) map_clip_off ();
  373.  
  374.     if (frame_info.plot) {
  375.         ps_setpaint (gmtdefs.basemap_frame_rgb[0], gmtdefs.basemap_frame_rgb[1], gmtdefs.basemap_frame_rgb[2]);
  376.         map_basemap ();
  377.         ps_setpaint (gmtdefs.background_rgb[0], gmtdefs.background_rgb[1], gmtdefs.background_rgb[2]);
  378.     }
  379.     ps_plotend (gmtdefs.last_page);
  380.  
  381.     free ((char *)r);
  382.     free ((char *)theta);
  383.     
  384.     gmt_end (argc, argv);
  385. }
  386.