home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Science / Science.zip / gmt_os2.zip / src / xyz2grd.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-28  |  8.5 KB  |  293 lines
  1. /*--------------------------------------------------------------------
  2.  *    The GMT-system:    @(#)xyz2grd.c    2.23  4/28/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.  * xyz2grd reads a xyz file from standard input and creates the
  9.  * corresponding grd-file. The input file does not have to be
  10.  * sorted. xyz2grd will report if some nodes are missing. nx and
  11.  * ny are computed as:
  12.  *    nx = rint((east-west)/dx) + 1, ny = rint((north-south)/dy) + 1
  13.  *
  14.  * Author:    Paul Wessel
  15.  * Date:    3-JAN-1991-1995
  16.  * Version:    2.0    Based on old v1.x
  17.  */
  18.  
  19. #include "gmt.h"
  20.  
  21. int *flag;
  22. float *a;
  23.  
  24. main (argc, argv)
  25. int argc;
  26. char **argv; {
  27.     int error = FALSE, got_input = FALSE, pixel = FALSE, z_only = FALSE, more, binary = FALSE, skip;
  28.     int single_precision = FALSE, xy_units = TRUE;
  29.     
  30.     int i, ij, x, y, z, nm, n_read = 0, n_filled = 0, n_used = 0, one_or_zero, dummy[4];
  31.     int n_empty = 0, n_stuffed = 0, n_bad = 0, entry, n_in, n_fields;
  32.     
  33.     double w, e, s, n, dx = 0.0, dy = 0.0, *in, no_data;
  34.     double off, idx, idy;
  35.     
  36.     char *grdfile, line[512], input[512], *ptr;
  37.     
  38.     FILE *fpi = NULL;
  39.     
  40.     struct GRD_HEADER grd;
  41.     
  42.     grdfile = CNULL;
  43.     input[0] = 0;
  44.     w = e = s = n = 0.0;
  45.     dummy[0] = dummy[1] = dummy[2] = dummy[3] = 0;
  46.     
  47.     argc = gmt_begin (argc, argv);
  48.  
  49.     no_data = gmt_NaN;
  50.     
  51.     for (i = 1; i < argc; i++) {
  52.         if (argv[i][0] == '-') {
  53.             switch (argv[i][1]) {
  54.                 /* Common parameters */
  55.             
  56.                 case 'H':
  57.                 case 'R':
  58.                 case 'V':
  59.                 case ':':
  60.                 case '\0':
  61.                     error += get_common_args (argv[i], &w, &e, &s, &n);
  62.                     break;
  63.  
  64.                 case 'b':
  65.                     single_precision = !(argv[i][strlen(argv[i])-1] == 'd');
  66.                     binary = TRUE;
  67.                     break;
  68.                 case 'D':
  69.                     strcpy (input, &argv[i][2]);
  70.                     got_input = TRUE;
  71.                     break;
  72.                 case 'F':
  73.                     pixel = TRUE;
  74.                     break;
  75.                 case 'G':
  76.                     grdfile = &argv[i][2];
  77.                     break;
  78.                 case 'I':
  79.                     gmt_getinc (&argv[i][2], &dx, &dy);
  80.                     break;
  81.                 case 'L':
  82.                     xy_units = FALSE;
  83.                     break;
  84.                 case 'N':
  85.                     if (!argv[i][2]) {
  86.                         fprintf (stderr, "%s: GMT SYNTAX ERROR -N option:  Must specify value or NaN\n", gmt_program);
  87.                         error++;
  88.                     }
  89.                     else
  90.                         no_data = (!strcmp (&argv[i][2], "NaN")) ? gmt_NaN : atof (&argv[i][2]);
  91.                     break;
  92.                 case 'Z':
  93.                     z_only = TRUE;
  94.                     break;
  95.                 default:
  96.                     error = TRUE;
  97.                     gmt_default_error (argv[i][1]);
  98.                     break;
  99.             }
  100.         }
  101.         else if ((fpi = fopen (argv[i], "r")) == NULL) {
  102.             fprintf (stderr, "xyz2grd: Cannot find ascii file %s\n", argv[i]);
  103.             exit (-1);
  104.         }
  105.     }
  106.         
  107.     if (argc == 1 || gmt_quick) {
  108.         fprintf (stderr, "xyz2grd %s - Converting an ASCII xyz-file to a netCDF grdfile\n\n", GMT_VERSION);
  109.         fprintf (stderr, "usage: xyz2grd [<xyzfile>] -G<grdfile> -I<dx[m|c]>[/<dy[m|c]>] -R<west/east/south/north>\n");
  110.         fprintf (stderr, "    [-D<xunit/yunit/zunit/scale/offset/title/remark>] [-F] [-H] [-L] [-N<nodata>] [-V] [-Z] [-:] [-b[d]]\n");
  111.         
  112.         if (gmt_quick) exit (-1);
  113.         
  114.         fprintf (stderr, "    xyzfile is a 3-column ascii file [or standard input] with x,y,z values\n");
  115.         fprintf (stderr, "    -G to name the output grdfile.\n");
  116.         fprintf (stderr, "    -I specifies grid size(s).  Append m [or c] to <dx> and/or <dy> for minutes [or seconds]\n");
  117.         explain_option ('R');
  118.         fprintf (stderr, "\n\tOPTIONS:\n");
  119.         fprintf (stderr, "    -D to enter information.  Specify '=' to get default value\n");
  120.         fprintf (stderr, "    -F will force pixel registration [Default is grid registration]\n");
  121.         explain_option ('H');
  122.         fprintf (stderr, "    -L means that x is longitude\n");
  123.         fprintf (stderr, "    -N set value for nodes without data [Default is NaN]\n");
  124.         explain_option ('V');
  125.         fprintf (stderr, "    -Z Input data is z-array in scanline orientation (no x,y pair)\n");
  126.         fprintf (stderr, "       This assumes all nodes have data values.\n");
  127.         explain_option (':');
  128.         fprintf (stderr, "    -b Input data is binary [Default is ascii].\n");
  129.         fprintf (stderr, "\t   Append d for double precision [Default is single precision].\n");        explain_option ('.');
  130.         explain_option ('.');
  131.         exit (-1);
  132.     }
  133.     
  134.     if (!project_info.region_supplied) {
  135.         fprintf (stderr, "%s: GMT SYNTAX ERROR:  Must specify -R option\n", gmt_program);
  136.         error++;
  137.     }
  138.     if (dx <= 0.0 || dy <= 0.0) {
  139.         fprintf (stderr, "%s: GMT SYNTAX ERROR -I option.  Must specify positive increment(s)\n", gmt_program);
  140.         error++;
  141.     }
  142.     if (!grdfile) {
  143.         fprintf (stderr, "%s: GMT SYNTAX ERROR option -G:  Must specify output file\n", gmt_program);
  144.         error++;
  145.     }    
  146.     if (binary && gmtdefs.io_header) {
  147.         fprintf (stderr, "%s: GMT SYNTAX ERROR.  Binary input data cannot have header -H\n", gmt_program);
  148.         error++;
  149.     }
  150.                 
  151.     if (error) exit (-1);
  152.  
  153.     if (binary) gmt_input = (single_precision) ? bin_float_input : bin_double_input;
  154.  
  155.     if (fpi == NULL) fpi = stdin;
  156.  
  157.     grd_init (&grd, argc, argv, FALSE);
  158.     
  159.     /* Decode grd information given, if any */
  160.     
  161.     if (got_input) {
  162.         for (i = 0; input[i]; i++) if (input[i] == '_') input[i] = ' ';
  163.         ptr = strtok (input, "/");
  164.         entry = 0;
  165.         while (ptr) {
  166.             switch (entry) {
  167.                 case 0:
  168.                     if (ptr[0] != '=') strcpy (grd.x_units, ptr);
  169.                     break;
  170.                 case 1:
  171.                     if (ptr[0] != '=') strcpy (grd.y_units, ptr);
  172.                     break;
  173.                 case 2:
  174.                     if (ptr[0] != '=') strcpy (grd.z_units, ptr);
  175.                     break;
  176.                 case 3:
  177.                     if (ptr[0] != '=') grd.z_scale_factor = atof (ptr);
  178.                     break;
  179.                 case 4:
  180.                     if (ptr[0] != '=') grd.z_add_offset = atof (ptr);
  181.                     break;
  182.                 case 5:
  183.                     if (ptr[0] != '=') strcpy (grd.title, ptr);
  184.                     break;
  185.                 case 6:
  186.                     if (ptr[0] != '=') strcpy (grd.remark, ptr);
  187.                     break;
  188.                 default:
  189.                     break;
  190.             }
  191.             ptr = strtok (CNULL, "/");
  192.             entry++;
  193.         }
  194.     }
  195.     
  196.     if (pixel) {
  197.         grd.node_offset = 1;
  198.         one_or_zero = 0;
  199.         off = 0.0;
  200.     }
  201.     else {
  202.         grd.node_offset = 0;
  203.         one_or_zero = 1;
  204.         off = 0.5;
  205.     }
  206.  
  207.     grd.nx = rint ((e-w)/dx) + one_or_zero;
  208.     grd.ny = rint ((n-s)/dy) + one_or_zero;
  209.     grd.x_min = w;    grd.x_max = e;
  210.     grd.y_min = s;    grd.y_max = n;
  211.     grd.x_inc = dx;    grd.y_inc = dy;
  212.     
  213.     if (gmtdefs.verbose) fprintf (stderr, "xyz2grd: nx = %d  ny = %d\n", grd.nx, grd.ny);
  214.     
  215.     nm = grd.nx * grd.ny;
  216.     
  217.     a = (float *) memory (CNULL, nm, sizeof (float), "xyz2grd");
  218.     flag = (int *) memory (CNULL, nm, sizeof (int), "xyz2grd");
  219.  
  220.     x = (gmtdefs.xy_toggle) ? 1 : 0;    y = 1 - x;        /* Set up which columns have x and y */
  221.  
  222.     if (gmtdefs.io_header) for (i = 0; i < gmtdefs.n_header_recs; i++) fgets (line, 512, fpi);
  223.     
  224.     idx = 1.0 / dx;    idy = 1.0 / dy;
  225.     n_in = (z_only) ? 1 : 3;
  226.     z = (z_only) ? 0 : 2;
  227.     
  228.     more = ((n_fields = gmt_input (fpi,  n_in, &in)) == n_in);
  229.         
  230.     ij = -1;    /* Will be incremented to 0 or set first time around */
  231.     
  232.     while (more) {
  233.         skip = FALSE;
  234.         n_read++;
  235.         
  236.         if (!z_only) {
  237.             if (in[y] >= s && in[y] <= n) {    /* y ok, check x */
  238.                 if (xy_units && (in[x] < w || in[x] > e))
  239.                     skip = TRUE;
  240.                 else {
  241.                     while (in[x] < w) in[x] += 360.0;
  242.                     if (in[x] > e) skip = TRUE;
  243.                 }
  244.             }
  245.             else
  246.                 skip = TRUE;
  247.         }
  248.         
  249.         if (!skip) {
  250.             ij = (z_only) ? ij + 1 : floor ((n - in[y]) * idy + off) * grd.nx + floor ((in[x] - w) * idx + off);
  251.             a[ij] = in[z];
  252.             flag[ij]++;
  253.             n_used++;
  254.         }
  255.         
  256.         more = ((n_fields = gmt_input (fpi,  n_in, &in)) == n_in);
  257.     }
  258.     if (fpi != stdin) fclose (fpi);
  259.     if (n_fields == -1) n_fields = 0;    /* Ok to get EOF */
  260.     if (n_fields%n_in) {    /* Got garbage or multiple segment subheader */
  261.         fprintf (stderr, "%s:  Cannot read line %d, aborts\n", gmt_program, n_read);
  262.         exit (-1);
  263.     }
  264.     
  265.     for (ij = 0; ij < nm; ij++) {    /* Check if all nodes got one value only */
  266.         if (flag[ij] == 1)
  267.             n_filled++;
  268.         else if (flag[ij] == 0) {
  269.             n_empty++;
  270.             a[ij] = no_data;
  271.         }
  272.         else    /* More than 1 value went to this node */
  273.             n_stuffed++;
  274.     }
  275.     if (n_read != n_used || n_filled != nm || n_empty > 0) gmtdefs.verbose = TRUE;
  276.     if (gmtdefs.verbose) {
  277.         fprintf (stderr, "xyz2grd: %d  n_used: %d  n_filled: %d n_empty: %d set to ",
  278.             n_read, n_used, n_filled, n_empty);
  279.         (bad_float (no_data)) ? fprintf (stderr, "NaN\n") : fprintf (stderr, "%lg\n", no_data);
  280.     }
  281.     if (n_bad) fprintf (stderr, "xyz2grd: %d records unreadable\n", n_bad);
  282.     if (n_stuffed) fprintf (stderr, "xyz2grd: Warning - %d nodes had multiple entries\n", n_stuffed);
  283.     if (write_grd (grdfile, &grd, a, 0.0, 0.0, 0.0, 0.0, dummy, FALSE)) {
  284.         fprintf (stderr, "xyz2grd: Error writing file %s\n", grdfile);
  285.         exit (-1);
  286.     }
  287.     
  288.     free ((char *)a);
  289.     free ((char *)flag);
  290.     
  291.     gmt_end (argc, argv);
  292. }
  293.