home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grdtrack.c < prev    next >
C/C++ Source or Header  |  1995-03-13  |  6KB  |  201 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grdtrack.c    2.19  3/13/95
  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.  * grdtrack reads a xyfile, opens the 2d binary gridded grdfile, 
  9.  * and samples the dataset at the xy positions with a bilinear or bicubic
  10.  * interpolant.  This new data is added to the input as an extra column
  11.  * and printed to standard output.  In order to evaluate derivatives along
  12.  * the edges of the grdfile region, we assume natural bicubic spline
  13.  * boundary conditions (d2z/dn2 = 0, n being the normal to the edge;
  14.  * d2z/dxdy = 0 in the corners).  Rectangles of size x_inc by y_inc are 
  15.  * mapped to [0,1] x [0,1] by affine transformation, and the interpolation
  16.  * done on the normalized rectangle.
  17.  *
  18.  * Author:    Walter H F Smith
  19.  * Date:    23-SEP-1993
  20.  * 
  21.  * Based on the original grdtrack, which had this authorship/date/history:
  22.  *
  23.  * Author:    Paul Wessel
  24.  * Date:    29-JUN-1988
  25.  * Revised:    5-JAN-1990    PW: Updated to v.2.0
  26.  *        4-AUG-1993    PW: Added -Q
  27.  */
  28.  
  29. #include "gmt.h"
  30. #include "gmt_bcr.h"
  31.  
  32. float *f;
  33.  
  34. main (argc, argv)
  35. int argc;
  36. char **argv; {
  37.     int i, ix, iy, mx, my, n_read, n_points = 0, pad[4];
  38.     
  39.     BOOLEAN error = FALSE, multi_segments = FALSE, bilinear = FALSE, periodic = FALSE;
  40.     
  41.     double value, xy[2];
  42.     
  43.     char *grdfile, stuff[80], line[512], format1[512], format2[512], EOL_flag = '>';
  44.     
  45.     FILE *fp = NULL;
  46.     
  47.     struct GRD_HEADER grd;
  48.     
  49.     grdfile = CNULL;
  50.     
  51.     argc = gmt_begin (argc, argv);
  52.  
  53.     for (i = 1; i < argc; i++) {
  54.         if (argv[i][0] == '-') {
  55.             switch (argv[i][1]) {
  56.             
  57.                 /* Common parameters */
  58.             
  59.                 case 'H':
  60.                 case 'V':
  61.                 case ':':
  62.                 case '\0':
  63.                     error += get_common_args (argv[i], 0, 0, 0, 0);
  64.                     break;
  65.  
  66.                 /* Supplemental parameters */
  67.                 
  68.                 case 'G':
  69.                     grdfile = &argv[i][2];
  70.                     break;
  71.                 case 'L':
  72.                     periodic = TRUE;
  73.                     break;
  74.                 case 'M':
  75.                     multi_segments = TRUE;
  76.                     if (argv[i][2]) EOL_flag = argv[i][2];
  77.                     break;
  78.                 case 'Q':
  79.                     bilinear = TRUE;
  80.                     break;
  81.  
  82.                 default:
  83.                     error = TRUE;
  84.                     gmt_default_error (argv[i][1]);
  85.                     break;
  86.             }
  87.         }
  88.         else if ((fp = fopen (argv[i], "r")) == NULL) {
  89.             fprintf (stderr, "grdtrack: Cannot open file %s\n", argv[i]);
  90.             exit (-1);
  91.         }
  92.     }
  93.     
  94.     if (argc == 1 || gmt_quick) {
  95.         fprintf (stderr,"grdtrack %s - Sampling of a 2-D gridded netCDF grdfile along 1-D trackline\n\n", GMT_VERSION);
  96.         fprintf (stderr, "usage: grdtrack <xyfile> -G<grdfile> [-H] [-L] [-M[<flag>]] [-Q] [-V] [-:]\n");
  97.         
  98.         if (gmt_quick) exit(-1);
  99.         
  100.         fprintf (stderr, "    <xyfile> is an multicolumn ASCII file with (lon,lat) in the first two columns\n");
  101.         fprintf (stderr, "    -G <grdfile> is the name of the 2-D binary data set\n");
  102.         fprintf (stderr, "\n\tOPTIONS:\n");
  103.         explain_option ('H');
  104.         fprintf (stderr, "    -L means that x is longitude\n");
  105.         fprintf (stderr, "    -M Input file consist of multiple segments separated by one record\n");
  106.         fprintf (stderr, "       whose first character is <flag> [>].\n");
  107.         fprintf (stderr, "    -Q Quick mode, use bilinear rather than bicubic interpolation\n");
  108.         explain_option ('V');
  109.         explain_option (':');
  110.         explain_option ('.');
  111.         exit(-1);
  112.     }
  113.     
  114.     if (!grdfile) {
  115.         fprintf (stderr, "%s: GMT SYNTAX ERROR -G:  Must specify output file\n", gmt_program);
  116.         error++;
  117.     }
  118.     
  119.     if (error) exit (-1);
  120.  
  121.     if (fp == NULL) {
  122.         fp = stdin;
  123.         fprintf (stderr, "grdtrack: Reads from standard input\n");
  124.     }
  125.  
  126.     if (read_grd_info (grdfile, &grd)) {
  127.         fprintf (stderr, "grdtrack: Error opening file %s\n", grdfile);
  128.         exit (-1);
  129.     }
  130.     
  131.     mx = grd.nx + 2;
  132.     my = grd.ny + 2;
  133.     
  134.     f = (float *) memory (CNULL, mx * my, sizeof (float), "grdtrack");
  135.  
  136.     pad[0] = pad[1] = pad[2] = pad[3] = 1;
  137.     if (read_grd (grdfile, &grd, f, grd.x_min, grd.x_max, grd.y_min, grd.y_max, pad, FALSE)) {
  138.         fprintf (stderr, "grdtrack: Error reading file %s\n", grdfile);
  139.         exit (-1);
  140.     }
  141.  
  142.     /* Initialize bcr structure:  */
  143.  
  144.     bcr_init (&grd, pad, bilinear);
  145.  
  146.     /* Set boundary conditions  */
  147.     
  148.     set_bcr_boundaries (&grd, f);
  149.     
  150.     sprintf (format1, "%s\t%s\t%%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  151.     sprintf (format2, "%s\t%s\t%s\n\0", gmtdefs.d_format, gmtdefs.d_format, gmtdefs.d_format);
  152.     
  153.     if (gmtdefs.io_header) {    /* First echo headers, if any */
  154.         for (i = 0; i < gmtdefs.n_header_recs - 1; i++) {
  155.             fgets (line, 512, fp);
  156.             printf ("%s", line);
  157.         }
  158.         fgets (line, 512, fp);
  159.         line[strlen(line)-1] = 0;
  160.         printf ("%s\tsample\n", line);
  161.     }
  162.  
  163.     ix = (gmtdefs.xy_toggle) ? 1 : 0;    iy = 1 - ix;        /* Set up which columns have x and y */
  164.     while (fgets (line, 512, fp)) {
  165.  
  166.         if (multi_segments && line[0] == EOL_flag) {
  167.             printf ("%s", line);
  168.             continue;
  169.         }
  170.  
  171.         n_read = sscanf (line, "%lf %lf %[^\n]", &xy[ix], &xy[iy], stuff);
  172.  
  173.         /* Check that we are inside the grd area */
  174.         
  175.         if (xy[1] < grd.y_min || xy[1] > grd.y_max) continue;
  176.         if (periodic) {
  177.             xy[0] -= 360.0;
  178.             while (xy[0] < grd.x_min) xy[0] += 360.0;
  179.             if (xy[0] > grd.x_max) continue;
  180.         }
  181.         else if (xy[0] < grd.x_min || xy[0] > grd.x_max)
  182.             continue;
  183.         
  184.         value = get_bcr_z(&grd, xy[0], xy[1], f);
  185.  
  186.         if (n_read == 2)
  187.             printf (format2, xy[0], xy[1], value);
  188.         else
  189.             printf (format1, xy[0], xy[1], stuff, value);
  190.         n_points++;
  191.     }
  192.     fclose(fp);
  193.     
  194.     if (gmtdefs.verbose) fprintf (stderr, "grdtrack: Sampled %d points from grid %s (%d x %d)\n",
  195.         n_points, grdfile, grd.nx, grd.ny);
  196.     
  197.     free ((char *)f);
  198.     
  199.     gmt_end (argc, argv);
  200. }
  201.