home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / grd2xyz.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-03-13  |  4.9 KB  |  176 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)grd2xyz.c    2.16  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.  * grd2xyz.c reads a grd file and prints out the x,y,z values to
  9.  * standard output.
  10.  *
  11.  * Author:    Paul Wessel
  12.  * Date:    3-JAN-1991-1995
  13.  * Version:    2.0    Based on old v1.x
  14.  */
  15.  
  16. #include "gmt.h"
  17.  
  18. float *z;
  19.  
  20. main (argc, argv)
  21. int argc;
  22. char **argv; {
  23.     int i, j, ij, nm, nx, ny, dummy[4], one_or_zero, ix, iy;
  24.     int error = FALSE, global = FALSE, z_only = FALSE, binary = FALSE, single_precision = FALSE;
  25.     double w, e, s, n, d, *x, *y, out[3];
  26.     char *grdfile = NULL;
  27.     struct GRD_HEADER grd;
  28.     
  29.     argc = gmt_begin (argc, argv);
  30.  
  31.     w = e = s = n = 0.0;
  32.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  33.     
  34.     for (i = 1; i < argc; i++) {
  35.         if (argv[i][0] == '-') {
  36.             switch (argv[i][1]) {
  37.                 /* Common parameters */
  38.             
  39.                 case 'H':
  40.                 case 'R':
  41.                 case 'V':
  42.                 case ':':
  43.                 case '\0':
  44.                     error += get_common_args (argv[i], &w, &e, &s, &n);
  45.                     break;
  46.                 case 'b':
  47.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  48.                     binary = TRUE;
  49.                     break;
  50.                 case 'Z':
  51.                     z_only = TRUE;
  52.                     break;
  53.                     
  54.                 default:
  55.                     error = TRUE;
  56.                     gmt_default_error (argv[i][1]);
  57.                     break;
  58.             }
  59.         }
  60.         else
  61.             grdfile = argv[i];
  62.     }
  63.     
  64.     if (argc == 1 || gmt_quick) {
  65.         fprintf (stderr, "grd2xyz %s - Converting a netCDF grdfile to an ASCII xyz-file\n\n", GMT_VERSION);
  66.         fprintf( stderr, "usage: grd2xyz grdfile [-H] [-Rw/e/s/n] [-V] [-Z] [-:] [-b[d]]> xyzfile\n");
  67.  
  68.         if (gmt_quick) exit (-1);
  69.         
  70.         fprintf (stderr, "    grdfile is the grd file to convert\n");
  71.         fprintf (stderr, "\n\tOPTIONS:\n");
  72.         explain_option ('H');
  73.         explain_option ('R');
  74.         explain_option ('V');
  75.         fprintf (stderr, "    -Z Only output the z-array in scanline orientation\n");
  76.         explain_option (':');
  77.         fprintf (stderr, "    -b means binary output.  Append d for double precision [Default is single].\n");
  78.         exit (-1);
  79.     }
  80.     
  81.     if (!grdfile) {
  82.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify input file\n", gmt_program);
  83.         error++;
  84.     }
  85.     if (binary && gmtdefs.io_header) {
  86.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  87.         error++;
  88.     }
  89.     if (error) exit (-1);
  90.  
  91.     if (read_grd_info (grdfile, &grd)) {
  92.         fprintf (stderr, "grd2xyz: Error opening file %s\n", grdfile);
  93.         exit (-1);
  94.     }
  95.     
  96.     if (binary) gmt_output = (single_precision) ? bin_float_output : bin_double_output;
  97.  
  98.     ix = (gmtdefs.xy_toggle) ? 1 : 0;        iy = 1 - ix;              /* Set up which columns have x and y */
  99.  
  100.     nm = grd.nx * grd.ny;
  101.     
  102.     if (e > w && n > s) {
  103.         global = (fabs (grd.x_max - grd.x_min) == 360.0);
  104.         if (!global && (w < grd.x_min || e > grd.x_max)) error = TRUE;
  105.         if (s < grd.y_min || n > grd.y_max) error = TRUE;
  106.         if (error) {
  107.             fprintf (stderr, "%s: GMT ERROR: Subset exceeds data domain!\n", gmt_program);
  108.             exit (-1);
  109.         }
  110.         one_or_zero = (grd.node_offset) ? 0 : 1;
  111.         nx = rint ((e - w) / grd.x_inc) + one_or_zero;
  112.         ny = rint ((n - s) / grd.y_inc) + one_or_zero;
  113.         
  114.         z = (float *) memory (CNULL, nx * ny, sizeof (float), "grd2xyz");
  115.         
  116.         if (read_grd (grdfile, &grd, z, w, e, s, n, dummy, FALSE)) {
  117.             fprintf (stderr, "grd2xyz: Error reading file %s\n", grdfile);
  118.             exit (-1);
  119.         }
  120.     }
  121.     else {
  122.         z = (float *) memory (CNULL, nm, sizeof (float), "grd2xyz");
  123.  
  124.         if (read_grd (grdfile, &grd, z, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  125.             fprintf (stderr, "grd2xyz: Error reading file %s\n", grdfile);
  126.             exit (-1);
  127.         }
  128.     }
  129.  
  130.     if (!z_only) {
  131.     
  132.         x = (double *) memory (CNULL, grd.nx, sizeof (double), "grd2xyz");
  133.         y = (double *) memory (CNULL, grd.ny, sizeof (double), "grd2xyz");
  134.  
  135.         /* Compute grid node positions once only */
  136.  
  137.         d = (grd.node_offset) ? 0.5 : 0.0;
  138.  
  139.         for (j = 0; j < grd.ny; j++) y[j] = (j == (grd.ny-1)) ? grd.y_min + d * grd.y_inc : grd.y_max - (j + d) * grd.y_inc;
  140.  
  141.         for (i = 0; i < grd.nx; i++) x[i] = (i == (grd.nx-1)) ? grd.x_max - d * grd.x_inc: grd.x_min + (i + d) * grd.x_inc;
  142.     }
  143.  
  144.     if (gmtdefs.io_header) {
  145.         if (!grd.x_units[0]) strcpy (grd.x_units, "x");
  146.         if (!grd.y_units[0]) strcpy (grd.y_units, "y");
  147.         if (!grd.z_units[0]) strcpy (grd.z_units, "z");
  148.         if (gmtdefs.xy_toggle)
  149.             printf ("%s\t%s\t%s\n", grd.y_units, grd.x_units, grd.z_units);
  150.         else
  151.             printf ("%s\t%s\t%s\n", grd.x_units, grd.y_units, grd.z_units);
  152.     }
  153.     
  154.     if (z_only) {
  155.         for (ij = 0; ij < nm; ij++) {
  156.             out[0] = z[ij];
  157.             gmt_output (stdout, 1, out);
  158.         }
  159.     }
  160.     else  {
  161.         for (j = ij = 0; j < grd.ny; j++) for (i = 0; i < grd.nx; i++, ij++) {
  162.             out[ix] = x[i];    out[iy] = y[j];    out[2] = z[ij];
  163.             gmt_output (stdout, 3, out);
  164.         }
  165.     }
  166.     
  167.     free ((char *)z);
  168.     
  169.     if (!z_only) {
  170.         free ((char *)x);
  171.         free ((char *)y);
  172.     }
  173.     
  174.     gmt_end (argc, argv);
  175. }
  176.